教你一步一步用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、二零一二年十月十一日。