http://blog.csdn.net/chenyusiyuan/article/details/8710462

OpenCV学习笔记（27）KAZE 算法原理与源码分析（一）非线性扩散滤波

2013-03-23 17:44 16963人阅读 评论(28) 收藏 举报

分类：

机器视觉（34）

版权声明：本文为博主原创文章，未经博主允许不得转载。

目录(?)[+]

KAZE系列笔记：

1. OpenCV学习笔记（27）KAZE 算法原理与源码分析（一）非线性扩散滤波
2. OpenCV学习笔记（28）KAZE 算法原理与源码分析（二）非线性尺度空间构建
3. OpenCV学习笔记（29）KAZE 算法原理与源码分析（三）特征检测与描述
4. OpenCV学习笔记（30）KAZE 算法原理与源码分析（四）KAZE特征的性能分析与比较
5. OpenCV学习笔记（31）KAZE 算法原理与源码分析（五）KAZE的性能优化及与SIFT的比较

KAZE算法资源：

1. 论文：  http://www.robesafe.com/personal/pablo.alcantarilla/papers/Alcantarilla12eccv.pdf
2. 项目主页：http://www.robesafe.com/personal/pablo.alcantarilla/kaze.html
3. 作者代码：http://www.robesafe.com/personal/pablo.alcantarilla/code/kaze_features_1_4.tar 
   （需要boost库，另外其计时函数的使用比较复杂，可以用OpenCV的cv::getTickCount代替）
4. Computer Vision Talks的评测：http://computer-vision-talks.com/2013/03/porting-kaze-features-to-opencv/
5. Computer Vision Talks 博主Ievgen Khvedchenia将KAZE集成到OpenCV的cv::Feature2D类，但需要重新编译OpenCV，并且没有实现算法参数调整和按Mask过滤特征点的功能：https://github.com/BloodAxe/opencv/tree/kaze-features
6. 我在Ievgen的项目库中提取出KAZE，封装成继承cv::Feature2D的类，无需重新编译OpenCV，实现了参数调整和Mask过滤的功能: https://github.com/yuhuazou/kaze_opencv (2013-03-28更新，对KAZE代码进行了优化)
7. Matlab 版的接口程序，封装了1.0版的KAZE代码：https://github.com/vlfeat/vlbenchmarks/blob/unstable/%2BlocalFeatures/Kaze.m

简介

ECCV2012中出现了一种比SIFT更稳定的特征检测算法KAZE （[1]）。KAZE的取名是为了纪念尺度空间分析的开创者—日本学者Iijima。KAZE是日语‘风’的谐音，寓意是就像风的形成是空气在空间中非线性的流动过程一样，KAZE特征检测是在图像域中进行非线性扩散处理的过程。

传统的SIFT、SURF等特征检测算法都是基于线性的高斯金字塔进行多尺度分解来消除噪声和提取显著特征点。但高斯分解是牺牲了局部精度为代价的，容易造成边界模糊和细节丢失。非线性的尺度分解有望解决这种问题，但传统方法基于正向欧拉法（forward Euler scheme）求解非线性扩散（Non-linear diffusion）方程时迭代收敛的步长太短，耗时长、计算复杂度高。由此，KAZE算法的作者提出采用加性算子分裂算法(Additive Operator Splitting, AOS)来进行非线性扩散滤波，可以采用任意步长来构造稳定的非线性尺度空间。

注：KAZE算法的原理与SIFT和SURF有很多相似之处，在深入了解KAZE之前，可以参考以下的博客文章对SIFT和SURF的介绍：

1. 【OpenCV】SIFT原理与源码分析 - 小魏的修行路 （很不错的博客，用心、认真、图文并茂）
2. SIFT特征提取分析 - Rachel Zhang的专栏 (这里的分析也很详细)
3. Surf算法学习心得(一)——算法原理 - yang_yong’ blog
4. Opencv学习笔记（六）SURF学习笔记 - crazy_sparrow的专栏

1.1 非线性扩散滤波

1.1.1 Perona-Malik扩散方程

具体地，非线性扩散滤波方法是将图像亮度（L）在不同尺度上的变化视为某种形式的流动函数（flow function）的散度（divergence），可以通过非线性偏微分方程来描述：

通过设置合适的传导函数 c(x,y,t) ，可以使得扩散自适应于图像的局部结构。传导函数可以是标量、也可以是张量。时间t作为尺度参数，其值越大、则图像的表示形式越简单。Perona和Malik提出了传导函数的构造方式：

其中的▽L_σ是高斯平滑后的图像L_σ的梯度（gradient）。函数 g() 的形式有如下几种：

其中函数g₁优先保留高对比度的边缘，g₂优先保留宽度较大的区域，g₃能够有效平滑区域内部而保留边界信息（KAZE代码中默认采用函数g₂）。

函数g₁、g₂、g₃的实现代码如下（在文件 kaze_nldiffusion_functions.cpp 中）：

[cpp] view plain copy

1. //*************************************************************************************
2. //*************************************************************************************
3. /**
4. * @brief This function computes the Perona and Malik conductivity coefficient g1
5. * g1 = exp(-|dL|^2/k^2)
6. * @param src Input image
7. * @param dst Output image
8. * @param Lx First order image derivative in X-direction (horizontal)
9. * @param Ly First order image derivative in Y-direction (vertical)
10. * @param k Contrast factor parameter
11. */
12. void PM_G1(const cv::Mat &src, cv::Mat &dst, cv::Mat &Lx, cv::Mat &Ly, float k )
13. {
14. cv::exp(-(Lx.mul(Lx) + Ly.mul(Ly))/(k*k),dst);
15. }
16. //*************************************************************************************
17. //*************************************************************************************
18. /**
19. * @brief This function computes the Perona and Malik conductivity coefficient g2
20. * g2 = 1 / (1 + dL^2 / k^2)
21. * @param src Input image
22. * @param dst Output image
23. * @param Lx First order image derivative in X-direction (horizontal)
24. * @param Ly First order image derivative in Y-direction (vertical)
25. * @param k Contrast factor parameter
26. */
27. void PM_G2(const cv::Mat &src, cv::Mat &dst, cv::Mat &Lx, cv::Mat &Ly, float k )
28. {
29. dst = 1./(1. + (Lx.mul(Lx) + Ly.mul(Ly))/(k*k));
30. }
31. //*************************************************************************************
32. //*************************************************************************************
33. /**
34. * @brief This function computes Weickert conductivity coefficient g3
35. * @param src Input image
36. * @param dst Output image
37. * @param Lx First order image derivative in X-direction (horizontal)
38. * @param Ly First order image derivative in Y-direction (vertical)
39. * @param k Contrast factor parameter
40. * @note For more information check the following paper: J. Weickert
41. * Applications of nonlinear diffusion in image processing and computer vision,
42. * Proceedings of Algorithmy 2000
43. */
44. void Weickert_Diffusivity(const cv::Mat &src, cv::Mat &dst, cv::Mat &Lx, cv::Mat &Ly, float k )
45. {
46. cv::Mat modg;
47. cv::pow((Lx.mul(Lx) + Ly.mul(Ly))/(k*k),4,modg);
48. cv::exp(-3.315/modg, dst);
49. dst = 1.0 - dst;
50. }

参数k是控制扩散级别的对比度因子（contrast factor），能够决定保留多少边缘信息，其值越大，保留的边缘信息越少。在KAZE算法中，参数k的取值是梯度图像▽L_σ的直方图70% 百分位上的值：

计算参数k的实现源码如下（在文件 kaze_nldiffusion_functions.cpp 中）：

[cpp] view plain copy

1. //*************************************************************************************
2. //*************************************************************************************
3. /**
4. * @brief This function computes a good empirical value for the k contrast factor
5. * given an input image, the percentile (0-1), the gradient scale and the number of
6. * bins in the histogram
7. * @param img Input image
8. * @param perc Percentile of the image gradient histogram (0-1)
9. * @param gscale Scale for computing the image gradient histogram
10. * @param nbins Number of histogram bins
11. * @param ksize_x Kernel size in X-direction (horizontal) for the Gaussian smoothing kernel
12. * @param ksize_y Kernel size in Y-direction (vertical) for the Gaussian smoothing kernel
13. * @return k contrast factor
14. */
15. float Compute_K_Percentile(const cv::Mat &img, float perc, float gscale, unsigned int nbins, unsigned int ksize_x, unsigned int ksize_y)
16. {
17. float kperc = 0.0, modg = 0.0, lx = 0.0, ly = 0.0;
18. unsigned int nbin = 0, nelements = 0, nthreshold = 0, k = 0;
19. float npoints = 0.0;    // number of points of which gradient greater than zero
20. float hmax = 0.0;       // maximum gradient
21. // Create the array for the histogram
22. float *hist = new float[nbins];
23. // Create the matrices
24. cv::Mat gaussian = cv::Mat::zeros(img.rows,img.cols,CV_32F);
25. cv::Mat Lx = cv::Mat::zeros(img.rows,img.cols,CV_32F);
26. cv::Mat Ly = cv::Mat::zeros(img.rows,img.cols,CV_32F);
27. // Set the histogram to zero, just in case
28. for( unsigned int i = 0; i < nbins; i++ )
29. {
30. hist[i] = 0.0;
31. }
32. // Perform the Gaussian convolution
33. Gaussian_2D_Convolution(img,gaussian,ksize_x,ksize_y,gscale);
34. // Compute the Gaussian derivatives Lx and Ly
35. Image_Derivatives_Scharr(gaussian,Lx,1,0);
36. Image_Derivatives_Scharr(gaussian,Ly,0,1);
37. // Skip the borders for computing the histogram
38. for( int i = 1; i < gaussian.rows-1; i++ )
39. {
40. for( int j = 1; j < gaussian.cols-1; j++ )
41. {
42. lx = *(Lx.ptr<float>(i)+j);
43. ly = *(Ly.ptr<float>(i)+j);
44. modg = sqrt(lx*lx + ly*ly);
45. // Get the maximum
46. if( modg > hmax )
47. {
48. hmax = modg;
49. }
50. }
51. }
52. // Skip the borders for computing the histogram
53. for( int i = 1; i < gaussian.rows-1; i++ )
54. {
55. for( int j = 1; j < gaussian.cols-1; j++ )
56. {
57. lx = *(Lx.ptr<float>(i)+j);
58. ly = *(Ly.ptr<float>(i)+j);
59. modg = sqrt(lx*lx + ly*ly);     // modg can be stored in a matrix in last step in oder to avoid re-computation
60. // Find the correspondent bin
61. if( modg != 0.0 )
62. {
63. nbin = floor(nbins*(modg/hmax));
64. if( nbin == nbins )
65. {
66. nbin--;
67. }
68. hist[nbin]++;
69. npoints++;
70. }
71. }
72. }
73. // Now find the perc of the histogram percentile
74. nthreshold = (unsigned int)(npoints*perc);
75. // find the bin (k) in which accumulated points are greater than 70% (perc) of total valid points (npoints)
76. for( k = 0; nelements < nthreshold && k < nbins; k++)
77. {
78. nelements = nelements + hist[k];
79. }
80. if( nelements < nthreshold )
81. {
82. kperc = 0.03;
83. }
84. else
85. {
86. kperc = hmax*((float)(k)/(float)nbins);
87. }
88. delete hist;
89. return kperc;
90. }

注：有关非线性扩散滤波的应用，参见[2]。

1.1.2 AOS算法

由于非线性偏微分方程并没有解析解，一般通过数值分析的方法进行迭代求解。传统上采用显式差分格式的求解方法只能采用小步长，收敛缓慢。为此，将方程离散化为以下的隐式差分格式：

其中A_l是表示图像在各维度（l）上传导性的矩阵。该方程的解如下：

这种求解方法对任意时间步长（τ）都有效。上式中矩阵A_l是三对角矩阵并且对角占优（tridiagonal and diagonally dominant matrix），这样的线性系统可以通过Thomas算法快速求解。（有关AOS的应用，参见[3]）

该算法的实现源码如下（在文件 kaze.cpp 中）：

[cpp] view plain copy

1. //*************************************************************************************
2. //*************************************************************************************
3. /**
4. * @brief This method performs a scalar non-linear diffusion step using AOS schemes
5. * @param Ld Image at a given evolution step
6. * @param Ldprev Image at a previous evolution step
7. * @param c Conductivity image
8. * @param stepsize Stepsize for the nonlinear diffusion evolution
9. * @note If c is constant, the diffusion will be linear
10. * If c is a matrix of the same size as Ld, the diffusion will be nonlinear
11. * The stepsize can be arbitrarilly large
12. */
13. void KAZE::AOS_Step_Scalar(cv::Mat &Ld, const cv::Mat &Ldprev, const cv::Mat &c, const float stepsize)
14. {
15. AOS_Rows(Ldprev,c,stepsize);
16. AOS_Columns(Ldprev,c,stepsize);
17. Ld = 0.5*(Lty + Ltx.t());
18. }
19. //*************************************************************************************
20. //*************************************************************************************
21. /**
22. * @brief This method performs a scalar non-linear diffusion step using AOS schemes
23. * Diffusion in each dimension is computed independently in a different thread
24. * @param Ld Image at a given evolution step
25. * @param Ldprev Image at a previous evolution step
26. * @param c Conductivity image
27. * @param stepsize Stepsize for the nonlinear diffusion evolution
28. * @note If c is constant, the diffusion will be linear
29. * If c is a matrix of the same size as Ld, the diffusion will be nonlinear
30. * The stepsize can be arbitrarilly large
31. */
32. #if HAVE_THREADING_SUPPORT
33. void KAZE::AOS_Step_Scalar_Parallel(cv::Mat &Ld, const cv::Mat &Ldprev, const cv::Mat &c, const float stepsize)
34. {
35. boost::thread *AOSth1 = new boost::thread(&KAZE::AOS_Rows,this,Ldprev,c,stepsize);
36. boost::thread *AOSth2 = new boost::thread(&KAZE::AOS_Columns,this,Ldprev,c,stepsize);
37. AOSth1->join();
38. AOSth2->join();
39. Ld = 0.5*(Lty + Ltx.t());
40. delete AOSth1;
41. delete AOSth2;
42. }
43. #endif
44. //*************************************************************************************
45. //*************************************************************************************
46. /**
47. * @brief This method performs performs 1D-AOS for the image rows
48. * @param Ldprev Image at a previous evolution step
49. * @param c Conductivity image
50. * @param stepsize Stepsize for the nonlinear diffusion evolution
51. */
52. void KAZE::AOS_Rows(const cv::Mat &Ldprev, const cv::Mat &c, const float stepsize)
53. {
54. // Operate on rows
55. for( int i = 0; i < qr.rows; i++ )
56. {
57. for( int j = 0; j < qr.cols; j++ )
58. {
59. *(qr.ptr<float>(i)+j) = *(c.ptr<float>(i)+j) + *(c.ptr<float>(i+1)+j);
60. }
61. }
62. for( int j = 0; j < py.cols; j++ )
63. {
64. *(py.ptr<float>(0)+j) = *(qr.ptr<float>(0)+j);
65. }
66. for( int j = 0; j < py.cols; j++ )
67. {
68. *(py.ptr<float>(py.rows-1)+j) = *(qr.ptr<float>(qr.rows-1)+j);
69. }
70. for( int i = 1; i < py.rows-1; i++ )
71. {
72. for( int j = 0; j < py.cols; j++ )
73. {
74. *(py.ptr<float>(i)+j) = *(qr.ptr<float>(i-1)+j) + *(qr.ptr<float>(i)+j);
75. }
76. }
77. // a = 1 + t.*p; (p is -1*p)
78. // b = -t.*q;
79. ay = 1.0 + stepsize*py; // p is -1*p
80. by = -stepsize*qr;
81. // Call to Thomas algorithm now
82. Thomas(ay,by,Ldprev,Lty);
83. }
84. //*************************************************************************************
85. //*************************************************************************************
86. /**
87. * @brief This method performs performs 1D-AOS for the image columns
88. * @param Ldprev Image at a previous evolution step
89. * @param c Conductivity image
90. * @param stepsize Stepsize for the nonlinear diffusion evolution
91. */
92. void KAZE::AOS_Columns(const cv::Mat &Ldprev, const cv::Mat &c, const float stepsize)
93. {
94. // Operate on columns
95. for( int j = 0; j < qc.cols; j++ )
96. {
97. for( int i = 0; i < qc.rows; i++ )
98. {
99. *(qc.ptr<float>(i)+j) = *(c.ptr<float>(i)+j) + *(c.ptr<float>(i)+j+1);
100. }
101. }
102. for( int i = 0; i < px.rows; i++ )
103. {
104. *(px.ptr<float>(i)) = *(qc.ptr<float>(i));
105. }
106. for( int i = 0; i < px.rows; i++ )
107. {
108. *(px.ptr<float>(i)+px.cols-1) = *(qc.ptr<float>(i)+qc.cols-1);
109. }
110. for( int j = 1; j < px.cols-1; j++ )
111. {
112. for( int i = 0; i < px.rows; i++ )
113. {
114. *(px.ptr<float>(i)+j) = *(qc.ptr<float>(i)+j-1) + *(qc.ptr<float>(i)+j);
115. }
116. }
117. // a = 1 + t.*p‘;
118. ax = 1.0 + stepsize*px.t();
119. // b = -t.*q‘;
120. bx = -stepsize*qc.t();
121. // Call Thomas algorithm again
122. // But take care since we need to transpose the solution!!
123. Thomas(ax,bx,Ldprev.t(),Ltx);
124. }
125. //*************************************************************************************
126. //*************************************************************************************
127. /**
128. * @brief This method does the Thomas algorithm for solving a tridiagonal linear system
129. * @note The matrix A must be strictly diagonally dominant for a stable solution
130. */
131. void KAZE::Thomas(cv::Mat a, cv::Mat b, cv::Mat Ld, cv::Mat x)
132. {
133. // Auxiliary variables
134. int n = a.rows;
135. cv::Mat m = cv::Mat::zeros(a.rows,a.cols,CV_32F);
136. cv::Mat l = cv::Mat::zeros(b.rows,b.cols,CV_32F);
137. cv::Mat y = cv::Mat::zeros(Ld.rows,Ld.cols,CV_32F);
138. /** A*x = d;                                                                            */
139. /**  / a1 b1  0  0 0  ...    0 \  / x1 \ = / d1 \                                           */
140. /**  | c1 a2 b2  0 0  ...    0 |  | x2 | = | d2 |                                           */
141. /**  |  0 c2 a3 b3 0  ...    0 |  | x3 | = | d3 |                                           */
142. /**  |  :  :  :  : 0  ...    0 |  |  : | = |  : |                                           */
143. /**  |  :  :  :  : 0  cn-1  an |  | xn | = | dn |                                           */
144. /** 1. LU decomposition
145. / L = / 1                 \      U = / m1 r1            \
146. /     | l1 1              |          |    m2 r2         |
147. /     |    l2 1          |           |       m3 r3      |
148. /      |     : : :        |          |       :  :  :    |
149. /      \           ln-1 1 /          \               mn /    */
150. for( int j = 0; j < m.cols; j++ )
151. {
152. *(m.ptr<float>(0)+j) = *(a.ptr<float>(0)+j);
153. }
154. for( int j = 0; j < y.cols; j++ )
155. {
156. *(y.ptr<float>(0)+j) = *(Ld.ptr<float>(0)+j);
157. }
158. // 2. Forward substitution L*y = d for y
159. for( int k = 1; k < n; k++ )
160. {
161. for( int j=0; j < l.cols; j++ )
162. {
163. *(l.ptr<float>(k-1)+j) = *(b.ptr<float>(k-1)+j) / *(m.ptr<float>(k-1)+j);
164. }
165. for( int j=0; j < m.cols; j++ )
166. {
167. *(m.ptr<float>(k)+j) = *(a.ptr<float>(k)+j) - *(l.ptr<float>(k-1)+j)*(*(b.ptr<float>(k-1)+j));
168. }
169. for( int j=0; j < y.cols; j++ )
170. {
171. *(y.ptr<float>(k)+j) = *(Ld.ptr<float>(k)+j) - *(l.ptr<float>(k-1)+j)*(*(y.ptr<float>(k-1)+j));
172. }
173. }
174. // 3. Backward substitution U*x = y
175. for( int j=0; j < y.cols; j++ )
176. {
177. *(x.ptr<float>(n-1)+j) = (*(y.ptr<float>(n-1)+j))/(*(m.ptr<float>(n-1)+j));
178. }
179. for( int i = n-2; i >= 0; i-- )
180. {
181. for( int j = 0; j < x.cols; j++ )
182. {
183. *(x.ptr<float>(i)+j) = (*(y.ptr<float>(i)+j) - (*(b.ptr<float>(i)+j))*(*(x.ptr<float>(i+1)+j)))/(*(m.ptr<float>(i)+j));
184. }
185. }
186. }

上面介绍了非线性扩散滤波和AOS求解隐性差分方程的原理，是KAZE算法求解非线性尺度空间的基础，下一节我们将介绍KAZE算法的非线性尺度空间构建、特征检测与描述等内容。

待续...

Ref：

[1] http://www.robesafe.com/personal/pablo.alcantarilla/papers/Alcantarilla12eccv.pdf

[2] http://manu16.magtech.com.cn/geoprog/CN/article/downloadArticleFile.do?attachType=PDF&id=3146

[3] http://file.lw23.com/8/8e/8ec/8ecd21e4-b030-4e05-9333-40cc2d97bde4.pdf