SIFT特征提取和代码分析

本文只记录sift特征提取过程和sift的扩展应用，并分析了opensift的代码。如果想详细理解sift的理论知识请参见Rachel-Zhang的文章。这里没分析OpenCV的代码，是因为相比之下opensift代码结构更加清楚，可读性更好。

一、SIFT提取过程

对图像宽高放大1倍，并假定图像已被0.5高斯滤波过，为了达到初始为1.6高斯的效果，再用1.62−0.52−−−−−−−−−√\sqrt {{{1.6}^2} - {{0.5}^2}} 高斯滤波一次，公式参考wiki-高斯模糊；

创建高斯金字塔，塔数o=log2min(w,h)o = {\log _2}\min \left( {w,h} \right)，每塔含6层，高斯方差为σ,kσ,k2σ,k3σ,k4σ,k5σ\sigma ,k\sigma ,{k^2}\sigma ,{k^3}\sigma ,{k^4}\sigma ,{k^5}\sigma ，下一个塔时将图像宽高缩小1倍再高斯滤过；

高斯差分，对每个塔内各层做差，形成DoG图像；

对DoG图像求局部极值点并亚像素化，并用X^=−∂2D−1∂X2∂D∂X\hat X = - {{{\partial ^2}{D^{ - 1}}} \over {\partial {X^2}}}{{\partial D} \over {\partial X}}去迭代XX，这里的理论本质是利用牛顿迭代法来最优化DD，可参见李航的《统计学习方法》讲牛顿迭代法那部分；

过滤过DD值太小的点；

去除边缘响应的点，利用Hessian矩阵模和迹的关系设计滤波器，H=[DxxDxyDxyDyy]H = \left[ {\matrix{{{D_{xx}}} & {{D_{xy}}} \cr{{D_{xy}}} & {{D_{yy}}} \cr} } \right]，这里区分一下harris角点中的自相关矩阵。图像二阶导数为Dxx=Dx+1+Dx−1−2Dx{D_{xx}} = {D_{x + 1}} + {D_{x - 1}} - 2{D_x}，Dxy=(Dx+1,y+1−Dx+1,y−1−Dx−1,y+1+Dx−1,y−1)/4{D_{xy}} = ({D_{x + 1,y + 1}} - {D_{x + 1,y - 1}} - {D_{x - 1,y + 1}} + {D_{x - 1,y - 1}})/4；

计算关键点窗口下的梯度直方图，窗口大小为3×1.5σoct3 \times 1.5{\sigma _{oct}}，高斯加权中方差为1.5σoct1.5{\sigma _{oct}}，其中octoct为关键点所在塔内的尺度，注意这里的尺度已亚像素化了，并且累加梯度强度时未采用线性插值法，计算比较粗略，但在后面计算描述子时就非常精细了；

计算梯度直方图的峰值，取峰值的80%作为阈值，直方图中大于阈值的极值点会成为新关键点，并附上该方向加后到特征点列表中，所以一个关键点可以会多次加到关键点列表中，但关键点的附属方向不一样；

计算关键点旋转窗口下的梯度直方图，窗口大小为3σoct×2√×(d+1)+12{{3{\sigma _{oct}} \times \sqrt 2 \times \left( {d + 1} \right) + 1} \over 2}，其中dd表示描述子的宽度为4。计算时以关键点为中点，从(−r,−r)(-r, -r)、(−r+1,−r)(-r+1, -r)、…、(r,r)(r, r)坐标进行旋转[cosθsinθ−sinθcosθ]×[xy]\left[ {\matrix{{\cos \theta } & { - \sin \theta } \cr {\sin \theta } & {\cos \theta } \cr } } \right] \times \left[ {\matrix{x \cr y \cr } } \right]，然后[x′y′]\left[ {\matrix{{{x'}} \cr {{y'}} \cr } } \right]映射到4×44 \times 4的网格内，角度值为tan−1dydx−{\tan ^{ - 1}}{{dy} \over {dx}} - 主方向角，累加梯度强度时采用了线性插值法，最后形成4×4×8=1284 \times 4 \times 8 = 128维向量；

对特征向量归一化，可以用L2−HysL2-Hys归一化。此时所有关键点的描述子计算完成；

对所有关键点按照所在尺度（σ×2oct+kS\sigma \times {2^{oct + {k \over S}}}，其中kk已被亚像素化过）进行降序排序，这是为了方便以后关键点匹配。

二、扩展知识

SIFT算法的扩展

PCA-SIFT，其它步骤都一样，只在特征描述阶段有所区别，具体为：对每个特征点周围选择一个大小为41*41像素的像块（主方向旋转后的），计算垂直和水平的梯度，形成39*39*2=3042的向量，对所有的特征点做成一个矩阵k*3042做PCA，选取n个特征向量做成投影矩阵。一般为36维，匹配速度更快，但区分度下降，并且延长了特征的计算时间。

GLOH（Gradient Location-Orientation Histogram），也只是在特征描述阶段有所区别，具体为：把原来SIFT中4*4棋盘格的子块改成放射状的3个同心圆，对大的两个同心圆等分8份，共形成8+8+1=17的区域，直方图16维，再做PCA。区分度更高但是数据压缩花销时间太长。

SIFT的一些应用

sparse sift，做特征比对，有图像拼接的应用，图像拷贝检测；

dense sift，用BoW模型，做图像分类，目标识别（链接以后加上）；

三、代码分析

SIFT特征检测的数据结构

struct detection_data
{
int r;            // 所在图像中的行号
int c;            // 所在图像中的列号
int octv;         // 所在尺度空间的金字塔数，octv越大图像越小
int intvl;        // 所在金字塔中的层数
double subintvl;  // 亚像素化时金字塔中的层数的小数位
double scl_octv;  // sigma * pow(2, (intvl + subintval) / intvls)
};

SIFT关键点检测与描述的代码结构

int _sift_features(...)
{
init_img = create_init_img(...);    // 图像宽高放大1倍，假定图像初始化尺度为1.6
gauss_pyr = build_gauss_pyr(...);   // 建立高斯金字塔图像，存在gauss_pyr[octvs][intval+3]中，函数中为了减少使用pow()，做了优化，见补充1
dog_pyr = build_dog_pyr(...);       // 建立高斯差分金字塔
features = scale_space_extrema(...) // 寻找高斯差分空间极大值点，完成了亚像素化和边缘响应
calc_feature_scales(...);           // 计算关键点的各尺度，即关键点的尺度值=sigma*pow(2, octv + (intvl + subintvl) / intvls)
if(img_dbl)                         // 如果初始化时图像double过，则需要把关键点的坐标缩小
adjust_for_img_dbl( features );
calc_feature_oris(...);             // 计算关键点主方向，此时一个关键点有多个主方向，则会复制产生新的关键点，这里计算的直方图未插值，但最后对直方图做了平滑
compute_descriptors(...);           // 对关键点进行描述
cvSeqSort(...);                     // 对关键点按尺度降序排序
}

补充

1、计算高斯金字塔时的trick，普通情况下是利用σk=sigma×2ks\sigma _k =sigma\times 2 ^{k \over s}核去高斯模糊，而代码中用了递推式sig[i]=(kiσ)2−(ki−1σ)2−−−−−−−−−−−−−−√=ki−1σk2−1−−−−−√sig[i]=\sqrt {{(k^i\sigma)}^2-{(k^{i-1}\sigma)}^2}=k^{i-1}\sigma\sqrt{k^2-1}，所以这个样递推sig[0]=σsig[0]=\sigma、sig[1]=σk2−1−−−−−√sig[1]=\sigma\sqrt{k^2-1}、sig[i]=k×sig[i−1]sig[i]=k\times sig[i-1]