sift算法c语言实现
2024-06-08 20:58:12
前段时间在做三维测量方面的研究，需要得到物体表面三维数据，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 );//计算高斯金字塔的组数octvsoctvs = 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类型的数组featn = 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_diffcvSmooth( 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_diffcvSmooth( 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}elsefree( 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指定了关键点的最大插值次数，即最多修正多少次，默认是5while( 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 );//σ方向，即层方向//若坐标修正后超出范围，则结束插值，返回NULLif( 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 * Xcontr = interp_contr( dog_pyr, octv, intvl, r, c, xi, xr, xc );if( ABS( contr ) < contr_thr / intvls )//若该点对比度过小，舍弃，返回NULLreturn NULL;//为一个特征点feature结构分配空间并初始化，返回特征点指针feat = new_feature();//调用宏feat_detection_data来提取参数feat中的feature_data成员并转换为//detection_data类型的指针ddata = feat_detection_data( feat );//将修正后的坐标赋值给特征点feat//原图中特征点的x坐标，因为第octv组中的图的尺寸比原图小2^octv倍，//所以坐标值要乘以2^octvfeat->img_pt.x = feat->x = ( c + xc ) * pow( 2.0, octv );//原图中特征点的y坐标，因为第octv组中的图的尺寸比原图小2^octv倍，//所以坐标值要乘以2^octvfeat->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 * XcvGEMM( 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/rr越大，越可能是边缘点；伴随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 );//计算指定像素点的梯度方向直方图，返回存放直方图的数组给histhist = 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的值,返回给omaxomax = 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)的个数，默认是36
rad：区域半径，在此区域中计算梯度方向直方图
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;}//行列坐标值不合法，返回0elsereturn 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;//遍历直方图，找到最大的binfor( 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为4
n：每个直方图中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) / 2radius = 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_THRfor( 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来变为整型，并且最大值不能超过255for( 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;
}sift算法c语言实现相关推荐
- 教你一步一步用C语言实现sift算法、上
原文:http://blog.csdn.net/v_july_v/article/details/6245939 引言: 在我写的关于sift算法的前倆篇文章里头,已经对sift算法有了初步的 ...
- C语言实现寻找极值点,九之再续：教你一步一步用c语言实现sift算法、上
教你一步一步用c语言实现sift算法.上 作者:July.二零一一年三月十二日 出处:http://blog.csdn.net/v_JULY_v 参考:Rob Hess维护的sift 库 环境:win ...
- CS131-专题7：图像特征（SIFT算法）
速记要点: SIFT是什么:全称Scale Invariant Feature Transform尺度不变特征转换,2004年的论文.可以检测出图像中的局部特征点. SIFT算法特点: 稳定性:SIF ...
- 尺度不变特征变换（SIFT算法）Matlab程序代码测试例子的说明(Lowe的代码）
目前网络上可以找到的关于SIFT算法Matlab测试代码的资源就是: 1 加拿大University of British Columbia 大学计算机科学系教授 David G. Lowe发表于20 ...
- 【OpenCV】 ⚠️高手勿入! 半小时学会基本操作 24⚠️ SIFT 算法
[OpenCV] ⚠️高手勿入! 半小时学会基本操作 24⚠️ SIFT 算法 概述 图像尺度空间 多分辨率金字塔 高斯差分金字塔 计算极值点 SIFT 算法 函数 实战 概述 OpenCV 是一个跨 ...
- SIFT算法用VL_feat库实现（matlab）
sift算法是非常经典的特征提取算法,之后可以用于 对应特征匹配,从而进行图像拼接,求图像之间的转换矩阵,三维重建等工作.最近上课学习了这个算法,本打算能手敲源码,后来还是选择了调包,真香~ 毕竟前人 ...
- 特征提取与检测（二） --- SIFT算法
SIFT(Scale-invariant feature transform)是一种检测局部特征的算法,该算法通过求一幅图中的特征点(interest points,or corner points) ...
- SIFT算法原理详解
通过<图像局部不变性特征与描述>学习SIFT,遇到各种Issue,总结了这篇博客和另外九篇博客.感谢关注,希望可以互相学习,不断提升.转载请注明链接:https://www.cnblogs ...
- SURF算法与SIFT算法的性能比较——图像特征点检测与提取算法分析
图像特征点提取算法的算法研究(SURF和SIFT算法) 1. 摘要 计算机视觉中,很大一部分研究集中在图像特征提取和特征生成算法上.对图像的优化,不同于一般数学问题的优化方法,图像的优化是对像素点,在 ...
最新文章
- RDKit | 基于Lipinski规则过滤化合物库
- js文件里获取路由 vue_「如何优雅的使用Vue?」不可不知的Vue实战技巧
- Python入门100题 | 第069题
- 【收藏】deepin环境安装nodejs
- 对于ARM的启动，系统升级，烧写过程和文件系统等方面的总结分析
- 蓝桥杯--算法入门级题目及答案解析
- selenium-行为链-ActionChains-0223
- 我的个人博客live2d插件模型模块汇总（仅本人可使用，无需看）
- figma 导入导出 fig 文件
- Halcon 汉字识别
- 微信开发工具使用git
- java笔试中易考的概念
- 什么是创意啊？这才是创意
- PC 音频，视频硬件输出设置
- FBT熔融拉锥大芯径多模光纤耦合器简介
- recyclerview的条目添加点击事件
- 大数据医疗面临着哪些挑战？
- 【Spring实战】----Spring事务管理配置解析
- Ubuntu安装re2c和ninja
- 瑞萨单片机iap串口升级boot工程的构建-学习记录