ITdaan
搜索本文相关内容

sift算法研究(3)c语言实现sift算法、下

本文转载自  ruifangcui7758  查看原文  2011/10/11  1908    算法/ image/ sift/ 语言/ c语言/ 实现/ float           
 作者：July、二零一一年三月十二日
出处： http://blog.csdn.net/v_JULY_v。
参考：Rob Hess维护的sift 库
环境：windows xp+vc6.0
条件：c语言实现。
说明：本BLOG内会陆续一一实现所有经典算法。
------------------------


本文接上，教你一步一步用c语言实现sift算法、上，而来：
函数编写
    ok，接上文，咱们一个一个的来编写main函数中所涉及到所有函数，这也是本文的关键部分：

view plain copy to clipboard print ?
  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);  
  8.    
  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. }  
  16.   
  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);  
  24.    
  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)]  
  27.    
  28.  for ( j = 0; j < h; j++)   
  29.   for ( i = 0; i < w; i++)   
  30.    Imnew(j,i)=Im(j/2, i/2);  
  31.     
  32.   return imnew;  
  33. }  
  34.   
  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);  
  42.    
  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)]  
  45.    
  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. }  
  75.   
  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)]  
  85.    
  86.  irow = (int) row;  
  87.  icol = (int) col;  
  88.    
  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. }  
  118.   
  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;  
  124.    
  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. }  
  132.   
  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. }  
  143.   
  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. }  
  152.   
  153. //产生1D高斯核   
  154. float* GaussianKernel1D(float sigma, int dim)   
  155. {  
  156.    
  157.  unsigned int i;  
  158.  //printf("GaussianKernel1D(): Creating 1x%d vector for sigma=%.3f gaussian kernel/n", dim, sigma);  
  159.    
  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. }  
  177.   
  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. }  
  208.   
  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. }  
  231.   
  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;  
  237.    
  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. }  
  247.   
  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. }  
  269.   
  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. }  
  284.   
  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);  
  296.    
  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算法第一步：扩大图像，预滤波剔除噪声，得到金字塔的最底层-第一阶的第一层：

view plain copy to clipboard print ?
  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)]  
  16.     
  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);   
  20.     
  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决定，给定一个基本图像，计算它的高斯金字塔图像，返回外部向量是阶梯指针，内部向量是每一个阶梯内部的不同尺度图像。

view plain copy to clipboard print ?
  1. //SIFT算法第二步   
  2. ImageOctaves* BuildGaussianOctaves(CvMat * image)   
  3. {  
  4.     ImageOctaves *octaves;  
  5.  CvMat *tempMat;  
  6.     CvMat *dst;  
  7.  CvMat *temp;  
  8.    
  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) );  
  20.    
  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);  
  23.    
  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);  
  29.       
  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;  
  40.     
  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);  
  49.     
  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);  
  83.      
  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金字塔中的局部最大值，找到之后，还要经过两个检验才能确认为特征点：一是它必须有明显的差异，二是他不应该是边缘点，（也就是说，在极值点处的主曲率比应该小于某一个阈值）。

view plain copy to clipboard print ?
  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)]  
  8.    
  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.      {  
  23.         
  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  
  36.           
  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  
  45.           
  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))&&  
  65.           
  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))&&   
  74.           
  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. }  
  134.   
  135. //在图像中，显示SIFT特征点的位置   
  136. void DisplayKeypointLocation(IplImage* image, ImageOctaves *GaussianPyr)  
  137. {  
  138.    
  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. }  
  154.   
  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算法第四步
view plain copy to clipboard print ?
  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)]  
  14.    
  15.     int keypoint_count = 0;  
  16.  Keypoint p = keypoints; // p指向第一个结点   
  17.    
  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)]  
  33.      
  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");  
  74.       
  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;  
  115.        
  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.      }  
  123.        
  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. }  
  152.   
  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. }  
  165.   
  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. }  
  185.   
  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;  
  200.    
  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. }  
  212.   
  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  
  231.     
  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网格内梯度方向的描述，旋转网格到主方向，插值得到网格处梯度值。
一个特征点可以用2*2*8=32维的向量，也可以用4*4*8=128维的向量更精确的进行描述。

view plain copy to clipboard print ?
  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;  
  32.           
  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);   
  97.       
  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);  
  117.       
  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. }  
  147.   
  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. }  
  166.   
  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);  
  176.    
  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) ;  
  184.      
  185.    return mosaic;  
  186. }  

 

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

view plain copy to clipboard print ?
  1. int main( void )  
  2. {  
  3.  //声明当前帧IplImage指针   
  4.  IplImage* src = NULL;   
  5.  IplImage* image1 = NULL;   
  6.  IplImage* grey_im1 = NULL;   
  7.  IplImage* DoubleSizeImage = NULL;  
  8.    
  9.  IplImage* mosaic1 = NULL;   
  10.  IplImage* mosaic2 = NULL;   
  11.    
  12.  CvMat* mosaicHorizen1 = NULL;  
  13.  CvMat* mosaicHorizen2 = NULL;  
  14.  CvMat* mosaicVertical1 = NULL;  
  15.    
  16.  CvMat* image1Mat = NULL;  
  17.  CvMat* tempMat=NULL;  
  18.    
  19.  ImageOctaves *Gaussianpyr;  
  20.  int rows,cols;  
  21.   
  22. #define Im1Mat(ROW,COL) ((float *)(image1Mat->data.fl + image1Mat->step/sizeof(float) *(ROW)))[(COL)]  
  23.    
  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]  
  28.    
  29.  storage = cvCreateMemStorage(0);   
  30.    
  31.  //读取图片   
  32.  if( (src = cvLoadImage( "street1.jpg", 1)) == 0 )  // test1.jpg einstein.pgm back1.bmp  
  33.   return -1;  
  34.   
  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);  
  39.   
  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);  
  46.    
  47.  double t = (double)cvGetTickCount();  
  48.  //图像归一化   
  49.  cvConvertScale( image1Mat, image1Mat, 1.0/255, 0 );  
  50.    
  51.  int dim = min(image1Mat->rows, image1Mat->cols);  
  52.  numoctaves = (int) (log((double) dim) / log(2.0)) - 2;    //金字塔阶数  
  53.  numoctaves = min(numoctaves, MAXOCTAVES);  
  54.   
  55.  //SIFT算法第一步，预滤波除噪声，建立金字塔底层   
  56.  tempMat = ScaleInitImage(image1Mat) ;  
  57.  //SIFT算法第二步，建立Guassian金字塔和DOG金字塔  
  58.  Gaussianpyr = BuildGaussianOctaves(tempMat) ;  
  59.    
  60.  t = (double)cvGetTickCount() - t;  
  61.  printf( "the time of build Gaussian pyramid and DOG pyramid is %.1f/n", t/(cvGetTickFrequency()*1000.) );  
  62.    
  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 );  
  97.    
  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 );  
  142.    
  143.  //  cvSaveImage("DOGPyramid of me.jpg",mosaic2);  
  144.  cvNamedWindow("mosaic1",1);  
  145.  cvShowImage("mosaic1", mosaic2);  
  146.  cvWaitKey(0);  
  147.    
  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);  
  153.    
  154.  cvPyrUp( image1, DoubleSizeImage, CV_GAUSSIAN_5x5 );  
  155.  cvNamedWindow("image1",1);  
  156.  cvShowImage("image1", DoubleSizeImage);  
  157.  cvWaitKey(0);    
  158.  cvDestroyWindow("image1");  
  159.    
  160.  //SIFT算法第四步：计算高斯图像的梯度方向和幅值，计算各个特征点的主方向   
  161.  ComputeGrad_DirecandMag(numoctaves, Gaussianpyr);  
  162.  AssignTheMainOrientation( numoctaves, Gaussianpyr,mag_pyr,grad_pyr);  
  163.  cvCopy(src,image1,NULL);  
  164.  DisplayOrientation ( image1, Gaussianpyr);  
  165.    
  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);  
  171.    
  172.  //SIFT算法第五步：抽取各个特征点处的特征描述字   
  173.  ExtractFeatureDescriptors( numoctaves, Gaussianpyr);  
  174.  cvWaitKey(0);  
  175.    
  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年照）：

完。 

 

版权声明：
1、本文版权归本人和CSDN共同拥有。转载，请注明出处及作者本人。
2、版权侵犯者，无论任何人，任何网站，1、永久追踪，2、永久谴责，3、永久追究法律责任的权利。
       July、二零一一年三月十二日声明。

 

本文转载自：http://blog.csdn.net/v_JULY_v/article/details/6246213，如有转载请注明原出处。谢谢

智能推荐

注意！

本站转载的文章为个人学习借鉴使用，本站对版权不负任何法律责任。如果侵犯了您的隐私权益，请联系我们删除。



猜您在找
sift算法c语言实现 九之再续：教你一步一步用c语言实现sift算法、下 九之再续：教你一步一步用c语言实现sift算法、下 教你一步一步用c语言实现sift算法、下 编译Rob Hess基于OpenCV的SIFT算法的C语言实现
智能推荐
 
© 2014-2019 ITdaan.com 粤ICP备14056181号  

赞助商广告