SIFT原理

转载请注明出处：http://blog.csdn.NET/luoshixian099/article/details/47377611

相关： KD树+BBF算法解析

           SURF原理与源码解析

     SIFT的原理已经有很多大牛的博客上做了解析，本文重点将以Rob Hess等人用C实现的代码做解析，结合代码SIFT原理会更容易理解。一些难理解点的用了☆标注。

      欢迎大家批评指正！

     SIFT（Scale-invariant feature transform）即尺度不变特征转换，提取的局部特征点具有尺度不变性，且对于旋转，亮度，噪声等有很高的稳定性。

下图中，涉及到图像的旋转，仿射，光照等变化，SIFT算法依然有很好的匹配效果。

SIFT特征点提取

本文将以下函数为参照顺序介绍SIFT特征点提取与描述方法。

 1.图像预处理

 2.构建高斯金字塔（不同尺度下的图像）

 3.生成DOG尺度空间

 4.关键点搜索与定位

 5.计算特征点所在的尺度

 6.为特征点分配方向角

 7.构建特征描述子

[cpp]  view plain  copy

 print ?

1. /** 
2.    Finds SIFT features in an image using user-specified parameter values.  All 
3.    detected features are stored in the array pointed to by \a feat. 
4. */  
5. int _sift_features( IplImage* img, struct feature** feat, int intvls,  
6.             double sigma, double contr_thr, int curv_thr,  
7.             int img_dbl, int descr_width, int descr_hist_bins )  
8. {  
9.   IplImage* init_img;  
10.   IplImage*** gauss_pyr, *** dog_pyr;  
11.   CvMemStorage* storage;  
12.   CvSeq* features;  
13.   int octvs, i, n = 0;  
14.     
15.   /* check arguments */  
16.   if( ! img )  
17.     fatal_error( "NULL pointer error, %s, line %d",  __FILE__, __LINE__ );  
18.   if( ! feat )  
19.     fatal_error( "NULL pointer error, %s, line %d",  __FILE__, __LINE__ );  
20.   
21.   /* build scale space pyramid; smallest dimension of top level is ~4 pixels */  
22.   init_img = create_init_img( img, img_dbl, sigma );                            //对进行图片预处理         
23.   octvs = log( MIN( init_img->width, init_img->height ) ) / log(2) - 2;  //计算高斯金字塔的组数(octave),同时保证顶层至少有4个像素点  
24.   gauss_pyr = build_gauss_pyr( init_img, octvs, intvls, sigma );  //建立高斯金字塔  
25.   dog_pyr = build_dog_pyr( gauss_pyr, octvs, intvls );   //DOG尺度空间  
26.     
27.   storage = cvCreateMemStorage( 0 );  
28.   features = scale_space_extrema( dog_pyr, octvs, intvls, contr_thr,   //极值点检测，并去除不稳定特征点  
29.                   curv_thr, storage );  
30.   calc_feature_scales( features, sigma, intvls );                      //计算特征点所在的尺度  
31.   if( img_dbl )  
32.     adjust_for_img_dbl( features );                                       //如果图像初始被扩大了2倍，所有坐标与尺度要除以2  
33.   calc_feature_oris( features, gauss_pyr );                               //计算特征点所在尺度内的方向角  
34.   compute_descriptors( features, gauss_pyr, descr_width, descr_hist_bins );//计算特征描述子 128维向量  
35.   
36.   /* sort features by decreasing scale and move from CvSeq to array */  
37.   cvSeqSort( features, (CvCmpFunc)feature_cmp, NULL );   //对特征点按尺度排序  
38.   n = features->total;  
39.   *feat = calloc( n, sizeof(struct feature) );  
40.   *feat = cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ );         
41.   for( i = 0; i < n; i++ )  
42.     {  
43.       free( (*feat)[i].feature_data );  
44.       (*feat)[i].feature_data = NULL;  
45.     }  
46.     
47.   cvReleaseMemStorage( &storage );  
48.   cvReleaseImage( &init_img );  
49.   release_pyr( &gauss_pyr, octvs, intvls + 3 );  
50.   release_pyr( &dog_pyr, octvs, intvls + 2 );  
51.   return n;  
52. }  

  —————————————————————————————————————————————————————

 1.图像预处理

[cpp]  view plain  copy

 print ?

1. /************************ Functions prototyped here **************************/  
2.   
3. /* 
4.   Converts an image to 8-bit grayscale and Gaussian-smooths it.  The image is 
5.   optionally doubled in size prior to smoothing. 
6.  
7.   @param img input image 
8.   @param img_dbl if true, image is doubled in size prior to smoothing 
9.   @param sigma total std of Gaussian smoothing 
10. */  
11. static IplImage* create_init_img( IplImage* img, int img_dbl, double sigma )  
12. {  
13.   IplImage* gray, * dbl;  
14.   double sig_diff;  
15.   
16.   gray = convert_to_gray32( img );   //转换为32位灰度图  
17.   if( img_dbl )                                  // 图像被放大二倍  
18.     {  
19.       sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4 );   //  sigma = 1.6 , SIFT_INIT_SIGMA = 0.5  lowe认为图像在尺度0.5下最清晰  
20.       dbl = cvCreateImage( cvSize( img->width*2, img->height*2 ),  
21.                IPL_DEPTH_32F, 1 );  
22.       cvResize( gray, dbl, CV_INTER_CUBIC );  //双三次插值方法 放大图像  
23.       cvSmooth( dbl, dbl, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff );     //高斯平滑  
24.       cvReleaseImage( &gray );  
25.       return dbl;  
26.     }  
27.   else  
28.     {  
29.       sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA );  
30.       cvSmooth( gray, gray, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff ); // 高斯平滑  
31.       return gray;  
32.     }  
33. }  

lowe建议把初始图像放大二倍，可以得到更多的特征点，提取到更多细节，并且认为图像在尺度σ = 0.5时图像最清晰，初始高斯尺度为σ = 1.6。

☆第19行因为图像被放大二倍，此时σ = 1.0 。因为对二倍化后的图像平滑是在σ = 0.5 上叠加的高斯模糊，

  所以有模糊系数有sig_diff = sqrt (sigma *sigma - 0.5*0.5*4)=sqrt(1.6*1.6 -1) ；

 2.构建高斯金字塔

构建高斯金字塔过程即构建出图像在不同尺度上图像，提取到的特征点可有具有尺度不变性。

图像的尺度空间L(x,y,σ)可以用一个高斯函数G(x,y,σ)与图像I(x,y)卷积产生，即L(x,y,σ) = G(x,y,σ) * I(x,y) 

其中二维高斯核的计算为             。

☆不同的尺度空间即用不同的高斯核函数平滑图像， 平滑系数越大，图像越模糊。即模拟出动物的视觉效果，因为事先不知道物体的大小，在不同的尺度下，图像的细节会表现的不同。当尺度由小变大的过程中，是一个细节逐步简化的过程，图像中特征不够明显的物体，就模糊的多了，而有些物体还可以看得到大致的轮廓。所以要在不同尺度下，观察物体的尺度响应，提取到的特征才能具有尺度不变性。

SIFT算法采用高斯金字塔实现连续的尺度空间的图像。金字塔共分为O(octave)组，每组有S(intervals)层 ，下一组是由上一组隔点采样得到（即降2倍分辨率），这是为了减轻卷积运算的工作量。

构建高斯金字塔（octave = 5, intervals +3=6）：

                        

全部空间尺度为 ： 

                                         

☆1.这个尺度因子都是在原图上进行的，而在算法实现过程中，采用高斯平滑是在上一层图像上再叠加高斯平滑，即我们在程序中看到的平滑因子为   

            

Eg. 在第一层上为了得到kσ的高斯模糊图像，可以在原图上直接采用kσ平滑，也可以在上一层图像上（已被σ平滑）的图像上采用平滑因子为平滑图像，效果是一样的。

 ☆2.我们在源码上同时也没有看到组间的2倍的关系，实际在对每组的平滑因子都是一样的，2倍的关系是由于在降采样的过程中产生的，第二层的第一张图是由第一层的平滑因子为2σ的图像上（即倒数第三张）降采样得到，此时图像平滑因子为σ，所以继续采用以上的平滑因子。但在原图上看，形成了全部的空间尺度。

☆3. 每组（octave）有S+3层图像，是由于在DOG尺度空间上寻找极值点的方法是在一个立方体内进行，即上下层比较，所以不在DOG空间的第一层与最后一层寻找，即DOG需要S+2层图像，由于DOG尺度空间是由高斯金字塔相邻图像相减得到，即每组需要S+3层图像。

[cpp]  view plain  copy

 print ?

1. /* 
2.   Builds Gaussian scale space pyramid from an image 
3.   @param base base image of the pyramid 
4.   @param octvs number of octaves of scale space 
5.   @param intvls number of intervals per octave 
6.   @param sigma amount of Gaussian smoothing per octave 
7.  
8.   @return Returns a Gaussian scale space pyramid as an octvs x (intvls + 3) 
9.     array  
10.      
11. 给定组数(octave)和层数(intvls)，以及初始平滑系数sigma,构建高斯金字塔 
12. 返回的每组中层数为intvls+3 
13. */  
14. static IplImage*** build_gauss_pyr( IplImage* base, int octvs,  
15.                  int intvls, double sigma )  
16. {  
17.   IplImage*** gauss_pyr;  
18.   const int _intvls = intvls;             // lowe 采用了每组层数(intvls)为 3  
19.  // double  sig_total, sig_prev;  
20.      double  k;  
21.   int i, o;  
22.   double *sig = (double *)malloc(sizeof(int)*(_intvls+3));  //存储每组的高斯平滑因子，每组对应的平滑因子都相同  
23.   
24.   gauss_pyr = calloc( octvs, sizeof( IplImage** ) );                
25.   for( i = 0; i < octvs; i++ )  
26.     gauss_pyr[i] = calloc( intvls + 3, sizeof( IplImage *) );  
27.   
28.   /* 
29.     precompute Gaussian sigmas using the following formula: 
30.  
31.     \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2 
32.  
33.     sig[i] is the incremental sigma value needed to compute  
34.     the actual sigma of level i. Keeping track of incremental 
35.     sigmas vs. total sigmas keeps the gaussian kernel small. 
36.   */  
37.   k = pow( 2.0, 1.0 / intvls );                 // k = 2^(1/S)  
38.   sig[0] = sigma;  
39.   sig[1] = sigma * sqrt( k*k- 1 );  
40.   for (i = 2; i < intvls + 3; i++)  
41.       sig[i] = sig[i-1] * k;                       //每组对应的平滑因子为 σ ,  sqrt(k^2 -1)* σ, sqrt(k^2 -1)* kσ , ...  
42.   
43.   for( o = 0; o < octvs; o++ )  
44.     for( i = 0; i < intvls + 3; i++ )  
45.       {  
46.     if( o == 0  &&  i == 0 )  
47.       gauss_pyr[o][i] = cvCloneImage(base);                       //第一组，第一层为原图  
48.   
49.     /* base of new octvave is halved image from end of previous octave */  
50.     else if( i == 0 )  
51.       gauss_pyr[o][i] = downsample( gauss_pyr[o-1][intvls] );  //第一层图像由上一层倒数第三张隔点采样得到  
52.         
53.     /* blur the current octave's last image to create the next one */  
54.     else  
55.       {  
56.         gauss_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i-1]),  
57.                          IPL_DEPTH_32F, 1 );  
58.         cvSmooth( gauss_pyr[o][i-1], gauss_pyr[o][i],  
59.               CV_GAUSSIAN, 0, 0, sig[i], sig[i] );                       //高斯平滑  
60.       }  
61.       }  
62.        
63.      
64.   return gauss_pyr;  
65. }  

3.生成DOG尺度空间

Lindeberg发现高斯差分函数（Difference of Gaussian ，简称DOG算子）与尺度归一化的高斯拉普拉斯函数非常近似，且

 

差分近似：

lowe建议采用相邻尺度的图像相减来获得高斯差分图像D(x,y,σ)来近似LOG来进行极值检测。

D(x,y,σ) = G(x,y,kσ)*I(x,y)-G(x,y,σ)*I(x,y)

              =L(x,y,kσ) - L(x,y,σ)

对高斯金字塔的每组内相邻图像相减，形成DOG尺度空间，这时DOG中每组有S+2层图像

[cpp]  view plain  copy

 print ?

1. static IplImage*** build_dog_pyr( IplImage*** gauss_pyr, int octvs, int intvls )  
2. {  
3.   IplImage*** dog_pyr;  
4.   int i, o;  
5.   
6.   dog_pyr = calloc( octvs, sizeof( IplImage** ) );  
7.   for( i = 0; i < octvs; i++ )  
8.     dog_pyr[i] = calloc( intvls + 2, sizeof(IplImage*) );  
9.   
10.   for( o = 0; o < octvs; o++ )  
11.     for( i = 0; i < intvls + 2; i++ )  
12.       {  
13.     dog_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i]),  
14.                        IPL_DEPTH_32F, 1 );  
15.     cvSub( gauss_pyr[o][i+1], gauss_pyr[o][i], dog_pyr[o][i], NULL );   //相邻两层图像相减，结果放在dog_pyr数组内  
16.       }  
17.   
18.   return dog_pyr;  
19. }  

 4.关键点搜索与定位

   在DOG尺度空间上，首先寻找极值点，插值处理，找到准确的极值点坐标，再排除不稳定的特征点（边界点）

[cpp]  view plain  copy

 print ?

1. /* 
2.   Detects features at extrema in DoG scale space.  Bad features are discarded 
3.   based on contrast and ratio of principal curvatures. 
4.  
5.   @return Returns an array of detected features whose scales, orientations, 
6.     and descriptors are yet to be determined. 
7. */  
8. static CvSeq* scale_space_extrema( IplImage*** dog_pyr, int octvs, int intvls,  
9.                    double contr_thr, int curv_thr,  
10.                    CvMemStorage* storage )  
11. {  
12.   CvSeq* features;  
13.   double prelim_contr_thr = 0.5 * contr_thr / intvls; //极值比较前的阈值处理  
14.   struct feature* feat;  
15.   struct detection_data* ddata;  
16.   int o, i, r, c;  
17.   
18.   features = cvCreateSeq( 0, sizeof(CvSeq), sizeof(struct feature), storage );  
19.   for( o = 0; o < octvs; o++ )                     //对DOG尺度空间上，遍历从第二层图像开始到倒数第二层图像上，每个像素点  
20.     for( i = 1; i <= intvls; i++ )  
21.       for(r = SIFT_IMG_BORDER; r < dog_pyr[o][0]->height-SIFT_IMG_BORDER; r++)          
22.     for(c = SIFT_IMG_BORDER; c < dog_pyr[o][0]->width-SIFT_IMG_BORDER; c++)  
23.       /* perform preliminary check on contrast */  
24.       if( ABS( pixval32f( dog_pyr[o][i], r, c ) ) > prelim_contr_thr )    // 排除像素值小于阈值prelim_contr_thr的点，提高稳定性  
25.         if( is_extremum( dog_pyr, o, i, r, c ) )             //与周围26个像素值比较，是否极大值或者极小值点  
26.           {  
27.         feat = interp_extremum(dog_pyr, o, i, r, c, intvls, contr_thr); //插值处理，找到准确的特征点坐标  
28.         if( feat )  
29.           {  
30.             ddata = feat_detection_data( feat );  
31.             if( ! is_too_edge_like( dog_pyr[ddata->octv][ddata->intvl],   //根据Hessian矩阵 判断是否为边缘上的点  
32.                         ddata->r, ddata->c, curv_thr ) )  
33.               {  
34.             cvSeqPush( features, feat );          //是特征点进入特征点序列  
35.               }  
36.             else  
37.               free( ddata );  
38.             free( feat );  
39.           }  
40.           }  
41.     
42.   return features;  
43. }  

4.1

寻找极值点

在DOG尺度空间上，每组有S+2层图像，每一组都从第二层开始每一个像素点都要与它相邻的像素点比较，看是否比它在图像域或尺度域的所有点的值大或者小。与它同尺度的相邻像素点有8个，上下相邻尺度的点共有2×9=18，共有26个像素点。也就在一个3×3的立方体内进行。搜索的过程是第二层开始到倒数第二层结束，共检测了octave组，每组S层。

  

[cpp]  view plain  copy

 print ?

1. /* 
2.   Determines whether a pixel is a scale-space extremum by comparing it to it's 
3.   3x3x3 pixel neighborhood. 
4. */  
5. static int is_extremum( IplImage*** dog_pyr, int octv, int intvl, int r, int c )  
6. {  
7.   double val = pixval32f( dog_pyr[octv][intvl], r, c );  
8.   int i, j, k;  
9.   
10.   /* check for maximum */  
11.   if( val > 0 )  
12.     {  
13.       for( i = -1; i <= 1; i++ )  
14.     for( j = -1; j <= 1; j++ )  
15.       for( k = -1; k <= 1; k++ )  
16.         if( val < pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )  
17.           return 0;  
18.     }  
19.   
20.   /* check for minimum */  
21.   else  
22.     {  
23.       for( i = -1; i <= 1; i++ )  
24.     for( j = -1; j <= 1; j++ )  
25.       for( k = -1; k <= 1; k++ )  
26.         if( val > pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )  
27.           return 0;  
28.     }  
29.   
30.   return 1;  
31. }  

4.2

准确定位特征点

      以上的极值点搜索是在离散空间进行的，极值点不真正意义上的极值点。通过对空间尺度函数拟合，可以得到亚像素级像素点坐标。

尺度空间的Taylor展开式：

                                       ，其中

求导并令其为0，得到亚像素级：

                                           

对应的函数值为：

                                        

 

是一个三维矢量，矢量在任何一个方向上的偏移量大于0.5时，意味着已经偏离了原像素点，这样的特征坐标位置需要更新或者继续插值计算。算法实现过程中，为了保证插值能够收敛，设置了最大插值次数（lowe 设置了5次）。同时当 时（本文阈值采用了0.04/S） ，特征点才被保留，因为响应值过小的点，容易受噪声的干扰而不稳定。

对离散空间进行函数拟合（插值）：

[cpp]  view plain  copy

 print ?

1. /* 
2.   Performs one step of extremum interpolation.  Based on Eqn. (3) in Lowe's 
3.   paper. 
4.  
5.   r,c 为特征点位置，xi,xr,xc,保存三个方向的偏移量 
6. */  
7.   
8. static void interp_step( IplImage*** dog_pyr, int octv, int intvl, int r, int c,  
9.              double* xi, double* xr, double* xc )  
10. {  
11.   CvMat* dD, * H, * H_inv, X;  
12.   double x[3] = { 0 };  
13.     
14.   dD = deriv_3D( dog_pyr, octv, intvl, r, c );      //计算三个方向的梯度  
15.   H = hessian_3D( dog_pyr, octv, intvl, r, c );    // 计算3维空间的hessian矩阵  
16.   H_inv = cvCreateMat( 3, 3, CV_64FC1 );  
17.   cvInvert( H, H_inv, CV_SVD );           //计算逆矩阵  
18.   cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );  
19.   cvGEMM( H_inv, dD, -1, NULL, 0, &X, 0 );   //广义乘法   
20.     
21.   cvReleaseMat( &dD );  
22.   cvReleaseMat( &H );  
23.   cvReleaseMat( &H_inv );  
24.   
25.   *xi = x[2];  
26.   *xr = x[1];  
27.   *xc = x[0];  
28. }  

[cpp]  view plain  copy

 print ?

1. /* 
2.   Interpolates a scale-space extremum's location and scale to subpixel 
3.   accuracy to form an image feature.  
4. */  
5. static struct feature* interp_extremum( IplImage*** dog_pyr, int octv,          //通过拟合求取准确的特征点位置  
6.                     int intvl, int r, int c, int intvls,  
7.                     double contr_thr )  
8. {  
9.   struct feature* feat;  
10.   struct detection_data* ddata;  
11.   double xi, xr, xc, contr;  
12.   int i = 0;  
13.     
14.   while( i < SIFT_MAX_INTERP_STEPS )   //在最大迭代次数范围内进行  
15.     {  
16.       interp_step( dog_pyr, octv, intvl, r, c, &xi, &xr, &xc );          //插值后得到的三个方向的偏移量(xi,xr,xc)  
17.       if( ABS( xi ) < 0.5  &&  ABS( xr ) < 0.5  &&  ABS( xc ) < 0.5 )  
18.     break;  
19.         
20.       c += cvRound( xc );    //更新位置  
21.       r += cvRound( xr );  
22.       intvl += cvRound( xi );  
23.         
24.       if( intvl < 1  ||                            
25.       intvl > intvls  ||  
26.       c < SIFT_IMG_BORDER  ||  
27.       r < SIFT_IMG_BORDER  ||  
28.       c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER  ||  
29.       r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER )  
30.     {  
31.       return NULL;  
32.     }  
33.         
34.       i++;  
35.     }  
36.     
37.   /* ensure convergence of interpolation */  
38.   if( i >= SIFT_MAX_INTERP_STEPS )     
39.     return NULL;  
40.     
41.   contr = interp_contr( dog_pyr, octv, intvl, r, c, xi, xr, xc );     //计算插值后对应的函数值  
42.   if( ABS( contr ) < contr_thr / intvls )   //小于阈值（0.04/S）的点，则丢弃  
43.     return NULL;  
44.   
45.   feat = new_feature();  
46.   ddata = feat_detection_data( feat );  
47.   feat->img_pt.x = feat->x = ( c + xc ) * pow( 2.0, octv );       // 计算特征点根据降采样的次数对应于原图中位置  
48.   feat->img_pt.y = feat->y = ( r + xr ) * pow( 2.0, octv );  
49.   ddata->r = r;                  // 在本尺度内的坐标位置  
50.   ddata->c = c;  
51.   ddata->octv = octv;                 //组信息  
52.   ddata->intvl = intvl;                 // 层信息  
53.   ddata->subintvl = xi;              // 层方向的偏移量  
54.   
55.   return feat;  
56. }  

4.3

删除边缘效应

为了得到稳定的特征点，要删除掉落在图像边缘上的点。一个落在边缘上的点，可以根据主曲率计算判断。主曲率可以通过2维的 Hessian矩阵求出；

在边缘上的点，必定使得Hessian矩阵的两个特征值相差比较大，而特征值与矩阵元素有以下关系；

令α=rβ ， 所以有：

我们可以判断上述公式的比值大小，大于阈值（lowe采用 r =10）的点排除。

[cpp]  view plain  copy

 print ?

1. static int is_too_edge_like( IplImage* dog_img, int r, int c, int curv_thr )  
2. {  
3.   double d, dxx, dyy, dxy, tr, det;  
4.   
5.   /* principal curvatures are computed using the trace and det of Hessian */              
6.   d = pixval32f(dog_img, r, c);                                                                             //计算Hessian 矩阵内的4个元素值  
7.   dxx = pixval32f( dog_img, r, c+1 ) + pixval32f( dog_img, r, c-1 ) - 2 * d;  
8.   dyy = pixval32f( dog_img, r+1, c ) + pixval32f( dog_img, r-1, c ) - 2 * d;  
9.   dxy = ( pixval32f(dog_img, r+1, c+1) - pixval32f(dog_img, r+1, c-1) -  
10.       pixval32f(dog_img, r-1, c+1) + pixval32f(dog_img, r-1, c-1) ) / 4.0;  
11.   tr = dxx + dyy;                          //矩阵的迹  
12.   det = dxx * dyy - dxy * dxy;     //矩阵的值  
13.   
14.   /* negative determinant -> curvatures have different signs; reject feature */  
15.   if( det <= 0 )     // 矩阵值为负值，说明曲率有不同符号，丢弃  
16.     return 1;  
17.   
18.   if( tr * tr / det < ( curv_thr + 1.0 )*( curv_thr + 1.0 ) / curv_thr )   //比值小于阈值的特征点被保留  curv_thr = 10  
19.     return 0;  
20.   return 1;  
21. }  

 5.计算特征点对应的尺度

[cpp]  view plain  copy

 print ?

1. static void calc_feature_scales( CvSeq* features, double sigma, int intvls )  
2. {  
3.   struct feature* feat;  
4.   struct detection_data* ddata;  
5.   double intvl;  
6.   int i, n;  
7.   
8.   n = features->total;  
9.   for( i = 0; i < n; i++ )  
10.     {  
11.       feat = CV_GET_SEQ_ELEM( struct feature, features, i );  
12.       ddata = feat_detection_data( feat );  
13.       intvl = ddata->intvl + ddata->subintvl;                          
14.       feat->scl = sigma * pow( 2.0, ddata->octv + intvl / intvls );      // feat->scl 保存特征点在总体上尺度  
15.       ddata->scl_octv = sigma * pow( 2.0, intvl / intvls );     // feat->feature_data->scl__octv 保存特征点在组内的尺度，用来下面计算方向角  
16.     }  
17. }  

 6.为特征点分配方向角

这部分包括：计算邻域内梯度直方图，平滑直方图，复制特征点（有辅方向的特征点）

[cpp]  view plain  copy

 print ?

1. static void calc_feature_oris( CvSeq* features, IplImage*** gauss_pyr )    
2. {  
3.   struct feature* feat;  
4.   struct detection_data* ddata;  
5.   double* hist;  
6.   double omax;  
7.   int i, j, n = features->total;  
8.   
9.   for( i = 0; i < n; i++ )  
10.     {  
11.       feat = malloc( sizeof( struct feature ) );  
12.       cvSeqPopFront( features, feat );  
13.       ddata = feat_detection_data( feat );  
14.       hist = ori_hist( gauss_pyr[ddata->octv][ddata->intvl],     // 计算邻域内的梯度直方图,邻域半径radius = 3*1.5*sigma;  高斯加权系数= 1.5 *sigma   
15.                ddata->r, ddata->c, SIFT_ORI_HIST_BINS,  
16.                cvRound( SIFT_ORI_RADIUS * ddata->scl_octv ),  
17.                SIFT_ORI_SIG_FCTR * ddata->scl_octv );  
18.       for( j = 0; j < SIFT_ORI_SMOOTH_PASSES; j++ )  
19.     smooth_ori_hist( hist, SIFT_ORI_HIST_BINS );    // 对直方图平滑  
20.       omax = dominant_ori( hist, SIFT_ORI_HIST_BINS ); // 返回直方图的主方向  
21.       add_good_ori_features( features, hist, SIFT_ORI_HIST_BINS,//大于主方向的80%为辅方向  
22.                  omax * SIFT_ORI_PEAK_RATIO, feat );  
23.       free( ddata );  
24.       free( feat );  
25.       free( hist );  
26.     }  
27. }  

6.1

建立特征点邻域内的直方图

上一步scl_octv保存了每个特征点所在的组内尺度，下面计算特征点所在尺度内的高斯图像，以3×1.5×scl_octv为半径的区域内的所有像素点的梯度幅值与幅角；

计算公式：

在计算完所有特征点的幅值与幅角后，使用直方图统计。直方图横轴为梯度方向角，纵轴为对应幅值的累加值（与权重），梯度方向范围为0~360度，划分为36个bin,每个bin的宽度为10。

下图描述的划分为8个bin,每个bin的宽度为45的效果图：

其次，每个被加入直方图的幅值，要进行权重处理，权重也是采用高斯加权函数，其中高斯系数为1.5×scl_octv。通过高斯加权使特征点附近的点有较大的权重，可以弥补部分因没有仿射不变性而产生的不稳定问题；

即每个bin值按下面的公式累加，mag是幅值，后面为权重；i，j,为偏离特征点距离：

☆程序上可以帮助你理解上面的概念：

[cpp]  view plain  copy

 print ?

1. static double* ori_hist( IplImage* img, int r, int c, int n, int rad,  
2.              double sigma )  
3. {  
4.   double* hist;  
5.   double mag, ori, w, exp_denom, PI2 = CV_PI * 2.0;  
6.   int bin, i, j;  
7.   
8.   hist = calloc( n, sizeof( double ) );  
9.   exp_denom = 2.0 * sigma * sigma;  
10.   for( i = -rad; i <= rad; i++ )  
11.     for( j = -rad; j <= rad; j++ )  
12.       if( calc_grad_mag_ori( img, r + i, c + j, &mag, &ori ) )    //计算领域像素点的梯度幅值mag与方向ori  
13.     {  
14.       w = exp( -( i*i + j*j ) / exp_denom );                    //高斯权重  
15.       bin = cvRound( n * ( ori + CV_PI ) / PI2 );  
16.       bin = ( bin < n )? bin : 0;  
17.       hist[bin] += w * mag;                             //直方图上累加  
18.     }  
19.   
20.   return hist;  //返回累加完成的直方图  
21. }  

6.2

平滑直方图

lowe建议对直方图进行平滑，减少突变的影响。

[cpp]  view plain  copy

 print ?

1. static void smooth_ori_hist( double* hist, int n )  
2. {  
3.   double prev, tmp, h0 = hist[0];  
4.   int i;  
5.   
6.   prev = hist[n-1];  
7.   for( i = 0; i < n; i++ )  
8.     {  
9.       tmp = hist[i];  
10.       hist[i] = 0.25 * prev + 0.5 * hist[i] +   
11.     0.25 * ( ( i+1 == n )? h0 : hist[i+1] );  
12.       prev = tmp;  
13.     }  
14. }  

6.3

复制特征点

在上面的直方图上，我们已经找到了特征点主方向的峰值omax,当存在另一个大于主峰值80%的峰值时，将这个方向作为特征点的辅方向，即一个特征点有多个方向，这可以增强匹配的鲁棒性。在算法上，即把该特征点复制多份作为新的特征点，新特征点的方向为这些辅方向，其他属性保持一致。

[cpp]  view plain  copy

 print ?

1. static void add_good_ori_features( CvSeq* features, double* hist, int n,  
2.                    double mag_thr, struct feature* feat )  
3. {  
4.   struct feature* new_feat;  
5.   double bin, PI2 = CV_PI * 2.0;  
6.   int l, r, i;  
7.   
8.   for( i = 0; i < n; i++ )                   //检查所有的方向  
9.     {  
10.       l = ( i == 0 )? n - 1 : i-1;  
11.       r = ( i + 1 ) % n;  
12.         
13.       if( hist[i] > hist[l]  &&  hist[i] > hist[r]  &&  hist[i] >= mag_thr ) //判断是不是幅方向  
14.     {  
15.       bin = i + interp_hist_peak( hist[l], hist[i], hist[r] );     //插值离散处理，取得更精确的方向  
16.       bin = ( bin < 0 )? n + bin : ( bin >= n )? bin - n : bin;  
17.       new_feat = clone_feature( feat );       //复制特征点  
18.       new_feat->ori = ( ( PI2 * bin ) / n ) - CV_PI;//  为特征点方向赋值[-180,180]  
19.       cvSeqPush( features, new_feat );  //  
20.       free( new_feat );  
21.     }  
22.     }  
23. }  

 7.构建特征描述子

目前每个特征点具有属性有位置、方向、尺度三个信息，现在要用一个向量去描述这个特征点，使其具有高度的唯一特征性。

1.lowe采用了把特征点邻域划分成 d×d (lowe建议d=4) 个子区域，然后再统计每个子区域的方向直方图（8个方向），直方图横轴有8个bin,纵轴为梯度幅值（×权重）的累加。这样描述这个特征点的向量为4×4×8=128维。每个子区域的宽度建议为3×octv，octv为组内的尺度。考虑到插值问题，特征点的邻域范围边长为3×octv×(d+1)，考虑到旋转问题,邻域的范围边长为3×octv×(d+1)×sqrt(2)，最后半径为：

2.把坐标系旋转到主方向位置，再次统计邻域内所有像素点的梯度幅值与方向，计算所在子区域，并把幅值×权重累加到这个子区域的直方图上。

算法上即统计每个邻域的方向直方图时，全部是相对于这个特征点的主方向的方向。如果主方向为30度,某个像素点的梯度方向为50度，这时统计到该子区域直方图上就成了20度。同时由于旋转，这时权重也必须是按旋转后的坐标。

计算所在的子区域的位置：

  

权重按高斯加权函数，系数为描述子宽度的一半，即0.5d：

[cpp]  view plain  copy

 print ?

1. static double*** descr_hist( IplImage* img, int r, int c, double ori,  
2.                  double scl, int d, int n )  
3. {  
4.   double*** hist;  
5.   double cos_t, sin_t, hist_width, exp_denom, r_rot, c_rot, grad_mag,  
6.     grad_ori, w, rbin, cbin, obin, bins_per_rad, PI2 = 2.0 * CV_PI;  
7.   int radius, i, j;  
8.   
9.   hist = calloc( d, sizeof( double** ) );  
10.   for( i = 0; i < d; i++ )  
11.     {  
12.       hist[i] = calloc( d, sizeof( double* ) );  
13.       for( j = 0; j < d; j++ )  
14.     hist[i][j] = calloc( n, sizeof( double ) );  
15.     }  
16.     
17.   cos_t = cos( ori );  
18.   sin_t = sin( ori );  
19.   bins_per_rad = n / PI2;  
20.   exp_denom = d * d * 0.5;  
21.   hist_width = SIFT_DESCR_SCL_FCTR * scl;  
22.   radius = hist_width * sqrt(2) * ( d + 1.0 ) * 0.5 + 0.5;   //计算邻域范围半径，+0.5为了取得更多元素  
23.   for( i = -radius; i <= radius; i++ )  
24.     for( j = -radius; j <= radius; j++ )  
25.       {  
26.     /* 
27.       Calculate sample's histogram array coords rotated relative to ori. 
28.       Subtract 0.5 so samples that fall e.g. in the center of row 1 (i.e. 
29.       r_rot = 1.5) have full weight placed in row 1 after interpolation. 
30.     */  
31.     c_rot = ( j * cos_t - i * sin_t ) / hist_width;              //  
32.     r_rot = ( j * sin_t + i * cos_t ) / hist_width;  
33.     rbin = r_rot + d / 2 - 0.5;                                        //旋转后对应的x``及y''  
34.     cbin = c_rot + d / 2 - 0.5;  
35.       
36.     if( rbin > -1.0  &&  rbin < d  &&  cbin > -1.0  &&  cbin < d )  
37.       if( calc_grad_mag_ori( img, r + i, c + j, &grad_mag, &grad_ori ))       //计算每一个像素点的梯度方向、幅值、  
38.         {  
39.           grad_ori -= ori;              //每个像素点相对于特征点的梯度方向  
40.           while( grad_ori < 0.0 )  
41.         grad_ori += PI2;  
42.           while( grad_ori >= PI2 )  
43.         grad_ori -= PI2;  
44.             
45.           obin = grad_ori * bins_per_rad;     //像素梯度方向转换成8个方向  
46.           w = exp( -(c_rot * c_rot + r_rot * r_rot) / exp_denom );     //每个子区域内像素点，对应的权重  
47.           interp_hist_entry( hist, rbin, cbin, obin, grad_mag * w, d, n );   //为了减小突变影响对附近三个bin值三线性插值处理  
48.         }  
49.       }  
50.   
51.   return hist;  
52. }  

每个维度上bin值累加方法，即计算一个像素的幅值对于相邻的方向，以及位置的贡献，dr，dc为相邻位置，do为相邻方向

， 这就是128维向量的数据，计算方法

[cpp]  view plain  copy

 print ?

1. static void interp_hist_entry( double*** hist, double rbin, double cbin,  
2.                    double obin, double mag, int d, int n )  
3. {  
4.   double d_r, d_c, d_o, v_r, v_c, v_o;  
5.   double** row, * h;  
6.   int r0, c0, o0, rb, cb, ob, r, c, o;  
7.   
8.   r0 = cvFloor( rbin );  
9.   c0 = cvFloor( cbin );  
10.   o0 = cvFloor( obin );  
11.   d_r = rbin - r0;  
12.   d_c = cbin - c0;  
13.   d_o = obin - o0;  
14.   
15.   /* 
16.     The entry is distributed into up to 8 bins.  Each entry into a bin 
17.     is multiplied by a weight of 1 - d for each dimension, where d is the 
18.     distance from the center value of the bin measured in bin units. 
19.   */  
20.   for( r = 0; r <= 1; r++ )  
21.     {  
22.       rb = r0 + r;  
23.       if( rb >= 0  &&  rb < d )  
24.     {  
25.       v_r = mag * ( ( r == 0 )? 1.0 - d_r : d_r );  
26.       row = hist[rb];  
27.       for( c = 0; c <= 1; c++ )  
28.         {  
29.           cb = c0 + c;  
30.           if( cb >= 0  &&  cb < d )  
31.         {  
32.           v_c = v_r * ( ( c == 0 )? 1.0 - d_c : d_c );  
33.           h = row[cb];  
34.           for( o = 0; o <= 1; o++ )  
35.             {  
36.               ob = ( o0 + o ) % n;  
37.               v_o = v_c * ( ( o == 0 )? 1.0 - d_o : d_o );  
38.               h[ob] += v_o;  
39.             }  
40.         }  
41.         }  
42.     }  
43.     }  
44. }  

最后为了去除光照的影响，对128维向量进行归一化处理，同时设置门限，大于0.2的梯度幅值截断

[cpp]  view plain  copy

 print ?

1. static void hist_to_descr( double*** hist, int d, int n, struct feature* feat )  
2. {  
3.   int int_val, i, r, c, o, k = 0;  
4.   
5.   for( r = 0; r < d; r++ )  
6.     for( c = 0; c < d; c++ )  
7.       for( o = 0; o < n; o++ )  
8.     feat->descr[k++] = hist[r][c][o];  
9.   
10.   feat->d = k;  
11.   normalize_descr( feat );          //向量归一化  
12.   for( i = 0; i < k; i++ )  
13.     if( feat->descr[i] > SIFT_DESCR_MAG_THR )   //设置门限，门限为0.2  
14.       feat->descr[i] = SIFT_DESCR_MAG_THR;  
15.   normalize_descr( feat );      //向量归一化  
16.   
17.   /* convert floating-point descriptor to integer valued descriptor */  
18.   for( i = 0; i < k; i++ )              //换成整形值  
19.     {  
20.       int_val = SIFT_INT_DESCR_FCTR * feat->descr[i];      
21.       feat->descr[i] = MIN( 255, int_val );  
22.     }  
23. }  

最后对特征点按尺度大小进行排序，强特征点放在前面；

这样每个特征点就对应一个128维的向量，接下来可以用可以用向量做以后的匹配工作了。

特征点匹配原理后序文章会更新~

------------------------------------------------------------------------------------

在此非常感谢优快云上几位图像上的大牛，我也是通过他们的文章去学习研究的，本文也是参考了他们的文章才写成！

推荐看大牛们的文章，原理写的很好！

http://blog.csdn.net/abcjennifer/article/details/7639681

http://blog.csdn.Net/zddblog/article/details/7521424

http://blog.youkuaiyun.com/chen825919148/article/details/7685952

http://blog.youkuaiyun.com/xiaowei_cqu/article/details/8069548