• 教你一步一步用c语言实现sift算法、下
    • 函数编写
    • 五个步骤
      • SIFT算法第一步
      • SIFT算法第二步
      • SIFT算法第三步
      • SIFT算法第四步
      • SIFT算法第五步

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

    本文接上，教你一步一步用c语言实现sift算法、上而来：

    函数编写

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

    1. //下采样原来的图像，返回缩小2倍尺寸的图像  
    2. CvMat * halfSizeImage(CvMat * im)   
    3. {  
    4.     unsigned int i,j;  
    5.     int w = im->cols/2;  
    6.     int h = im->rows/2;   
    7.     CvMat *imnew = cvCreateMat(h, w, CV_32FC1);  
    9. #define Im(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]  
    10. #define Imnew(ROW,COL) ((float *)(imnew->data.fl + imnew->step/sizeof(float) *(ROW)))[(COL)]  
    11.     for ( j = 0; j < h; j++)   
    12.         for ( i = 0; i < w; i++)   
    13.             Imnew(j,i)=Im(j*2, i*2);  
    14.     return imnew;  
    15. }  
    17. //上采样原来的图像，返回放大2倍尺寸的图像  
    18. CvMat * doubleSizeImage(CvMat * im)   
    19. {  
    20.     unsigned int i,j;  
    21.     int w = im->cols*2;  
    22.     int h = im->rows*2;   
    23.     CvMat *imnew = cvCreateMat(h, w, CV_32FC1);  
    25. #define Im(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]  
    26. #define Imnew(ROW,COL) ((float *)(imnew->data.fl + imnew->step/sizeof(float) *(ROW)))[(COL)]  
    28.     for ( j = 0; j < h; j++)   
    29.         for ( i = 0; i < w; i++)   
    30.             Imnew(j,i)=Im(j/2, i/2);  
    32.     return imnew;  
    33. }  
    35. //上采样原来的图像，返回放大2倍尺寸的线性插值图像  
    36. CvMat * doubleSizeImage2(CvMat * im)   
    37. {  
    38.     unsigned int i,j;  
    39.     int w = im->cols*2;  
    40.     int h = im->rows*2;   
    41.     CvMat *imnew = cvCreateMat(h, w, CV_32FC1);  
    43. #define Im(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]  
    44. #define Imnew(ROW,COL) ((float *)(imnew->data.fl + imnew->step/sizeof(float) *(ROW)))[(COL)]  
    46.     // fill every pixel so we don't have to worry about skipping pixels later  
    47.     for ( j = 0; j < h; j++)   
    48.     {  
    49.         for ( i = 0; i < w; i++)   
    50.         {  
    51.             Imnew(j,i)=Im(j/2, i/2);  
    52.         }  
    53.     }  
    54.     /* 
    55.     A B C 
    56.     E F G 
    57.     H I J 
    58.     pixels A C H J are pixels from original image 
    59.     pixels B E G I F are interpolated pixels 
    60.     */  
    61.     // interpolate pixels B and I  
    62.     for ( j = 0; j < h; j += 2)  
    63.         for ( i = 1; i < w - 1; i += 2)  
    64.             Imnew(j,i)=0.5*(Im(j/2, i/2)+Im(j/2, i/2+1));  
    65.     // interpolate pixels E and G  
    66.     for ( j = 1; j < h - 1; j += 2)  
    67.         for ( i = 0; i < w; i += 2)  
    68.             Imnew(j,i)=0.5*(Im(j/2, i/2)+Im(j/2+1, i/2));  
    69.     // interpolate pixel F  
    70.     for ( j = 1; j < h - 1; j += 2)  
    71.         for ( i = 1; i < w - 1; i += 2)  
    72.             Imnew(j,i)=0.25*(Im(j/2, i/2)+Im(j/2+1, i/2)+Im(j/2, i/2+1)+Im(j/2+1, i/2+1));  
    73.     return imnew;  
    74. }  
    76. //双线性插值，返回像素间的灰度值  
    77. float getPixelBI(CvMat * im, float col, float row)   
    78. {  
    79.     int irow, icol;  
    80.     float rfrac, cfrac;  
    81.     float row1 = 0, row2 = 0;  
    82.     int width=im->cols;  
    83.     int height=im->rows;  
    84. #define ImMat(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]  
    86.     irow = (int) row;  
    87.     icol = (int) col;  
    89.     if (irow < 0 || irow >= height  
    90.         || icol < 0 || icol >= width)  
    91.         return 0;  
    92.     if (row > height - 1)  
    93.         row = height - 1;  
    94.     if (col > width - 1)  
    95.         col = width - 1;  
    96.     rfrac = 1.0 - (row - (float) irow);  
    97.     cfrac = 1.0 - (col - (float) icol);  
    98.     if (cfrac < 1)   
    99.     {  
    100.         row1 = cfrac * ImMat(irow,icol) + (1.0 - cfrac) * ImMat(irow,icol+1);  
    101.     }   
    102.     else   
    103.     {  
    104.         row1 = ImMat(irow,icol);  
    105.     }  
    106.     if (rfrac < 1)   
    107.     {  
    108.         if (cfrac < 1)   
    109.         {  
    110.             row2 = cfrac * ImMat(irow+1,icol) + (1.0 - cfrac) * ImMat(irow+1,icol+1);  
    111.         } else   
    112.         {  
    113.             row2 = ImMat(irow+1,icol);  
    114.         }  
    115.     }  
    116.     return rfrac * row1 + (1.0 - rfrac) * row2;  
    117. }  
    119. //矩阵归一化  
    120. void normalizeMat(CvMat* mat)   
    121. {  
    122. #define Mat(ROW,COL) ((float *)(mat->data.fl + mat->step/sizeof(float) *(ROW)))[(COL)]  
    123.     float sum = 0;  
    125.     for (unsigned int j = 0; j < mat->rows; j++)   
    126.         for (unsigned int i = 0; i < mat->cols; i++)   
    127.             sum += Mat(j,i);  
    128.     for ( j = 0; j < mat->rows; j++)   
    129.         for (unsigned int i = 0; i < mat->rows; i++)   
    130.             Mat(j,i) /= sum;  
    131. }  
    133. //向量归一化  
    134. void normalizeVec(float* vec, int dim)   
    135. {  
    136.     unsigned int i;  
    137.     float sum = 0;  
    138.     for ( i = 0; i < dim; i++)  
    139.         sum += vec[i];  
    140.     for ( i = 0; i < dim; i++)  
    141.         vec[i] /= sum;  
    142. }  
    144. //得到向量的欧式长度，2-范数  
    145. float GetVecNorm( float* vec, int dim )  
    146. {  
    147.     float sum=0.0;  
    148.     for (unsigned int i=0;i<dim;i++)  
    149.         sum+=vec[i]*vec[i];  
    150.     return sqrt(sum);  
    151. }  
    153. //产生1D高斯核  
    154. float* GaussianKernel1D(float sigma, int dim)   
    155. {  
    157.     unsigned int i;  
    158.     //printf("GaussianKernel1D(): Creating 1x%d vector for sigma=%.3f gaussian kernel/n", dim, sigma);  
    160.     float *kern=(float*)malloc( dim*sizeof(float) );  
    161.     float s2 = sigma * sigma;  
    162.     int c = dim / 2;  
    163.     float m= 1.0/(sqrt(2.0 * CV_PI) * sigma);  
    164.     double v;   
    165.     for ( i = 0; i < (dim + 1) / 2; i++)   
    166.     {  
    167.         v = m * exp(-(1.0*i*i)/(2.0 * s2)) ;  
    168.         kern[c+i] = v;  
    169.         kern[c-i] = v;  
    170.     }  
    171.     //   normalizeVec(kern, dim);  
    172.     // for ( i = 0; i < dim; i++)  
    173.     //  printf("%f  ", kern[i]);  
    174.     //  printf("/n");  
    175.     return kern;  
    176. }  
    178. //产生2D高斯核矩阵  
    179. CvMat* GaussianKernel2D(float sigma)   
    180. {  
    181.     // int dim = (int) max(3.0f, GAUSSKERN * sigma);  
    182.     int dim = (int) max(3.0f, 2.0 * GAUSSKERN *sigma + 1.0f);  
    183.     // make dim odd  
    184.     if (dim % 2 == 0)  
    185.         dim++;  
    186.     //printf("GaussianKernel(): Creating %dx%d matrix for sigma=%.3f gaussian/n", dim, dim, sigma);  
    187.     CvMat* mat=cvCreateMat(dim, dim, CV_32FC1);  
    188. #define Mat(ROW,COL) ((float *)(mat->data.fl + mat->step/sizeof(float) *(ROW)))[(COL)]  
    189.     float s2 = sigma * sigma;  
    190.     int c = dim / 2;  
    191.     //printf("%d %d/n", mat.size(), mat[0].size());  
    192.     float m= 1.0/(sqrt(2.0 * CV_PI) * sigma);  
    193.     for (int i = 0; i < (dim + 1) / 2; i++)   
    194.     {  
    195.         for (int j = 0; j < (dim + 1) / 2; j++)   
    196.         {  
    197.             //printf("%d %d %d/n", c, i, j);  
    198.             float v = m * exp(-(1.0*i*i + 1.0*j*j) / (2.0 * s2));  
    199.             Mat(c+i,c+j) =v;  
    200.             Mat(c-i,c+j) =v;  
    201.             Mat(c+i,c-j) =v;  
    202.             Mat(c-i,c-j) =v;  
    203.         }  
    204.     }  
    205.     // normalizeMat(mat);  
    206.     return mat;  
    207. }  
    209. //x方向像素处作卷积  
    210. float ConvolveLocWidth(float* kernel, int dim, CvMat * src, int x, int y)   
    211. {  
    212. #define Src(ROW,COL) ((float *)(src->data.fl + src->step/sizeof(float) *(ROW)))[(COL)]  
    213.     unsigned int i;  
    214.     float pixel = 0;  
    215.     int col;  
    216.     int cen = dim / 2;  
    217.     //printf("ConvolveLoc(): Applying convoluation at location (%d, %d)/n", x, y);  
    218.     for ( i = 0; i < dim; i++)   
    219.     {  
    220.         col = x + (i - cen);  
    221.         if (col < 0)  
    222.             col = 0;  
    223.         if (col >= src->cols)  
    224.             col = src->cols - 1;  
    225.         pixel += kernel[i] * Src(y,col);  
    226.     }  
    227.     if (pixel > 1)  
    228.         pixel = 1;  
    229.     return pixel;  
    230. }  
    232. //x方向作卷积  
    233. void Convolve1DWidth(float* kern, int dim, CvMat * src, CvMat * dst)   
    234. {  
    235. #define DST(ROW,COL) ((float *)(dst->data.fl + dst->step/sizeof(float) *(ROW)))[(COL)]  
    236.     unsigned int i,j;  
    238.     for ( j = 0; j < src->rows; j++)   
    239.     {  
    240.         for ( i = 0; i < src->cols; i++)   
    241.         {  
    242.             //printf("%d, %d/n", i, j);  
    243.             DST(j,i) = ConvolveLocWidth(kern, dim, src, i, j);  
    244.         }  
    245.     }  
    246. }  
    248. //y方向像素处作卷积  
    249. float ConvolveLocHeight(float* kernel, int dim, CvMat * src, int x, int y)   
    250. {  
    251. #define Src(ROW,COL) ((float *)(src->data.fl + src->step/sizeof(float) *(ROW)))[(COL)]  
    252.     unsigned int j;  
    253.     float pixel = 0;  
    254.     int cen = dim / 2;  
    255.     //printf("ConvolveLoc(): Applying convoluation at location (%d, %d)/n", x, y);  
    256.     for ( j = 0; j < dim; j++)   
    257.     {  
    258.         int row = y + (j - cen);  
    259.         if (row < 0)  
    260.             row = 0;  
    261.         if (row >= src->rows)  
    262.             row = src->rows - 1;  
    263.         pixel += kernel[j] * Src(row,x);  
    264.     }  
    265.     if (pixel > 1)  
    266.         pixel = 1;  
    267.     return pixel;  
    268. }  
    270. //y方向作卷积  
    271. void Convolve1DHeight(float* kern, int dim, CvMat * src, CvMat * dst)   
    272. {  
    273. #define Dst(ROW,COL) ((float *)(dst->data.fl + dst->step/sizeof(float) *(ROW)))[(COL)]  
    274.     unsigned int i,j;  
    275.     for ( j = 0; j < src->rows; j++)   
    276.     {  
    277.         for ( i = 0; i < src->cols; i++)   
    278.         {  
    279.             //printf("%d, %d/n", i, j);  
    280.             Dst(j,i) = ConvolveLocHeight(kern, dim, src, i, j);  
    281.         }  
    282.     }  
    283. }  
    285. //卷积模糊图像  
    286. int BlurImage(CvMat * src, CvMat * dst, float sigma)   
    287. {  
    288.     float* convkernel;  
    289.     int dim = (int) max(3.0f, 2.0 * GAUSSKERN * sigma + 1.0f);  
    290.     CvMat *tempMat;  
    291.     // make dim odd  
    292.     if (dim % 2 == 0)  
    293.         dim++;  
    294.     tempMat = cvCreateMat(src->rows, src->cols, CV_32FC1);  
    295.     convkernel = GaussianKernel1D(sigma, dim);  
    297.     Convolve1DWidth(convkernel, dim, src, tempMat);  
    298.     Convolve1DHeight(convkernel, dim, tempMat, dst);  
    299.     cvReleaseMat(&tempMat);  
    300.     return dim;  
    301. }

    五个步骤

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

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

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

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

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

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

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

    int DetectKeypoint(int numoctaves, ImageOctaves *GaussianPyr);

    4、 SIFT算法第四步：计算高斯图像的梯度方向和幅值，计算各个特征点的主方向

    void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr);

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

    void ExtractFeatureDescriptors(int numoctaves, ImageOctaves *GaussianPyr);

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

    SIFT算法第一步

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

    1. CvMat *ScaleInitImage(CvMat * im)   
    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决定，给定一个基本图像，计算它的高斯金字塔图像，返回外部向量是阶梯指针，内部向量是每一个阶梯内部的不同尺度图像。

    1. //SIFT算法第二步  
    2. ImageOctaves* BuildGaussianOctaves(CvMat * image)   
    3. {  
    4.     ImageOctaves *octaves;  
    5.     CvMat *tempMat;  
    6.     CvMat *dst;  
    7.     CvMat *temp;  
    9.     int i,j;  
    10.     double k = pow(2, 1.0/((float)SCALESPEROCTAVE));  //方差倍数  
    11.     float preblur_sigma, initial_sigma , sigma1,sigma2,sigma,absolute_sigma,sigma_f;  
    12.     //计算金字塔的阶梯数目  
    13.     int dim = min(image->rows, image->cols);  
    14.     int numoctaves = (int) (log((double) dim) / log(2.0)) - 2;    //金字塔阶数  
    15.     //限定金字塔的阶梯数  
    16.     numoctaves = min(numoctaves, MAXOCTAVES);  
    17.     //为高斯金塔和DOG金字塔分配内存  
    18.     octaves=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );  
    19.     DOGoctaves=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );  
    21.     printf("BuildGaussianOctaves(): Base image dimension is %dx%d/n", (int)(0.5*(image->cols)), (int)(0.5*(image->rows)) );  
    22.     printf("BuildGaussianOctaves(): Building %d octaves/n", numoctaves);  
    24.     // start with initial source image  
    25.     tempMat=cvCloneMat( image );  
    26.     // preblur_sigma = 1.0;//sqrt(2 - 4*INITSIGMA*INITSIGMA);  
    27.     initial_sigma = sqrt(2);//sqrt( (4*INITSIGMA*INITSIGMA) + preblur_sigma * preblur_sigma );  
    28.     //   initial_sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA * 4);  
    30.     //在每一阶金字塔图像中建立不同的尺度图像  
    31.     for ( i = 0; i < numoctaves; i++)   
    32.     {     
    33.         //首先建立金字塔每一阶梯的最底层，其中0阶梯的最底层已经建立好  
    34.         printf("Building octave %d of dimesion (%d, %d)/n", i, tempMat->cols,tempMat->rows);  
    35.         //为各个阶梯分配内存  
    36.         octaves[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE + 3) * sizeof(ImageLevels) );  
    37.         DOGoctaves[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE + 2) * sizeof(ImageLevels) );  
    38.         //存储各个阶梯的最底层  
    39.         (octaves[i].Octave)[0].Level=tempMat;  
    41.         octaves[i].col=tempMat->cols;  
    42.         octaves[i].row=tempMat->rows;  
    43.         DOGoctaves[i].col=tempMat->cols;  
    44.         DOGoctaves[i].row=tempMat->rows;  
    45.         if (DOUBLE_BASE_IMAGE_SIZE)  
    46.             octaves[i].subsample=pow(2,i)*0.5;  
    47.         else  
    48.             octaves[i].subsample=pow(2,i);  
    50.         if(i==0)       
    51.         {  
    52.             (octaves[0].Octave)[0].levelsigma = initial_sigma;  
    53.             (octaves[0].Octave)[0].absolute_sigma = initial_sigma;  
    54.             printf("0 scale and blur sigma : %f /n", (octaves[0].subsample) * ((octaves[0].Octave)[0].absolute_sigma));  
    55.         }  
    56.         else  
    57.         {  
    58.             (octaves[i].Octave)[0].levelsigma = (octaves[i-1].Octave)[SCALESPEROCTAVE].levelsigma;  
    59.             (octaves[i].Octave)[0].absolute_sigma = (octaves[i-1].Octave)[SCALESPEROCTAVE].absolute_sigma;  
    60.             printf( "0 scale and blur sigma : %f /n", ((octaves[i].Octave)[0].absolute_sigma) );  
    61.         }  
    62.         sigma = initial_sigma;  
    63.         //建立本阶梯其他层的图像  
    64.         for ( j =  1; j < SCALESPEROCTAVE + 3; j++)   
    65.         {  
    66.             dst = cvCreateMat(tempMat->rows, tempMat->cols, CV_32FC1);//用于存储高斯层  
    67.             temp = cvCreateMat(tempMat->rows, tempMat->cols, CV_32FC1);//用于存储DOG层  
    68.             // 2 passes of 1D on original  
    69.             //   if(i!=0)  
    70.             //   {  
    71.             //       sigma1 = pow(k, j - 1) * ((octaves[i-1].Octave)[j-1].levelsigma);  
    72.             //          sigma2 = pow(k, j) * ((octaves[i].Octave)[j-1].levelsigma);  
    73.             //       sigma = sqrt(sigma2*sigma2 - sigma1*sigma1);  
    74.             sigma_f= sqrt(k*k-1)*sigma;  
    75.             //   }  
    76.             //   else  
    77.             //   {  
    78.             //       sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA * 4)*pow(k,j);  
    79.             //   }    
    80.             sigma = k*sigma;  
    81.             absolute_sigma = sigma * (octaves[i].subsample);  
    82.             printf("%d scale and Blur sigma: %f  /n", j, absolute_sigma);  
    84.             (octaves[i].Octave)[j].levelsigma = sigma;  
    85.             (octaves[i].Octave)[j].absolute_sigma = absolute_sigma;  
    86.             //产生高斯层  
    87.             int length=BlurImage((octaves[i].Octave)[j-1].Level, dst, sigma_f);//相应尺度  
    88.             (octaves[i].Octave)[j].levelsigmalength = length;  
    89.             (octaves[i].Octave)[j].Level=dst;  
    90.             //产生DOG层  
    91.             cvSub( ((octaves[i].Octave)[j]).Level, ((octaves[i].Octave)[j-1]).Level, temp, 0 );  
    92.             //         cvAbsDiff( ((octaves[i].Octave)[j]).Level, ((octaves[i].Octave)[j-1]).Level, temp );  
    93.             ((DOGoctaves[i].Octave)[j-1]).Level=temp;  
    94.         }  
    95.         // halve the image size for next iteration  
    96.         tempMat  = halfSizeImage( ( (octaves[i].Octave)[SCALESPEROCTAVE].Level ) );  
    97.     }  
    98.     return octaves;  
    99. }

    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.                                     //Tr_H = Dxx + Dyy;  
    99.                                     //Det_H = Dxx*Dyy - Dxy^2;  
    100.                                     float Dxx,Dyy,Dxy,Tr_H,Det_H,curvature_ratio;  
    101.                                     Dxx = ImLevels(i,j,m,n-1) + ImLevels(i,j,m,n+1)-2.0*ImLevels(i,j,m,n);  
    102.                                     Dyy = ImLevels(i,j,m-1,n) + ImLevels(i,j,m+1,n)-2.0*ImLevels(i,j,m,n);  
    103.                                     Dxy = ImLevels(i,j,m-1,n-1) + ImLevels(i,j,m+1,n+1) - ImLevels(i,j,m+1,n-1) - ImLevels(i,j,m-1,n+1);  
    104.                                     Tr_H = Dxx + Dyy;  
    105.                                     Det_H = Dxx*Dyy - Dxy*Dxy;  
    106.                                     // Compute the ratio of the principal curvatures.  
    107.                                     curvature_ratio = (1.0*Tr_H*Tr_H)/Det_H;  
    108.                                     if ( (Det_H>=0.0) && (curvature_ratio <= curvature_threshold) )  //最后得到最具有显著性特征的特征点  
    109.                                     {  
    110.                                         //将其存储起来，以计算后面的特征描述字  
    111.                                         keypoint_count++;  
    112.                                         Keypoint k;  
    113.                                         /* Allocate memory for the keypoint. */  
    114.                                         k = (Keypoint) malloc(sizeof(struct KeypointSt));  
    115.                                         k->next = keypoints;  
    116.                                         keypoints = k;  
    117.                                         k->row = m*(GaussianPyr[i].subsample);  
    118.                                         k->col =n*(GaussianPyr[i].subsample);  
    119.                                         k->sy = m;    //行  
    120.                                         k->sx = n;    //列  
    121.                                         k->octave=i;  
    122.                                         k->level=j;  
    123.                                         k->scale = (GaussianPyr[i].Octave)[j].absolute_sigma;        
    124.                                     }//if >curvature_thresh  
    125.                                 }//if >contrast  
    126.                             }//if inf value  
    127.                         }//if non zero  
    128.                     }//if >contrast  
    129.                 }  //for concrete image level col  
    130.         }//for levels  
    131.     }//for octaves  
    132.     return keypoint_count;  
    133. }  
    135. //在图像中，显示SIFT特征点的位置  
    136. void DisplayKeypointLocation(IplImage* image, ImageOctaves *GaussianPyr)  
    137. {  
    139.     Keypoint p = keypoints; // p指向第一个结点  
    140.     while(p) // 没到表尾  
    141.     {     
    142.         cvLine( image, cvPoint((int)((p->col)-3),(int)(p->row)),   
    143.             cvPoint((int)((p->col)+3),(int)(p->row)), CV_RGB(255,255,0),  
    144.             1, 8, 0 );  
    145.         cvLine( image, cvPoint((int)(p->col),(int)((p->row)-3)),   
    146.             cvPoint((int)(p->col),(int)((p->row)+3)), CV_RGB(255,255,0),  
    147.             1, 8, 0 );  
    148.         //  cvCircle(image,cvPoint((uchar)(p->col),(uchar)(p->row)),  
    149.         //   (int)((GaussianPyr[p->octave].Octave)[p->level].absolute_sigma),  
    150.         //   CV_RGB(255,0,0),1,8,0);  
    151.         p=p->next;  
    152.     }   
    153. }  
    155. // Compute the gradient direction and magnitude of the gaussian pyramid images  
    156. void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr)  
    157. {  
    158.     // ImageOctaves *mag_thresh ;  
    159.     mag_pyr=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );  
    160.     grad_pyr=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );  
    161.     // float sigma=( (GaussianPyr[0].Octave)[SCALESPEROCTAVE+2].absolute_sigma ) / GaussianPyr[0].subsample;  
    162.     // int dim = (int) (max(3.0f, 2 * GAUSSKERN *sigma + 1.0f)*0.5+0.5);  
    163. #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(GaussianPyr[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + GaussianPyr[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]  
    164.     for (int i=0; i<numoctaves; i++)    
    165.     {          
    166.         mag_pyr[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE) * sizeof(ImageLevels) );  
    167.         grad_pyr[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE) * sizeof(ImageLevels) );  
    168.         for(int j=1;j<SCALESPEROCTAVE+1;j++)//取中间的scaleperoctave个层  
    169.         {    
    170.             CvMat *Mag = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);  
    171.             CvMat *Ori = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);  
    172.             CvMat *tempMat1 = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);  
    173.             CvMat *tempMat2 = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);  
    174.             cvZero(Mag);  
    175.             cvZero(Ori);  
    176.             cvZero(tempMat1);  
    177.             cvZero(tempMat2);   
    178. #define MAG(ROW,COL) ((float *)(Mag->data.fl + Mag->step/sizeof(float) *(ROW)))[(COL)]     
    179. #define ORI(ROW,COL) ((float *)(Ori->data.fl + Ori->step/sizeof(float) *(ROW)))[(COL)]    
    180. #define TEMPMAT1(ROW,COL) ((float *)(tempMat1->data.fl + tempMat1->step/sizeof(float) *(ROW)))[(COL)]  
    181. #define TEMPMAT2(ROW,COL) ((float *)(tempMat2->data.fl + tempMat2->step/sizeof(float) *(ROW)))[(COL)]  
    182.             for (int m=1;m<(GaussianPyr[i].row-1);m++)   
    183.                 for(int n=1;n<(GaussianPyr[i].col-1);n++)  
    184.                 {  
    185.                     //计算幅值  
    186.                     TEMPMAT1(m,n) = 0.5*( ImLevels(i,j,m,n+1)-ImLevels(i,j,m,n-1) );  //dx  
    187.                     TEMPMAT2(m,n) = 0.5*( ImLevels(i,j,m+1,n)-ImLevels(i,j,m-1,n) );  //dy  
    188.                     MAG(m,n) = sqrt(TEMPMAT1(m,n)*TEMPMAT1(m,n)+TEMPMAT2(m,n)*TEMPMAT2(m,n));  //mag  
    189.                     //计算方向  
    190.                     ORI(m,n) =atan( TEMPMAT2(m,n)/TEMPMAT1(m,n) );  
    191.                     if (ORI(m,n)==CV_PI)  
    192.                         ORI(m,n)=-CV_PI;  
    193.                 }  
    194.                 ((mag_pyr[i].Octave)[j-1]).Level=Mag;  
    195.                 ((grad_pyr[i].Octave)[j-1]).Level=Ori;  
    196.                 cvReleaseMat(&tempMat1);  
    197.                 cvReleaseMat(&tempMat2);  
    198.         }//for levels  
    199.     }//for octaves  
    200. }

    SIFT算法第四步

    1. //SIFT算法第四步：计算各个特征点的主方向，确定主方向  
    2. void AssignTheMainOrientation(int numoctaves, ImageOctaves *GaussianPyr,ImageOctaves *mag_pyr,ImageOctaves *grad_pyr)  
    3. {  
    4.     // Set up the histogram bin centers for a 36 bin histogram.  
    5.     int num_bins = 36;  
    6.     float hist_step = 2.0*PI/num_bins;  
    7.     float hist_orient[36];  
    8.     for (int i=0;i<36;i++)  
    9.         hist_orient[i]=-PI+i*hist_step;  
    10.     float sigma1=( ((GaussianPyr[0].Octave)[SCALESPEROCTAVE].absolute_sigma) ) / (GaussianPyr[0].subsample);//SCALESPEROCTAVE+2  
    11.     int zero_pad = (int) (max(3.0f, 2 * GAUSSKERN *sigma1 + 1.0f)*0.5+0.5);  
    12.     //Assign orientations to the keypoints.  
    13. #define ImLevels(OCTAVES,LEVELS,ROW,COL) ((float *)((GaussianPyr[(OCTAVES)].Octave[(LEVELS)].Level)->data.fl + (GaussianPyr[(OCTAVES)].Octave[(LEVELS)].Level)->step/sizeof(float) *(ROW)))[(COL)]  
    15.     int keypoint_count = 0;  
    16.     Keypoint p = keypoints; // p指向第一个结点  
    18.     while(p) // 没到表尾  
    19.     {  
    20.         int i=p->octave;  
    21.         int j=p->level;  
    22.         int m=p->sy;   //行  
    23.         int n=p->sx;   //列  
    24.         if ((m>=zero_pad)&&(m<GaussianPyr[i].row-zero_pad)&&  
    25.             (n>=zero_pad)&&(n<GaussianPyr[i].col-zero_pad) )  
    26.         {  
    27.             float sigma=( ((GaussianPyr[i].Octave)[j].absolute_sigma) ) / (GaussianPyr[i].subsample);  
    28.             //产生二维高斯模板  
    29.             CvMat* mat = GaussianKernel2D( sigma );           
    30.             int dim=(int)(0.5 * (mat->rows));  
    31.             //分配用于存储Patch幅值和方向的空间  
    32. #define MAT(ROW,COL) ((float *)(mat->data.fl + mat->step/sizeof(float) *(ROW)))[(COL)]  
    34.             //声明方向直方图变量  
    35.             double* orienthist = (double *) malloc(36 * sizeof(double));  
    36.             for ( int sw = 0 ; sw < 36 ; ++sw)   
    37.             {  
    38.                 orienthist[sw]=0.0;    
    39.             }  
    40.             //在特征点的周围统计梯度方向  
    41.             for (int x=m-dim,mm=0;x<=(m+dim);x++,mm++)   
    42.                 for(int y=n-dim,nn=0;y<=(n+dim);y++,nn++)  
    43.                 {       
    44.                     //计算特征点处的幅值  
    45.                     double dx = 0.5*(ImLevels(i,j,x,y+1)-ImLevels(i,j,x,y-1));  //dx  
    46.                     double dy = 0.5*(ImLevels(i,j,x+1,y)-ImLevels(i,j,x-1,y));  //dy  
    47.                     double mag = sqrt(dx*dx+dy*dy);  //mag  
    48.                     //计算方向  
    49.                     double Ori =atan( 1.0*dy/dx );  
    50.                     int binIdx = FindClosestRotationBin(36, Ori);                   //得到离现有方向最近的直方块  
    51.                     orienthist[binIdx] = orienthist[binIdx] + 1.0* mag * MAT(mm,nn);//利用高斯加权累加进直方图相应的块  
    52.                 }  
    53.                 // Find peaks in the orientation histogram using nonmax suppression.  
    54.                 AverageWeakBins (orienthist, 36);  
    55.                 // find the maximum peak in gradient orientation  
    56.                 double maxGrad = 0.0;  
    57.                 int maxBin = 0;  
    58.                 for (int b = 0 ; b < 36 ; ++b)   
    59.                 {  
    60.                     if (orienthist[b] > maxGrad)   
    61.                     {  
    62.                         maxGrad = orienthist[b];  
    63.                         maxBin = b;  
    64.                     }  
    65.                 }  
    66.                 // First determine the real interpolated peak high at the maximum bin  
    67.                 // position, which is guaranteed to be an absolute peak.  
    68.                 double maxPeakValue=0.0;  
    69.                 double maxDegreeCorrection=0.0;  
    70.                 if ( (InterpolateOrientation ( orienthist[maxBin == 0 ? (36 - 1) : (maxBin - 1)],  
    71.                     orienthist[maxBin], orienthist[(maxBin + 1) % 36],  
    72.                     &maxDegreeCorrection, &maxPeakValue)) == false)  
    73.                     printf("BUG: Parabola fitting broken");  
    75.                 // Now that we know the maximum peak value, we can find other keypoint  
    76.                 // orientations, which have to fulfill two criterias:  
    77.                 //  
    78.                 //  1. They must be a local peak themselves. Else we might add a very  
    79.                 //     similar keypoint orientation twice (imagine for example the  
    80.                 //     values: 0.4 1.0 0.8, if 1.0 is maximum peak, 0.8 is still added  
    81.                 //     with the default threshhold, but the maximum peak orientation  
    82.                 //     was already added).  
    83.                 //  2. They must have at least peakRelThresh times the maximum peak  
    84.                 //     value.  
    85.                 bool binIsKeypoint[36];  
    86.                 for ( b = 0 ; b < 36 ; ++b)   
    87.                 {  
    88.                     binIsKeypoint[b] = false;  
    89.                     // The maximum peak of course is  
    90.                     if (b == maxBin)   
    91.                     {  
    92.                         binIsKeypoint[b] = true;  
    93.                         continue;  
    94.                     }  
    95.                     // Local peaks are, too, in case they fulfill the threshhold  
    96.                     if (orienthist[b] < (peakRelThresh * maxPeakValue))  
    97.                         continue;  
    98.                     int leftI = (b == 0) ? (36 - 1) : (b - 1);  
    99.                     int rightI = (b + 1) % 36;  
    100.                     if (orienthist[b] <= orienthist[leftI] || orienthist[b] <= orienthist[rightI])  
    101.                         continue; // no local peak  
    102.                     binIsKeypoint[b] = true;  
    103.                 }  
    104.                 // find other possible locations  
    105.                 double oneBinRad = (2.0 * PI) / 36;  
    106.                 for ( b = 0 ; b < 36 ; ++b)   
    107.                 {  
    108.                     if (binIsKeypoint[b] == false)  
    109.                         continue;  
    110.                     int bLeft = (b == 0) ? (36 - 1) : (b - 1);  
    111.                     int bRight = (b + 1) % 36;  
    112.                     // Get an interpolated peak direction and value guess.  
    113.                     double peakValue;  
    114.                     double degreeCorrection;  
    116.                     double maxPeakValue, maxDegreeCorrection;                
    117.                     if (InterpolateOrientation ( orienthist[maxBin == 0 ? (36 - 1) : (maxBin - 1)],  
    118.                         orienthist[maxBin], orienthist[(maxBin + 1) % 36],  
    119.                         °reeCorrection, &peakValue) == false)  
    120.                     {  
    121.                         printf("BUG: Parabola fitting broken");  
    122.                     }  
    124.                     double degree = (b + degreeCorrection) * oneBinRad - PI;  
    125.                     if (degree < -PI)  
    126.                         degree += 2.0 * PI;  
    127.                     else if (degree > PI)  
    128.                         degree -= 2.0 * PI;  
    129.                     //存储方向，可以直接利用检测到的链表进行该步主方向的指定;  
    130.                     //分配内存重新存储特征点  
    131.                     Keypoint k;  
    132.                     /* Allocate memory for the keypoint Descriptor. */  
    133.                     k = (Keypoint) malloc(sizeof(struct KeypointSt));  
    134.                     k->next = keyDescriptors;  
    135.                     keyDescriptors = k;  
    136.                     k->descrip = (float*)malloc(LEN * sizeof(float));  
    137.                     k->row = p->row;  
    138.                     k->col = p->col;  
    139.                     k->sy = p->sy;    //行  
    140.                     k->sx = p->sx;    //列  
    141.                     k->octave = p->octave;  
    142.                     k->level = p->level;  
    143.                     k->scale = p->scale;        
    144.                     k->ori = degree;  
    145.                     k->mag = peakValue;    
    146.                 }//for  
    147.                 free(orienthist);  
    148.         }  
    149.         p=p->next;  
    150.     }   
    151. }  
    153. //寻找与方向直方图最近的柱，确定其index   
    154. int FindClosestRotationBin (int binCount, float angle)  
    155. {  
    156.     angle += CV_PI;  
    157.     angle /= 2.0 * CV_PI;  
    158.     // calculate the aligned bin  
    159.     angle *= binCount;  
    160.     int idx = (int) angle;  
    161.     if (idx == binCount)  
    162.         idx = 0;  
    163.     return (idx);  
    164. }  
    166. // Average the content of the direction bins.  
    167. void AverageWeakBins (double* hist, int binCount)  
    168. {  
    169.     // TODO: make some tests what number of passes is the best. (its clear  
    170.     // one is not enough, as we may have something like  
    171.     // ( 0.4, 0.4, 0.3, 0.4, 0.4 ))  
    172.     for (int sn = 0 ; sn < 2 ; ++sn)   
    173.     {  
    174.         double firstE = hist[0];  
    175.         double last = hist[binCount-1];  
    176.         for (int sw = 0 ; sw < binCount ; ++sw)   
    177.         {  
    178.             double cur = hist[sw];  
    179.             double next = (sw == (binCount - 1)) ? firstE : hist[(sw + 1) % binCount];  
    180.             hist[sw] = (last + cur + next) / 3.0;  
    181.             last = cur;  
    182.         }  
    183.     }  
    184. }  
    186. // Fit a parabol to the three points (-1.0 ; left), (0.0 ; middle) and  
    187. // (1.0 ; right).  
    188. // Formulas:  
    189. // f(x) = a (x - c)^2 + b  
    190. // c is the peak offset (where f'(x) is zero), b is the peak value.  
    191. // In case there is an error false is returned, otherwise a correction  
    192. // value between [-1 ; 1] is returned in 'degreeCorrection', where -1  
    193. // means the peak is located completely at the left vector, and -0.5 just  
    194. // in the middle between left and middle and > 0 to the right side. In  
    195. // 'peakValue' the maximum estimated peak value is stored.  
    196. bool InterpolateOrientation (double left, double middle,double right, double *degreeCorrection, double *peakValue)  
    197. {  
    198.     double a = ((left + right) - 2.0 * middle) / 2.0;   //抛物线捏合系数a  
    199.     // degreeCorrection = peakValue = Double.NaN;  
    201.     // Not a parabol  
    202.     if (a == 0.0)  
    203.         return false;  
    204.     double c = (((left - middle) / a) - 1.0) / 2.0;  
    205.     double b = middle - c * c * a;  
    206.     if (c < -0.5 || c > 0.5)  
    207.         return false;  
    208.     *degreeCorrection = c;  
    209.     *peakValue = b;  
    210.     return true;  
    211. }  
    213. //显示特征点处的主方向  
    214. void DisplayOrientation (IplImage* image, ImageOctaves *GaussianPyr)  
    215. {  
    216.     Keypoint p = keyDescriptors; // p指向第一个结点  
    217.     while(p) // 没到表尾  
    218.     {  
    219.         float scale=(GaussianPyr[p->octave].Octave)[p->level].absolute_sigma;  
    220.         float autoscale = 3.0;   
    221.         float uu=autoscale*scale*cos(p->ori);  
    222.         float vv=autoscale*scale*sin(p->ori);  
    223.         float x=(p->col)+uu;  
    224.         float y=(p->row)+vv;  
    225.         cvLine( image, cvPoint((int)(p->col),(int)(p->row)),   
    226.             cvPoint((int)x,(int)y), CV_RGB(255,255,0),  
    227.             1, 8, 0 );  
    228.         // Arrow head parameters  
    229.         float alpha = 0.33; // Size of arrow head relative to the length of the vector  
    230.         float beta = 0.33;  // Width of the base of the arrow head relative to the length  
    232.         float xx0= (p->col)+uu-alpha*(uu+beta*vv);  
    233.         float yy0= (p->row)+vv-alpha*(vv-beta*uu);  
    234.         float xx1= (p->col)+uu-alpha*(uu-beta*vv);  
    235.         float yy1= (p->row)+vv-alpha*(vv+beta*uu);  
    236.         cvLine( image, cvPoint((int)xx0,(int)yy0),   
    237.             cvPoint((int)x,(int)y), CV_RGB(255,255,0),  
    238.             1, 8, 0 );  
    239.         cvLine( image, cvPoint((int)xx1,(int)yy1),   
    240.             cvPoint((int)x,(int)y), CV_RGB(255,255,0),  
    241.             1, 8, 0 );  
    242.         p=p->next;  
    243.     }   
    244. }

    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.                 float y=feat_grid[i*2*GridSpacing+j+1];  
    43.                 featcenter[i*2*GridSpacing+j]=((cosine * x + sine * y) + p->sx);  
    44.                 featcenter[i*2*GridSpacing+j+1]=((-sine * x + cosine * y) + p->sy);  
    45.             }  
    46.         }  
    47.         // calculate sample window coordinates (rotated along keypoint)  
    48.         float *feat=(float *) malloc( 2*256 * sizeof(float));  
    49.         for ( i=0;i<64*GridSpacing;i++,i++)  
    50.         {  
    51.             float x=feat_samples[i];  
    52.             float y=feat_samples[i+1];  
    53.             feat[i]=((cosine * x + sine * y) + p->sx);  
    54.             feat[i+1]=((-sine * x + cosine * y) + p->sy);  
    55.         }  
    56.         //Initialize the feature descriptor.  
    57.         float *feat_desc = (float *) malloc( 128 * sizeof(float));  
    58.         for (i=0;i<128;i++)  
    59.         {  
    60.             feat_desc[i]=0.0;  
    61.             // printf("%f  ",feat_desc[i]);    
    62.         }  
    63.         //printf("/n");  
    64.         for ( i=0;i<512;++i,++i)  
    65.         {  
    66.             float x_sample = feat[i];  
    67.             float y_sample = feat[i+1];  
    68.             // Interpolate the gradient at the sample position  
    69.             /* 
    70.             0   1   0 
    71.             1   *   1 
    72.             0   1   0   具体插值策略如图示 
    73.             */  
    74.             float sample12=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample, y_sample-1);  
    75.             float sample21=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample-1, y_sample);   
    76.             float sample22=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample, y_sample);   
    77.             float sample23=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample+1, y_sample);   
    78.             float sample32=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample, y_sample+1);   
    79.             //float diff_x = 0.5*(sample23 - sample21);  
    80.             //float diff_y = 0.5*(sample32 - sample12);  
    81.             float diff_x = sample23 - sample21;  
    82.             float diff_y = sample32 - sample12;  
    83.             float mag_sample = sqrt( diff_x*diff_x + diff_y*diff_y );  
    84.             float grad_sample = atan( diff_y / diff_x );  
    85.             if(grad_sample == CV_PI)  
    86.                 grad_sample = -CV_PI;  
    87.             // Compute the weighting for the x and y dimensions.  
    88.             float *x_wght=(float *) malloc( GridSpacing * GridSpacing * sizeof(float));  
    89.             float *y_wght=(float *) malloc( GridSpacing * GridSpacing * sizeof(float));  
    90.             float *pos_wght=(float *) malloc( 8*GridSpacing * GridSpacing * sizeof(float));;  
    91.             for (int m=0;m<32;++m,++m)  
    92.             {  
    93.                 float x=featcenter[m];  
    94.                 float y=featcenter[m+1];  
    95.                 x_wght[m/2] = max(1 - (fabs(x - x_sample)*1.0/GridSpacing), 0);  
    96.                 y_wght[m/2] = max(1 - (fabs(y - y_sample)*1.0/GridSpacing), 0);   
    98.             }  
    99.             for ( m=0;m<16;++m)  
    100.                 for (int n=0;n<8;++n)  
    101.                     pos_wght[m*8+n]=x_wght[m]*y_wght[m];  
    102.             free(x_wght);  
    103.             free(y_wght);  
    104.             //计算方向的加权，首先旋转梯度场到主方向，然后计算差异   
    105.             float diff[8],orient_wght[128];  
    106.             for ( m=0;m<8;++m)  
    107.             {   
    108.                 float angle = grad_sample-(p->ori)-orient_angles[m]+CV_PI;  
    109.                 float temp = angle / (2.0 * CV_PI);  
    110.                 angle -= (int)(temp) * (2.0 * CV_PI);  
    111.                 diff[m]= angle - CV_PI;  
    112.             }  
    113.             // Compute the gaussian weighting.  
    114.             float x=p->sx;  
    115.             float y=p->sy;  
    116.             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);  
    118.             for ( m=0;m<128;++m)  
    119.             {  
    120.                 orient_wght[m] = max((1.0 - 1.0*fabs(diff[m%8])/orient_bin_spacing),0);  
    121.                 feat_desc[m] = feat_desc[m] + orient_wght[m]*pos_wght[m]*g*mag_sample;  
    122.             }  
    123.             free(pos_wght);     
    124.         }  
    125.         free(feat);  
    126.         free(featcenter);  
    127.         float norm=GetVecNorm( feat_desc, 128);  
    128.         for (int m=0;m<128;m++)  
    129.         {  
    130.             feat_desc[m]/=norm;  
    131.             if (feat_desc[m]>0.2)  
    132.                 feat_desc[m]=0.2;  
    133.         }  
    134.         norm=GetVecNorm( feat_desc, 128);  
    135.         for ( m=0;m<128;m++)  
    136.         {  
    137.             feat_desc[m]/=norm;  
    138.             printf("%f  ",feat_desc[m]);    
    139.         }  
    140.         printf("/n");  
    141.         p->descrip = feat_desc;  
    142.         p=p->next;  
    143.     }  
    144.     free(feat_grid);  
    145.     free(feat_samples);  
    146. }  
    148. //为了显示图象金字塔，而作的图像水平拼接  
    149. CvMat* MosaicHorizen( CvMat* im1, CvMat* im2 )  
    150. {  
    151.     int row,col;  
    152.     CvMat *mosaic = cvCreateMat( max(im1->rows,im2->rows),(im1->cols+im2->cols),CV_32FC1);  
    153. #define Mosaic(ROW,COL) ((float*)(mosaic->data.fl + mosaic->step/sizeof(float)*(ROW)))[(COL)]  
    154. #define Im11Mat(ROW,COL) ((float *)(im1->data.fl + im1->step/sizeof(float) *(ROW)))[(COL)]  
    155. #define Im22Mat(ROW,COL) ((float *)(im2->data.fl + im2->step/sizeof(float) *(ROW)))[(COL)]  
    156.     cvZero(mosaic);  
    157.     /* Copy images into mosaic1. */  
    158.     for ( row = 0; row < im1->rows; row++)  
    159.         for ( col = 0; col < im1->cols; col++)  
    160.             Mosaic(row,col)=Im11Mat(row,col) ;  
    161.     for (  row = 0; row < im2->rows; row++)  
    162.         for (  col = 0; col < im2->cols; col++)  
    163.             Mosaic(row, (col+im1->cols) )= Im22Mat(row,col) ;  
    164.     return mosaic;  
    165. }  
    167. //为了显示图象金字塔，而作的图像垂直拼接  
    168. CvMat* MosaicVertical( CvMat* im1, CvMat* im2 )  
    169. {  
    170.     int row,col;  
    171.     CvMat *mosaic = cvCreateMat(im1->rows+im2->rows,max(im1->cols,im2->cols), CV_32FC1);  
    172. #define Mosaic(ROW,COL) ((float*)(mosaic->data.fl + mosaic->step/sizeof(float)*(ROW)))[(COL)]  
    173. #define Im11Mat(ROW,COL) ((float *)(im1->data.fl + im1->step/sizeof(float) *(ROW)))[(COL)]  
    174. #define Im22Mat(ROW,COL) ((float *)(im2->data.fl + im2->step/sizeof(float) *(ROW)))[(COL)]  
    175.     cvZero(mosaic);  
    177.     /* Copy images into mosaic1. */  
    178.     for ( row = 0; row < im1->rows; row++)  
    179.         for ( col = 0; col < im1->cols; col++)  
    180.             Mosaic(row,col)= Im11Mat(row,col) ;  
    181.     for ( row = 0; row < im2->rows; row++)  
    182.         for ( col = 0; col < im2->cols; col++)  
    183.             Mosaic((row+im1->rows),col)=Im22Mat(row,col) ;  
    185.     return mosaic;  
    186. }

    ok，为了版述清晰，再贴一下上文所述的主函数（注，上文已贴出，此是为了版述清晰，重复造轮）：

    1. int main( void )  
    2. {  
    3.     //声明当前帧IplImage指针  
    4.     IplImage* src = NULL;   
    5.     IplImage* image1 = NULL;   
    6.     IplImage* grey_im1 = NULL;   
    7.     IplImage* DoubleSizeImage = NULL;  
    9.     IplImage* mosaic1 = NULL;   
    10.     IplImage* mosaic2 = NULL;   
    12.     CvMat* mosaicHorizen1 = NULL;  
    13.     CvMat* mosaicHorizen2 = NULL;  
    14.     CvMat* mosaicVertical1 = NULL;  
    16.     CvMat* image1Mat = NULL;  
    17.     CvMat* tempMat=NULL;  
    19.     ImageOctaves *Gaussianpyr;  
    20.     int rows,cols;  
    22. #define Im1Mat(ROW,COL) ((float *)(image1Mat->data.fl + image1Mat->step/sizeof(float) *(ROW)))[(COL)]  
    24.     //灰度图象像素的数据结构  
    25. #define Im1B(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3]  
    26. #define Im1G(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3+1]  
    27. #define Im1R(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3+2]  
    29.     storage = cvCreateMemStorage(0);   
    31.     //读取图片  
    32.     if( (src = cvLoadImage( "street1.jpg", 1)) == 0 )  // test1.jpg einstein.pgm back1.bmp  
    33.         return -1;  
    35.     //为图像分配内存   
    36.     image1 = cvCreateImage(cvSize(src->width, src->height),  IPL_DEPTH_8U,3);  
    37.     grey_im1 = cvCreateImage(cvSize(src->width, src->height),  IPL_DEPTH_8U,1);  
    38.     DoubleSizeImage = cvCreateImage(cvSize(2*(src->width), 2*(src->height)),  IPL_DEPTH_8U,3);  
    40.     //为图像阵列分配内存，假设两幅图像的大小相同，tempMat跟随image1的大小  
    41.     image1Mat = cvCreateMat(src->height, src->width, CV_32FC1);  
    42.     //转化成单通道图像再处理  
    43.     cvCvtColor(src, grey_im1, CV_BGR2GRAY);  
    44.     //转换进入Mat数据结构,图像操作使用的是浮点型操作  
    45.     cvConvert(grey_im1, image1Mat);  
    47.     double t = (double)cvGetTickCount();  
    48.     //图像归一化  
    49.     cvConvertScale( image1Mat, image1Mat, 1.0/255, 0 );  
    51.     int dim = min(image1Mat->rows, image1Mat->cols);  
    52.     numoctaves = (int) (log((double) dim) / log(2.0)) - 2;    //金字塔阶数  
    53.     numoctaves = min(numoctaves, MAXOCTAVES);  
    55.     //SIFT算法第一步，预滤波除噪声，建立金字塔底层  
    56.     tempMat = ScaleInitImage(image1Mat) ;  
    57.     //SIFT算法第二步，建立Guassian金字塔和DOG金字塔  
    58.     Gaussianpyr = BuildGaussianOctaves(tempMat) ;  
    60.     t = (double)cvGetTickCount() - t;  
    61.     printf( "the time of build Gaussian pyramid and DOG pyramid is %.1f/n", t/(cvGetTickFrequency()*1000.) );  
    63. #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(Gaussianpyr[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + Gaussianpyr[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]  
    64.     //显示高斯金字塔  
    65.     for (int i=0; i<numoctaves;i++)  
    66.     {  
    67.         if (i==0)  
    68.         {  
    69.             mosaicHorizen1=MosaicHorizen( (Gaussianpyr[0].Octave)[0].Level, (Gaussianpyr[0].Octave)[1].Level );  
    70.             for (int j=2;j<SCALESPEROCTAVE+3;j++)  
    71.                 mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (Gaussianpyr[0].Octave)[j].Level );  
    72.             for ( j=0;j<NUMSIZE;j++)  
    73.                 mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
    74.         }  
    75.         else if (i==1)  
    76.         {  
    77.             mosaicHorizen2=MosaicHorizen( (Gaussianpyr[1].Octave)[0].Level, (Gaussianpyr[1].Octave)[1].Level );  
    78.             for (int j=2;j<SCALESPEROCTAVE+3;j++)  
    79.                 mosaicHorizen2=MosaicHorizen( mosaicHorizen2, (Gaussianpyr[1].Octave)[j].Level );  
    80.             for ( j=0;j<NUMSIZE;j++)  
    81.                 mosaicHorizen2=halfSizeImage(mosaicHorizen2);  
    82.             mosaicVertical1=MosaicVertical( mosaicHorizen1, mosaicHorizen2 );  
    83.         }  
    84.         else  
    85.         {  
    86.             mosaicHorizen1=MosaicHorizen( (Gaussianpyr[i].Octave)[0].Level, (Gaussianpyr[i].Octave)[1].Level );  
    87.             for (int j=2;j<SCALESPEROCTAVE+3;j++)  
    88.                 mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (Gaussianpyr[i].Octave)[j].Level );  
    89.             for ( j=0;j<NUMSIZE;j++)  
    90.                 mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
    91.             mosaicVertical1=MosaicVertical( mosaicVertical1, mosaicHorizen1 );  
    92.         }  
    93.     }  
    94.     mosaic1 = cvCreateImage(cvSize(mosaicVertical1->width, mosaicVertical1->height),  IPL_DEPTH_8U,1);  
    95.     cvConvertScale( mosaicVertical1, mosaicVertical1, 255.0, 0 );  
    96.     cvConvertScaleAbs( mosaicVertical1, mosaic1, 1, 0 );  
    98.     //  cvSaveImage("GaussianPyramid of me.jpg",mosaic1);  
    99.     cvNamedWindow("mosaic1",1);  
    100.     cvShowImage("mosaic1", mosaic1);  
    101.     cvWaitKey(0);  
    102.     cvDestroyWindow("mosaic1");  
    103.     //显示DOG金字塔  
    104.     for ( i=0; i<numoctaves;i++)  
    105.     {  
    106.         if (i==0)  
    107.         {  
    108.             mosaicHorizen1=MosaicHorizen( (DOGoctaves[0].Octave)[0].Level, (DOGoctaves[0].Octave)[1].Level );  
    109.             for (int j=2;j<SCALESPEROCTAVE+2;j++)  
    110.                 mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (DOGoctaves[0].Octave)[j].Level );  
    111.             for ( j=0;j<NUMSIZE;j++)  
    112.                 mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
    113.         }  
    114.         else if (i==1)  
    115.         {  
    116.             mosaicHorizen2=MosaicHorizen( (DOGoctaves[1].Octave)[0].Level, (DOGoctaves[1].Octave)[1].Level );  
    117.             for (int j=2;j<SCALESPEROCTAVE+2;j++)  
    118.                 mosaicHorizen2=MosaicHorizen( mosaicHorizen2, (DOGoctaves[1].Octave)[j].Level );  
    119.             for ( j=0;j<NUMSIZE;j++)  
    120.                 mosaicHorizen2=halfSizeImage(mosaicHorizen2);  
    121.             mosaicVertical1=MosaicVertical( mosaicHorizen1, mosaicHorizen2 );  
    122.         }  
    123.         else  
    124.         {  
    125.             mosaicHorizen1=MosaicHorizen( (DOGoctaves[i].Octave)[0].Level, (DOGoctaves[i].Octave)[1].Level );  
    126.             for (int j=2;j<SCALESPEROCTAVE+2;j++)  
    127.                 mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (DOGoctaves[i].Octave)[j].Level );  
    128.             for ( j=0;j<NUMSIZE;j++)  
    129.                 mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
    130.             mosaicVertical1=MosaicVertical( mosaicVertical1, mosaicHorizen1 );  
    131.         }  
    132.     }  
    133.     //考虑到DOG金字塔各层图像都会有正负，所以，必须寻找最负的，以将所有图像抬高一个台阶去显示  
    134.     double min_val=0;  
    135.     double max_val=0;  
    136.     cvMinMaxLoc( mosaicVertical1, &min_val, &max_val,NULL, NULL, NULL );  
    137.     if ( min_val<0.0 )  
    138.         cvAddS( mosaicVertical1, cvScalarAll( (-1.0)*min_val ), mosaicVertical1, NULL );  
    139.     mosaic2 = cvCreateImage(cvSize(mosaicVertical1->width, mosaicVertical1->height),  IPL_DEPTH_8U,1);  
    140.     cvConvertScale( mosaicVertical1, mosaicVertical1, 255.0/(max_val-min_val), 0 );  
    141.     cvConvertScaleAbs( mosaicVertical1, mosaic2, 1, 0 );  
    143.     //  cvSaveImage("DOGPyramid of me.jpg",mosaic2);  
    144.     cvNamedWindow("mosaic1",1);  
    145.     cvShowImage("mosaic1", mosaic2);  
    146.     cvWaitKey(0);  
    148.     //SIFT算法第三步：特征点位置检测，最后确定特征点的位置  
    149.     int keycount=DetectKeypoint(numoctaves, Gaussianpyr);  
    150.     printf("the keypoints number are %d ;/n", keycount);  
    151.     cvCopy(src,image1,NULL);  
    152.     DisplayKeypointLocation( image1 ,Gaussianpyr);  
    154.     cvPyrUp( image1, DoubleSizeImage, CV_GAUSSIAN_5x5 );  
    155.     cvNamedWindow("image1",1);  
    156.     cvShowImage("image1", DoubleSizeImage);  
    157.     cvWaitKey(0);    
    158.     cvDestroyWindow("image1");  
    160.     //SIFT算法第四步：计算高斯图像的梯度方向和幅值，计算各个特征点的主方向  
    161.     ComputeGrad_DirecandMag(numoctaves, Gaussianpyr);  
    162.     AssignTheMainOrientation( numoctaves, Gaussianpyr,mag_pyr,grad_pyr);  
    163.     cvCopy(src,image1,NULL);  
    164.     DisplayOrientation ( image1, Gaussianpyr);  
    166.     //  cvPyrUp( image1, DoubleSizeImage, CV_GAUSSIAN_5x5 );  
    167.     cvNamedWindow("image1",1);  
    168.     //  cvResizeWindow("image1", 2*(image1->width), 2*(image1->height) );  
    169.     cvShowImage("image1", image1);  
    170.     cvWaitKey(0);  
    172.     //SIFT算法第五步：抽取各个特征点处的特征描述字  
    173.     ExtractFeatureDescriptors( numoctaves, Gaussianpyr);  
    174.     cvWaitKey(0);  
    176.     //销毁窗口  
    177.     cvDestroyWindow("image1");  
    178.     cvDestroyWindow("mosaic1");  
    179.     //释放图像  
    180.     cvReleaseImage(&image1);  
    181.     cvReleaseImage(&grey_im1);  
    182.     cvReleaseImage(&mosaic1);  
    183.     cvReleaseImage(&mosaic2);  
    184.     return 0;  
    185. }

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

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

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

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

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

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

    完。

    updated

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