RobHess的SIFT源码分析：sift.h和sift.c文件(转)

SIFT源码分析系列文章的索引在这里：

RobHess的SIFT源码分析：综述

这两个文件是RobHess的SIFT库中最重要的两个文件，里面包括用SIFT算法进行特征点检测的函数。

文件中的内容说白了很简单，就是两个特征点检测函数sift_features()和 _sift_features()，

sift_features() 是用默认参数进行特征点检测， _sift_features()允许用户输入各种检测参数，其实 sift_features() 中也是再次调用_sift_features()函数。

SIFT算法大致分为四个步骤：

██步骤一：██：建立尺度空间，即建立高斯差分(DoG)金字塔dog_pyr

██步骤二：██：在尺度空间中检测极值点，并进行精确定位和筛选

██步骤三：██：特征点方向赋值，完成此步骤后，每个特征点有三个信息：位置、尺度、方向

██步骤四：██：计算特征描述子

函数_sift_features()中的主要部分正好对应这几个步骤：

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;//原图经初始化后的图像，归一化的32位灰度图

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 );

//若设置了将图像放大为原图的2倍

if( img_dbl )//将特征点序列中每个特征点的坐标减半(当设置了将图像放大为原图的2倍时，特征点检测完之后调用)

adjust_for_img_dbl( features );

//██步骤三：██：特征点方向赋值，完成此步骤后，每个特征点有三个信息：位置、尺度、方向

//计算每个特征点的梯度直方图，找出其主方向，若一个特征点有不止一个主方向，将其分为两个特征点

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;//返回检测到的特征点的个数

}

下面是这两个文件的详细分析：

sift.h文件

#ifndef SIFT_H

#define SIFT_H

#include "cxcore.h"

//极值点检测中用到的结构

//在SIFT特征提取过程中，此类型数据会被赋值给feature结构的feature_data成员

struct detection_data

{

int r; //特征点所在的行

int c; //特征点所在的列

int octv; //高斯差分金字塔中，特征点所在的组

int intvl; //高斯差分金字塔中，特征点所在的组中的层

double subintvl; //特征点在层方向(σ方向,intvl方向)上的亚像素偏移量

double scl_octv; //特征点所在的组的尺度

};

struct feature;

//高斯金字塔每组内的层数

#define SIFT_INTVLS 3

//第0层的初始尺度，即第0层高斯模糊所使用的参数

#define SIFT_SIGMA 1.6

//对比度阈值，针对归一化后的图像，用来去除不稳定特征

#define SIFT_CONTR_THR 0.04

//主曲率比值的阈值，用来去除边缘特征

#define SIFT_CURV_THR 10

//是否将图像放大为之前的两倍

#define SIFT_IMG_DBL 1

//输入图像的尺度为0.5

#define SIFT_INIT_SIGMA 0.5

//边界的像素宽度，检测过程中将忽略边界线中的极值点，即只检测边界线以内是否存在极值点

#define SIFT_IMG_BORDER 5

//通过插值进行极值点精确定位时，最大差值次数，即关键点修正次数

#define SIFT_MAX_INTERP_STEPS 5

//特征点方向赋值过程中，梯度方向直方图中柱子(bin)的个数

#define SIFT_ORI_HIST_BINS 36

//特征点方向赋值过程中，搜索邻域的半径为：3 * 1.5 * σ

#define SIFT_ORI_SIG_FCTR 1.5

//特征点方向赋值过程中，搜索邻域的半径为：3 * 1.5 * σ

#define SIFT_ORI_RADIUS 3.0 * SIFT_ORI_SIG_FCTR

//特征点方向赋值过程中，梯度方向直方图的平滑次数，计算出梯度直方图后还要进行高斯平滑

#define SIFT_ORI_SMOOTH_PASSES 2

//特征点方向赋值过程中，梯度幅值达到最大值的80%则分裂为两个特征点

#define SIFT_ORI_PEAK_RATIO 0.8

//计算特征描述子过程中，计算方向直方图时，将特征点附近划分为d*d个区域，每个区域生成一个直方图，SIFT_DESCR_WIDTH即d的默认值

#define SIFT_DESCR_WIDTH 4

//计算特征描述子过程中，每个方向直方图的bin个数

#define SIFT_DESCR_HIST_BINS 8

//计算特征描述子过程中，特征点周围的d*d个区域中，每个区域的宽度为m*σ个像素，SIFT_DESCR_SCL_FCTR即m的默认值，σ为特征点的尺度

#define SIFT_DESCR_SCL_FCTR 3.0

//计算特征描述子过程中，特征描述子向量中元素的阈值(最大值，并且是针对归一化后的特征描述子)，超过此阈值的元素被强行赋值为此阈值

#define SIFT_DESCR_MAG_THR 0.2

//计算特征描述子过程中，将浮点型的特征描述子变为整型时乘以的系数

#define SIFT_INT_DESCR_FCTR 512.0

//定义了一个带参数的函数宏，用来提取参数f中的feature_data成员并转换为detection_data格式的指针

#define feat_detection_data(f) ( (struct detection_data*)(f->feature_data) )

extern int sift_features( IplImage* img, struct feature** feat );

extern 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 );

#endif

sift.c文件

#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 * 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 )

{

//调用_sift_features()函数进行特征点检测

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 );

}

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;//原图经初始化后的图像，归一化的32位灰度图

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 );

//若设置了将图像放大为原图的2倍

if( img_dbl )//将特征点序列中每个特征点的坐标减半(当设置了将图像放大为原图的2倍时，特征点检测完之后调用)

adjust_for_img_dbl( features );

//██步骤三：██：特征点方向赋值，完成此步骤后，每个特征点有三个信息：位置、尺度、方向

//计算每个特征点的梯度直方图，找出其主方向，若一个特征点有不止一个主方向，将其分为两个特征点

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;//返回检测到的特征点的个数

}

static IplImage* create_init_img( IplImage* img, int img_dbl, double sigma )

{

IplImage* gray, * dbl;

float sig_diff;

//调用函数，将输入图像转换为32位灰度图，并归一化

gray = convert_to_gray32( img );

//若设置了将图像放大为原图的2倍

if( img_dbl )

{

//将图像长宽扩展一倍时，便有了底-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;

}

}

static IplImage* convert_to_gray32( IplImage* img )

{

IplImage* gray8, * gray32;

gray32 = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );//创建32位单通道图像

//首先将原图转换为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 );//释放临时图像

return gray32;//返回32位单通道图像

}

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;//返回高斯金字塔

}

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;

}

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;//返回高斯差分金字塔

}

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;//返回特征点序列

}

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;

}

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;//返回特征点指针

}

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方向上偏移量

}

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;

}

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;

}

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;

}

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;

}

static int is_too_edge_like( IplImage* dog_img, int r, int c, int curv_thr )

{

double d, dxx, dyy, dxy, tr, det;

d = pixval32f(dog_img, r, c);//调用函数pixval32f获取图像dog_img的第r行第c列的点的坐标值

//用差分近似代替偏导，求出海森矩阵的几个元素值

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;//海森矩阵的行列式

//若行列式为负，表明曲率有不同的符号，去除此点

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;//是边缘点

}

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 );

}

}

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;

}

}

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 );

add_good_ori_features( features, hist, SIFT_ORI_HIST_BINS,

omax * SIFT_ORI_PEAK_RATIO, feat );

//释放内存

free( ddata );

free( feat );

free( hist );

}

}

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;//返回直方图数组

}

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;

}

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;

}

}

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)) )

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 );

}

}

}

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;//返回克隆生成的特征点的指针

}

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 );

}

}

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;

}

}

}

}

}

}

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 );

//遍历特征描述子向量，每个元素乘以系数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;

}

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;

}

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;

}

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;

}

http://www.tuicool.com/articles/jEJVja