鲁棒张量主成分分析的非凸框架
The Nonconvex Framework for Robust Tensor Principal Component Analysis
摘要: 非凸鲁棒张量主成分分析问题包括从被噪声破坏的张量中恢复低秩和稀疏部分,这在广泛的实际应用中引起了极大的关注。然而,现有的非凸方法面临许多问题,其中最重要的两个问题是对特定非凸函数的限制和低秩部分的信息损失。在本文,我们提出了一种广义非凸鲁棒张量主成分分析模型(N-RTPCA),其中包括一些最常用的非凸函数。并且提出了一个非凸ADMM算法来求解广义非凸鲁棒张量主成分分析模型(N-RTPCA)。最后,实验部分通过模拟实验和真实图片的实验验证了所提方法的优越性。
Abstract: The problem of nonconvex robust tensor principal component analysis, which consists of recovering the low-rank part and sparse part of a tensor corrupted by noise, has attracted great attention in a wide range of practical applications. However, existing nonconvex methods face many problems, the two most important of which are the restriction to specific nonconvex functions and the loss of information in the low-rank part. In this paper, we propose a generalized nonconvex robust tensor principal component analysis model (N-RTPCA) which includes some of the most commonly used nonconvex functions. And a nonconvex ADMM algorithm is proposed to solve the generalized nonconvex robust tensor principal component analysis model (N-RTPCA). Finally, the experimental part verifies the superiority of the proposed method by simulation experiments and experiments on real pictures.
文章引用:唐开煜, 樊亚莉. 鲁棒张量主成分分析的非凸框架[J]. 建模与仿真, 2024, 13(4): 4171-4179. https://doi.org/10.12677/mos.2024.134378

1. 引言

张量是多维数组或矩阵的扩展,常用于表示和处理多维数据,目前已被广泛应用于实践中[1]。然而,在实际应用中,张量的获取往往并不完美,通常会受到异常值或非高斯噪声(如稀疏噪声)的严重破坏,因此,找到能处理带有噪声的张量数据的去噪方法至关重要。

鲁棒张量主成分分析(RTPCA),目的是从被稀疏噪声污染的张量中精确地恢复低秩部分,已被广泛用于实践中,见[2][3]。具体来说,假设 X 是一个被稀疏噪声破坏的张量,RTPCA指出该张量可以被分解为两部分,即 X=+ ,其中 是低秩张量, 是稀释张量。RTPCA的最小化问题表述如下:

min , rank( )+λ 1 ,s.t.X=+, (1)

其中 rank( ) 代表张量秩, 1 代表 1 -范数。要实现这种想法,张量秩的定义及其严格松弛就成了一个不可忽视的问题。一种被称为张量乘积(t-product) [4]的新方法被引入,它被作为矩阵乘法到张量的自然扩展。这一进步促成了一种新颖的张量分解框架—张量奇异值分解(t-SVD)的发展。与其他分解不同,t-SVD无需对张量进行矩阵化即可进行计算。在t-SVD的基础上,[5]提出了张量管道秩的概念,以解决张量恢复问题。从本质上讲,张量管道秩侧重于利用奇异值信息来揭示张量的内在结构。为了解决RTPCA问题,文献[6]引入了张量平均秩的概念,证明了它与张量管道秩的密切联系。此外,他们还定义了一种新的张量核范数,并在张量谱范数的单位球内建立了其凸包特性。那么,(1)中的RTPCA问题可以放宽为

min , +λ 1 ,s.tX=+, (2)

其中 表示张量核范数。模型(2)可用交替方向乘法(ADMM) [7]求解。该研究首次将RPCA直接扩展到三向张量,并包含了矩阵情况下的所有特性。为了和统一的RTPCA问题区分,将模型(2)记为TRPCA。此外,张量管道秩和张量核范数已被广泛应用于张量恢复[8]和图像聚类[9]等众多领域。同时,一些基于TNN的改进方法也被提出来提高TRPCA的性能。文献[10]提出了一种基于低秩核矩阵的张量管道秩的近似凸替代(IRTPCA)。然而,尽管张量核范数在RTPCA问题上取得了很好的效果,TRPCA的这些改进方法在张量互补以及张量去噪方面也取得了很好的性能,但它们仍然存在一些不可避免的缺点。

首先,需要注意的是,张量核范数是张量平均秩的紧凸包络,而不是张量管道秩的凸包。因此,使用这种松弛方法可能会导致(1)中原始RTPCA问题的次优解,从而与最小化管道秩产生差异。为了解决这一局限性,人们在张量恢复问题中广泛探索了非凸松弛技术[11]-[13]。与经典的凸方法相比,这些非凸方法表现出更高的精确度,突出了它们实现更精确解的潜力。例如,文献[14]提出了一个专门针对管道秩的广义非凸框架,旨在解决张量补全问题。他们为这一框架内的精确和鲁棒恢复提供了理论保证。然而,他们并没有考虑张量主成分分析问题(RTPCA问题)。此外,现有的非凸算法主要适用于特定的非凸函数[15] [16]。因此,有必要为RTPCA问题建立一个通用、统一的非凸框架。

本文重点研究了非凸框架下的鲁棒张量主成分分析问题。主要贡献有三个方面。

首先,本文提出了广义非凸鲁棒张量主成分分析模型(N-RTPCA),它保证了满足特定条件的非凸惩罚函数的普遍性,比经典TRPCA获得了更好的低阶解。

其次,文本提出了一个新的非凸ADMM算法来求解N-RTPCA,并给出了算法的计算复杂度分析。

最后,实验验证了N-RTPCA模型具有良好的张量恢复效果。同时,我们将这模型应用于实际图片数据,结果证实了所提出的N-RTPCA模型优于经典TRPCA模型和较先进的IRTPCA方法。

2. 相关知识

为了更清晰的表述本文的工作,需要引入一些记号。对于 X 1 , X 2 n 1 × n 2 × n 3 ,张量内积被定义为 X 1 , X 2 = i X 1 (i) , X 2 (i) 。然后,沿着张量 X 第三维度的离散傅里叶变换(DFT)记为 X ¯ =fft( X,[],3 ) ,同时逆DFT记为 X=ifft( X ¯ ,[],3 ) 。接下来,沿着张量第三维度矩阵化的算子被定义为 unfold( X )=[ X (1) ; X (2) ;; X ( n 3 ) ] ,该算子的逆算子被定义为 fold( unfold( X ) )=X

定义1 (T-product [4]) 假设 X 1 n 1 × n 2 × n 3 and X 2 n 2 ×l× n 3 ,它们的t-product被定义为一个大小为 n 1 ×l× n 3 的张量,记为 X 1 X 2 :=fold( bcirc( X 1 )unfold( X 1 ) )

引理1 (T-SVD [6]) 对于任意张量 X n 1 × n 2 × n 3 ,它存在t-SVD如下式:

X=US V

其中 U n 1 × n 1 × n 3 V n 2 × n 2 × n 3 是正交张量, S n 1 × n 2 × n 3 是一个F-对角张量。

为了介绍所提出的非凸最小化模型,需要介绍张量管道秩的统一非凸替代、几个相关定义以及张量管道核范数的部分和。

Table 1. Some of the most commonly used non-convex penalty functions

1. 一些最常用的非凸惩罚函数

非凸函数

表达式( x0,β>0,γ>0 )

Laplace

β( 1exp( x γ ) )

Logarithm

β log( γ+1 ) log( γx+1 )

L p

β x p ,0<p<1

MCP

{ βx x 2 2γ , xγβ, 1 2 γ β 2 , x>γβ.

Capped L 1

{ βx, x<γ, βγ, xγ.

SCAD

{ βx, x<β, x 2 +2γβx β 2 2( γ1 ) , βx<γβ, β 2 ( γ+1 ) 2 , xγβ.

定义2 (统一非凸替代[14]) 假设 X=US V 是张量 X n 1 × n 2 × n 3 的t-SVD,其张量管道秩的统一非凸替代定义为

Ψ( X ):= 1 n 3 i=1 m j=1 n 3 ψ ( σ ij ) (3)

其中 ψ: + + 是一个连续、单调非递减的凹函数,且满足

ψ( 0 )=0, lim x+ ψ( x )/x =0

表1总结了文献[14]中提到的一系列常用的非凸惩罚函数。

定义3 (加权t-TNN [14]) 对于张量 X n 1 × n 2 × n 3 ,它的加权t-TNN被定义为 X 在傅里叶域上所有前向切片的所有奇异值加权和,表示如下

X W, := 1 n 3 i=1 m j=1 n 3 ω ij σ ij (4)

其中 W= ( ω ij ) m× n 3

定义4 (加权T-SVT [14]) 对于任意常数 τ>0 以及张量 Y n 1 × n 2 × n 3 ,假设其t-SVD为 X=US V ,其加权T-SVT (WTSVT)算子被定义为

D W,τ ( Y ):=U S W,τ V (5)

其中 S W,τ =ifft( ( S ¯ τW ),[],3 ) W n 1 × n 2 × n 3 是一个F-对角张量且它的第j个前向切片的对角元素等于权重矩阵 W 的第j列元素。

3. 基于广义非凸框架的鲁棒张量主成分分析模型

本节通过[14]中的广义非凸替代给出RTPCA的广义非凸框架。根据经典的基于TNN的RTPCA问题,对于原始张量 X n 1 × n 2 × n 3 来说,它可以自然地由低秩张量 和稀疏张量 组成。要从 X 中恢复理想的 ,需要求解最小化问题(2),其中张量核范数被用作张量管秩的凸替代。

3.1. 广义非凸鲁棒主成分分析模型

然而,与TNN相比,非凸框架能带来更严密的管道秩替代。考虑到非凸替代的通用性,我们将广义非凸鲁棒张量主成分分析模型称为N-RTPCA,并将该模型表述为

min , Ψ( )+λ 1 ,s.tX=+ (6)

从(6)可以看出,N-RTPCA作为一个广义的非凸框架,包含一系列张量核范数的非凸替代。例如,文献[15] [16]就是(6)的特例。

尽管文献[15]提出了一种加权张量SVT算子来求解广义非凸管道秩最小化,但其重点在于张量补全而非去噪,这导致其无法探究噪声的组成和分布。简而言之,与[15]中的管道秩最小化模型相比,模型(6)多了一个 1 -范数约束,其目的是从原始张量中恢复低秩部分和稀疏部分。为此,我们提出了一种使用WTSVT算子的ADMM算法。下面,我们将详细介绍ADMM求解(6)的方法。

3.2. 模型的优化算法

首先(6)的增广拉格朗日函数写法如下:

L μ ( ,, )= Ψ( )+λ 1 + ,+X + μ 2 +X F 2 (7)

其中 μ 是增广拉格朗日惩罚参数, 是对偶变量。在(7)的ADMM算法中,我们选择0.015作为惩罚参数的初始值,并在每一步迭代中更新它。

更新 :关于变量 的子问题表述为:

k+1 =arg min Ψ( )+ μ k 2 + k X+ k μ k F 2 (8)

直接解决该子问题比较困难,因此需要对其进行松弛。考虑到非凸函数的凹性,上述子问题可以放宽如下:

k+1 =arg min 1 n 3 i=1 m j=1 n 3 ( ψ( σ ij k )+ ω ij k ( σ ij σ ij k ) ) + μ k 2 + k X+ k μ k F 2 (9)

然后,去掉常数项,根据方程(4)可以得出方程(9)与下面的问题等价:

k+1 =arg min 1 μ k W, + 1 2 + k X+ k μ k F 2 (10)

根据相关理论[14],我们可以通过以下WTSVT算子求解方程(10):

k+1 = D W,τ ( X k k μ k ) (11)

更新 :关于变量 的子问题表述为:

k+1 =arg min λ 1 + μ k 2 + k X+ k μ k F 2 = S λ/ μ k ( X k k μ k ) (12)

其中 S τ ( ) 是软阈值算子,被定义为

S τ ( X )=sgn( X )max{ | X |τ,0 } (13)

对于参数 μ 的选择和更新,本文设定初始值 μ 0 =0.015 ,算法更新的第k步中的参数 μ k 通过下式计算:

μ k =ρ μ k1

其中 ρ=1.1 。最后,算法1总结了N-RTPCA的整个优化过程。

算法1 N-RTPCA的ADMM算法

输入:非凸函数 Ψ 的公式, ρ>1 X n 1 × n 2 × n 3 λ

初始化: 0 = 0 = 0 = W 0 =0 , μ 0 =0.015 , μ max =1e8 , ε=1e8 , k=0

While不收敛:

更新 k+1 :通过算子(11)求解问题(10)得到 k+1

更新 k+1 :通过算子(13)求解问题(12)得到 k+1

更新 W k+1 :通过 ω ij ψ( σ ij ( k+1 ) ) 得到 W k+1 .

k+1 = k + μ k ( k+1 + k+1 X )

续表

μ k+1 =min( ρ μ k , μ max )

检查收敛条件:

k+1 k ε, k+1 k ε, k+1 + k X ε

更新 k=k+1

end while

输出:低秩张量 ,稀疏张量

3.3. 计算复杂度分析

现在,给定一个原始张量 X n 1 × n 2 × n 3 ,我们分别给出算法1的计算复杂度分析。算法1的每一个迭代步需要分别更新 W 。根据公式(13),每一步中更新 的计算复杂度为 O( n 1 n 2 n 3 ) 。然后,由于更新 需要计算WTSVT算子,而WTSVT算子已经包含了 W 的更新,更新 的计算复杂度与

同,因此算法1每个迭代步骤的主要计算复杂度来自于更新 ,这需要计算DFT和 n 3 +1 2 n 1 × n 2

阵的SVDs。总的来说,算法1每一个迭代步中的计算复杂度为 O( n 1 n 2 n 3 log n 3 + n (1) n (2) 2 n 3 ) ,其中 n (1) =max( n 1 , n 2 ) n (2) =min( n 1 , n 2 )

4. 实验

在这些实验中,我们将实验结果与目前最先进的基于t-SVD方向的张量主成分分析方法进行了比较,即TRPCA [6]、IRTPCA [10]。上述方法的超参数都是根据理论值或实验最优值设定的,而我们提出的方法的超参数设置为理论值 λ=1/ max( n 1 , n 2 ) n 3 。在所有实验中,我们从常用的满足属性[14]的非凸惩罚函数中选择了性能非最佳的SCAD函数( γ=100 β=0.95 )。

所有实验均在装有Windows 10和MATLAB (R2019b)的电脑上进行,CPU为英特尔酷睿i5-7300HQ 2.50 GHz,内存为8 GB。

4.1. 生成实验

我们首先通过对生成的张量数据进行数值实验验证了所提方法的合理性和收敛性能,然后验证了所提N-RTPCA对具有不同管道秩和不同稀疏噪声的张量的恢复能力。我们采用与文献[6]相同的生成方法,随机生成一个大小为 n×n×n 、管道秩为r的张量 0 ,并分别取 n=50 n=100 0 的支持集 Ω (大小为 m )是均匀随机选择的。为简单起见,张量 0 的稀疏度(即 0 上的非零元素个数)为 m

表2列出了两种拟议方法在不同管道等级和稀疏性条件下的恢复结果。具体而言,N-RTPCA的结果用 ^ ^ 表示。可以看出,N-RTPCA在所有情况下都能提供正确的秩估计,而且相对误差都小于10e−7。同时,N-RTPCA恢复稀疏部分的相对误差小于10e−8,它在估计 0 方面也较为成功(在噪声稀疏时,能近似到 m )。

4.2. 在Berkeley Segmentation数据集上的恢复实验

在本节中,为了评估所提出的方法在张量恢复方面的性能,我们将N-RTPCA应用于恢复被随机噪声严重破坏的真实图像。由于彩色图像的R、G、B三通道可以定义为张量的三个前向切片,一个公认的事实是彩色图像可以通过低秩矩阵或张量获得。实验表明,在恢复真实彩色图像数据方面,我们提出的方法优于目前最先进的几种方法。

Table 2. Correct recovery for randomly generated tensors of different sizes

2. 对不同大小的随机生成张量的正确恢复

n

m

r t ( ^ )

^ 0

^ 0 F 0 F

^ 0 F 0 F

r t ( 0 )=0.1n,m= 0 0 =0.05 n 3

50

6250

5

6431

3.08e−7

2.44e−9

100

5e4

10

53,214

5.14e−7

2.52e−9

r t ( 0 )=0.1n,m= 0 0 =0.1 n 3

50

12,500

5

13,873

3.42e−7

3.13e−9

100

1e5

10

130,449

5.98e−7

3.39e−9

r t ( 0 )=0.2n,m= 0 0 =0.05 n 3

50

6250

10

13,168

3.01e−7

8.95e−9

100

5e4

20

84,998

5.16e−7

7.15e−9

r t ( 0 )=0.2n,m= 0 0 =0.1 n 3

50

12,500

10

16,219

3.06e−7

5.89e−9

100

1e5

20

126,464

5.35e−7

4.99e−9

我们从伯克利分割数据集(BSD) [17]中随机选取了20幅大小为321 × 481 × 3的彩色图像进行恢复实验。首先,我们用“CR”表示损坏率,它是张量中被稀疏噪声损坏的条目的比率。其次,对于每幅图像,我们分别将20%、25%和30%的像素随机损坏为[0, 255]中的值。然后,针对带有不同比例噪声的图像,我们使用TRPCA、IRTPCA和提出的新方法进行了恢复实验。最后,我们使用峰值信噪比(PSNR)和结构相似性(SSIM)评估了可恢复性。众所周知,PSNR和SSIM值越高,相应方法的恢复能力就越强。

Table 3. Comparison of image “Lake” recovery when CR is 20%, 25% and 30% respectively

3. CR分别为20%、25%和30%时,图像“Lake”的恢复情况对比

方法

CR = 20%

CR = 25%

CR = 30%

PSNR

SSIM

PSNR

SSIM

PSNR

SSIM

TRPCA

27.12

0.9152

25.99

0.8875

24.85

0.8480

IRTPCA

27.95

0.9225

25.88

0.8823

24.46

0.8253

N-RTPCA

28.23

0.9278

26.93

0.8974

25.42

0.8472

图1显示了CR = 25%时三幅图像的恢复结果。为了突出我们提出的方法的优越性,表3详细介绍了“湖泊”图像在三种不同噪声率下的恢复结果。最大的PSNR和SSIM值以粗体标出。首先,图1中的可视化示例表明,与TRPCA相比,基于管道秩的非凸方法恢复的图像具有更清晰的结构和纹理。其次,从表3所示的PSNR和SSIM值中,我们可以得出结论:在大多数真实图像数据上,所提出的方法在恢复性能上优于其他方法。这要归功于提出的方法对张量管道秩进行了更严格的替换,从而保留了张量低秩部分的更多信息。因此,我们提出的方法具有更强的抗噪声能力。

Figure 1. Comparison of the recovery performance of three example images corrupted by 25% random noise. From top to bottom, geese, lakes and starfish. (a) Original image. (b) Observed image. (c)~(e) Images recovered with TRPCA, IRTPCA and N-RTPCA, respectively

1. 被25%随机噪音破坏的三幅示例图像的恢复性能比较。从上到下依次为鹅、湖泊和海星。(a) 原始图像。(b) 观察到的图像。(c)~(e) 分别用TRPCA、IRTPCA和N-RTPCA恢复的图像

5. 总结

本章采用广义非凸方法解决RTPCA问题,这项工作不仅将广义非凸框架很好地应用于RTPCA问题,包括算法和理论的扩展,而且在开发改进模型的同时,还实现一些实验证明了模型在实际应用中的优越性。具体地,非凸RTPCA对张量管道秩有更优的近似,通过求解N-RTPCA问题,能得到更接近于原始张量的低秩张量。与此同时,通过我们的非凸ADMM算法求解N-RTPCA问题,能得到N-RTPCA问题的最优解。不过,虽然这项工作成功地促进了t-SVD框架下RTPCA问题的发展,但这项工作的理论方面研究甚少。探索基于张量非相干性条件的非凸理论恢复保证也是一个极其复杂但可行的方向。

参考文献

[1] Kolda, T.G. and Bader, B.W. (2009) Tensor Decompositions and Applications. SIAM Review, 51, 455-500.
https://doi.org/10.1137/07070111xs
[2] Lu, C., Feng, J., Chen, Y., Liu, W., Lin, Z. and Yan, S. (2016). Tensor Robust Principal Component Analysis: Exact Recovery of Corrupted Low-Rank Tensors via Convex Optimization. 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), Las Vegas, 27-30 June 2016, 5249-5257.
https://doi.org/10.1109/cvpr.2016.567
[3] Chen, L., Liu, Y. and Zhu, C. (2017). Iterative Block Tensor Singular Value Thresholding for Extraction of Lowrank Component of Image Data. 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), New Orleans, 5-9 March 2017, 1862-1866.
https://doi.org/10.1109/icassp.2017.7952479
[4] Kilmer, M.E. and Martin, C.D. (2011) Factorization Strategies for Third-Order Tensors. Linear Algebra and Its Applications, 435, 641-658.
https://doi.org/10.1016/j.laa.2010.09.020
[5] Kilmer, M.E., Braman, K., Hao, N. and Hoover, R.C. (2013) Third-Order Tensors as Operators on Matrices: A Theoretical and Computational Framework with Applications in Imaging. SIAM Journal on Matrix Analysis and Applications, 34, 148-172.
https://doi.org/10.1137/110837711
[6] Lu, C., Feng, J., Chen, Y., Liu, W., Lin, Z. and Yan, S. (2020) Tensor Robust Principal Component Analysis with a New Tensor Nuclear Norm. IEEE Transactions on Pattern Analysis and Machine Intelligence, 42, 925-938.
https://doi.org/10.1109/tpami.2019.2891760
[7] Lu, C., Feng, J., Yan, S. and Lin, Z. (2018) A Unified Alternating Direction Method of Multipliers by Majorization Minimization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 40, 527-541.
https://doi.org/10.1109/tpami.2017.2689021
[8] Zhang, F., Wang, J., Wang, W. and Xu, C. (2021) Low-Tubal-Rank Plus Sparse Tensor Recovery with Prior Subspace Information. IEEE Transactions on Pattern Analysis and Machine Intelligence, 43, 3492-3507.
https://doi.org/10.1109/tpami.2020.2986773
[9] Zhou, P., Lu, C., Feng, J., Lin, Z. and Yan, S. (2021) Tensor Low-Rank Representation for Data Recovery and Clustering. IEEE Transactions on Pattern Analysis and Machine Intelligence, 43, 1718-1732.
https://doi.org/10.1109/tpami.2019.2954874
[10] Liu, Y., Chen, L. and Zhu, C. (2018) Improved Robust Tensor Principal Component Analysis via Low-Rank Core Matrix. IEEE Journal of Selected Topics in Signal Processing, 12, 1378-1389.
https://doi.org/10.1109/jstsp.2018.2873142
[11] Shi, Q., Cheung, Y. and Lou, J. (2022) Robust Tensor SVD and Recovery with Rank Estimation. IEEE Transactions on Cybernetics, 52, 10667-10682.
https://doi.org/10.1109/tcyb.2021.3067676
[12] Xu, W., Zhao, X., Ji, T., Miao, J., Ma, T., Wang, S., et al. (2019) Laplace Function Based Nonconvex Surrogate for Low-Rank Tensor Completion. Signal Processing: Image Communication, 73, 62-69.
https://doi.org/10.1016/j.image.2018.11.007
[13] Sun, X., Zhang, X., Xu, C., Xiao, M. and Tang, Y. (2023) Tensorial Multiview Representation for Saliency Detection via Nonconvex Approach. IEEE Transactions on Cybernetics, 53, 1816-1829.
https://doi.org/10.1109/tcyb.2021.3139037
[14] Wang, H., Zhang, F., Wang, J., Huang, T., Huang, J. and Liu, X. (2022) Generalized Nonconvex Approach for Low-Tubal-Rank Tensor Recovery. IEEE Transactions on Neural Networks and Learning Systems, 33, 3305-3319.
https://doi.org/10.1109/tnnls.2021.3051650
[15] Cai, S., Luo, Q., Yang, M., Li, W. and Xiao, M. (2019) Tensor Robust Principal Component Analysis via Non-Convex Low Rank Approximation. Applied Sciences, 9, Article 1411.
https://doi.org/10.3390/app9071411
[16] Li, T. and Ma, J. (2021). T-SVD Based Non-Convex Tensor Completion and Robust Principal Component Analysis. 2020 25th International Conference on Pattern Recognition (ICPR), Milan, 10-15 January 2021, 6980-6987.
https://doi.org/10.1109/icpr48806.2021.9412248
[17] Martin D, Fowlkes C, Tal D, et al. (2001) A Database of Human Segmented Natural Images and Its Application to Evaluating Segmentation Algorithms and Measuring Ecological Statistics. Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001, Vancouver, 7-14 July 2001, 416-423.
https://doi.org/10.1109/ICCV.2001.937655