前段时间在做三维測量方面的研究。须要得到物体表面三维数据。sift算法是立体匹配中的经典算法。以下是对RobHess的SIFT源码的凝视。部分内容參考网上，在这里向各位大神表示感谢。
http://blog.****.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, 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)的个数，默认是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;
}
//行列坐标值不合法，返回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为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) / 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;
}