前段时间在做三维測量方面的研究。须要得到物体表面三维数据。sift算法是立体匹配中的经典算法。以下是对RobHess的SIFT源码的凝视。部分内容參考网上，在这里向各位大神表示感谢。

头文件及函数声明

#include "sift.h"#include "imgfeatures.h"#include "utils.h"#include 
      
       #include 
       
        //将原图转换为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 * Xstatic 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和方向orientationstatic 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, 1 );        //放大原图的尺寸        cvResize( gray, dbl, CV_INTER_CUBIC );        //高斯平滑。高斯核在x,y方向上的标准差都是sig_diff        cvSmooth( dbl, dbl, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff );        cvReleaseImage( &gray );        return dbl;    }    else//不用放大为原图的2倍    {        //计算第0层的尺度        sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA );        //高斯平滑，高斯核在x,y方向上的标准差都是sig_diff        cvSmooth( gray, gray, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff );        return gray;    }}/*将输入图像转换为32位灰度图,并进行归一化參数：img：输入图像，3通道8位彩色图(BGR)或8位灰度图返回值：32位灰度图*/static IplImage* convert_to_gray32( IplImage* img ){    IplImage* gray8, * gray32;    gray32 = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );    //首先将原图转换为8位单通道图像    if( img->nChannels == 1 )//若原图本身就是单通道，直接克隆原图        gray8 = cvClone( img );    else//若原图是3通道图像    {        gray8 = cvCreateImage( cvGetSize(img), IPL_DEPTH_8U, 1 );//创建8位单通道图像        cvCvtColor( img, gray8, CV_BGR2GRAY );//将原图转换为8为单通道图像    }    //然后将8为单通道图像gray8转换为32位单通道图像，并进行归一化处理(除以255)    cvConvertScale( gray8, gray32, 1.0 / 255.0, 0 );    //释放暂时图像    cvReleaseImage( &gray8 );    //返回32位单通道图像    return gray32;}/*依据输入參数建立高斯金字塔參数：base：输入图像，作为高斯金字塔的基图像octvs：高斯金字塔的组数intvls：每组的层数sigma：初始尺度返回值：高斯金字塔，是一个octvs*(intvls+3)的图像数组*/static IplImage*** build_gauss_pyr( IplImage* base, int octvs,                                   int intvls, double sigma ){    IplImage*** gauss_pyr;    //为了保证连续性。在每一层的顶层继续用高斯模糊生成3幅图像，所以高斯金字塔每组有    //intvls+3层。DOG金字塔每组有intvls+2层    double* sig = calloc( intvls + 3, sizeof(double));//每层的sigma參数数组    double sig_total, sig_prev, k;    int i, o;    //为高斯金字塔gauss_pyr分配空间。共octvs个元素。每一个元素是一组图像的首指针    gauss_pyr = calloc( octvs, sizeof( IplImage** ) );    //为第i组图像gauss_pyr[i]分配空间。共intvls+3个元素，每一个元素是一个图像指针    for( i = 0; i < octvs; i++ )        gauss_pyr[i] = calloc( intvls + 3, sizeof( IplImage* ) );    //计算每次高斯模糊的sigma參数    sig[0] = sigma;//初始尺度    k = pow( 2.0, 1.0 / intvls );    for( i = 1; i < intvls + 3; i++ )    {        sig_prev = pow( k, i - 1 ) * sigma;//i-1层的尺度        sig_total = sig_prev * k;//i层的尺度        sig[i] = sqrt( sig_total * sig_total - sig_prev * sig_prev );    }    //逐组逐层生成高斯金字塔    for( o = 0; o < octvs; o++ )//遍历组        for( i = 0; i < intvls + 3; i++ )//遍历层        {            if( o == 0  &&  i == 0 )//第0组。第0层。就是原图像                gauss_pyr[o][i] = cvCloneImage(base);             else if( i == 0 )//新的一组的首层图像是由上一组最后一层图像下採样得到                gauss_pyr[o][i] = downsample( gauss_pyr[o-1][intvls] );             else//对上一层图像进行高斯平滑得到当前层图像            {   //创建图像                gauss_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i-1]),IPL_DEPTH_32F, 1 );                //对上一层图像gauss_pyr[o][i-1]进行參数为sig[i]的高斯平滑，得到当前层图像                //gauss_pyr[o][i]                cvSmooth( gauss_pyr[o][i-1], gauss_pyr[o][i], CV_GAUSSIAN, 0, 0, sig[i], sig[i] );            }        }        free( sig );//释放sigma參数数组        return gauss_pyr;//返回高斯金字塔}/*对输入图像做下採样生成其四分之中的一个大小的图像(每一个维度上减半)，使用近期邻差值方法參数：img：输入图像返回值：下採样后的图像*/static IplImage* downsample( IplImage* img ){    //下採样图像    IplImage* smaller = cvCreateImage( cvSize(img->width / 2, img->height / 2), img->depth, img->nChannels );    cvResize( img, smaller, CV_INTER_NN );//尺寸变换    return smaller;}/*通过对高斯金字塔中每相邻两层图像相减来建立高斯差分金字塔參数：gauss_pyr：高斯金字塔octvs：组数intvls：每组的层数返回值：高斯差分金字塔，是一个octvs*(intvls+2)的图像数组*/static IplImage*** build_dog_pyr( IplImage*** gauss_pyr, int octvs, int intvls ){    IplImage*** dog_pyr;    int i, o;    //为高斯差分金字塔分配空间，共octvs个元素，每一个元素是一组图像的首指针    dog_pyr = calloc( octvs, sizeof( IplImage** ) );    //为第i组图像dog_pyr[i]分配空间，共(intvls+2)个元素。每一个元素是一个图像指针    for( i = 0; i < octvs; i++ )        dog_pyr[i] = calloc( intvls + 2, sizeof(IplImage*) );    //逐组逐层计算差分图像    for( o = 0; o < octvs; o++ )//遍历组        for( i = 0; i < intvls + 2; i++ )//遍历层        {   //创建DoG金字塔的第o组第i层的差分图像            dog_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i]), IPL_DEPTH_32F, 1 );            //将高斯金字塔的第o组第i+1层图像和第i层图像相减来得到DoG金字塔的第o组第i层            //图像            cvSub( gauss_pyr[o][i+1], gauss_pyr[o][i], dog_pyr[o][i], NULL );        }        return dog_pyr;//返回高斯差分金字塔}

局部空间极值点检測

/*在尺度空间中检測极值点，通过插值精确定位。去除低对照度的点。去除边缘点。返回检測到的特征点序列參数：dog_pyr：高斯差分金字塔octvs：高斯差分金字塔的组数intvls：每组的层数contr_thr：对照度阈值，针对归一化后的图像。用来去除不稳定特征cur_thr：主曲率比值的阈值，用来去除边缘特征storage：存储器返回值：返回检測到的特征点的序列*/static CvSeq* scale_space_extrema( IplImage*** dog_pyr, int octvs, int intvls,                                  double contr_thr, int curv_thr, CvMemStorage* storage ){    CvSeq* features;//特征点序列    double prelim_contr_thr = 0.5 * contr_thr / intvls;//像素的对照度阈值    struct feature* feat;    struct detection_data* ddata;    int o, i, r, c;    //在存储器storage上创建存储极值点的序列，当中存储feature结构类型的数据    features = cvCreateSeq( 0, sizeof(CvSeq), sizeof(struct feature), storage );    /*遍历高斯差分金字塔，检測极值点*/    //SIFT_IMG_BORDER指明边界宽度，仅仅检測边界线以内的极值点    for( o = 0; o < octvs; o++ )//第o组        for( i = 1; i <= intvls; i++ )//遍i层            for(r = SIFT_IMG_BORDER; r < dog_pyr[o][0]->height-SIFT_IMG_BORDER; r++)//第r行                for(c = SIFT_IMG_BORDER; c < dog_pyr[o][0]->width-SIFT_IMG_BORDER; c++)//第c列                    //进行初步的对照度检查，仅仅有当归一化后的像素值大于对照度                    //阈值prelim_contr_thr时才继续检測此像素点是否可能是极值                    //调用函数pixval32f获取图像dog_pyr[o][i]的第r行第c列的点的坐标值。                    //然后调用ABS宏求其绝对值                    if( ABS( pixval32f( dog_pyr[o][i], r, c ) ) > prelim_contr_thr )                        //通过在尺度空间中将一个像素点的值与其周围3*3*3邻域内的点比較来                        //决定此点是否极值点(极大值或极小都行)                        if( is_extremum( dog_pyr, o, i, r, c ) )//若是极值点                        {                            //由于极值点的检測是在离散空间中进行的。所以检測到的极值点并                            //不一定是真正意义上的极值点                            //由于真正的极值点可能位于两个像素之间，而在离散空间中仅仅能精                            //确到坐标点精度上                            //通过亚像素级插值进行极值点精确定位(修正极值点坐标)，并去除                            //低对照度的极值点，将修正后的特征点组成feature结构返回                            feat = interp_extremum(dog_pyr, o, i, r, c, intvls, contr_thr);                            //返回值非空，表明此点已被成功修正                            if( feat )                            {                                //调用宏feat_detection_data来提取參数feat中的feature_data成员                                //并转换为detection_data类型的指针                                ddata = feat_detection_data( feat );                                //去除边缘响应，即通过计算主曲率比值推断某点是否边缘点,                                //返回值为0表示不是边缘点，可做特征点                                if( ! is_too_edge_like( dog_pyr[ddata->octv][ddata->intvl],                                     ddata->r, ddata->c, curv_thr ) )                                {                                    cvSeqPush( features, feat );//向特征点序列features末尾插                                    //入新检測到的特征点feat                                }                                else                                    free( ddata );                                free( feat );                            }                        }                        return features;//返回特征点序列}/*通过在尺度空间中将一个像素点的值与其周围3*3*3邻域内的点比較来决定此点是否极值点(极大值或极小都行)參数：dog_pyr：高斯差分金字塔octv：像素点所在的组intvl：像素点所在的层r：像素点所在的行c：像素点所在的列返回值：若指定的像素点是极值点(极大值或极小值)，返回1。否则返回0*/static int is_extremum( IplImage*** dog_pyr, int octv, int intvl, int r, int c ){    //调用函数pixval32f获取图像dog_pyr[octv][intvl]的第r行第c列的点的坐标值    float val = pixval32f( dog_pyr[octv][intvl], r, c );    int i, j, k;    //检查是否最大值    if( val > 0 )    {        for( i = -1; i <= 1; i++ )//层            for( j = -1; j <= 1; j++ )//行                for( k = -1; k <= 1; k++ )//列                    if( val < pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )                        return 0;    }    //检查是否最小值    else    {        for( i = -1; i <= 1; i++ )//层            for( j = -1; j <= 1; j++ )//行                for( k = -1; k <= 1; k++ )//列                    if( val > pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )                        return 0;    }    return 1;}

剔除不稳定点。精确定位关键点位置

/*通过亚像素级插值进行极值点精确定位(修正极值点坐标)，并去除低对照度的极值点，将修正后的特征点组成feature结构返回參数：dog_pyr：高斯差分金字塔octv：像素点所在的组intvl：像素点所在的层r：像素点所在的行c：像素点所在的列intvls：每组的层数contr_thr：对照度阈值，针对归一化后的图像，用来去除不稳定特征返回值：返回经插值修正后的特征点(feature类型)；若经有限次插值依旧无法精确到理想情况或者该点对照度过低，返回NULL*/static struct feature* interp_extremum( IplImage*** dog_pyr, int octv, int intvl,                                       int r, int c, int intvls, double contr_thr ){    struct feature* feat;//修正后的特征点    struct detection_data* ddata;//与特征检測有关的结构，存在feature结构的feature_data成员中    double xi, xr, xc, contr;//xi,xr,xc分别为亚像素的intvl(层),row(y),col(x)方向上的    //增量(偏移量)    int i = 0;//插值次数    //SIFT_MAX_INTERP_STEPS指定了关键点的最大插值次数，即最多修正多少次，默认是5    while( i < SIFT_MAX_INTERP_STEPS )    {        //进行一次极值点差值，计算σ(层方向,intvl方向)，y。x方向上的子像素偏移量(增量)        interp_step( dog_pyr, octv, intvl, r, c, &xi, &xr, &xc );        //若在随意方向上的偏移量大于0.5时，意味着差值中心已经偏移到它的临近点上。        //所以必须改变当前关键点的位置坐标        if( ABS( xi ) < 0.5  &&  ABS( xr ) < 0.5  &&  ABS( xc ) < 0.5 )//若三方向上偏移量            //都小于0.5。表示已经够精确。则不用继续插值            break;        //修正关键点的坐标。x,y,σ三方向上的原坐标加上偏移量取整(四舍五入)        c += cvRound( xc );//x坐标修正        r += cvRound( xr );//y坐标修正        intvl += cvRound( xi );//σ方向，即层方向        //若坐标修正后超出范围，则结束插值。返回NULL        if( intvl < 1  ||           //层坐标插之后越界            intvl > intvls  ||            c < SIFT_IMG_BORDER  ||   //行列坐标插之后到边界线内            r < SIFT_IMG_BORDER  ||            c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER  ||            r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER )        {            return NULL;        }        i++;    }    //若经过SIFT_MAX_INTERP_STEPS次插值后还没有修正到理想的精确位置。则返回NULL。    //即舍弃此极值点    if( i >= SIFT_MAX_INTERP_STEPS )        return NULL;    //计算被插值点的对照度：D + 0.5 * dD^T * X    contr = interp_contr( dog_pyr, octv, intvl, r, c, xi, xr, xc );    if( ABS( contr ) < contr_thr / intvls )//若该点对照度过小。舍弃，返回NULL        return NULL;    //为一个特征点feature结构分配空间并初始化。返回特征点指针    feat = new_feature();    //调用宏feat_detection_data来提取參数feat中的feature_data成员并转换为    //detection_data类型的指针    ddata = feat_detection_data( feat );    //将修正后的坐标赋值给特征点feat    //原图中特征点的x坐标，由于第octv组中的图的尺寸比原图小2^octv倍。    //所以坐标值要乘以2^octv    feat->img_pt.x = feat->x = ( c + xc ) * pow( 2.0, octv );    //原图中特征点的y坐标。由于第octv组中的图的尺寸比原图小2^octv倍，    //所以坐标值要乘以2^octv    feat->img_pt.y = feat->y = ( r + xr ) * pow( 2.0, octv );    ddata->r = r;//特征点所在的行    ddata->c = c;//特征点所在的列    ddata->octv = octv;//高斯差分金字塔中，特征点所在的组    ddata->intvl = intvl;//高斯差分金字塔中。特征点所在的组中的层    ddata->subintvl = xi;//特征点在层方向(σ方向,intvl方向)上的亚像素偏移量    return feat;//返回特征点指针}/*进行一次极值点差值，计算x，y，σ方向(层方向)上的子像素偏移量(增量)參数：dog_pyr：高斯差分金字塔octv：像素点所在的组intvl：像素点所在的层r：像素点所在的行c：像素点所在的列xi：输出參数，层方向上的子像素增量(偏移)xr：输出參数，y方向上的子像素增量(偏移)xc：输出參数，x方向上的子像素增量(偏移)*/static void interp_step( IplImage*** dog_pyr, int octv, int intvl, int r, int c,                        double* xi, double* xr, double* xc ){    CvMat* dD, * H, * H_inv, X;    double x[3] = { 0 };    //在DoG金字塔中计算某点的x方向、y方向以及尺度方向上的偏导数，结果存放在列向量dD中    dD = deriv_3D( dog_pyr, octv, intvl, r, c );    //在DoG金字塔中计算某点的3*3海森矩阵    H = hessian_3D( dog_pyr, octv, intvl, r, c );    H_inv = cvCreateMat( 3, 3, CV_64FC1 );//海森矩阵的逆阵    cvInvert( H, H_inv, CV_SVD );    cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );    //X = - H^(-1) * dD，H的三个元素各自是x,y,σ方向上的偏移量(详细见SIFT算法说明)    cvGEMM( H_inv, dD, -1, NULL, 0, &X, 0 );    cvReleaseMat( &dD );    cvReleaseMat( &H );    cvReleaseMat( &H_inv );    *xi = x[2];//σ方向(层方向)偏移量    *xr = x[1];//y方向上偏移量    *xc = x[0];//x方向上偏移量}/*在DoG金字塔中计算某点的x方向、y方向以及尺度方向上的偏导数參数：dog_pyr：高斯差分金字塔octv：像素点所在的组intvl：像素点所在的层r：像素点所在的行c：像素点所在的列返回值：返回3个偏导数组成的列向量{ dI/dx, dI/dy, dI/ds }^T*/static CvMat* deriv_3D( IplImage*** dog_pyr, int octv, int intvl, int r, int c ){    CvMat* dI;    double dx, dy, ds;    //求差分来取代偏导，这里是用的隔行求差取中值的梯度计算方法    //求x方向上的差分来近似取代偏导数    dx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) -        pixval32f( dog_pyr[octv][intvl], r, c-1 ) ) / 2.0;    //求y方向上的差分来近似取代偏导数    dy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) -        pixval32f( dog_pyr[octv][intvl], r-1, c ) ) / 2.0;    //求层间的差分来近似取代尺度方向上的偏导数    ds = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) -        pixval32f( dog_pyr[octv][intvl-1], r, c ) ) / 2.0;    //组成列向量    dI = cvCreateMat( 3, 1, CV_64FC1 );    cvmSet( dI, 0, 0, dx );    cvmSet( dI, 1, 0, dy );    cvmSet( dI, 2, 0, ds );    return dI;}/*在DoG金字塔中计算某点的3*3海森矩阵/ Ixx  Ixy  Ixs \| Ixy  Iyy  Iys |\ Ixs  Iys  Iss /參数：dog_pyr：高斯差分金字塔octv：像素点所在的组intvl：像素点所在的层r：像素点所在的行c：像素点所在的列返回值：返回3*3的海森矩阵*/static CvMat* hessian_3D( IplImage*** dog_pyr, int octv, int intvl, int r, int c ){    CvMat* H;    double v, dxx, dyy, dss, dxy, dxs, dys;    v = pixval32f( dog_pyr[octv][intvl], r, c );//该点的像素值    //用差分近似取代倒数(详细公式见各种梯度的求法)    //dxx = f(i+1,j) - 2f(i,j) + f(i-1,j)    //dyy = f(i,j+1) - 2f(i,j) + f(i,j-1)    dxx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) +         pixval32f( dog_pyr[octv][intvl], r, c-1 ) - 2 * v );    dyy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) +        pixval32f( dog_pyr[octv][intvl], r-1, c ) - 2 * v );    dss = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) +        pixval32f( dog_pyr[octv][intvl-1], r, c ) - 2 * v );    dxy = ( pixval32f( dog_pyr[octv][intvl], r+1, c+1 ) -        pixval32f( dog_pyr[octv][intvl], r+1, c-1 ) -        pixval32f( dog_pyr[octv][intvl], r-1, c+1 ) +        pixval32f( dog_pyr[octv][intvl], r-1, c-1 ) ) / 4.0;    dxs = ( pixval32f( dog_pyr[octv][intvl+1], r, c+1 ) -        pixval32f( dog_pyr[octv][intvl+1], r, c-1 ) -        pixval32f( dog_pyr[octv][intvl-1], r, c+1 ) +        pixval32f( dog_pyr[octv][intvl-1], r, c-1 ) ) / 4.0;    dys = ( pixval32f( dog_pyr[octv][intvl+1], r+1, c ) -        pixval32f( dog_pyr[octv][intvl+1], r-1, c ) -        pixval32f( dog_pyr[octv][intvl-1], r+1, c ) +        pixval32f( dog_pyr[octv][intvl-1], r-1, c ) ) / 4.0;    //组成海森矩阵    H = cvCreateMat( 3, 3, CV_64FC1 );    cvmSet( H, 0, 0, dxx );    cvmSet( H, 0, 1, dxy );    cvmSet( H, 0, 2, dxs );    cvmSet( H, 1, 0, dxy );    cvmSet( H, 1, 1, dyy );    cvmSet( H, 1, 2, dys );    cvmSet( H, 2, 0, dxs );    cvmSet( H, 2, 1, dys );    cvmSet( H, 2, 2, dss );    return H;}/*计算被插值点的对照度：D + 0.5 * dD^T * X參数：dog_pyr：高斯差分金字塔octv：像素点所在的组intvl：像素点所在的层r：像素点所在的行c：像素点所在的列xi：层方向上的子像素增量xr：y方向上的子像素增量xc：x方向上的子像素增量返回值：插值点的对照度*/static double interp_contr( IplImage*** dog_pyr, int octv, int intvl, int r,                           int c, double xi, double xr, double xc ){    CvMat* dD, X, T;    double t[1], x[3] = { xc, xr, xi };    //偏移量组成的列向量X，当中是x。y，σ三方向上的偏移量    cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );    //矩阵乘法的结果T，是一个数值    cvInitMatHeader( &T, 1, 1, CV_64FC1, t, CV_AUTOSTEP );    //在DoG金字塔中计算某点的x方向、y方向以及尺度方向上的偏导数，结果存放在列向量dD中    dD = deriv_3D( dog_pyr, octv, intvl, r, c );    //矩阵乘法：T = dD^T * X    cvGEMM( dD, &X, 1, NULL, 0, &T,  CV_GEMM_A_T );    cvReleaseMat( &dD );    //返回计算出的对照度值：D + 0.5 * dD^T * X (详细公式推导见SIFT算法说明)    return pixval32f( dog_pyr[octv][intvl], r, c ) + t[0] * 0.5;}/*为一个feature结构分配空间并初始化返回值：初始化完毕的feature结构的指针*/static struct feature* new_feature( void ){    struct feature* feat;//特征点指针    struct detection_data* ddata;//与特征检測相关的结构    feat = malloc( sizeof( struct feature ) );//分配空间    memset( feat, 0, sizeof( struct feature ) );//清零    ddata = malloc( sizeof( struct detection_data ) );    memset( ddata, 0, sizeof( struct detection_data ) );    feat->feature_data = ddata;//将特征检測相关的结构指针赋值给特征点的feature_data成员    feat->type = FEATURE_LOWE;//默认是LOWE类型的特征点    return feat;}/*去除边缘响应，即通过计算主曲率比值推断某点是否边缘点參数：dog_img：此特征点所在的DoG图像r：特征点所在的行c：特征点所在的列cur_thr：主曲率比值的阈值。用来去除边缘特征返回值：0：此点是非边缘点；1：此点是边缘点*/static int is_too_edge_like( IplImage* dog_img, int r, int c, int curv_thr ){    double d, dxx, dyy, dxy, tr, det;    /*某点的主曲率与其海森矩阵的特征值成正比，为了避免直接计算特征值，这里仅仅考虑    特征值的比值可通过计算海森矩阵的迹tr(H)和行列式det(H)来计算特征值的比值    设a是海森矩阵的较大特征值。b是较小的特征值。有a = r*b。r是大小特征值的比值    tr(H) = a + b; det(H) = a*b;    tr(H)^2 / det(H) = (a+b)^2 / ab = (r+1)^2/r    r越大，越可能是边缘点。伴随r的增大，(r+1)^2/r 的值也增大。所以可通过(r+1)^2/r 推断    主曲率比值是否满足条件*/    /* principal curvatures are computed using the trace and det of Hessian */    d = pixval32f(dog_img, r, c);//调用函数pixval32f获取图像dog_img的第r行第c列的点的坐标值    //用差分近似取代偏导，求出海森矩阵的几个元素值    /*  / dxx  dxy \    \ dxy  dyy /   */    dxx = pixval32f( dog_img, r, c+1 ) + pixval32f( dog_img, r, c-1 ) - 2 * d;    dyy = pixval32f( dog_img, r+1, c ) + pixval32f( dog_img, r-1, c ) - 2 * d;    dxy = ( pixval32f(dog_img, r+1, c+1) - pixval32f(dog_img, r+1, c-1) -        pixval32f(dog_img, r-1, c+1) + pixval32f(dog_img, r-1, c-1) ) / 4.0;    tr = dxx + dyy;//海森矩阵的迹    det = dxx * dyy - dxy * dxy;//海森矩阵的行列式    //若行列式为负，表明曲率有不同的符号，去除此点    /* negative determinant -> curvatures have different signs; reject feature */    if( det <= 0 )        return 1;//返回1表明此点是边缘点    //通过式子：(r+1)^2/r 推断主曲率的比值是否满足条件，若小于阈值，表明不是边缘点    if( tr * tr / det < ( curv_thr + 1.0 )*( curv_thr + 1.0 ) / curv_thr )        return 0;//不是边缘点    return 1;//是边缘点}/*计算特征点序列中每一个特征点的尺度參数：features：特征点序列sigma：初始高斯平滑參数，即初始尺度intvls：尺度空间中每组的层数*/static void calc_feature_scales( CvSeq* features, double sigma, int intvls ){    struct feature* feat;    struct detection_data* ddata;    double intvl;    int i, n;    n = features->total;//总的特征点个数    //遍历特征点    for( i = 0; i < n; i++ )    {        //调用宏，获取序列features中的第i个元素。并强制转换为struct feature类型        feat = CV_GET_SEQ_ELEM( struct feature, features, i );        //调用宏feat_detection_data来提取參数feat中的feature_data成员并转换为        //detection_data类型的指针        ddata = feat_detection_data( feat );        //特征点所在的层数ddata->intvl加上特征点在层方向上的亚像素偏移量。得到        //特征点的较为精确的层数        intvl = ddata->intvl + ddata->subintvl;        //计算特征点的尺度(公式见SIFT算法说明)，并赋值给scl成员        feat->scl = sigma * pow( 2.0, ddata->octv + intvl / intvls );        //计算特征点所在的组的尺度。给detection_data的scl_octv成员赋值        ddata->scl_octv = sigma * pow( 2.0, intvl / intvls );    }}/*将特征点序列中每一个特征点的坐标减半(当设置了将图像放大为原图的2倍时，特征点检測完之后调用)參数：features：特征点序列*/static void adjust_for_img_dbl( CvSeq* features ){    struct feature* feat;    int i, n;    n = features->total;//总的特征点个数    //遍历特征点    for( i = 0; i < n; i++ )    {        //调用宏，获取序列features中的第i个元素。并强制转换为struct feature类型        feat = CV_GET_SEQ_ELEM( struct feature, features, i );        //将特征点的x,y坐标和尺度都减半        feat->x /= 2.0;        feat->y /= 2.0;        feat->scl /= 2.0;        feat->img_pt.x /= 2.0;        feat->img_pt.y /= 2.0;    }}

确定关键点的大小和方向

/*计算每一个特征点的梯度直方图，找出其主方向，若一个特征点有不止一个主方向，将其分为两个特征点參数：features：特征点序列gauss_pyr：高斯金字塔*/static void calc_feature_oris( CvSeq* features, IplImage*** gauss_pyr ){    struct feature* feat;    struct detection_data* ddata;    double* hist;//存放梯度直方图的数组    double omax;    int i, j, n = features->total;//特征点个数    //遍历特征点序列    for( i = 0; i < n; i++ )    {        //给每一个特征点分配feature结构大小的内存        feat = malloc( sizeof( struct feature ) );        //移除列首元素，放到feat中        cvSeqPopFront( features, feat );        //调用宏feat_detection_data来提取參数feat中的feature_data成员并转换为        //detection_data类型的指针        //detection_data数据中存放有此特征点的行列坐标和尺度。以及所在的层和组        ddata = feat_detection_data( feat );        //计算指定像素点的梯度方向直方图，返回存放直方图的数组给hist        hist = ori_hist( gauss_pyr[ddata->octv][ddata->intvl],       //特征点所在的图像            ddata->r, ddata->c,                          //特征点的行列坐标            SIFT_ORI_HIST_BINS,                          //默认的梯度直方图的bin(柱子)个数            cvRound( SIFT_ORI_RADIUS * ddata->scl_octv ),//特征点方向赋值过程中，搜索邻域                                                                  //的半径为：3 * 1.5 * σ            SIFT_ORI_SIG_FCTR * ddata->scl_octv );       //计算直翻图时梯度幅值的高斯权重                                                                  //的初始值        //对梯度直方图进行高斯平滑。弥补因没有仿射不变性而产生的特征点不稳定的问题,        //SIFT_ORI_SMOOTH_PASSES指定了平滑次数        for( j = 0; j < SIFT_ORI_SMOOTH_PASSES; j++ )            smooth_ori_hist( hist, SIFT_ORI_HIST_BINS );        //查找梯度直方图中主方向的梯度幅值，即查找直方图中最大bin的值,返回给omax        omax = dominant_ori( hist, SIFT_ORI_HIST_BINS );        /*若当前特征点的直方图中某个bin的值大于给定的阈值，则新生成一个特征点并添        加到特征点序列末尾传入的特征点指针feat是已经从特征点序列features中移除的，        所以即使此特征点没有辅方向(第二个大于幅值阈值的方向),在函数        add_good_ori_features中也会运行一次克隆feat，对其方向进行插值修正。并插入特征        点序列的操作幅值阈值一般设置为当前特征点的梯度直方图的最大bin值的80%                   */        add_good_ori_features( features, hist, SIFT_ORI_HIST_BINS,            omax * SIFT_ORI_PEAK_RATIO, feat );        //释放内存        free( ddata );        free( feat );        free( hist );    }}/*计算指定像素点的梯度方向直方图。返回存放直方图的数组參数：img：图像指针r：特征点所在的行c：特征点所在的列n：直方图中柱(bin)的个数，默认是36rad：区域半径，在此区域中计算梯度方向直方图sigma：计算直翻图时梯度幅值的高斯权重的初始值返回值：返回一个n元数组，当中是方向直方图的统计数据*/static double* ori_hist( IplImage* img, int r, int c, int n, int rad, double sigma){    double* hist;//直方图数组    double mag, ori, w, exp_denom, PI2 = CV_PI * 2.0;    int bin, i, j;    //为直方图数组分配空间，共n个元素，n是柱的个数    hist = calloc( n, sizeof( double ) );    exp_denom = 2.0 * sigma * sigma;    //遍历以指定点为中心的搜索区域    for( i = -rad; i <= rad; i++ )        for( j = -rad; j <= rad; j++ )            //计算指定点的梯度的幅值mag和方向ori，返回值为1表示计算成功            if( calc_grad_mag_ori( img, r + i, c + j, &mag, &ori ) )            {                w = exp( -( i*i + j*j ) / exp_denom );//该点的梯度幅值权重                bin = cvRound( n * ( ori + CV_PI ) / PI2 );//计算梯度的方向相应的直方图中                                                            //的bin下标                bin = ( bin < n )?

bin : 0; hist[bin] += w * mag;//在直方图的某个bin中累加加权后的幅值 } return hist;//返回直方图数组 } /*计算指定点的梯度的幅值magnitude和方向orientation 參数： img：图像指针 r：特征点所在的行 c：特征点所在的列 img：输出參数，此点的梯度幅值 ori：输出參数，此点的梯度方向 返回值：假设指定的点是合法点并已计算出幅值和方向。返回1；否则返回0*/ static int calc_grad_mag_ori( IplImage* img, int r, int c, double* mag, double* ori ) { double dx, dy; //对输入的坐标值进行检查 if( r > 0 && r < img->height - 1 && c > 0 && c < img->width - 1 ) { //用差分近似取代偏导，来求梯度的幅值和方向 dx = pixval32f( img, r, c+1 ) - pixval32f( img, r, c-1 );//x方向偏导 dy = pixval32f( img, r-1, c ) - pixval32f( img, r+1, c );//y方向偏导 *mag = sqrt( dx*dx + dy*dy );//梯度的幅值，即梯度的模 *ori = atan2( dy, dx );//梯度的方向 return 1; } //行列坐标值不合法，返回0 else return 0; } /*对梯度方向直方图进行高斯平滑，弥补因没有仿射不变性而产生的特征点不稳定的问题 參数： hist：存放梯度直方图的数组 n：梯度直方图中bin的个数*/ static void smooth_ori_hist( double* hist, int n ) { double prev, tmp, h0 = hist[0]; int i; prev = hist[n-1]; //相似均值漂移的一种邻域平滑，降低突变的影响 for( i = 0; i < n; i++ ) { tmp = hist[i]; hist[i] = 0.25 * prev + 0.5 * hist[i] + 0.25 * ( ( i+1 == n )? h0 : hist[i+1] ); prev = tmp; } } /*查找梯度直方图中主方向的梯度幅值，即查找直方图中最大bin的值 參数： hist：存放直方图的数组 n：直方图中bin的个数 返回值：返回直方图中最大的bin的值*/ static double dominant_ori( double* hist, int n ) { double omax; int maxbin, i; omax = hist[0]; maxbin = 0; //遍历直方图。找到最大的bin for( i = 1; i < n; i++ ) if( hist[i] > omax ) { omax = hist[i]; maxbin = i; } return omax;//返回最大的bin的值 } //依据左、中、右三个bin的值对当前bin进行直方图插值，以求取更精确的方向角度值 #define interp_hist_peak( l, c, r ) ( 0.5 * ((l)-(r)) / ((l) - 2.0*(c) + (r)) ) /*若当前特征点的直方图中某个bin的值大于给定的阈值。则新生成一个特征点并加入到特征点序列末尾 传入的特征点指针feat是已经从特征点序列features中移除的。所以即使此特征点没有辅方向(第二个大于幅值阈值的方向) 也会运行一次克隆feat。对其方向进行插值修正。并插入特征点序列的操作 參数： features：特征点序列 hist：梯度直方图 n：直方图中bin的个数 mag_thr：幅值阈值。若直方图中有bin的值大于此阈值，则添加新特征点 feat：一个特征点指针，新的特征点克隆自feat。但方向不同 */ static void add_good_ori_features( CvSeq* features, double* hist, int n, double mag_thr, struct feature* feat ) { struct feature* new_feat; double bin, PI2 = CV_PI * 2.0; int l, r, i; //遍历直方图 for( i = 0; i < n; i++ ) { l = ( i == 0 )? n - 1 : i-1;//前一个(左边的)bin的下标 r = ( i + 1 ) % n;//后一个(右边的)bin的下标 //若当前的bin是局部极值(比前一个和后一个bin都大)，而且值大于给定的幅值阈值，则新生成一个特征点并加入到特征点序列末尾 if( hist[i] > hist[l] && hist[i] > hist[r] && hist[i] >= mag_thr ) { //依据左、中、右三个bin的值对当前bin进行直方图插值 bin = i + interp_hist_peak( hist[l], hist[i], hist[r] ); bin = ( bin < 0 )? n + bin : ( bin >= n )?

bin - n : bin;//将插值结果规范到[0,n]内 new_feat = clone_feature( feat );//克隆当前特征点为新特征点 new_feat->ori = ( ( PI2 * bin ) / n ) - CV_PI;//新特征点的方向 cvSeqPush( features, new_feat );//插入到特征点序列末尾 free( new_feat ); } } } /*对输入的feature结构特征点做深拷贝，返回克隆生成的特征点的指针 參数： feat：将要被克隆的特征点的指针 返回值：拷贝生成的特征点的指针 */ static struct feature* clone_feature( struct feature* feat ) { struct feature* new_feat; struct detection_data* ddata; //为一个feature结构分配空间并初始化 new_feat = new_feature(); //调用宏feat_detection_data来提取參数feat中的feature_data成员并转换为detection_data类型的指针 ddata = feat_detection_data( new_feat ); //对内存空间进行赋值 memcpy( new_feat, feat, sizeof( struct feature ) ); memcpy( ddata, feat_detection_data(feat), sizeof( struct detection_data ) ); new_feat->feature_data = ddata; return new_feat;//返回克隆生成的特征点的指针 }

生成特征描写叙述子

/*计算特征点序列中每一个特征点的特征描写叙述子向量參数：features：特征点序列gauss_pyr：高斯金字塔图像组d：计算方向直方图时。将特征点附近划分为d*d个区域，每一个区域生成一个直方图n：每一个方向直方图的bin个数*/static void compute_descriptors( CvSeq* features, IplImage*** gauss_pyr, int d, int n){    struct feature* feat;    struct detection_data* ddata;    double*** hist;//d*d*n的三维直方图数组    int i, k = features->total;//特征点的个数    //遍历特征点序列中的特征点    for( i = 0; i < k; i++ )    {        //调用宏，获取序列features中的第i个元素，并强制转换为struct feature类型        feat = CV_GET_SEQ_ELEM( struct feature, features, i );        //调用宏feat_detection_data来提取參数feat中的feature_data成员并转换为detection_data类型的指针        ddata = feat_detection_data( feat );        //计算特征点附近区域的方向直方图，此直方图在计算特征描写叙述子中要用到，返回值是一个d*d*n的三维数组        hist = descr_hist( gauss_pyr[ddata->octv][ddata->intvl], ddata->r,            ddata->c, feat->ori, ddata->scl_octv, d, n );        //将某特征点的方向直方图转换为特征描写叙述子向量，对特征描写叙述子归一化并将全部元素转化为整型。存入特征点feat中        hist_to_descr( hist, d, n, feat );        //释放特征点的方向直方图        release_descr_hist( &hist, d );    }}/*计算特征点附近区域的方向直方图，此直方图在计算特征描写叙述子中要用到，返回值是一个d*d*n的三维数组參数：img：图像指针r：特征点所在的行c：特征点所在的列ori：特征点的主方向scl：特征点的尺度d：计算方向直方图时，将特征点附近划分为d*d个区域，每一个区域生成一个直方图，默认d为4n：每一个直方图中bin的个数返回值：double类型的三维数组，即一个d*d的二维数组，数组中每一个元素是一个有n个bin的直方图数组*/static double*** descr_hist( IplImage* img, int r, int c, double ori,                            double scl, int d, int n ){    double*** hist;//d*d*n的三维直方图数组    double cos_t, sin_t, hist_width, exp_denom, r_rot, c_rot, grad_mag,        grad_ori, w, rbin, cbin, obin, bins_per_rad, PI2 = 2.0 * CV_PI;    int radius, i, j;    //为直方图数组分配空间    hist = calloc( d, sizeof( double** ) );//为第一维分配空间    for( i = 0; i < d; i++ )    {        hist[i] = calloc( d, sizeof( double* ) );//为第二维分配空间        for( j = 0; j < d; j++ )            hist[i][j] = calloc( n, sizeof( double ) );//为第三维分配空间    }    //为了保证特征描写叙述子具有旋转不变性，要以特征点为中心，在附近邻域内旋转θ角，即旋转为特征点的方向    cos_t = cos( ori );    sin_t = sin( ori );    bins_per_rad = n / PI2;    exp_denom = d * d * 0.5;    //计算特征描写叙述子过程中，特征点周围的d*d个区域中，每一个区域的宽度为m*σ个像素，        SIFT_DESCR_SCL_FCTR即m的默认值，σ为特征点的尺度    hist_width = SIFT_DESCR_SCL_FCTR * scl;    //考虑到要进行双线性插值，每一个区域的宽度应为:SIFT_DESCR_SCL_FCTR * scl * ( d + 1.0 )    //在考虑到旋转因素，每一个区域的宽度应为：SIFT_DESCR_SCL_FCTR * scl * ( d + 1.0 ) * sqrt(2)    //所以搜索的半径是：SIFT_DESCR_SCL_FCTR * scl * ( d + 1.0 ) * sqrt(2) / 2    radius = hist_width * sqrt(2) * ( d + 1.0 ) * 0.5 + 0.5;    //遍历每一个区域的像素    for( i = -radius; i <= radius; i++ )        for( j = -radius; j <= radius; j++ )        {            //坐标旋转为主方向            c_rot = ( j * cos_t - i * sin_t ) / hist_width;            r_rot = ( j * sin_t + i * cos_t ) / hist_width;            rbin = r_rot + d / 2 - 0.5;            cbin = c_rot + d / 2 - 0.5;            if( rbin > -1.0  &&  rbin < d  &&  cbin > -1.0  &&  cbin < d )                if( calc_grad_mag_ori( img, r + i, c + j, &grad_mag, &grad_ori ))                {                    grad_ori -= ori;                    while( grad_ori < 0.0 )                        grad_ori += PI2;                    while( grad_ori >= PI2 )                        grad_ori -= PI2;                    obin = grad_ori * bins_per_rad;                    w = exp( -(c_rot * c_rot + r_rot * r_rot) / exp_denom );                    interp_hist_entry( hist, rbin, cbin, obin, grad_mag * w, d, n );                }        }        return hist;}//双线性插值static void interp_hist_entry( double*** hist, double rbin, double cbin,                              double obin, double mag, int d, int n ){    double d_r, d_c, d_o, v_r, v_c, v_o;    double** row, * h;    int r0, c0, o0, rb, cb, ob, r, c, o;    r0 = cvFloor( rbin );    c0 = cvFloor( cbin );    o0 = cvFloor( obin );    d_r = rbin - r0;    d_c = cbin - c0;    d_o = obin - o0;    for( r = 0; r <= 1; r++ )    {        rb = r0 + r;        if( rb >= 0  &&  rb < d )        {            v_r = mag * ( ( r == 0 )? 1.0 - d_r : d_r );            row = hist[rb];            for( c = 0; c <= 1; c++ )            {                cb = c0 + c;                if( cb >= 0  &&  cb < d )                {                    v_c = v_r * ( ( c == 0 )? 1.0 - d_c : d_c );                    h = row[cb];                    for( o = 0; o <= 1; o++ )                    {                        ob = ( o0 + o ) % n;                        v_o = v_c * ( ( o == 0 )? 1.0 - d_o : d_o );                        h[ob] += v_o;                    }                }            }        }    }}/*将某特征点的方向直方图转换为特征描写叙述子向量，对特征描写叙述子归一化并将全部元素转化为整型，存入指定特征点中參数：hist：d*d*n的三维直方图数组d：计算方向直方图时。将特征点附近划分为d*d个区域。每一个区域生成一个直方图n：每一个直方图的bin个数feat：特征点指针，将计算好的特征描写叙述子存入当中*/static void hist_to_descr( double*** hist, int d, int n, struct feature* feat ){    int int_val, i, r, c, o, k = 0;    //遍历d*d*n的三维直方图数组，将当中的全部数据(通常是128个)都存入feat结构的descr成员中    for( r = 0; r < d; r++ )        for( c = 0; c < d; c++ )            for( o = 0; o < n; o++ )                feat->descr[k++] = hist[r][c][o];    feat->d = k;//特征描写叙述子的维数。通常是128    //归一化特征点的特征描写叙述子，即将特征描写叙述子数组中每一个元素除以特征描写叙述子的模    normalize_descr( feat );    //遍历特征描写叙述子向量。将超过阈值SIFT_DESCR_MAG_THR的元素强行赋值为SIFT_DESCR_MAG_THR    for( i = 0; i < k; i++ )        if( feat->descr[i] > SIFT_DESCR_MAG_THR )            feat->descr[i] = SIFT_DESCR_MAG_THR;    //再次归一化特征描写叙述子向量    normalize_descr( feat );    /* convert floating-point descriptor to integer valued descriptor */    //遍历特征描写叙述子向量。每一个元素乘以系数SIFT_INT_DESCR_FCTR来变为整型，而且最大值不能超过255    for( i = 0; i < k; i++ )    {        int_val = SIFT_INT_DESCR_FCTR * feat->descr[i];        feat->descr[i] = MIN( 255, int_val );    }}/*归一化特征点的特征描写叙述子，即将特征描写叙述子数组中每一个元素除以特征描写叙述子的模*/static void normalize_descr( struct feature* feat ){    double cur, len_inv, len_sq = 0.0;    int i, d = feat->d;//特征描写叙述子的维数    //求特征描写叙述子的模    for( i = 0; i < d; i++ )    {        cur = feat->descr[i];        len_sq += cur*cur;    }    len_inv = 1.0 / sqrt( len_sq );    //特征描写叙述子中每一个元素除以特征描写叙述子的模，完毕归一化    for( i = 0; i < d; i++ )        feat->descr[i] *= len_inv;}/*比較函数，将特征点按尺度的降序排列，用在序列排序函数CvSeqSort中參数：feat1：第一个特征点的指针feat2：第二个特征点的指针param：用户自己定义參数，这里不使用返回值：假设feat1的尺度大于feat2的尺度，返回1；否则返回-1；若相等返回0(好像反了)*/static int feature_cmp( void* feat1, void* feat2, void* param ){    //将输入的參数强制转换为struct feature类型的指针    struct feature* f1 = (struct feature*) feat1;    struct feature* f2 = (struct feature*) feat2;    //比較两个特征点的尺度值    if( f1->scl < f2->scl )        return 1;    if( f1->scl > f2->scl )        return -1;    return 0;}/*释放计算特征描写叙述子过程中用到的方向直方图的内存空间參数：hist：方向直方图的指针。是一个d*d*n的三维直方图数组d：直方图数组前两维的维数*/static void release_descr_hist( double**** hist, int d ){    int i, j;    for( i = 0; i < d; i++)    {        for( j = 0; j < d; j++ )            free( (*hist)[i][j] );//释放第三维的内存        free( (*hist)[i] );//释放第二维的内存    }    free( *hist );//释放第一维的内存    *hist = NULL;}/*释放金字塔图像组的存储空间參数:pyr：金字塔图像组的指针octvs：金字塔的组数n：每一组的图像数*/static void release_pyr( IplImage**** pyr, int octvs, int n ){    int i, j;    for( i = 0; i < octvs; i++ )    {        for( j = 0; j < n; j++ )            cvReleaseImage( &(*pyr)[i][j] );//释放每一个图像        free( (*pyr)[i] );//释放每一个组    }    free( *pyr );//释放金字塔    *pyr = NULL;}