SIFT算法实现理解及注释详解(基于RobHess源码)

1、Rob Hess的SIFT算法实现理解及注释 SIFT算法不用我多解释了，这是一个很强大的算法，主要用于图像配准 和物体识别等领域，但是其计算量相比也比较大，性价比比较高的算法包括 PCA-SIFT和SURF其中OpenCV提供了 SURF算法，但是为了方便理解。这 里给出了 Rob Hess所实现的SIFT算法的实现以及注释，结合我自己的理解， 如果，您有关于SIFT算法不理解的地方咱们可以一起交流一下。或者您认为不 详细的地方提出来。 SIFT算法的主要实现在sift.c这个文件，其主要流程为: (1)首先创建初始图像，即通过将图像转换为 32位的灰度图，然后将图像使用三 次插值来方大，之

2、后通过高斯模糊处理 (2) 在此基础上进行高斯金字塔的构建以及高斯差分金字塔的构建 (3) 对图像进行极值点检测 (4) 计算特征向量的尺度 调整图像大小 (6)计算特征的方向 (7)计算描述子，其中包括计算二维方向直方图并转换直方图为特征描述子 首先给出sift算法的整体框架代码: 输入参数: img 为输入图像; feat 为所要提取的特征指针; intvl 指的是高斯金字塔和差分金字塔的层数; sigma指的是图像初始化过程中高斯模糊所使用的参数; contr_thr是归一化之后的去除不稳定特征的阈值; curv_thr指的是去除边缘的特征的主曲率阈值; img_dbl是是否将图像放大为

3、之前的两倍; descr_with用来计算特征描述子的方向直方图的宽度; descr_hist_bi ns 是直方图中的条数 初始化图像 输入参数: 这里不需要解释了 该函数主要用来初始化图像，转换图像为 32位灰度图以及进行高斯模糊。 (2)构建高斯金字塔 输入参数: octvs是高斯金字塔的组 invls是高斯金字塔的层数 sigma是初始的高斯模糊参数，后续也通过它计算每一层所使用的sigma 7. 9. for ( i = 0; i ( calloc( intvls + 3, sizeof ( Ip lima ge* ); 11. 12. /* 13. p reco mpute Gau

4、ssian sigmas using the following formula: 14. 预计算每次高斯模糊的sigma 15. 16. sigma_totalA2 = sigma_i2 + sigma_i-12 17. */ 18. sig0 = sigma; 19. k = pow( 2.0, 1.0 / intvls ); 20. for ( i = 1; i intvls + 3; i+ ) 21. 22. sig_ prev = pow( k, i - 1 ) * sigma; 23. sig_total = sig_ prev * k; 24. sigi = sqrt( sig_

5、total * sig_total - sig_ prev * sig_ prev ); 25. 26. 27. 28. for ( o = 0; o octvs; o+ ) 29. for ( i = 0; i intvls + 3; i+ ) 30. 31. /对每一层进行降采样，形成高斯金字塔的每一层 32. if ( o = 0 34. 35. /* base of new octvave is halved image from end of p revious octave */ 36. /每一组的第一层都是通过对前面一组的最上面一层的降采样实现的 37. else if ( i

6、= 0 ) 38. gauss_ pyroi = downsa mple( gauss_ pyro-1intvls); 39. 40. /* blur the current octaves last image to create the next one */ 41. /每一组的其他层则使通过使用不同sigma的高斯模糊来进行处理 42. else 43. 44. gauss_ pyroi = cvCreateImage( cvGetSize(gauss_ pyroi-1). 45. IP L_DE PTH_32F, 1 ); 46. cvSmooth( gauss_ pyroi-1, g

7、auss_ pyroi, 47. CV_GAUSSIAN, 0, 0, sigi, sigi); 48. 降采样处理 输入参数: 不解释 这就是降采样，其实就是将图像通过最近邻算法缩小为原来的一半 (3)构建高斯差分金字塔 输入参数: 不解释了参见上面的说明即可 实际上差分金字塔的构成是通过对相邻层的图像进行相减获得的 (4)极值点检测 输入参数: contr thr是去除对比度低的点所采用的阈值 curv thr是去除边缘特征的阈值 SIFT IMG BORDE是预定义的图像边缘; 通过和对比度阈值比较去掉低对比度的点; 而通过is_extremum来判断是否为极值点，如果是则通过极值点插值

8、的方式获取 亚像素的极值点的位置。 然后通过is_too_eage_like和所给的主曲率阈值判断是否为边缘点 *判断是否为极值点 其原理为：通过和高斯金字塔的上一层的9个像素+本层的除了本像素自 己的其他的8个像素和下一层的9个像素进行比较看是否为这 26个像 素中最小的一个或者是否为最大的一个，如果是则为极值点。 *获取亚像素的极值点的位置 5. struct detection_data* ddata; 6. double xi, xr, xc, contr;/ 分另U为亚像素的 intval,row,col 的偏移offset ，和对比度 7. 9. while ( i SIFT_MA

9、X_INTERP_STEPS )/重新确定极值点并重新定位的操作只能循环5次 10. 11. interp_ste p( dog_ pyr, octv, intvl, r, c, 12. if ( ABS( xi ) 0.5 14. /否则继续寻找极值点 15. c += cvRound( xc ); 16. r += cvRound( xr ); 17. intvl += cvRound( xi ); 18. 19. if ( intvl intvls | 21. c SIFT_IMG_BORDER | 22. r = dog_ pyroctv0-width - SIFT_IMG_BORDE

10、R | 24. r = dog_ pyroctvO-height - SIFT_IMG_BORDER ) 25. 26. return NULL; 27. 28. 29. i+; 30. 31. 32. /确保极值点是经过最大5步找到的 33. /* ensure convergence of inter po lation */ 34. if ( i = SIFT_MAX_INTERP_STEPS ) 35. return NULL; 36. 37. /获取找到的极值点的对比度 38. contr = interp_contr( dog_ pyr, octv, intvl, r, c, xi

11、, xr, xc ); 39. /判断极值点是否小于某一个阈值 40. if ( ABS( contr ) img _p t.x = feat-x = ( c + xc ) * pow( 2.0, octv ); 46. feat-img_ pt.y = feat-y = ( r + xr ) * pow( 2.0, octv ); 47. ddata-r = r; *获取亚像素位置中所用到的函数 *计算三维偏导数 计算在x和y方向上的偏导数，高斯差分尺度空间金字塔中像素的尺度 实际上在离散数据中计算偏导数是通过相邻像素的相减来计算的 比如说计算x方向的偏导数dx，则通过该向所的x方向的后一个

12、减去 前一个然后除以2即可求的dx *计算三维海森矩阵 不需要讲什么，其实就是计算二次导数，计算方法也和一次导数的计算 如出一辙。 然后将结果放入到一个矩阵中去。 *计算插入像素的对比度 其中cvGEMM 是矩阵的通用计算函数，至于CV_GEMM_A_T 是计算dD的转置矩阵放入 *去除边缘相应 通过计算所在特征向量的主曲率半径来判断特征是边缘的从而导致不 稳定 即去除边缘响应 (4)计算特征向量的尺度 实际上是通过最初的sigma来获得每一层每一组的尺度 (5)调整图像特征坐标、尺度、点的坐标大小为原来的一半 (6)给每一个图像特征向量计算规范化的方向 *对所给像素计算灰度方向直方图 以关键

13、点为中心的邻域窗口内采样，并用直方图统计邻域像素的梯度 方向。梯度直方图的范围是 0360度，其中每10度一个柱，总共36 个柱 *计算所给像素的梯度大小和方向 每一个小格都代表了特征点邻域所在的尺度空间的一个像素 ，箭头方 向代表了像素梯 度方向，箭头长度代表该像素的幅值也就是梯度的值 *对方向直方图进行高斯模糊 使用高斯函数对直方图进行平滑，减少突变的影响。 *在直方图中找到主方向的梯度 利用关键点邻域像素的梯度方向分布特性为每个关键点指定方向参数, 使算子具备 旋转不变性。 *将大于某一个梯度大小阈值的特征向量加入到直方图中去 n为方向的个数 (7)计算特征描述子 *计算二维方向直方图

14、7. 8 int radius, i, j; 8. 9. hist =static_castvdouble *( calloc( d,sizeof ( double * ); 10. for ( i = 0; i v d; i+ ) 11. 12. histi =static_cast vdouble *( calloc( d,sizeof ( double * ); 13. for ( j = 0; j v d; j+ ) 14. histij =static_castvdouble *( calloc( n,sizeof ( double ); 15. 16. 17. cos_t = co

15、s( ori ); 18. sin_t = sin( ori ); 19. bins_ per_rad = n / P12; 20. exp_denom = d * d * 0.5; 21. hist_width = SIFT_DESCR_SCL_FCTR * scl; 22. radius = hist_width * sqrt(2.0) * ( d + 1.0 ) * 0.5 + 0.5; 23. for ( i = -radius; i v= radius; i+ ) 24. for ( j = -radius; j v= radius; j+ ) 25. 26. /* 27. 即将坐标

16、移至关键点主方向 28. 计算采用的直方图数组中相对于方向旋转的坐标 29. Calculate sampi es histogram array coords rotated relative to ori. 30. Subtract 0.5 so samples that fall e.g. in the center of row 1 (i.e. 31. r_rot = 1.5) have full weight pl aced in row 1 after inter po lation. 32. */ 33. c_rot = ( j * cos_t - i * sin_t ) / h

17、ist_width; 34. r_rot = ( j * sin_t + i * cos_t ) / hist_width; 35. rbin = r_rot + d / 2 - 0.5; 36. cbin = c_rot + d / 2 - 0.5; 37. 38. if ( rbin -1.0 42. while ( grad_ori v 0.0 ) 43. grad_ori += PI2; 44. while ( grad_ori = PI2 ) 45. grad_ori -= P12; 46. 47. obin = grad_ori * bins_ per_rad; 48. w = exp( -(c_rot * c_rot + r_rot * r_rot) / exp_denom ); 49. interp_hist_entry( hist, rbin, c