sift算法c语言实现

前段时间在做三维测量方面的研究，需要得到物体表面三维数据，sift算法是立体匹配中的经典算法，下面是对RobHess的SIFT源代码的注释。部分内容参考网上，在这里向各位大神表示感谢！
http://blog.csdn.net/lsh_2013/article/details/46826141

头文件及函数声明

#include "sift.h"
#include "imgfeatures.h"
#include "utils.h"
#include <cxcore.h>
#include <cv.h>
//将原图转换为32位灰度图并归一化，然后进行一次高斯平滑，并根据参数img_dbl决定
//是否将图像尺寸放大为原图的2倍
static IplImage* create_init_img( IplImage*, int, double );
//将输入图像转换为32位灰度图,并进行归一化
static IplImage* convert_to_gray32( IplImage* );
//根据输入参数建立高斯金字塔
static IplImage*** build_gauss_pyr( IplImage*, int, int, double );
//对输入图像做下采样生成其四分之一大小的图像(每个维度上减半)，使用最近邻差值方法
static IplImage* downsample( IplImage* );
//通过对高斯金字塔中每相邻两层图像相减来建立高斯差分金字塔
static IplImage*** build_dog_pyr( IplImage***, int, int );
//在尺度空间中检测极值点，通过插值精确定位，去除低对比度的点，去除边缘点，返回检测
//到的特征点序列
static CvSeq* scale_space_extrema( IplImage***, int, int, double, int, CvMemStorage*);
//通过在尺度空间中将一个像素点的值与其周围3*3*3邻域内的点比较来决定此点是否极值
//点(极大值或极小都行)
static int is_extremum( IplImage***, int, int, int, int );
//通过亚像素级插值进行极值点精确定位(修正极值点坐标)，并去除低对比度的极值点，
//将修正后的特征点组成feature结构返回
static struct feature* interp_extremum( IplImage***, int, int, int, int, int, double);
//进行一次极值点差值，计算x，y，σ方向(层方向)上的子像素偏移量(增量)
static void interp_step( IplImage***, int, int, int, int, double*, double*, double* );
//在DoG金字塔中计算某点的x方向、y方向以及尺度方向上的偏导数
static CvMat* deriv_3D( IplImage***, int, int, int, int );
//在DoG金字塔中计算某点的3*3海森矩阵
static CvMat* hessian_3D( IplImage***, int, int, int, int );
//计算被插值点的对比度：D + 0.5 * dD^T * X
static double interp_contr( IplImage***, int, int, int, int, double, double, double );
//为一个feature结构分配空间并初始化
static struct feature* new_feature( void );
//去除边缘响应，即通过计算主曲率比值判断某点是否边缘点
static int is_too_edge_like( IplImage*, int, int, int );
//计算特征点序列中每个特征点的尺度
static void calc_feature_scales( CvSeq*, double, int );
//将特征点序列中每个特征点的坐标减半(当设置了将图像放大为原图的2倍时，特征点检测完之后
//调用)
static void adjust_for_img_dbl( CvSeq* );
//计算每个特征点的梯度直方图，找出其主方向，若一个特征点有不止一个主方向，将其分为
//两个特征点
static void calc_feature_oris( CvSeq*, IplImage*** );
//计算指定像素点的梯度方向直方图，返回存放直方图的数组
static double* ori_hist( IplImage*, int, int, int, int, double );
//计算指定点的梯度的幅值magnitude和方向orientation
static int calc_grad_mag_ori( IplImage*, int, int, double*, double* );
//对梯度方向直方图进行高斯平滑，弥补因没有仿射不变性而产生的特征点不稳定的问题
static void smooth_ori_hist( double*, int );
//查找梯度直方图中主方向的梯度幅值，即查找直方图中最大bin的值
static double dominant_ori( double*, int );
//若当前特征点的直方图中某个bin的值大于给定的阈值，则新生成一个特征点并添加到特征点
//序列末尾
static void add_good_ori_features( CvSeq*, double*, int, double, struct feature* );
//对输入的feature结构特征点做深拷贝，返回克隆生成的特征点的指针
static struct feature* clone_feature( struct feature* );
//计算特征点序列中每个特征点的特征描述子向量
static void compute_descriptors( CvSeq*, IplImage***, int, int );
//计算特征点附近区域的方向直方图，此直方图在计算特征描述子中要用到，返回值是一个
//d*d*n的三维数组
static double*** descr_hist( IplImage*, int, int, double, double, int, int );
static void interp_hist_entry( double***, double, double, double, double, int, int);
//将某特征点的方向直方图转换为特征描述子向量，对特征描述子归一化并将所有元素转化为整型，
//存入指定特征点中
static void hist_to_descr( double***, int, int, struct feature* );
//归一化特征点的特征描述子，即将特征描述子数组中每个元素除以特征描述子的模
static void normalize_descr( struct feature* );
//比较函数，将特征点按尺度的降序排列，用在序列排序函数CvSeqSort中
static int feature_cmp( void*, void*, void* );
//释放计算特征描述子过程中用到的方向直方图的内存空间
static void release_descr_hist( double****, int );
//释放金字塔图像组的存储空间
static void release_pyr( IplImage****, int, int );
int sift_features( IplImage* img, struct feature** feat )
{
    return _sift_features( img, feat, SIFT_INTVLS, SIFT_SIGMA, SIFT_CONTR_THR,
                            SIFT_CURV_THR, SIFT_IMG_DBL, SIFT_DESCR_WIDTH,
                            SIFT_DESCR_HIST_BINS );
}

函数实现

主函数

/*使用用户指定的参数在图像中提取SIFT特征点
参数：
img：输入图像
feat：存储特征点的数组的指针，此数组的内存将在本函数中被分配，使用完后必须在调用出释放：free(*feat)
intvls：每组的层数
sigma：初始高斯平滑参数σ
contr_thr：对比度阈值，针对归一化后的图像，用来去除不稳定特征
curv_thr：去除边缘的特征的主曲率阈值
img_dbl：是否将图像放大为之前的两倍
descr_width：特征描述过程中，计算方向直方图时，将特征点附近划分为descr_width*descr_width个区域，每个区域生成一个直方图
descr_hist_bins：特征描述过程中，每个直方图中bin的个数
返回值：提取的特征点个数，若返回-1表明提取失败*/
int _sift_features( IplImage* img, struct feature** feat, int intvls,
                   double sigma, double contr_thr, int curv_thr,
                   int img_dbl, int descr_width, int descr_hist_bins )
{
    IplImage* init_img;//原图经初始化后的图像
    IplImage*** gauss_pyr, *** dog_pyr;//三级指针，高斯金字塔图像组，DoG金字塔图像组
    CvMemStorage* storage;//存储器
    CvSeq* features;//存储特征点的序列，序列中存放的是struct feature类型的指针
    int octvs, i, n = 0;

    //输入参数检查
    if( ! img )
        fatal_error( "NULL pointer error, %s, line %d",  __FILE__, __LINE__ );

    if( ! feat )
        fatal_error( "NULL pointer error, %s, line %d",  __FILE__, __LINE__ );
    /*步骤一：建立尺度空间，即建立高斯差分(DoG)金字塔dog_pyr
    将原图转换为32位灰度图并归一化，然后进行一次高斯平滑，
    并根据参数img_dbl决定是否将图像尺寸放大为原图的2倍*/
    init_img = create_init_img( img, img_dbl, sigma );
    //计算高斯金字塔的组数octvs
    octvs = log( MIN( init_img->width, init_img->height ) ) / log(2) - 2;
    //为了保证连续性，在每一层的顶层继续用高斯模糊生成3幅图像，所以高斯金字塔每组
    //有intvls+3层，DOG金字塔每组有intvls+2层

    //建立高斯金字塔gauss_pyr，是一个octvs*(intvls+3)的图像数组
    gauss_pyr = build_gauss_pyr( init_img, octvs, intvls, sigma );
    //建立高斯差分(DoG)金字塔dog_pyr，是一个octvs*(intvls+2)的图像数组
    dog_pyr = build_dog_pyr( gauss_pyr, octvs, intvls );
    /*步骤二：在尺度空间中检测极值点，并进行精确定位和筛选创建默认大小的内存存储器*/
    storage = cvCreateMemStorage( 0 );
       //在尺度空间中检测极值点，通过插值精确定位，去除低对比度的点，去除边缘点，
    //返回检测到的特征点序列
    features = scale_space_extrema( dog_pyr, octvs, intvls, contr_thr,
        curv_thr, storage );
     //计算特征点序列features中每个特征点的尺度
    calc_feature_scales( features, sigma, intvls );
    if( img_dbl )  //若设置了将图像放大为原图的2倍 
        adjust_for_img_dbl( features );//将特征点序列中每个特征点的坐标减半
    //(当设置了将图像放大为原图的2倍时，特征点检测完之后调用)

    /*步骤三：特征点方向赋值，完成此步骤后，每个特征点有三个信息：位置、尺度、方向*/
    //计算每个特征点的梯度直方图，找出其主方向，若一个特征点有不止一个主方向，将其分为
    //两个特征点
    calc_feature_oris( features, gauss_pyr );
    /*步骤四：计算特征描述子*/
    //计算特征点序列中每个特征点的特征描述子向量 
    compute_descriptors( features, gauss_pyr, descr_width, descr_hist_bins );
       //按特征点尺度的降序排列序列中的元素的顺序，feature_cmp是自定义的比较函数
    cvSeqSort( features, (CvCmpFunc)feature_cmp, NULL );
     //将CvSeq类型的特征点序列features转换为通用的struct feature类型的数组feat
    n = features->total;//特征点个数
    *feat = calloc( n, sizeof(struct feature) );//分配控件
    //将序列features中的元素拷贝到数组feat中，返回数组指针给feat
    *feat = cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ );
    //释放特征点数组feat中所有特征点的feature_data成员，因为此成员中的数据在检测完特征
    //点后就没用了
    for( i = 0; i < n; i++ )
    {
        free( (*feat)[i].feature_data );
        (*feat)[i].feature_data = NULL;
    }
        //释放各种临时数据的存储空间
    cvReleaseMemStorage( &storage );
    cvReleaseImage( &init_img );
    release_pyr( &gauss_pyr, octvs, intvls + 3 );
    release_pyr( &dog_pyr, octvs, intvls + 2 );
    //返回检测到的特征点的个数
    return n;
}

尺度空间构造

/*将原图转换为32位灰度图并归一化，然后进行一次高斯平滑，并根据参数img_dbl决定是否将
图像尺寸放大为原图的2倍
参数：
img：输入的原图像
img_dbl：是否将图像放大为之前的两倍
sigma：初始高斯平滑参数σ
返回值：初始化完成的图像*/
static IplImage* create_init_img( IplImage* img, int img_dbl, double sigma )
{
    IplImage* gray, * dbl;
    float sig_diff;
    //调用函数，将输入图像转换为32位灰度图，并归一化
    gray = convert_to_gray32( img );
    if( img_dbl ) //若设置了将图像放大为原图的2倍
    {
        //将图像长宽扩展一倍时，便有了底-1层，该层尺度为：
        sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4 );
        //创建放大图像
        dbl = cvCreateImage( cvSize( img->width*2, img->height*2 ),
            IPL_DEPTH_32F,