教你一步一步用c语言实现sift算法、下

    ok，接上文，咱们一个一个的来编写main函数中所涉及到所有函数，这也是本文的关键部分：

    五个步骤

    ok，接下来，进入重点部分，咱们依据上文介绍的sift算法的几个步骤，来一一实现这些函数。

    为了版述清晰，再贴一下，主函数，顺便再加强下对sift 算法的五个步骤的认识：

    1、 SIFT算法第一步：图像预处理

    CvMat *ScaleInitImage(CvMat * im) ; //金字塔初始化

    2、 SIFT算法第二步：建立高斯金字塔函数

    ImageOctaves* BuildGaussianOctaves(CvMat * image) ; //建立高斯金字塔

    3、 SIFT算法第三步：特征点位置检测，最后确定特征点的位置

    int DetectKeypoint(int numoctaves, ImageOctaves *GaussianPyr);

    void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr);

    5、 SIFT算法第五步：抽取各个特征点处的特征描述字

    void ExtractFeatureDescriptors(int numoctaves, ImageOctaves *GaussianPyr);

    ok，接下来一一具体实现这几个函数：

    SIFT算法第一步

    SIFT算法第一步：扩大图像，预滤波剔除噪声，得到金字塔的最底层-第一阶的第一层：

    2. {  
    3.     double sigma,preblur_sigma;  
    4.     CvMat *imMat;  
    5.     CvMat * dst;  
    6.     CvMat *tempMat;  
    7.     //首先对图像进行平滑滤波，抑制噪声  
    8.     imMat = cvCreateMat(im->rows, im->cols, CV_32FC1);  
    9.     BlurImage(im, imMat, INITSIGMA);  
    10.     //针对两种情况分别进行处理：初始化放大原始图像或者在原图像基础上进行后续操作  
    11.     //建立金字塔的最底层  
    12.     if (DOUBLE_BASE_IMAGE_SIZE)   
    13.     {  
    14.         tempMat = doubleSizeImage2(imMat);//对扩大两倍的图像进行二次采样，采样率为0.5，采用线性插值  
    15. #define TEMPMAT(ROW,COL) ((float *)(tempMat->data.fl + tempMat->step/sizeof(float) * (ROW)))[(COL)]  
    17.         dst = cvCreateMat(tempMat->rows, tempMat->cols, CV_32FC1);  
    18.         preblur_sigma = 1.0;//sqrt(2 - 4*INITSIGMA*INITSIGMA);  
    19.         BlurImage(tempMat, dst, preblur_sigma);   
    21.         // The initial blurring for the first image of the first octave of the pyramid.  
    22.         sigma = sqrt( (4*INITSIGMA*INITSIGMA) + preblur_sigma * preblur_sigma );  
    23.         //  sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA * 4);  
    24.         //printf("Init Sigma: %f/n", sigma);  
    25.         BlurImage(dst, tempMat, sigma);       //得到金字塔的最底层-放大2倍的图像  
    26.         cvReleaseMat( &dst );   
    27.         return tempMat;  
    28.     }   
    29.     else   
    30.     {  
    31.         dst = cvCreateMat(im->rows, im->cols, CV_32FC1);  
    32.         //sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA);  
    33.         preblur_sigma = 1.0;//sqrt(2 - 4*INITSIGMA*INITSIGMA);  
    34.         sigma = sqrt( (4*INITSIGMA*INITSIGMA) + preblur_sigma * preblur_sigma );  
    35.         //printf("Init Sigma: %f/n", sigma);  
    36.         BlurImage(imMat, dst, sigma);        //得到金字塔的最底层：原始图像大小  
    37.         return dst;  
    38.     }   
    39. }

    SIFT算法第二步

    SIFT第二步，建立Gaussian金字塔，给定金字塔第一阶第一层图像后，计算高斯金字塔其他尺度图像，
    每一阶的数目由变量SCALESPEROCTAVE决定，给定一个基本图像，计算它的高斯金字塔图像，返回外部向量是阶梯指针，内部向量是每一个阶梯内部的不同尺度图像。

    SIFT算法第三步

    SIFT算法第三步，特征点位置检测，最后确定特征点的位置检测DOG金字塔中的局部最大值，找到之后，还要经过两个检验才能确认为特征点：一是它必须有明显的差异，二是他不应该是边缘点，（也就是说，在极值点处的主曲率比应该小于某一个阈值）。

    1. //SIFT算法第三步，特征点位置检测，  
    2. int DetectKeypoint(int numoctaves, ImageOctaves *GaussianPyr)  
    3. {  
    4.     //计算用于DOG极值点检测的主曲率比的阈值  
    5.     double curvature_threshold;  
    6.     curvature_threshold= ((CURVATURE_THRESHOLD + 1)*(CURVATURE_THRESHOLD + 1))/CURVATURE_THRESHOLD;  
    7. #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(DOGoctaves[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + DOGoctaves[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]  
    9.     int   keypoint_count = 0;     
    10.     for (int i=0; i<numoctaves; i++)    
    11.     {          
    12.         for(int j=1;j<SCALESPEROCTAVE+1;j++)//取中间的scaleperoctave个层  
    13.         {    
    14.             //在图像的有效区域内寻找具有显著性特征的局部最大值  
    15.             //float sigma=(GaussianPyr[i].Octave)[j].levelsigma;  
    16.             //int dim = (int) (max(3.0f, 2.0*GAUSSKERN *sigma + 1.0f)*0.5);  
    17.             int dim = (int)(0.5*((GaussianPyr[i].Octave)[j].levelsigmalength)+0.5);  
    18.             for (int m=dim;m<((DOGoctaves[i].row)-dim);m++)   
    19.                 for(int n=dim;n<((DOGoctaves[i].col)-dim);n++)  
    20.                 {       
    21.                     if ( fabs(ImLevels(i,j,m,n))>= CONTRAST_THRESHOLD )  
    22.                     {  
    24.                         if ( ImLevels(i,j,m,n)!=0.0 )  //1、首先是非零  
    25.                         {  
    26.                             float inf_val=ImLevels(i,j,m,n);  
    27.                             if(( (inf_val <= ImLevels(i,j-1,m-1,n-1))&&  
    28.                                 (inf_val <= ImLevels(i,j-1,m  ,n-1))&&  
    29.                                 (inf_val <= ImLevels(i,j-1,m+1,n-1))&&  
    30.                                 (inf_val <= ImLevels(i,j-1,m-1,n  ))&&  
    31.                                 (inf_val <= ImLevels(i,j-1,m  ,n  ))&&  
    32.                                 (inf_val <= ImLevels(i,j-1,m+1,n  ))&&  
    33.                                 (inf_val <= ImLevels(i,j-1,m-1,n+1))&&  
    34.                                 (inf_val <= ImLevels(i,j-1,m  ,n+1))&&  
    35.                                 (inf_val <= ImLevels(i,j-1,m+1,n+1))&&    //底层的小尺度9  
    37.                                 (inf_val <= ImLevels(i,j,m-1,n-1))&&  
    38.                                 (inf_val <= ImLevels(i,j,m  ,n-1))&&  
    39.                                 (inf_val <= ImLevels(i,j,m+1,n-1))&&  
    40.                                 (inf_val <= ImLevels(i,j,m-1,n  ))&&  
    41.                                 (inf_val <= ImLevels(i,j,m+1,n  ))&&  
    42.                                 (inf_val <= ImLevels(i,j,m-1,n+1))&&  
    43.                                 (inf_val <= ImLevels(i,j,m  ,n+1))&&  
    44.                                 (inf_val <= ImLevels(i,j,m+1,n+1))&&     //当前层8  
    46.                                 (inf_val <= ImLevels(i,j+1,m-1,n-1))&&  
    47.                                 (inf_val <= ImLevels(i,j+1,m  ,n-1))&&  
    48.                                 (inf_val <= ImLevels(i,j+1,m+1,n-1))&&  
    49.                                 (inf_val <= ImLevels(i,j+1,m-1,n  ))&&  
    50.                                 (inf_val <= ImLevels(i,j+1,m  ,n  ))&&  
    51.                                 (inf_val <= ImLevels(i,j+1,m+1,n  ))&&  
    52.                                 (inf_val <= ImLevels(i,j+1,m-1,n+1))&&  
    53.                                 (inf_val <= ImLevels(i,j+1,m  ,n+1))&&  
    54.                                 (inf_val <= ImLevels(i,j+1,m+1,n+1))     //下一层大尺度9          
    55.                                 ) ||   
    56.                                 ( (inf_val >= ImLevels(i,j-1,m-1,n-1))&&  
    57.                                 (inf_val >= ImLevels(i,j-1,m  ,n-1))&&  
    58.                                 (inf_val >= ImLevels(i,j-1,m+1,n-1))&&  
    59.                                 (inf_val >= ImLevels(i,j-1,m-1,n  ))&&  
    60.                                 (inf_val >= ImLevels(i,j-1,m  ,n  ))&&  
    61.                                 (inf_val >= ImLevels(i,j-1,m+1,n  ))&&  
    62.                                 (inf_val >= ImLevels(i,j-1,m-1,n+1))&&  
    63.                                 (inf_val >= ImLevels(i,j-1,m  ,n+1))&&  
    64.                                 (inf_val >= ImLevels(i,j-1,m+1,n+1))&&  
    66.                                 (inf_val >= ImLevels(i,j,m-1,n-1))&&  
    67.                                 (inf_val >= ImLevels(i,j,m  ,n-1))&&  
    68.                                 (inf_val >= ImLevels(i,j,m+1,n-1))&&  
    69.                                 (inf_val >= ImLevels(i,j,m-1,n  ))&&  
    70.                                 (inf_val >= ImLevels(i,j,m+1,n  ))&&  
    71.                                 (inf_val >= ImLevels(i,j,m-1,n+1))&&  
    72.                                 (inf_val >= ImLevels(i,j,m  ,n+1))&&  
    73.                                 (inf_val >= ImLevels(i,j,m+1,n+1))&&   
    75.                                 (inf_val >= ImLevels(i,j+1,m-1,n-1))&&  
    76.                                 (inf_val >= ImLevels(i,j+1,m  ,n-1))&&  
    77.                                 (inf_val >= ImLevels(i,j+1,m+1,n-1))&&  
    78.                                 (inf_val >= ImLevels(i,j+1,m-1,n  ))&&  
    79.                                 (inf_val >= ImLevels(i,j+1,m  ,n  ))&&  
    80.                                 (inf_val >= ImLevels(i,j+1,m+1,n  ))&&  
    81.                                 (inf_val >= ImLevels(i,j+1,m-1,n+1))&&  
    82.                                 (inf_val >= ImLevels(i,j+1,m  ,n+1))&&  
    83.                                 (inf_val >= ImLevels(i,j+1,m+1,n+1))   
    84.                                 ) )      //2、满足26个中极值点  
    85.                             {     
    86.                                 //此处可存储  
    87.                                 //然后必须具有明显的显著性，即必须大于CONTRAST_THRESHOLD=0.02  
    88.                                 if ( fabs(ImLevels(i,j,m,n))>= CONTRAST_THRESHOLD )  
    89.                                 {  
    90.                                     //最后显著处的特征点必须具有足够的曲率比，CURVATURE_THRESHOLD=10.0，首先计算Hessian矩阵  
    91.                                     // Compute the entries of the Hessian matrix at the extrema location.  
    92.                                     /* 
    93.                                     1   0   -1 
    94.                                     0   0   0 
    95.                                     -1   0   1         *0.25 
    96.                                     */  
    97.                                     // Compute the trace and the determinant of the Hessian.  
    98.                                     //Det_H = Dxx*Dyy - Dxy^2;  
    99.                                     float Dxx,Dyy,Dxy,Tr_H,Det_H,curvature_ratio;  
    100.                                     Dxx = ImLevels(i,j,m,n-1) + ImLevels(i,j,m,n+1)-2.0*ImLevels(i,j,m,n);  
    101.                                     Dyy = ImLevels(i,j,m-1,n) + ImLevels(i,j,m+1,n)-2.0*ImLevels(i,j,m,n);  
    103.                                     Tr_H = Dxx + Dyy;  
    104.                                     Det_H = Dxx*Dyy - Dxy*Dxy;  
    105.                                     // Compute the ratio of the principal curvatures.  
    106.                                     curvature_ratio = (1.0*Tr_H*Tr_H)/Det_H;  
    107.                                     if ( (Det_H>=0.0) && (curvature_ratio <= curvature_threshold) )  //最后得到最具有显著性特征的特征点  
    108.                                     {  
    109.                                         //将其存储起来，以计算后面的特征描述字  
    110.                                         keypoint_count++;  
    111.                                         Keypoint k;  
    112.                                         /* Allocate memory for the keypoint. */  
    113.                                         k = (Keypoint) malloc(sizeof(struct KeypointSt));  
    114.                                         k->next = keypoints;  
    115.                                         keypoints = k;  
    116.                                         k->row = m*(GaussianPyr[i].subsample);  
    117.                                         k->col =n*(GaussianPyr[i].subsample);  
    118.                                         k->sy = m;    //行  
    119.                                         k->sx = n;    //列  
    120.                                         k->octave=i;  
    121.                                         k->level=j;  
    122.                                         k->scale = (GaussianPyr[i].Octave)[j].absolute_sigma;        
    123.                                     }//if >curvature_thresh  
    124.                                 }//if >contrast  
    125.                             }//if inf value  
    126.                         }//if non zero  
    127.                     }//if >contrast  
    128.                 }  //for concrete image level col  
    129.         }//for levels  
    130.     }//for octaves  
    131.     return keypoint_count;  
    132. }  
    134. //在图像中，显示SIFT特征点的位置  
    135. void DisplayKeypointLocation(IplImage* image, ImageOctaves *GaussianPyr)  
    136. {  
    138.     Keypoint p = keypoints; // p指向第一个结点  
    139.     while(p) // 没到表尾  
    140.     {     
    141.         cvLine( image, cvPoint((int)((p->col)-3),(int)(p->row)),   
    142.             cvPoint((int)((p->col)+3),(int)(p->row)), CV_RGB(255,255,0),  
    143.             1, 8, 0 );  
    144.         cvLine( image, cvPoint((int)(p->col),(int)((p->row)-3)),   
    145.             cvPoint((int)(p->col),(int)((p->row)+3)), CV_RGB(255,255,0),  
    146.             1, 8, 0 );  
    147.         //  cvCircle(image,cvPoint((uchar)(p->col),(uchar)(p->row)),  
    148.         //   (int)((GaussianPyr[p->octave].Octave)[p->level].absolute_sigma),  
    149.         //   CV_RGB(255,0,0),1,8,0);  
    150.         p=p->next;  
    151.     }   
    152. }  
    154. // Compute the gradient direction and magnitude of the gaussian pyramid images  
    155. void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr)  
    156. {  
    157.     // ImageOctaves *mag_thresh ;  
    158.     mag_pyr=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );  
    159.     grad_pyr=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );  
    160.     // float sigma=( (GaussianPyr[0].Octave)[SCALESPEROCTAVE+2].absolute_sigma ) / GaussianPyr[0].subsample;  
    161.     // int dim = (int) (max(3.0f, 2 * GAUSSKERN *sigma + 1.0f)*0.5+0.5);  
    162. #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(GaussianPyr[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + GaussianPyr[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]  
    163.     for (int i=0; i<numoctaves; i++)    
    164.     {          
    165.         mag_pyr[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE) * sizeof(ImageLevels) );  
    166.         grad_pyr[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE) * sizeof(ImageLevels) );  
    167.         for(int j=1;j<SCALESPEROCTAVE+1;j++)//取中间的scaleperoctave个层  
    168.         {    
    169.             CvMat *Mag = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);  
    170.             CvMat *Ori = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);  
    171.             CvMat *tempMat1 = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);  
    172.             CvMat *tempMat2 = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);  
    173.             cvZero(Mag);  
    174.             cvZero(Ori);  
    175.             cvZero(tempMat1);  
    176.             cvZero(tempMat2);   
    177. #define MAG(ROW,COL) ((float *)(Mag->data.fl + Mag->step/sizeof(float) *(ROW)))[(COL)]     
    178. #define ORI(ROW,COL) ((float *)(Ori->data.fl + Ori->step/sizeof(float) *(ROW)))[(COL)]    
    179. #define TEMPMAT1(ROW,COL) ((float *)(tempMat1->data.fl + tempMat1->step/sizeof(float) *(ROW)))[(COL)]  
    180. #define TEMPMAT2(ROW,COL) ((float *)(tempMat2->data.fl + tempMat2->step/sizeof(float) *(ROW)))[(COL)]  
    181.             for (int m=1;m<(GaussianPyr[i].row-1);m++)   
    182.                 for(int n=1;n<(GaussianPyr[i].col-1);n++)  
    183.                 {  
    184.                     //计算幅值  
    185.                     TEMPMAT1(m,n) = 0.5*( ImLevels(i,j,m,n+1)-ImLevels(i,j,m,n-1) );  //dx  
    186.                     TEMPMAT2(m,n) = 0.5*( ImLevels(i,j,m+1,n)-ImLevels(i,j,m-1,n) );  //dy  
    187.                     MAG(m,n) = sqrt(TEMPMAT1(m,n)*TEMPMAT1(m,n)+TEMPMAT2(m,n)*TEMPMAT2(m,n));  //mag  
    188.                     //计算方向  
    189.                     ORI(m,n) =atan( TEMPMAT2(m,n)/TEMPMAT1(m,n) );  
    190.                     if (ORI(m,n)==CV_PI)  
    191.                         ORI(m,n)=-CV_PI;  
    192.                 }  
    193.                 ((mag_pyr[i].Octave)[j-1]).Level=Mag;  
    194.                 ((grad_pyr[i].Octave)[j-1]).Level=Ori;  
    195.                 cvReleaseMat(&tempMat1);  
    196.                 cvReleaseMat(&tempMat2);  
    197.         }//for levels  
    198.     }//for octaves  
    199. }

    SIFT算法第四步

    SIFT算法第五步

    SIFT算法第五步：抽取各个特征点处的特征描述字，确定特征点的描述字。描述字是Patch网格内梯度方向的描述，旋转网格到主方向，插值得到网格处梯度值。

    一个特征点可以用228=32维的向量，也可以用448=128维的向量更精确的进行描述。

    1. void ExtractFeatureDescriptors(int numoctaves, ImageOctaves *GaussianPyr)  
    2. {  
    3.     // The orientation histograms have 8 bins  
    4.     float orient_bin_spacing = PI/4;  
    5.     float orient_angles[8]={-PI,-PI+orient_bin_spacing,-PI*0.5, -orient_bin_spacing,  
    6.         0.0, orient_bin_spacing, PI*0.5,  PI+orient_bin_spacing};  
    7.     //产生描述字中心各点坐标  
    8.     float *feat_grid=(float *) malloc( 2*16 * sizeof(float));  
    9.     for (int i=0;i<GridSpacing;i++)  
    10.     {  
    11.         for (int j=0;j<2*GridSpacing;++j,++j)  
    12.         {  
    13.             feat_grid[i*2*GridSpacing+j]=-6.0+i*GridSpacing;  
    14.             feat_grid[i*2*GridSpacing+j+1]=-6.0+0.5*j*GridSpacing;  
    15.         }  
    16.     }  
    17.     //产生网格  
    18.     float *feat_samples=(float *) malloc( 2*256 * sizeof(float));  
    19.     for ( i=0;i<4*GridSpacing;i++)  
    20.     {  
    21.         for (int j=0;j<8*GridSpacing;j+=2)  
    22.         {  
    23.             feat_samples[i*8*GridSpacing+j]=-(2*GridSpacing-0.5)+i;  
    24.             feat_samples[i*8*GridSpacing+j+1]=-(2*GridSpacing-0.5)+0.5*j;  
    25.         }  
    26.     }  
    27.     float feat_window = 2*GridSpacing;  
    28.     Keypoint p = keyDescriptors; // p指向第一个结点  
    29.     while(p) // 没到表尾  
    30.     {  
    31.         float scale=(GaussianPyr[p->octave].Octave)[p->level].absolute_sigma;  
    33.         float sine = sin(p->ori);  
    34.         float cosine = cos(p->ori);    
    35.         //计算中心点坐标旋转之后的位置  
    36.         float *featcenter=(float *) malloc( 2*16 * sizeof(float));  
    37.         for (int i=0;i<GridSpacing;i++)  
    38.         {  
    39.             for (int j=0;j<2*GridSpacing;j+=2)  
    40.             {  
    41.                 float x=feat_grid[i*2*GridSpacing+j];  
    42.                 featcenter[i*2*GridSpacing+j]=((cosine * x + sine * y) + p->sx);  
    44.             }  
    45.         }  
    46.         // calculate sample window coordinates (rotated along keypoint)  
    47.         float *feat=(float *) malloc( 2*256 * sizeof(float));  
    48.         for ( i=0;i<64*GridSpacing;i++,i++)  
    49.         {  
    50.             float x=feat_samples[i];  
    51.             float y=feat_samples[i+1];  
    52.             feat[i]=((cosine * x + sine * y) + p->sx);  
    53.             feat[i+1]=((-sine * x + cosine * y) + p->sy);  
    54.         }  
    55.         //Initialize the feature descriptor.  
    56.         float *feat_desc = (float *) malloc( 128 * sizeof(float));  
    57.         for (i=0;i<128;i++)  
    58.         {  
    59.             feat_desc[i]=0.0;  
    60.             // printf("%f  ",feat_desc[i]);    
    61.         }  
    62.         //printf("/n");  
    63.         for ( i=0;i<512;++i,++i)  
    64.         {  
    65.             float x_sample = feat[i];  
    66.             float y_sample = feat[i+1];  
    67.             // Interpolate the gradient at the sample position  
    68.             /* 
    69.             0   1   0 
    70.             1   *   1 
    71.             0   1   0   具体插值策略如图示 
    72.             */  
    73.             float sample12=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample, y_sample-1);  
    74.             float sample21=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample-1, y_sample);   
    75.             float sample22=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample, y_sample);   
    76.             float sample23=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample+1, y_sample);   
    77.             float sample32=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample, y_sample+1);   
    78.             //float diff_x = 0.5*(sample23 - sample21);  
    79.             //float diff_y = 0.5*(sample32 - sample12);  
    80.             float diff_x = sample23 - sample21;  
    81.             float diff_y = sample32 - sample12;  
    82.             float mag_sample = sqrt( diff_x*diff_x + diff_y*diff_y );  
    83.             float grad_sample = atan( diff_y / diff_x );  
    84.             if(grad_sample == CV_PI)  
    85.                 grad_sample = -CV_PI;  
    86.             // Compute the weighting for the x and y dimensions.  
    87.             float *x_wght=(float *) malloc( GridSpacing * GridSpacing * sizeof(float));  
    88.             float *y_wght=(float *) malloc( GridSpacing * GridSpacing * sizeof(float));  
    89.             float *pos_wght=(float *) malloc( 8*GridSpacing * GridSpacing * sizeof(float));;  
    90.             for (int m=0;m<32;++m,++m)  
    91.             {  
    92.                 float x=featcenter[m];  
    93.                 float y=featcenter[m+1];  
    94.                 x_wght[m/2] = max(1 - (fabs(x - x_sample)*1.0/GridSpacing), 0);  
    95.                 y_wght[m/2] = max(1 - (fabs(y - y_sample)*1.0/GridSpacing), 0);   
    97.             }  
    98.             for ( m=0;m<16;++m)  
    99.                 for (int n=0;n<8;++n)  
    100.                     pos_wght[m*8+n]=x_wght[m]*y_wght[m];  
    101.             free(x_wght);  
    102.             free(y_wght);  
    103.             //计算方向的加权，首先旋转梯度场到主方向，然后计算差异   
    104.             float diff[8],orient_wght[128];  
    105.             for ( m=0;m<8;++m)  
    106.             {   
    107.                 float angle = grad_sample-(p->ori)-orient_angles[m]+CV_PI;  
    108.                 float temp = angle / (2.0 * CV_PI);  
    109.                 angle -= (int)(temp) * (2.0 * CV_PI);  
    110.                 diff[m]= angle - CV_PI;  
    111.             }  
    112.             // Compute the gaussian weighting.  
    113.             float x=p->sx;  
    114.             float y=p->sy;  
    115.             float g = exp(-((x_sample-x)*(x_sample-x)+(y_sample-y)*(y_sample-y))/(2*feat_window*feat_window))/(2*CV_PI*feat_window*feat_window);  
    117.             for ( m=0;m<128;++m)  
    118.             {  
    119.                 orient_wght[m] = max((1.0 - 1.0*fabs(diff[m%8])/orient_bin_spacing),0);  
    120.                 feat_desc[m] = feat_desc[m] + orient_wght[m]*pos_wght[m]*g*mag_sample;  
    121.             }  
    122.             free(pos_wght);     
    123.         }  
    124.         free(feat);  
    125.         free(featcenter);  
    126.         float norm=GetVecNorm( feat_desc, 128);  
    127.         for (int m=0;m<128;m++)  
    128.         {  
    129.             feat_desc[m]/=norm;  
    130.             if (feat_desc[m]>0.2)  
    131.                 feat_desc[m]=0.2;  
    132.         }  
    133.         norm=GetVecNorm( feat_desc, 128);  
    134.         for ( m=0;m<128;m++)  
    135.         {  
    136.             feat_desc[m]/=norm;  
    137.             printf("%f  ",feat_desc[m]);    
    138.         }  
    139.         printf("/n");  
    140.         p->descrip = feat_desc;  
    141.         p=p->next;  
    142.     }  
    143.     free(feat_grid);  
    144.     free(feat_samples);  
    145. }  
    147. //为了显示图象金字塔，而作的图像水平拼接  
    148. CvMat* MosaicHorizen( CvMat* im1, CvMat* im2 )  
    149. {  
    150.     int row,col;  
    151.     CvMat *mosaic = cvCreateMat( max(im1->rows,im2->rows),(im1->cols+im2->cols),CV_32FC1);  
    152. #define Mosaic(ROW,COL) ((float*)(mosaic->data.fl + mosaic->step/sizeof(float)*(ROW)))[(COL)]  
    153. #define Im11Mat(ROW,COL) ((float *)(im1->data.fl + im1->step/sizeof(float) *(ROW)))[(COL)]  
    154. #define Im22Mat(ROW,COL) ((float *)(im2->data.fl + im2->step/sizeof(float) *(ROW)))[(COL)]  
    155.     cvZero(mosaic);  
    156.     /* Copy images into mosaic1. */  
    157.     for ( row = 0; row < im1->rows; row++)  
    158.         for ( col = 0; col < im1->cols; col++)  
    159.             Mosaic(row,col)=Im11Mat(row,col) ;  
    160.     for (  row = 0; row < im2->rows; row++)  
    161.         for (  col = 0; col < im2->cols; col++)  
    162.             Mosaic(row, (col+im1->cols) )= Im22Mat(row,col) ;  
    163.     return mosaic;  
    164. }  
    166. //为了显示图象金字塔，而作的图像垂直拼接  
    167. CvMat* MosaicVertical( CvMat* im1, CvMat* im2 )  
    168. {  
    169.     int row,col;  
    170.     CvMat *mosaic = cvCreateMat(im1->rows+im2->rows,max(im1->cols,im2->cols), CV_32FC1);  
    171. #define Mosaic(ROW,COL) ((float*)(mosaic->data.fl + mosaic->step/sizeof(float)*(ROW)))[(COL)]  
    172. #define Im11Mat(ROW,COL) ((float *)(im1->data.fl + im1->step/sizeof(float) *(ROW)))[(COL)]  
    173. #define Im22Mat(ROW,COL) ((float *)(im2->data.fl + im2->step/sizeof(float) *(ROW)))[(COL)]  
    174.     cvZero(mosaic);  
    176.     /* Copy images into mosaic1. */  
    177.     for ( row = 0; row < im1->rows; row++)  
    178.         for ( col = 0; col < im1->cols; col++)  
    179.             Mosaic(row,col)= Im11Mat(row,col) ;  
    180.     for ( row = 0; row < im2->rows; row++)  
    181.         for ( col = 0; col < im2->cols; col++)  
    182.             Mosaic((row+im1->rows),col)=Im22Mat(row,col) ;

    最后，再看一下，运行效果（图中美女为老乡+朋友，何姐08年照）：

    教你一步一步用c语言实现sift算法、下 - 图1

    教你一步一步用c语言实现sift算法、下 - 图2

    完。

    updated

    有很多朋友都在本文评论下要求要本程序的完整源码包（注：本文代码未贴全，复制粘贴编译肯定诸多错误），但由于时隔太久，这份代码我自己也找不到了，不过，我可以提供一份sift + KD + BBF，且可以编译正确的代码供大家参考学习，有pudn帐号的朋友可以前去下载：（没有pudn账号的同学请加群：169056165，验证信息：sift，至群共享下载），然后用两幅不同的图片做了下匹配（当然，运行结果显示是不匹配的），效果还不错：http://weibo.com/1580904460/yDmzAEwcV#1348475194313! July、二零一二年十月十一日。