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

Rob Hess的SIFT算法实现理解及注释

        SIFT算法不用我多解释了，这是一个很强大的算法，主要用于图像配准和物体识别等领域，但是其计算量相比也比较大，性价比比较高的算法包括PCA-SIFT和SURF其中OpenCV提供了SURF算法，但是为了方便理解。这里给出了Rob Hess所实现的SIFT算法的实现以及注释，结合我自己的理解，如果，您有关于SIFT算法不理解的地方咱们可以一起交流一下。或者您认为不详细的地方提出来。

       SIFT算法的主要实现在sift.c这个文件，其主要流程为：

(1)首先创建初始图像，即通过将图像转换为32位的灰度图，然后将图像使用三次插值来方大，之后通过高斯模糊处理

(2)在此基础上进行高斯金字塔的构建以及高斯差分金字塔的构建

(3)对图像进行极值点检测

(4)计算特征向量的尺度

(5)调整图像大小

(6)计算特征的方向

(7)计算描述子，其中包括计算二维方向直方图并转换直方图为特征描述子

首先给出sift算法的整体框架代码：

输入参数：

img为输入图像；

feat为所要提取的特征指针；

intvl指的是高斯金字塔和差分金字塔的层数；

sigma指的是图像初始化过程中高斯模糊所使用的参数；

  contr_thr是归一化之后的去除不稳定特征的阈值；

curv_thr指的是去除边缘的特征的主曲率阈值；

img_dbl是是否将图像放大为之前的两倍；

descr_with用来计算特征描述子的方向直方图的宽度；

descr_hist_bins是直方图中的条数

1. int _sift_features( IplImage* img, struct feature** feat, int intvls,  
2.                    double sigma, double contr_thr, int curv_thr,  
3.                    int img_dbl, int descr_width, int descr_hist_bins )  
4. {  
5.     IplImage* init_img;  
6.     IplImage*** gauss_pyr, *** dog_pyr;  
7.     CvMemStorage* storage;  
8.     CvSeq* features;  
9.     int octvs, i, n = 0;  
10.   
11.     /* check arguments */  
12.     if( ! img )  
13.         fatal_error( "NULL pointer error, %s, line %d",  __FILE__, __LINE__ );  
14.   
15.     if( ! feat )  
16.         fatal_error( "NULL pointer error, %s, line %d",  __FILE__, __LINE__ );  
17.   
18.     /* build scale space pyramid; smallest dimension of top level is ~4 pixels */  
19.     /* 构建高斯尺度空间金字塔，顶层最小的为4像素 */  
20.     init_img = create_init_img( img, img_dbl, sigma );  
21.     octvs = log( double MIN( init_img->width, init_img->height ) ) / log(2.0) - 2;  
22.     //构建高斯金字塔和高斯差分金字塔  
23.     gauss_pyr = build_gauss_pyr( init_img, octvs, intvls, sigma );  
24.     dog_pyr = build_dog_pyr( gauss_pyr, octvs, intvls );  
25.   
26.     storage = cvCreateMemStorage( 0 );  
27.   
28.     //尺度空间极值点检测  
29.     features = scale_space_extrema( dog_pyr, octvs, intvls, contr_thr,  
30.         curv_thr, storage );  
31.       
32.     //画出去除低对比度的极值点  
33.     //draw_extrempoint(img , features);  
34.   
35.   
36.   
37.   
38.     //计算特征向量的尺度  
39.     calc_feature_scales( features, sigma, intvls );  
40.     if( img_dbl )  
41.         adjust_for_img_dbl( features );  
42.     //计算特征的方向  
43.     calc_feature_oris( features, gauss_pyr );  
44.     //计算描述子，包括计算二维方向直方图和转换其为特征描述子  
45.     compute_descriptors( features, gauss_pyr, descr_width, descr_hist_bins );  
46.   
47.     /* sort features by decreasing scale and move from CvSeq to array */  
48.     cvSeqSort( features, (CvCmpFunc)feature_cmp, NULL );  
49.     n = features->total;  
50.     *feat = static_cast<feature *>( calloc( n, sizeof(struct feature) ) );  
51.     *feat = static_cast<feature *>( cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ ) );  
52.   
53.   
54.   
55.   
56.     for( i = 0; i < n; i++ )  
57.     {  
58.         free( (*feat)[i].feature_data );  
59.         (*feat)[i].feature_data = NULL;  
60.     }  
61.   
62.     cvReleaseMemStorage( &storage );  
63.     cvReleaseImage( &init_img );  
64.     release_pyr( &gauss_pyr, octvs, intvls + 3 );  
65.     release_pyr( &dog_pyr, octvs, intvls + 2 );  
66.     return n;  
67. }  

(1)初始化图像

输入参数：

这里不需要解释了

该函数主要用来初始化图像，转换图像为32位灰度图以及进行高斯模糊。

1. static IplImage* create_init_img( IplImage* img, int img_dbl, double sigma )  
2. {  
3.     IplImage* gray, * dbl;  
4.     float sig_diff;  
5.   
6.     gray = convert_to_gray32( img );  
7.     if( img_dbl )  
8.     {  
9.         sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4 );  
10.         dbl = cvCreateImage( cvSize( img->width*2, img->height*2 ),  
11.             IPL_DEPTH_32F, 1 );  
12.         cvResize( gray, dbl, CV_INTER_CUBIC );  
13.         cvSmooth( dbl, dbl, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff );  
14.         cvReleaseImage( &gray );  
15.         return dbl;  
16.     }  
17.     else  
18.     {  
19.         sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA );  
20.         cvSmooth( gray, gray, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff );  
21.         return gray;  
22.     }  
23. }  

(2)构建高斯金字塔

输入参数：

octvs是高斯金字塔的组

invls是高斯金字塔的层数

sigma是初始的高斯模糊参数，后续也通过它计算每一层所使用的sigma

1. <span style="font-size:13px;">static IplImage*** build_gauss_pyr( IplImage* base, int octvs,int intvls, double sigma )  
2. {  
3.     IplImage*** gauss_pyr;  
4.     double* sig = static_cast<double *>( calloc( intvls + 3, sizeof(double)) );  
5.     double sig_total, sig_prev, k;  
6.     int i, o;  
7.   
8.     gauss_pyr = static_cast<IplImage ***>( calloc( octvs, sizeof( IplImage** ) ) );  
9.     for( i = 0; i < octvs; i++ )  
10.         gauss_pyr[i] = static_cast<IplImage **>( calloc( intvls + 3, sizeof( IplImage* ) ) );  
11.   
12.     /* 
13.         precompute Gaussian sigmas using the following formula: 
14.         预计算每次高斯模糊的sigma 
15.  
16.         \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2 
17.     */  
18.     sig[0] = 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.         sig[i] = sqrt( sig_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  &&  i == 0 )  
33.                 gauss_pyr[o][i] = cvCloneImage(base);  
34.   
35.             /* base of new octvave is halved image from end of previous octave */  
36.             //每一组的第一层都是通过对前面一组的最上面一层的降采样实现的  
37.             else if( i == 0 )  
38.                 gauss_pyr[o][i] = downsample( gauss_pyr[o-1][intvls] );  
39.   
40.             /* blur the current octave's last image to create the next one */  
41.             //每一组的其他层则使通过使用不同sigma的高斯模糊来进行处理  
42.             else  
43.             {  
44.                 gauss_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i-1]),  
45.                     IPL_DEPTH_32F, 1 );  
46.                 cvSmooth( gauss_pyr[o][i-1], gauss_pyr[o][i],  
47.                     CV_GAUSSIAN, 0, 0, sig[i], sig[i] );  
48.             }  
49.         }  
50.   
51.     free( sig );  
52.     return gauss_pyr;  
53. }</span>  

降采样处理

输入参数：

不解释

这就是降采样，其实就是将图像通过最近邻算法缩小为原来的一半

1. static IplImage* downsample( IplImage* img )  
2. {  
3.     IplImage* smaller = cvCreateImage( cvSize(img->width / 2, img->height / 2),  
4.         img->depth, img->nChannels );  
5.     cvResize( img, smaller, CV_INTER_NN );  
6.   
7.     return smaller;  
8. }  

(3)构建高斯差分金字塔

输入参数：

不解释了参见上面的说明即可

实际上差分金字塔的构成是通过对相邻层的图像进行相减获得的

1. <span style="font-size:16px;">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 = static_cast<IplImage ***>( calloc( octvs, sizeof( IplImage** ) ) );  
7.     for( i = 0; i < octvs; i++ )  
8.         dog_pyr[i] = static_cast<IplImage **>( 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 );  
16.         }  
17.   
18.     return dog_pyr;  
19. }</span>  

(4)极值点检测

输入参数：

contr_thr是去除对比度低的点所采用的阈值

curv_thr是去除边缘特征的阈值

1. static CvSeq* scale_space_extrema( IplImage*** dog_pyr, int octvs, int intvls,  
2.                                    double contr_thr, int curv_thr,  
3.                                    CvMemStorage* storage )  
4. {  
5.     CvSeq* features;  
6.     double prelim_contr_thr = 0.5 * contr_thr / intvls;  
7.     struct feature* feat;  
8.     struct detection_data* ddata;  
9.     int o, i, r, c;  
10.   
11.     features = cvCreateSeq( 0, sizeof(CvSeq), sizeof(struct feature), storage );  
12.     for( o = 0; o < octvs; o++ )  
13.         for( i = 1; i <= intvls; i++ )  
14.             for(r = SIFT_IMG_BORDER; r < dog_pyr[o][0]->height-SIFT_IMG_BORDER; r++)  
15.                 for(c = SIFT_IMG_BORDER; c < dog_pyr[o][0]->width-SIFT_IMG_BORDER; c++)  
16.                     /* perform preliminary check on contrast */  
17.                     if( ABS( pixval32f( dog_pyr[o][i], r, c ) ) > prelim_contr_thr )  
18.                         if( is_extremum( dog_pyr, o, i, r, c ) )  
19.                         {  
20.                             feat = interp_extremum(dog_pyr, o, i, r, c, intvls, contr_thr);  
21.                             if( feat )  
22.                             {  
23.                                 ddata = feat_detection_data( feat );  
24.                                 if( ! is_too_edge_like( dog_pyr[ddata->octv][ddata->intvl],  
25.                                     ddata->r, ddata->c, curv_thr ) )  
26.                                 {  
27.                                     cvSeqPush( features, feat );  
28.                                 }  
29.                                 else  
30.                                     free( ddata );  
31.                                 free( feat );  
32.                             }  
33.                         }  
34.   
35.     return features;  
36. }  

SIFT_IMG_BORDER是预定义的图像边缘；

通过和对比度阈值比较去掉低对比度的点；

而通过is_extremum来判断是否为极值点，如果是则通过极值点插值的方式获取亚像素的极值点的位置。

然后通过is_too_eage_like和所给的主曲率阈值判断是否为边缘点

*判断是否为极值点

其原理为:通过和高斯金字塔的上一层的9个像素+本层的除了本像素自己的其他的8个像素和下一层的9个像素进行比较看是否为这26个像素中最小的一个或者是否为最大的一个，如果是则为极值点。

1. static int is_extremum( IplImage*** dog_pyr, int octv, int intvl, int r, int c )  
2. {  
3.     float val = pixval32f( dog_pyr[octv][intvl], r, c );  
4.     int i, j, k;  
5.   
6.     /* check for maximum */  
7.     if( val > 0 )  
8.     {  
9.         for( i = -1; i <= 1; i++ )  
10.             for( j = -1; j <= 1; j++ )  
11.                 for( k = -1; k <= 1; k++ )  
12.                     if( val < pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )  
13.                         return 0;  
14.     }  
15.   
16.     /* check for minimum */  
17.     else  
18.     {  
19.         for( i = -1; i <= 1; i++ )  
20.             for( j = -1; j <= 1; j++ )  
21.                 for( k = -1; k <= 1; k++ )  
22.                     if( val > pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )  
23.                         return 0;  
24.     }  
25.   
26.     return 1;  
27. }  

*获取亚像素的极值点的位置

1. static struct feature* interp_extremum( IplImage*** dog_pyr, int octv, int intvl,  
2.                                         int r, int c, int intvls, double contr_thr )  
3. {  
4.     struct feature* feat;  
5.     struct detection_data* ddata;  
6.     double xi, xr, xc, contr;//分别为亚像素的intval,row,col的偏移offset，和对比度  
7.     int i = 0;  
8.   
9.     while( i < SIFT_MAX_INTERP_STEPS )//重新确定极值点并重新定位的操作只能循环 5次  
10.     {  
11.         interp_step( dog_pyr, octv, intvl, r, c, &xi, &xr, &xc );  
12.         if( ABS( xi ) < 0.5  &&  ABS( xr ) < 0.5  &&  ABS( xc ) < 0.5 )//如果满足条件就停止寻找  
13.             break;  
14.         //否则继续寻找极值点  
15.         c += cvRound( xc );  
16.         r += cvRound( xr );  
17.         intvl += cvRound( xi );  
18.   
19.         if( intvl < 1  ||  
20.             intvl > intvls  ||  
21.             c < SIFT_IMG_BORDER  ||  
22.             r < SIFT_IMG_BORDER  ||  
23.             c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER  ||  
24.             r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER )  
25.         {  
26.             return NULL;  
27.         }  
28.   
29.         i++;  
30.     }  
31.       
32.     //确保极值点是经过最大5步找到的  
33.     /* ensure convergence of interpolation */  
34.     if( i >= SIFT_MAX_INTERP_STEPS )  
35.         return NULL;  
36.       
37.     //获取找到的极值点的对比度  
38.     contr = interp_contr( dog_pyr, octv, intvl, r, c, xi, xr, xc );  
39.     //判断极值点是否小于某一个阈值  
40.     if( ABS( contr ) < contr_thr / intvls )  
41.         return NULL;  
42.     //若小于，则认为是极值点  
43.     feat = new_feature();  
44.     ddata = feat_detection_data( feat );  
45.     feat->img_pt.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;  
48.     ddata->c = c;  
49.     ddata->octv = octv;  
50.     ddata->intvl = intvl;  
51.     ddata->subintvl = xi;  
52.   
53.     return feat;  
54. }  

*获取亚像素位置中所用到的函数

1. static void interp_step( IplImage*** dog_pyr, int octv, int intvl, int r, int c,  
2.                          double* xi, double* xr, double* xc )  
3. {  
4.     CvMat* dD, * H, * H_inv, X;  
5.     double x[3] = { 0 };  
6.       
7.   
8.     //计算三维偏导数  
9.     dD = deriv_3D( dog_pyr, octv, intvl, r, c );  
10.     //计算三维海森矩阵  
11.     H = hessian_3D( dog_pyr, octv, intvl, r, c );  
12.     H_inv = cvCreateMat( 3, 3, CV_64FC1 );  
13.     cvInvert( H, H_inv, CV_SVD );  
14.     cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );  
15.   
16.     cvGEMM( H_inv, dD, -1, NULL, 0, &X, 0 );  
17.   
18.     cvReleaseMat( &dD );  
19.     cvReleaseMat( &H );  
20.     cvReleaseMat( &H_inv );  
21.   
22.     *xi = x[2];  
23.     *xr = x[1];  
24.     *xc = x[0];  
25. }  

*计算三维偏导数

计算在x和y方向上的偏导数，高斯差分尺度空间金字塔中像素的尺度

实际上在离散数据中计算偏导数是通过相邻像素的相减来计算的

比如说计算x方向的偏导数dx，则通过该向所的x方向的后一个减去前一个然后除以2即可求的dx

1. static CvMat* deriv_3D( IplImage*** dog_pyr, int octv, int intvl, int r, int c )  
2. {  
3.     CvMat* dI;  
4.     double dx, dy, ds;  
5.   
6.     dx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) -  
7.         pixval32f( dog_pyr[octv][intvl], r, c-1 ) ) / 2.0;  
8.     dy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) -  
9.         pixval32f( dog_pyr[octv][intvl], r-1, c ) ) / 2.0;  
10.     ds = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) -  
11.         pixval32f( dog_pyr[octv][intvl-1], r, c ) ) / 2.0;  
12.   
13.     dI = cvCreateMat( 3, 1, CV_64FC1 );  
14.     cvmSet( dI, 0, 0, dx );  
15.     cvmSet( dI, 1, 0, dy );  
16.     cvmSet( dI, 2, 0, ds );  
17.   
18.     return dI;  
19. }  

*计算三维海森矩阵

不需要讲什么，其实就是计算二次导数，计算方法也和一次导数的计算如出一辙。

然后将结果放入到一个矩阵中去。

1. static CvMat* hessian_3D( IplImage*** dog_pyr, int octv, int intvl, int r, int c )  
2. {  
3.     CvMat* H;  
4.     double v, dxx, dyy, dss, dxy, dxs, dys;  
5.   
6.     v = pixval32f( dog_pyr[octv][intvl], r, c );  
7.     dxx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) +   
8.             pixval32f( dog_pyr[octv][intvl], r, c-1 ) - 2 * v );  
9.     dyy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) +  
10.             pixval32f( dog_pyr[octv][intvl], r-1, c ) - 2 * v );  
11.     dss = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) +  
12.             pixval32f( dog_pyr[octv][intvl-1], r, c ) - 2 * v );  
13.     dxy = ( pixval32f( dog_pyr[octv][intvl], r+1, c+1 ) -  
14.             pixval32f( dog_pyr[octv][intvl], r+1, c-1 ) -  
15.             pixval32f( dog_pyr[octv][intvl], r-1, c+1 ) +  
16.             pixval32f( dog_pyr[octv][intvl], r-1, c-1 ) ) / 4.0;  
17.     dxs = ( pixval32f( dog_pyr[octv][intvl+1], r, c+1 ) -  
18.             pixval32f( dog_pyr[octv][intvl+1], r, c-1 ) -  
19.             pixval32f( dog_pyr[octv][intvl-1], r, c+1 ) +  
20.             pixval32f( dog_pyr[octv][intvl-1], r, c-1 ) ) / 4.0;  
21.     dys = ( pixval32f( dog_pyr[octv][intvl+1], r+1, c ) -  
22.             pixval32f( dog_pyr[octv][intvl+1], r-1, c ) -  
23.             pixval32f( dog_pyr[octv][intvl-1], r+1, c ) +  
24.             pixval32f( dog_pyr[octv][intvl-1], r-1, c ) ) / 4.0;  
25.   
26.     H = cvCreateMat( 3, 3, CV_64FC1 );  
27.     cvmSet( H, 0, 0, dxx );  
28.     cvmSet( H, 0, 1, dxy );  
29.     cvmSet( H, 0, 2, dxs );  
30.     cvmSet( H, 1, 0, dxy );  
31.     cvmSet( H, 1, 1, dyy );  
32.     cvmSet( H, 1, 2, dys );  
33.     cvmSet( H, 2, 0, dxs );  
34.     cvmSet( H, 2, 1, dys );  
35.     cvmSet( H, 2, 2, dss );  
36.   
37.     return H;  
38. }  

*计算插入像素的对比度

1. static double interp_contr( IplImage*** dog_pyr, int octv, int intvl, int r,  
2.                             int c, double xi, double xr, double xc )  
3. {  
4.     CvMat* dD, X, T;  
5.     double t[1], x[3] = { xc, xr, xi };  
6.   
7.     cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );  
8.     cvInitMatHeader( &T, 1, 1, CV_64FC1, t, CV_AUTOSTEP );  
9.     dD = deriv_3D( dog_pyr, octv, intvl, r, c );  
10.     cvGEMM( dD, &X, 1, NULL, 0, &T,  CV_GEMM_A_T );  
11.     cvReleaseMat( &dD );  
12.   
13.     return pixval32f( dog_pyr[octv][intvl], r, c ) + t[0] * 0.5;  
14. }  

其中cvGEMM是矩阵的通用计算函数，至于CV_GEMM_A_T是计算dD的转置矩阵放入T中

*去除边缘相应
通过计算所在特征向量的主曲率半径来判断特征是边缘的从而导致不稳定
即去除边缘响应

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);  
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 )  
19.         return 0;  
20.     return 1;  
21. }  

(4)计算特征向量的尺度

实际上是通过最初的sigma来获得每一层每一组的尺度

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 );  
15.         ddata->scl_octv = sigma * pow( 2.0, intvl / intvls );  
16.     }  
17. }  

(5)调整图像特征坐标、尺度、点的坐标大小为原来的一半

1. static void adjust_for_img_dbl( CvSeq* features )  
2. {  
3.     struct feature* feat;  
4.     int i, n;  
5.   
6.     n = features->total;  
7.     for( i = 0; i < n; i++ )  
8.     {  
9.         feat = CV_GET_SEQ_ELEM( struct feature, features, i );  
10.         feat->x /= 2.0;  
11.         feat->y /= 2.0;  
12.         feat->scl /= 2.0;  
13.         feat->img_pt.x /= 2.0;  
14.         feat->img_pt.y /= 2.0;  
15.     }  
16. }  

(6)给每一个图像特征向量计算规范化的方向

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.   
10.     //遍历整个检测出来的特征点，计算每个特征点的直方图，然后平滑直方图去除突变，然后找到每一个特征点的主方向，并加入到好的方向特征数组中去  
11.     for( i = 0; i < n; i++ )  
12.     {  
13.         feat = static_cast<feature *>( malloc( sizeof( struct feature ) ) );  
14.         cvSeqPopFront( features, feat );  
15.         ddata = feat_detection_data( feat );  
16.         //计算给定的某个像素的灰度方向直方图  
17.         hist = ori_hist( gauss_pyr[ddata->octv][ddata->intvl],  
18.                         ddata->r, ddata->c, SIFT_ORI_HIST_BINS,  
19.                         cvRound( SIFT_ORI_RADIUS * ddata->scl_octv ),  
20.                         SIFT_ORI_SIG_FCTR * ddata->scl_octv );  
21.         for( j = 0; j < SIFT_ORI_SMOOTH_PASSES; j++ )  
22.             smooth_ori_hist( hist, SIFT_ORI_HIST_BINS );  
23.         omax = dominant_ori( hist, SIFT_ORI_HIST_BINS );  
24.   
25.         //描述子向量元素门限化  
26.         add_good_ori_features( features, hist, SIFT_ORI_HIST_BINS,  
27.                                 omax * SIFT_ORI_PEAK_RATIO, feat );  
28.         free( ddata );  
29.         free( feat );  
30.         free( hist );  
31.     }  
32. }  

*对所给像素计算灰度方向直方图
以关键点为中心的邻域窗口内采样，并用直方图统计邻域像素的梯度
方向。梯度直方图的范围是0～360度，其中每10度一个柱，总共36个柱

1. static double* ori_hist( IplImage* img, int r, int c, int n, int rad, double sigma)  
2. {  
3.     double* hist;  
4.     double mag, ori, w, exp_denom, PI2 = CV_PI * 2.0;  
5.     int bin, i, j;  
6.   
7.     hist = static_cast<double *>( calloc( n, sizeof( double ) ) );  
8.     exp_denom = 2.0 * sigma * sigma;  
9.     for( i = -rad; i <= rad; i++ )  
10.         for( j = -rad; j <= rad; j++ )  
11.             if( calc_grad_mag_ori( img, r + i, c + j, &mag, &ori ) )  
12.             {  
13.                 w = exp( -( i*i + j*j ) / exp_denom );  
14.                 bin = cvRound( n * ( ori + CV_PI ) / PI2 );  
15.                 bin = ( bin < n )? bin : 0;  
16.                 hist[bin] += w * mag;  
17.             }  
18.   
19.     return hist;  
20. }  

*计算所给像素的梯度大小和方向
每一个小格都代表了特征点邻域所在的尺度空间的一个像素 ，箭头方向代表了像素梯
度方向，箭头长度代表该像素的幅值也就是梯度的值

1. static int calc_grad_mag_ori( IplImage* img, int r, int c, double* mag, double* ori )  
2. {  
3.     double dx, dy;  
4.   
5.     if( r > 0  &&  r < img->height - 1  &&  c > 0  &&  c < img->width - 1 )  
6.     {  
7.         dx = pixval32f( img, r, c+1 ) - pixval32f( img, r, c-1 );  
8.         dy = pixval32f( img, r-1, c ) - pixval32f( img, r+1, c );  
9.         *mag = sqrt( dx*dx + dy*dy );  
10.         *ori = atan2( dy, dx );  
11.         return 1;  
12.     }  
13.   
14.     else  
15.         return 0;  
16. }  

*对方向直方图进行高斯模糊
使用高斯函数对直方图进行平滑，减少突变的影响。

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. }  

*在直方图中找到主方向的梯度
利用关键点邻域像素的梯度方向分布特性为每个关键点指定方向参数，使算子具备
旋转不变性。

1. static double dominant_ori( double* hist, int n )  
2. {  
3.     double omax;  
4.     int maxbin, i;  
5.   
6.     omax = hist[0];  
7.     maxbin = 0;  
8.     for( i = 1; i < n; i++ )  
9.         if( hist[i] > omax )  
10.         {  
11.             omax = hist[i];  
12.             maxbin = i;  
13.         }  
14.     return omax;  
15. }  

*将大于某一个梯度大小阈值的特征向量加入到直方图中去

n为方向的个数

1. <span style="font-size:18px;">mag_thr描述子向量门限一般取0.2</span>  

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.   
14.         //描述子向量门限化，一般门限取0.2  
15.         if( hist[i] > hist[l]  &&  hist[i] > hist[r]  &&  hist[i] >= mag_thr )  
16.         {  
17.             bin = i + interp_hist_peak( hist[l], hist[i], hist[r] );  
18.             bin = ( bin < 0 )? n + bin : ( bin >= n )? bin - n : bin;  
19.             new_feat = clone_feature( feat );  
20.             new_feat->ori = ( ( PI2 * bin ) / n ) - CV_PI;  
21.             cvSeqPush( features, new_feat );  
22.             free( new_feat );  
23.         }  
24.     }  
25. }  

(7)计算特征描述子

1. static void compute_descriptors( CvSeq* features, IplImage*** gauss_pyr, int d, int n)  
2. {  
3.     struct feature* feat;  
4.     struct detection_data* ddata;  
5.     double*** hist;  
6.     int i, k = features->total;  
7.   
8.     for( i = 0; i < k; i++ )  
9.     {  
10.         feat = CV_GET_SEQ_ELEM( struct feature, features, i );  
11.         ddata = feat_detection_data( feat );  
12.         //计算二维方向直方图  
13.         hist = descr_hist( gauss_pyr[ddata->octv][ddata->intvl], ddata->r,  
14.             ddata->c, feat->ori, ddata->scl_octv, d, n );  
15.         //将二维方向直方图转换为特征描述子  
16.         hist_to_descr( hist, d, n, feat );  
17.         release_descr_hist( &hist, d );  
18.     }  
19. }  

*计算二维方向直方图

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 = static_cast<double ***>( calloc( d, sizeof( double** ) ) );  
10.     for( i = 0; i < d; i++ )  
11.     {  
12.         hist[i] =static_cast<double **>( calloc( d, sizeof( double* ) ) );  
13.         for( j = 0; j < d; j++ )  
14.             hist[i][j] = static_cast<double *>( 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.0) * ( d + 1.0 ) * 0.5 + 0.5;  
23.     for( i = -radius; i <= radius; i++ )  
24.         for( j = -radius; j <= radius; j++ )  
25.         {  
26.             /* 
27.             即将坐标移至关键点主方向 
28.             计算采用的直方图数组中相对于方向旋转的坐标 
29.             Calculate sample's 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 placed in row 1 after interpolation. 
32.             */  
33.             c_rot = ( j * cos_t - i * sin_t ) / hist_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  &&  rbin < d  &&  cbin > -1.0  &&  cbin < d )  
39.                 if( calc_grad_mag_ori( img, r + i, c + j, &grad_mag, &grad_ori ))  
40.                 {  
41.                     grad_ori -= ori;  
42.                     while( grad_ori < 0.0 )  
43.                         grad_ori += PI2;  
44.                     while( grad_ori >= PI2 )  
45.                         grad_ori -= PI2;  
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, cbin, obin, grad_mag * w, d, n );  
50.                 }  
51.         }  
52.   
53.     return hist;  
54. }  

*插入一个entry进入到方向直方图中从而形成特征描述子

这个，我也不怎么明白。。。

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. }  

*将二维直方图转换为特征描述子

实际上是归一化描述子和转换为整数

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 )  
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. }  

*归一化描述子

1. static void normalize_descr( struct feature* feat )  
2. {  
3.     double cur, len_inv, len_sq = 0.0;  
4.     int i, d = feat->d;//为描述子长度128维  
5.       
6.   
7.     //如何进行归一化特征描述子来降低对光照的影响  
8.     //主要就是将每一个特征的平方求和，然后开方，然后去其倒数，然后乘以每一个特征描述子的梯度值，这样就得到了归一化的特征描述子  
9.     for( i = 0; i < d; i++ )  
10.     {  
11.         cur = feat->descr[i];  
12.         len_sq += cur*cur;  
13.     }  
14.     len_inv = 1.0 / sqrt( len_sq );  
15.     for( i = 0; i < d; i++ )  
16.         feat->descr[i] *= len_inv;  
17. }  

下面给出存储在文件中的SIFT特征分析：

下面给出特征文件的部分特征

1. 114  128    
2. 101.350424   136.130888   40.169873   0.771085  orientation  
3.  0 0 0 0 3 1 0 0 2 23 46 15 18 3 0 0 6 20 13 1  
4.  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 88 36 0 0  
5.  81 95 57 47 185 114 2 7 185 155 19 6 19 6 1 22 22 0 0 0  
6.  0 0 0 1 0 0 0 0 37 8 0 0 91 12 0 1 185 144 11 35  
7.  185 50 0 0 23 28 8 95 40 1 0 0 0 0 0 4 0 0 0 0  
8.  0 0 0 0 11 5 0 0 4 2 0 0 49 20 0 0 1 0 0 1  
9.  0 0 0 0 0 0 0 0  
10. 127.871534 71.100559 15.768594 -2.024589  
11.  1 2 2 72 63 12 1 1 133 93 1 4 2 7 4 44 133 115 0 0  
12.  0 0 0 20 9 4 0 0 0 0 0 0 23 0 1 9 107 20 1 8  
13.  133 5 0 0 0 1 5 133 132 14 0 0 0 0 8 133 14 1 0 0  
14.  0 0 0 8 26 0 0 0 126 37 8 22 133 47 0 0 0 0 3 52  
15.  131 41 0 0 0 0 2 36 1 0 0 0 0 0 0 2 2 0 0 0  
16.  34 105 80 24 111 15 0 0 0 1 55 66 79 21 0 0 0 0 0 5  
17.  0 0 0 0 0 0 0 0  

下面给出说明：

1. 114 特征数目 128  向量维度  
2.   
3. 关键点坐标  
4. 101.350424  y坐标  
5. 136.130888  x坐标  
6. 40.169873  scale 尺度  
7. 0.771085  orientation  关键点的梯度方向  
8.   
9.   
10. 16个种子点的8个方向向量的信息共128个信息  
11.   
12.  0 0 0 0 3 1 0 0 2 23 46 15 18 3 0 0 6 20 13 1  
13.  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 88 36 0 0  
14.  81 95 57 47 185 114 2 7 185 155 19 6 19 6 1 22 22 0 0 0  
15.  0 0 0 1 0 0 0 0 37 8 0 0 91 12 0 1 185 144 11 35  
16.  185 50 0 0 23 28 8 95 40 1 0 0 0 0 0 4 0 0 0 0  
17.  0 0 0 0 11 5 0 0 4 2 0 0 49 20 0 0 1 0 0 1  
18.  0 0 0 0 0 0 0 0  
19.   
20. 下一组关键点向量  
21. 127.871534 y坐标  
22. 71.100559 x坐标  
23. 15.768594 尺度=sigma*2^(高斯模糊)  
24. -2.024589 梯度方向m(x,y)。。。。等等  

最后附上一个Rob Hess的源码下载地址吧

http://blogs.oregonstate.edu/hess/code/sift/

下面给出我的毕业设计：基于SIFT算法的图像伪造盲检测的初步实现效果，后续还有待改进。。。。