RobHess的SIFT代码解析步骤三

原创

@蓑衣客ʷᵖ 2022-12-26 15:52:31 ©著作权

文章标签 特征点直方图邻域 文章分类 虚拟化云计算

©著作权归作者所有：来自51CTO博客作者@蓑衣客ʷᵖ的原创作品，请联系作者获取转载授权，否则将追究法律责任

平台：win10 x64 +VS 2015专业版 +opencv-2.4.11 + gtk_-bundle_2.24.10_win32

主要参考：1.代码：RobHess的SIFT源码

2.书：王永明王贵锦《图像局部不变性特征与描述》

SIFT四步骤和特征匹配及筛选：

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

步骤二：在尺度空间中检测极值点，并进行精确定位和筛选创建默认大小的内存存储器

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

步骤四：计算特征描述子

SIFT后特征匹配：KD树+BBF算法

SIFT后特征匹配后错误点筛选：RANSAC算法

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

问题及解答：

（1）问题描述：怎么计算特征点的方向和幅值？思路是什么？

答：为了实现图像旋转的不变性，需要根据检测到的特征点的局部图像结构求得一个方向基准。我们使用图像梯度的方法求取该局部结构的稳定方向。对于己经检测到特征点，我们知道该特征点的尺度值σ，因此根据这一尺度值，在GSS中得到最接近这一尺度值的高斯图像。然后使用有限差分，计算以特征点为中心，以3X1.5σ为半径的区域内图像梯度的幅角和幅值，如下图所示。幅角和幅值计算公式加下：

RobHess的SIFT代码解析步骤三_邻域

参考书P130及图6-1

（2）问题描述：为什么要计算特征梯度直方图？怎么计算？

答：在完成特征点邻域的高斯图像的梯度计算后，使用直方图统计邻域内像素的梯度方向和幅值。梯度方向直方图的横轴是梯度方向角，纵轴是梯度方向角对应的（带高斯权重）梯度幅值累加值。梯度方向直方图将。0°~360°的范围，分为36个柱，每10°为一个柱。直方图的峰值代表了该特征点处邻域内图像梯度的主方向，也即该特征点的主方向，如下图所示：

RobHess的SIFT代码解析步骤三_邻域_02

绿色格点代表邻域范围，蓝色圆圈代表格点的高斯权重（稍后介绍），黑色箭头指向代表梯度方向，箭头长度代表梯度幅值。右边为梯度方向直方图（36柱，每柱代表10°，上图只显示了8柱）。获得梯度方向直方图的步骤如下：

1）生成领域各像元的高斯权重。其中高斯函数方差为该特征点的特征尺度σ的1.5倍。形式如下，其中（i,j）为该点距离特征点的相对位置，如上图，左上角点像元距离特征点（0,0）（即中心点）的相对位置坐标为（-4,-4），同理，右下角像元为（4,4）。

RobHess的SIFT代码解析步骤三_直方图_03

2）遍历邻域（绿色）中每个点，判断其梯度方向，将其加入相应的梯度方向直方图中，加入量为其梯度幅值 * wi,j ，例如左上角(-4,-4)的点，其梯度为方向为25°，梯度幅值为mag，我们将其加入到hist[2]中（假设hist[0]为0°~10°的直方柱，hist[1]为10°~20°的直方柱，以此类推至hist[35]为350°~360°）。加入的量为mag* w(-4,-4)，即hist[2] = hist[2] + mag* w(-4,-4)。直至遍历整个邻域，统计出该特征点出的梯度方向直方图。

3）平滑直方图。对上一步得出的直方图进行平滑，得到最终的梯度方向直方图。OpenCV中使用的 (1/16) * [1,4,6,4,1] 的高斯卷积和对直方图进行平滑处理，而vlfeat中使用了6次，邻域大小为3的平均处理，即hist[i] = (hist[i-1]+hist[i]+hist[i+1])/3。

参考书P131及图6-2

（3）问题描述：为什么每个点梯度幅值要使用高斯权重？
答：由于SIFT算法只考虑了尺度和旋转的不变性，并没有考虑仿射不变性。通过对各点梯度幅值进行高斯加权，使特征点附近的梯度幅值有较大的权重，这样可以部分弥补因没有仿射不变性而产生的特征点不稳定的问题。

（4）问题描述：如何确定特征点方向？

答：

有了梯度方向直方图之后，找到直方图中最大的值，则认为该方向为该特征点的主方向，如存在另一个方向大于最大值的80％，则认为该方向为该特征点的辅方向。一个特征点可能会有多个方向（一个主方向，一个以上的辅方向），这可以增强匹配的鲁棒性。具体而言，就是将该特征点复制成多份特征点（除了方向θ不同外，x,y,σ都相同）。

注意：在OpenCV中，若辅方向除了满足大于最大值80％外，还必须是局部最大值，即 hist[i] > hist[i-1] && hist[i] > hist[i+1]。

RobHess的SIFT代码解析步骤三_直方图_04

　　通常离散的梯度方向直方图，可以通过插值拟合处理，这样可以得到更精确的方向角度值。

RobHess的SIFT代码解析步骤三_特征点_05

　　经过上述过程，我们特征点的所有量（x,y,σ,θ）都已经求得，其中位置（x,y）、尺度σ都是在上一节中求得，而特征点方向θ是通过特征点邻域直方图求得。下一节，将介绍SIFT描述子的形成方式。

参考书P132-P133及图6-5和图6-6

（5）问题描述：RobHess的SIFT源码如何实现步骤三，大概思路是这样的？

答：（5.1）代码及说明：

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

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

calc_feature_oris( features, gauss_pyr );

（5.2）calc_feature_oris代码及说明：

确定关键点的大小和方向

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

为两个特征点

参数：

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

}

（5.2.1）cvSeqPopFront函数说明：

答：动态结构序列CvSeq是所有OpenCV动态数据结构的基础。
分为两类：
稠密序列
稀疏序列
（1）稠密序列都派生自CvSeq，他们用来代表可扩展的一维数组 — 向量、栈、队列和双端队列。数据间不存在空隙（连续存储）。如果元素元素从序列中间被删除或插入新的元素到序列，那么此元素后边的相关元素全部被移动。
（2）稀疏序列派生自CvSet，CvSet也是基于CvSeq的，他们都是由节点所组成，每一个节点要么被占用，那么为空，由标志位flag决定。这些序列作为无序数据结构被使用，如点集合、图、Hash表等。

CvSeq的定义如下：

CV_TREE_NODE_FIELDS(CvSeq);
    int       total;          //total表示稠密序列的元素个数，或者稀疏序列被分配的节点数。
    int       elem_size;      //elem_size表示序列中每个元素占用的字节数。
    schar*    block_max;      //block_max是最近一个内存的最大边界指针。
    schar*    ptr;            //ptr表示当写指针。
    int       delta_elems;    //delta_elems表示序列间隔尺寸。
    CvMemStorage* storage;    //storage指向序列存储的内存块的指针
    CvSeqBlock* free_blocks; //free_blocks表示空的块列表。
    CvSeqBlock* first; //first指向第一个序列块。h_next和h_prev并不是指向CvSeq内部元素的指针，它们是指向其它CvSeq的。

//头部添加
Char* cvSeqPushFront(CvSeq* seq,void* element=NULL);

//头部删除 Void cvSeqPopFront(CvSeq* seq,void* element=NULL);

参看：OpenCV——CvSeq动态结构序列—

（5.2.2）ori_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;

//遍历以指定点为中心的搜索区域

范围为2*3*1.5σ，参看书P131图6-2

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下标

　　　　　//ori范围(-Π，Π]，+Π后范围(0，2Π]

bin = ( bin < n )? bin : 0; //把bin的下标归到0~36

hist[bin] += w * mag;//在直方图的某个bin中累加加权后的幅值

}

return hist;//返回直方图数组

}

（5.2.2.1）calc_grad_mag_ori代码及说明：

答：

/*计算指定点的梯度的幅值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;

}

问题1：具体如何求梯度？为什么一阶偏导没有二分之一的系数？

答：差分代替偏导，一阶情况，具体参考：RobHess的SIFT代码解析步骤二（5.2.6.2.1.1）

RobHess的SIFT代码解析步骤三_邻域_06

问题2：为什么计算梯度的方向使用的是atan2，而不是（1）中介绍的arctan（atan）？

答：

说明：C 语言里 double atan2(double y,double x) 返回的是原点至点(x,y)的方位角，即与 x 轴的夹角。返回值的单位为弧度，取值范围为（-π, π]。结果为正表示从 X 轴逆时针旋转的角度，结果为负表示从 X 轴顺时针旋转的角度。

atan2与atan的区别：
atan接受的是一个正切值（直线的斜率）得到夹角，但是由于正切的规律性本可以有两个角度的但它却只返回一个，因为atan的值域是从-90~90 也就是它只处理一四象限，所以一般不用它。
atan2(double y,double x) 其中y代表已知点的Y坐标同理x ,返回值是此点与远点连线与x轴正方向的夹角，这样它就可以处理四个象限的任意情况了，它的值域相应的也就是-180~180了。
例如：
例1：斜率是1的直线的夹角
cout<<atan(1.0)*180/PI;//45°
cout<<atan2(1.0,1.0)*180/PI;//45° 第一象限
cout<<atan2(-1.0,-1.0)*180/PI;//-135°第三象限
后两个斜率都是1 但是atan只能求出一个45°

例2：斜率是-1的直线的角度
cout<<atan(-1.0)*180/PI;//-45°
cout<<atan2(-1.0,1.0)*180/PI;//-45° y为负在第四象限
cout<<atan2(1.0,-1.0)*180/PI;//135° x为负在第二象限
后两个斜率都是-1 但是atan只能求出一个-45°

此处与书P130页不同，atan的值域[ -Π/2，Π/2]，而atan2的值域[ -Π，Π]，涵盖了一周360度（2Π）的范围。

参看博客：atan2相关知识汇总

（5.2.2.2）梯度直方图某点梯度的幅值怎么计算？不是直接累加的，而是乘以一个高斯权重1.5sigma。

答：L为关键点所在的尺度空间值，按Lowe的建议，梯度的模值m(x,y)按 σ=1.5σ_oct 的高斯分布加成，按尺度采样的3σ原则，领域窗口半径为 3x1.5σ_oct。

在完成关键点的梯度计算后，使用直方图统计领域内像素的梯度和方向。梯度直方图将0~360度的方向范围分为36个柱(bins)，其中每柱10度。如图5.1所示，直方图的峰值方向代表了关键点的主方向，(为简化，图中只画了八个方向的直方图)。

RobHess的SIFT代码解析步骤三_直方图_07

参看博客

（5.2.3）smooth_ori_hist代码及说明：

答：

/*对梯度方向直方图进行高斯平滑，弥补因没有仿射不变性而产生的特征点不稳定的问题

参数：

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;

}

问题1：为什么要进行高斯平滑？

答：为了防止某个梯度方向角度因受到噪声的干扰而突变，我们还需要对梯度方向直方图进行平滑处理。Opencv 所使用的平滑公式为：

RobHess的SIFT代码解析步骤三_邻域_08

1）今天的内容是是平滑，先介绍均值和高斯平滑，高斯平滑是均值平滑的一个扩展，或者是一个进化版本。均值的原理是，一个规定的邻域内，所有像素的平局值作为最终计算的结果，每个像素的权值相同，为总像素的倒数，而高斯平滑是上述的升级版本，邻域内每个像素的权值不同，而权值是由高斯函数确定的。

均值平滑和高斯平滑都是线性的，也就是，一旦参数给定，模板就确定下来，不会因为位置和像素分布不同而改变，而线性模板的基本运算是卷积。

线性模板的另一个性质就是可以进行频域分析，比如高斯平滑的频域仍然是高斯的，而且是低通的，这与前面讲到的平滑是消除尖锐的噪声（高频）的操作相互证明了其正确性，均值平滑是一个盒状滤波器，其频域为sinc函数，也是一个近似的低通，但sinc函数有旁瓣，所以，模板宽度选择不好可能会有一些不好的效果，比如有些高频会被保留，这一点也比较好理解。

比较下两种均值（加权和不加权）。比如一维的一个序列{0，0，0，0，0，1000，0，0，0，0}，明显1000是个边缘，如果使用3个宽度的均值平滑，结果是{0，0，0，0，333，333，333，0，0，0}，边缘被完全模糊掉了。但如果使用{1，2，1}的近似高斯平滑模板，结果是{0，0，0，0，250，500，250，0，0，0}，边缘被保留。所以，加权平均（高斯）可以保留一定的细节。

参看博客：灰度图像--图像增强平滑之均值滤波、高斯滤波——http://www.mamicode.com/info-detail-447178.html

2）LOWE建议对直方图进行平滑，降低突变的影响，弥补因没有仿射不变性而产生的特征点不稳定的问题

问题2：对梯度直方图进行高斯平滑，为什么用的是1/4{1，2，1}？盒子滤波器的模板吗?

答：用的是1/4{1，2，1}不是盒子滤波器的模板，1/4{1，2，1}是高斯加权平均。.2 平滑直方图 lowe建议对直方图进行平滑，降低突变的影响。

（5.2.4）dominant_ori代码及说明：

答：

/*查找梯度直方图中主方向的梯度幅值，即查找直方图中最大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的值

}

（5.2.5）add_good_ori_features代码及说明：

答：

/*若当前特征点的直方图中某个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;//新特征点的方向

方向还原为（Π，Π]，bin最大值36，2Π*36/36-Π=Π，bin最大值0，2Π*0/36-Π=-Π

cvSeqPush( features, new_feat );//插入到特征点序列末尾

free( new_feat );

}

问题1：插值时小柱子到边界怎么处理，怎么插值？因为抛物线插值需要利用三个点的坐标，如0之前和35之后怎么处理？

答：由于角度是循环的，即00=3600，如果出现h(j)，j超出了(0,…,35)的范围，那么可以通过圆周循环的方法找到它所对应的、在00=3600之间的值，如h(-1) = h(35)。把小柱子看成是循环的，角度的取值为0~360度是一个圆周，0之前利用标号为35的小柱子和35之后利用标号为0的小柱子。

3.2.1、梯度图像的平滑处理

问题2：对直方图采用什么插值？插值需要利用几个点的坐标？

答：采用抛物线插值，需要三个点的坐标，hist[l], hist[i], hist[r]，传入interp_hist_peak函数。

（5.2.5.1）interp_hist_peak代码及说明：

答：

//根据左、中、右三个bin的值对当前bin进行直方图插值，以求取更精确的方向角度值

#define interp_hist_peak( l, c, r ) ( 0.5 * ((l)-(r)) / ((l) - 2.0*(c) + (r)) )

问题：梯度直方图抛物线插值小柱子在直方图中的索引号的公式怎么计算？

RobHess的SIFT代码解析步骤三_特征点_09