SIFT特征提取和代码分析

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

一、SIFT提取过程

1. 对图像宽高放大1倍，并假定图像已被0.5高斯滤波过，为了达到初始为1.6高斯的效果，再用$\sqrt {{{1.6}^2} - {{0.5}^2}}$高斯滤波一次，公式参考wiki-高斯模糊；
2. 创建高斯金字塔，塔数$o = {\log _2}\min \left( {w,h} \right)$，每塔含6层，高斯方差为$\sigma ,k\sigma ,{k^2}\sigma ,{k^3}\sigma ,{k^4}\sigma ,{k^5}\sigma$，下一个塔时将图像宽高缩小1倍再高斯滤过；
3. 高斯差分，对每个塔内各层做差，形成DoG图像；
4. 对DoG图像求局部极值点并亚像素化，并用$\hat X = - {{{\partial ^2}{D^{ - 1}}} \over {\partial {X^2}}}{{\partial D} \over {\partial X}}$去迭代$X$，这里的理论本质是利用牛顿迭代法来最优化$D$，可参见李航的《统计学习方法》讲牛顿迭代法那部分；
5. 过滤过$D$值太小的点；
6. 去除边缘响应的点，利用Hessian矩阵模和迹的关系设计滤波器，$H = \left[ {\matrix{{{D_{xx}}} & {{D_{xy}}} \cr{{D_{xy}}} & {{D_{yy}}} \cr} } \right]$，这里区分一下harris角点中的自相关矩阵。图像二阶导数为${D_{xx}} = {D_{x + 1}} + {D_{x - 1}} - 2{D_x}$，${D_{xy}} = ({D_{x + 1,y + 1}} - {D_{x + 1,y - 1}} - {D_{x - 1,y + 1}} + {D_{x - 1,y - 1}})/4$；
7. 计算关键点窗口下的梯度直方图，窗口大小为$3 \times 1.5{\sigma _{oct}}$，高斯加权中方差为$1.5{\sigma _{oct}}$，其中$oct$为关键点所在塔内的尺度，注意这里的尺度已亚像素化了，并且累加梯度强度时未采用线性插值法，计算比较粗略，但在后面计算描述子时就非常精细了；
8. 计算梯度直方图的峰值，取峰值的80%作为阈值，直方图中大于阈值的极值点会成为新关键点，并附上该方向加后到特征点列表中，所以一个关键点可以会多次加到关键点列表中，但关键点的附属方向不一样；
9. 计算关键点旋转窗口下的梯度直方图，窗口大小为${{3{\sigma _{oct}} \times \sqrt 2 \times \left( {d + 1} \right) + 1} \over 2}$，其中$d$表示描述子的宽度为4。计算时以关键点为中点，从$(-r, -r)$、$(-r+1, -r)$、…、$(r, r)$坐标进行旋转$\left[ {\matrix{{\cos \theta } & { - \sin \theta } \cr {\sin \theta } & {\cos \theta } \cr } } \right] \times \left[ {\matrix{x \cr y \cr } } \right]$，然后$\left[ {\matrix{{{x'}} \cr {{y'}} \cr } } \right]$映射到$4 \times 4$的网格内，角度值为${\tan ^{ - 1}}{{dy} \over {dx}} -$主方向角，累加梯度强度时采用了线性插值法，最后形成$4 \times 4 \times 8 = 128$维向量；
10. 对特征向量归一化，可以用$L2-Hys$归一化。此时所有关键点的描述子计算完成；
11. 对所有关键点按照所在尺度（$\sigma \times {2^{oct + {k \over S}}}$，其中$k$已被亚像素化过）进行降序排序，这是为了方便以后关键点匹配。

二、扩展知识

• 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，普通情况下是利用$\sigma _k =sigma\times 2 ^{k \over s}$核去高斯模糊，而代码中用了递推式$sig[i]=\sqrt {{(k^i\sigma)}^2-{(k^{i-1}\sigma)}^2}=k^{i-1}\sigma\sqrt{k^2-1}$，所以这个样递推$sig[0]=\sigma$、$sig[1]=\sigma\sqrt{k^2-1}$、$sig[i]=k\times sig[i-1]$