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

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

作者：July、二零一一年三月十二日
出处：http://blog.youkuaiyun.com/v_JULY_v
参考：Rob Hess维护的sift 库
环境：windows xp+vc6.0
条件：c语言实现。
说明：本BLOG内会陆续一一实现所有经典算法。
------------------------

引言：
    在我写的关于sift算法的前倆篇文章里头，已经对sift算法有了初步的介绍：九、图像特征提取与匹配之SIFT算法，而后在：九（续）、sift算法的编译与实现里，我也简单记录下了如何利用opencv，gsl等库编译运行sift程序。
    但据一朋友表示，是否能用c语言实现sift算法，同时，尽量不用到opencv，gsl等第三方库之类的东西。而且，Rob Hess维护的sift 库，也不好懂，有的人根本搞不懂是怎么一回事。
    那么本文，就教你如何利用c语言一步一步实现sift算法，同时，你也就能真正明白sift算法到底是怎么一回事了。

    ok，先看一下，本程序最终运行的效果图，sift 算法分为五个步骤（下文详述），对应以下第二--第六幅图：

 

sift算法的步骤
    要实现一个算法，首先要完全理解这个算法的原理或思想。咱们先来简单了解下，什么叫sift算法：
    sift，尺度不变特征转换，是一种电脑视觉的算法用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量，此算法由 David Lowe 在1999年所发表，2004年完善总结。

    所谓，Sift算法就是用不同尺度（标准差）的高斯函数对图像进行平滑，然后比较平滑后图像的差别，
差别大的像素就是特征明显的点。

    以下是sift算法的五个步骤：
    一、建立图像尺度空间(或高斯金字塔)，并检测极值点

    首先建立尺度空间，要使得图像具有尺度空间不变形，就要建立尺度空间，sift算法采用了高斯函数来建立尺度空间，高斯函数公式为：
        
    上述公式G(x,y,e)，即为尺度可变高斯函数。

    而，一个图像的尺度空间L（x,y,e) ,定义为原始图像I(x,y)与上述的一个可变尺度的2维高斯函数G(x,y,e) 卷积运算。
    即，原始影像I(x，y)在不同的尺度e下，与高斯函数G(x,y,e)进行卷积，得到L(x,y,e)，如下：
        

    以上的（x，y）是空间坐标， e，是尺度坐标，或尺度空间因子，e的大小决定平滑程度，大尺度对应图像的概貌特征，小尺度对应图像的细节特征。大的e值对应粗糙尺度(低分辨率)，反之，对应精细尺度(高分辨率)。

    尺度，受e这个参数控制的表示。而不同的L（x,y,e）就构成了尺度空间，具体计算的时候，即使连续的高斯函数，都被离散为（一般为奇数大小）(2*k+1) *(2*k+1)矩阵，来和数字图像进行卷积运算。
    随着e的变化，建立起不同的尺度空间，或称之为建立起图像的高斯金字塔。

    但，像上述L(x,y,e) = G(x,y,e)*I(x,y)的操作，在进行高斯卷积时，整个图像就要遍历所有的像素进行卷积(边界点除外)，于此，就造成了时间和空间上的很大浪费。
    为了更有效的在尺度空间检测到稳定的关键点，也为了缩小时间和空间复杂度，对上述的操作作了一个改建：即，提出了高斯差分尺度空间（DOG scale-space）。利用不同尺度的高斯差分与原始图像I(x,y)相乘 ，卷积生成。
        
    DOG算子计算简单，是尺度归一化的LOG算子的近似。

    ok，耐心点，咱们再来总结一下上述内容：
    1、高斯卷积
    在组建一组尺度空间后，再组建下一组尺度空间，对上一组尺度空间的最后一幅图像进行二分之一采样，得到下一组尺度空间的第一幅图像，然后进行像建立第一组尺度空间那样的操作，得到第二组尺度空间，公式定义为
         L(x,y,e) = G(x,y,e)*I(x,y)

    图像金字塔的构建：图像金字塔共O组，每组有S层，下一组的图像由上一组图像降采样得到，效果图，图A如下(左为上一组，右为下一组)：

    
    2、高斯差分
    在尺度空间建立完毕后，为了能够找到稳定的关键点，采用高斯差分的方法来检测那些在局部位置的极值点，即采用俩个相邻的尺度中的图像相减，即公式定义为：
        D(x,y,e) = ((G(x,y,ke) - G(x,y,e)) * I(x,y) 
                 = L(x,y,ke) - L(x,y,e)

    效果图，图B：

    SIFT的精妙之处在于采用图像金字塔的方法解决这一问题，我们可以把两幅图像想象成是连续的，分别以它们作为底面作四棱锥，就像金字塔，那么每一个 截面与原图像相似，那么两个金字塔中必然会有包含大小一致的物体的无穷个截面，但应用只能是离散的，所以我们只能构造有限层，层数越多当然越好，但处理时 间会相应增加，层数太少不行，因为向下采样的截面中可能找不到尺寸大小一致的两个物体的图像。

    咱们再来具体阐述下构造D(x,y,e)的详细步骤：
    1、首先采用不同尺度因子的高斯核对图像进行卷积以得到图像的不同尺度空间，将这一组图像作为金子塔图像的第一层。
    2、接着对第一层图像中的2倍尺度图像（相对于该层第一幅图像的2倍尺度）以2倍像素距离进行下采样来得到金子塔图像的第二层中的第一幅图像，对该图像采用不同尺度因子的高斯核进行卷积，以获得金字塔图像中第二层的一组图像。
    3、再以金字塔图像中第二层中的2倍尺度图像（相对于该层第一幅图像的2倍尺度）以2倍像素距离进行下采样来得到金字塔图像的第三层中的第一幅图像，对该图像采用不同尺度因子的高斯核进行卷积，以获得金字塔图像中第三层的一组图像。这样依次类推，从而获得了金字塔图像的每一层中的一组图像，如下图所示：

    4、对上图得到的每一层相邻的高斯图像相减，就得到了高斯差分图像，如下述第一幅图所示。下述第二幅图中的右列显示了将每组中相邻图像相减所生成的高斯差分图像的结果，限于篇幅，图中只给出了第一层和第二层高斯差分图像的计算（下述俩幅图统称为图2）：

图像金字塔的建立：对于一幅图像I,建立其在不同尺度(scale)的图像，也成为子八度（octave），这是为了scale-invariant，也就是在任何尺度都能够有对应的特征点，第一个子八度的scale为原图大小，后面每个octave为上一个octave降采样的结果，即原图的1/4（长宽分别减半），构成下一个子八度（高一层金字塔）。

    5、因为高斯差分函数是归一化的高斯拉普拉斯函数的近似，所以可以从高斯差分金字塔分层结构提取出图像中的极值点作为候选的特征点。对DOG 尺度空间每个点与相邻尺度和相邻位置的点逐个进行比较，得到的局部极值位置即为特征点所处的位置和对应的尺度。

    二、检测关键点
    为了寻找尺度空间的极值点，每一个采样点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。如下图，图3所示，中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。

    因为需要同相邻尺度进行比较，所以在一组高斯差分图像中只能检测到两个尺度的极值点（如上述第二幅图中右图的五角星标识），而其它尺度的极值点检测则需要在图像金字塔的上一层高斯差分图像中进行。依次类推，最终在图像金字塔中不同层的高斯差分图像中完成不同尺度极值的检测。
    当然这样产生的极值点并不都是稳定的特征点，因为某些极值点响应较弱，而且DOG算子会产生较强的边缘响应。

    三、关键点方向的分配
    为了使描述符具有旋转不变性，需要利用图像的局部特征为给每一个关键点分配一个方向。利用关键点邻域像素的梯度及方向分布的特性，可以得到梯度模值和方向如下：

   
    其中，尺度为每个关键点各自所在的尺度。
    在以关键点为中心的邻域窗口内采样，并用直方图统计邻域像素的梯度方向。梯度直方图的范围是0～360度，其中每10度一个方向，总共36个方向。
直方图的峰值则代表了该关键点处邻域梯度的主方向，即作为该关键点的方向。

    在计算方向直方图时，需要用一个参数等于关键点所在尺度1.5倍的高斯权重窗对方向直方图进行加权，上图中用蓝色的圆形表示，中心处的蓝色较重，表示权值最大，边缘处颜色潜，表示权值小。如下图所示，该示例中为了简化给出了8方向的方向直方图计算结果，实际sift创始人David Lowe的原论文中采用36方向的直方图。

    方向直方图的峰值则代表了该特征点处邻域梯度的方向，以直方图中最大值作为该关键点的主方向。为了增强匹配的鲁棒性，只保留峰值大于主方向峰值80％的方向作为该关键点的辅方向。因此，对于同一梯度值的多个峰值的关键点位置，在相同位置和尺度将会有多个关键点被创建但方向不同。仅有15％的关键点被赋予多个方向，但可以明显的提高关键点匹配的稳定性。
        
    至此，图像的关键点已检测完毕，每个关键点有三个信息：位置、所处尺度、方向。由此可以确定一个SIFT特征区域。

    四、特征点描述符
    通过以上步骤，对于每一个关键点，拥有三个信息：位置、尺度以及方向。接下来就是为每个关键点建立一个描述符，使其不随各种变化而改变，比如光照变化、视角变化等等。并且描述符应该有较高的独特性，以便于提高特征点正确匹配的概率。  
首先将坐标轴旋转为关键点的方向，以确保旋转不变性。

    接下来以关键点为中心取8×8的窗口。
    上图，图5中左部分的中央黑点为当前关键点的位置，每个小格代表关键点邻域所在尺度空间的一个像素，箭头方向代表该像素的梯度方向，箭头长度代表梯度模值，图中蓝色的圈代表高斯加权的范围（越靠近关键点的像素梯度方向信息贡献越大）。
    然后在每4×4的小块上计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，即可形成一个种子点，如图5右部分所示。此图中一个关键点由2×2共4个种子点组成，每个种子点有8个方向向量信息。这种邻域方向性信息联合的思想增强了算法抗噪声的能力，同时对于含有定位误差的特征匹配也提供了较好的容错性。

    实际计算过程中，为了增强匹配的稳健性，Lowe建议对每个关键点使用4×4共16个种子点来描述，这样对于一个关键点就可以产生128个数据，即最终形成128维的SIFT特征向量。此时SIFT特征向量已经去除了尺度变化、旋转等几何变形因素的影响，再继续将特征向量的长度归一化，则可以进一步去除光照变化的影响。

    五、最后一步：当两幅图像的SIFT特征向量生成后，下一步我们采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量。取上图中，图像A中的某个关键点，并找出其与图像B中欧式距离最近的前两个关键点，在这两个关键点中，如果最近的距离除以次近的距离少于某个比例阈值，则接受这一对匹配点。降低这个比例阈值，SIFT匹配点数目会减少，但更加稳定。关于sift 算法的更多理论介绍请参看此文：http://blog.youkuaiyun.com/abcjennifer/article/details/7639681。

 

sift算法的逐步c实现
    ok，上文搅了那么多的理论，如果你没有看懂它，咋办列?没关系，下面，咱们来一步一步实现此sift算法，即使你没有看到上述的理论，慢慢的，你也会明白sift算法到底是怎么一回事，sift算法到底是怎么实现的...。
   yeah，请看：

前期工作：
在具体编写核心函数之前，得先做几个前期的准备工作：

0、头文件：

1. #ifdef _CH_  
2. #pragma package <opencv>  
3. #endif  
4.   
5. #ifndef _EiC  
6. #include <stdio.h>  
7.   
8. #include "stdlib.h"  
9. #include "string.h"   
10. #include "malloc.h"   
11. #include "math.h"   
12. #include <assert.h>  
13. #include <ctype.h>  
14. #include <time.h>  
15. #include <cv.h>  
16. #include <cxcore.h>  
17. #include <highgui.h>  
18. #include <vector>  
19. #endif  
20.   
21. #ifdef _EiC  
22. #define WIN32  
23. #endif  

1、定义几个宏，及变量，以免下文函数中，突然冒出一个变量，而您却不知道怎么一回事：

1. #define NUMSIZE 2  
2. #define GAUSSKERN 3.5  
3. #define PI 3.14159265358979323846  
4.   
5. //Sigma of base image -- See D.L.'s paper.  
6. #define INITSIGMA 0.5  
7. //Sigma of each octave -- See D.L.'s paper.  
8. #define SIGMA sqrt(3)//1.6//  
9.   
10. //Number of scales per octave.  See D.L.'s paper.  
11. #define SCALESPEROCTAVE 2  
12. #define MAXOCTAVES 4  
13. int     numoctaves;  
14.   
15. #define CONTRAST_THRESHOLD   0.02  
16. #define CURVATURE_THRESHOLD  10.0  
17. #define DOUBLE_BASE_IMAGE_SIZE 1  
18. #define peakRelThresh 0.8  
19. #define LEN 128  
20.   
21. // temporary storage  
22. CvMemStorage* storage = 0;   

2、然后，咱们还得，声明几个变量，以及建几个数据结构（数据结构是一切程序事物的基础麻，:D。）：

1. //Data structure for a float image.  
2. typedef struct ImageSt {        /*金字塔每一层*/  
3.    
4.  float levelsigma;  
5.  int levelsigmalength;  
6.  float absolute_sigma;  
7.  CvMat *Level;       //CvMat是OPENCV的矩阵类，其元素可以是图像的象素值         
8. } ImageLevels;  
9.   
10. typedef struct ImageSt1 {      /*金字塔每一阶梯*/  
11.  int row, col;          //Dimensions of image.   
12.  float subsample;  
13.  ImageLevels *Octave;                
14. } ImageOctaves;  
15.   
16. ImageOctaves *DOGoctaves;        
17. //DOG pyr，DOG算子计算简单，是尺度归一化的LoG算子的近似。  
18.     
19. ImageOctaves *mag_thresh ;  
20. ImageOctaves *mag_pyr ;  
21. ImageOctaves *grad_pyr ;  
22.   
23. //keypoint数据结构，Lists of keypoints are linked by the "next" field.  
24. typedef struct KeypointSt   
25. {  
26.  float row, col; /* 反馈回原图像大小，特征点的位置 */  
27.  float sx,sy;    /* 金字塔中特征点的位置*/  
28.  int octave,level;/*金字塔中，特征点所在的阶梯、层次*/  
29.    
30.  float scale, ori,mag; /*所在层的尺度sigma,主方向orientation (range [-PI,PI])，以及幅值*/  
31.  float *descrip;       /*特征描述字指针：128维或32维等*/  
32.  struct KeypointSt *next;/* Pointer to next keypoint in list. */  
33. } *Keypoint;  
34.   
35. //定义特征点具体变量  
36. Keypoint keypoints=NULL;      //用于临时存储特征点的位置等  
37. Keypoint keyDescriptors=NULL; //用于最后的确定特征点以及特征描述字  

3、声明几个图像的基本处理函数：

1. CvMat * halfSizeImage(CvMat * im);     //缩小图像：下采样  
2. CvMat * doubleSizeImage(CvMat * im);   //扩大图像：最近临方法  
3. CvMat * doubleSizeImage2(CvMat * im);  //扩大图像：线性插值  
4. float getPixelBI(CvMat * im, float col, float row);//双线性插值函数  
5. void normalizeVec(float* vec, int dim);//向量归一化    
6. CvMat* GaussianKernel2D(float sigma);  //得到2维高斯核  
7. void normalizeMat(CvMat* mat) ;        //矩阵归一化  
8. float* GaussianKernel1D(float sigma, int dim) ; //得到1维高斯核  
9.   
10. //在具体像素处宽度方向进行高斯卷积  
11. float ConvolveLocWidth(float* kernel, int dim, CvMat * src, int x, int y) ;    
12. //在整个图像宽度方向进行1D高斯卷积  
13. void Convolve1DWidth(float* kern, int dim, CvMat * src, CvMat * dst) ;         
14. //在具体像素处高度方向进行高斯卷积  
15. float ConvolveLocHeight(float* kernel, int dim, CvMat * src, int x, int y) ;   
16. //在整个图像高度方向进行1D高斯卷积  
17. void Convolve1DHeight(float* kern, int dim, CvMat * src, CvMat * dst);       
18. //用高斯函数模糊图像    
19. int BlurImage(CvMat * src, CvMat * dst, float sigma) ;            

             
算法核心
   本程序中，sift算法被分为以下五个步骤及其相对应的函数（可能表述与上，或与前俩篇文章有所偏差，但都一个意思）：

1. //SIFT算法第一步：图像预处理  
2. CvMat *ScaleInitImage(CvMat * im) ;                  //金字塔初始化  
3.   
4. //SIFT算法第二步：建立高斯金字塔函数  
5. ImageOctaves* BuildGaussianOctaves(CvMat * image) ;  //建立高斯金字塔  
6.   
7. //SIFT算法第三步：特征点位置检测，最后确定特征点的位置  
8. int DetectKeypoint(int numoctaves, ImageOctaves *GaussianPyr);  
9. void DisplayKeypointLocation(IplImage* image, ImageOctaves *GaussianPyr);  
10.   
11. //SIFT算法第四步：计算高斯图像的梯度方向和幅值，计算各个特征点的主方向  
12. void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr);  
13.   
14. int FindClosestRotationBin (int binCount, float angle);  //进行方向直方图统计  
15. void AverageWeakBins (double* bins, int binCount);       //对方向直方图滤波  
16. //确定真正的主方向  
17. bool InterpolateOrientation (double left, double middle,double right, double *degreeCorrection, double *peakValue);  
18. //确定各个特征点处的主方向函数  
19. void AssignTheMainOrientation(int numoctaves, ImageOctaves *GaussianPyr,ImageOctaves *mag_pyr,ImageOctaves *grad_pyr);  
20. //显示主方向  
21. void DisplayOrientation (IplImage* image, ImageOctaves *GaussianPyr);  
22.   
23. //SIFT算法第五步：抽取各个特征点处的特征描述字  
24. void ExtractFeatureDescriptors(int numoctaves, ImageOctaves *GaussianPyr);  
25.   
26. //为了显示图象金字塔，而作的图像水平、垂直拼接  
27. CvMat* MosaicHorizen( CvMat* im1, CvMat* im2 );  
28. CvMat* MosaicVertical( CvMat* im1, CvMat* im2 );  
29.   
30. //特征描述点，网格    
31. #define GridSpacing 4  

主体实现
    ok，以上所有的工作都就绪以后，那么接下来，咱们就先来编写main函数，因为你一看主函数之后，你就立马能发现sift算法的工作流程及其原理了。
    （主函数中涉及到的函数，下一篇文章：一、教你一步一步用c语言实现sift算法、下，咱们自会一个一个编写）：

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

更多见下文：一、教你一步一步用c语言实现sift算法、下。本文完。

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

作者：July、二零一一年三月十二日
出处：http://blog.youkuaiyun.com/v_JULY_v。
参考：Rob Hess维护的sift 库
环境：windows xp+vc6.0
条件：c语言实现。
说明：本BLOG内会陆续一一实现所有经典算法。
------------------------

本文接上，教你一步一步用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);  
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算法第一步：扩大图像，预滤波剔除噪声，得到金字塔的最底层-第一阶的第一层：

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

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金字塔中的局部最大值，找到之后，还要经过两个检验才能确认为特征点：一是它必须有明显的差异，二是他不应该是边缘点，（也就是说，在极值点处的主曲率比应该小于某一个阈值）。

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算法第四步

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维的向量更精确的进行描述。

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，为了版述清晰，再贴一下上文所述的主函数（注，上文已贴出，此是为了版述清晰，重复造轮）：

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年照）：

完。 

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