原创作品，允许转载，转载时请务必以超链接形式标明文章 原始出处 、作者信息和本声明。否则将追究法律责任。http://underthehood.blog.51cto.com/2531780/658350
By RaySaint 2011/09/05
1综述
结合论文[1]和Rob Hess的开源SIFT代码(发现OpenCV2.3的源码里也是用的Rob Hess的SIFT代码)对SIFT算法进行了研究，下面是小结：
在计算机视觉的领域中，图像匹配是很多问题最重要的一个方面，包括物体和场景识别，通过多幅图像进行3D重构，立体匹配和运动跟踪。SIFT特征对于图像的旋转和尺度变化具有不变性，对于光照改变和摄像机角度变化具有部分的不变性。SIFT算法生成图像特征的主要步骤有以下几个：
（1）尺度空间极值检测：搜索所有尺度上的图像位置。通过高斯微分函数来识别潜在的对于尺度和旋转不变的兴趣点。
（2）关键点的定位：在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。关键点的选择依据于它们的稳定程度。
（3）方向的确定：基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向。所有后面的对图像数据的操作都相对于关键点的方向、尺度和位置进行变换，从而提供对于这些变换的不变性。
（4）关键点描述：在每个关键点周围的邻域内，在选定的尺度上测量图像局部的梯度。这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变化。
2详细过程
这里先根据Rob Hess的代码简化一个描述SIFT特征的结构体：
struct feature
{
double x; /*!< 在输入图像上的x坐标 */
double y; /*!< 在输入图像上的y坐标 */
double scl; /*!< 关键点的尺度 */
double ori; /*!< 关键点的方向 */
int d; /*!< 描述器的长度 */
double descr[FEATURE_MAX_D]; /*!< 描述器 */
int r; /*!< 检测到关键点那一层图像上关键点所在的行 */
int c; /*!< 检测到关键点那一层图像上关键点所在的列 */
int octv; /*!< 检测到关键点那一层图像所在的组的索引*/
int intvl; /*!< 检测到关键点那一层图像所在的层的索引 */
double subintvl; /*!< 检测到关键点那一层图像所在的层的插值 */
double scl_octv; /*!< 检测到关键点那一层图像所在的组的尺度 */
};
2.1尺度空间表示
图像的尺度空间表示成一个函数L(x, y, σ)，它是由一个变化尺度的高斯函数G(x, y, σ)与图像I(x, y)卷积生成的，见Lindeberg.T的论文[2]。即
其中
代表卷积操作，而G(x, y, σ)为：
SIFT算法建议在某一个尺度上的对关键点的检测，可以通过对两个相邻高斯尺度空间的图像相减，得到一个DoG（高斯微分）的响应值图像D(x, y, σ)。然后，通过对响应值图像D(x, y, σ)进行局部极大搜索，在位置空间和尺度空间中定位关键点。其中：
式中，k为两相邻尺度空间倍数的常数。
事实上DoG是对LoG(高斯拉普拉斯)的一个近似。根据热传导方程（注：这里t=σ^2）：
对上式进行有限差分运算，得到：
由于常数k-1并不会影响极值点的位置，即关键点中心的位置。所以，D(x, y, σ) 是
的近似表示，也即DoG是LoG的近似表示。
接下来是在如何计算DoG。首先是构建图像的高斯金字塔。将图像金字塔共分O组，一组称为一个Octave，每组又分为多层，层间隔数为S，因此有S+3(S+1+2，2代表在上下再各添一层图像，搜索极值只在中间的S+1层图像上搜索)层图像，下一组的第一层图像由上一组的倒数第三层（如果层索引从0开始，则为第S层）图像按照隔点采样得到，这样做的目的是为了减少卷积运算的工作量。DoG是通过高斯金字塔中的每组上下相邻两层的高斯尺度空间图像相减得到。如下图：
在高斯金字塔中，σ和o, s的关系如下：
其中，σ0是基准层尺度，o为组Octave的索引，s为组里图像的层索引。
我们认为初始图像有一个初始的σ=0.5的高斯模糊（为了防止显著的混淆效应）。从上式可以看出，（2）式中的k=2^(1/S)。
在最开始建立高斯金字塔时，要预先模糊输入图像来作为第0个组的第0层的图像，这时相当于丢弃了最高的空域的采样率。因此通常的做法是先将图像的尺度扩大一倍来生成第0组。我们假定初始的输入图像为了抗击混淆现象，已经对其进行σn=0.5的高斯模糊，如果输入图像的尺寸用双线性插值扩大一倍，那么相当于σn=1.0。因此，实际
上面用到了高斯的半群性质得到，即：
其中t=σ2。
在实际的编程实现中（2）式如下实现：
在D.Lowe的论文和Rob Hess的代码中，使用的默认参数为：
σn=0.5, σ0=1.6, S=2, O的值通过对图像的长和宽中较小的一个长度取2的对数再减去2得到（保证最后一组的图像至少要有4个像素）。
下面是尺度空间表示的一个例子：
2.2尺度空间局部极值检测（关键点的初步探查）
兴趣点的初步探查是通过同一组内各DoG相邻层之间比较完成的。为了寻找尺度空间的极值点，每一个采样点要和它的相邻点进行比较，看其是否比它的图像域和尺度域的相邻点大或者小。如下图所示：
上图中表示×的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点进行比较，以确保在尺度空间和二维图像位置空间都检测到极值点。搜索过程从每组的第二层开始，也就是从层索引为1的那一层开始，到层索引为S的层结束。
下面是局部极值检测的结果：
2.3局部极值位置的插值（关键点的精确定位）
以上极值点的搜索是在离散空间中进行的，检测到的极值点并不是真正意义上的极值点。下图演示了二维函数离散空间得到的极值点与连续空间极值点的差别。利用已知的离散空间点插值得到的连续空间极值点的方法叫做子像素插值（Sub-pixel Interpolation）。
插值的方法是根据泰勒级数展开，具体过程参见[3]，这里只给出结论：
极值点位置的偏移量（位置矢量为）：
极值点的极值为
其中
代表相对于插值中心的偏移量，当它在任一维度上的偏移量大于0.5时（即x或y或σ），意味着插值中心已经偏移到它的邻近点上，所以必须改变当前关键点的位置（加上
,实际编成中对
进行四舍五入取整），同时在新的位置上反复插值直到收敛；也有可能超出所设定的迭代次数或者超出图像边界的范围，此时这样的点应该删除。另外，当
的绝对值小于0.03时（假设图像的灰度值是在0到1.0之间），其响应值过小，这样的点易受噪声的干扰而变得不稳定，所以也要删除。
下面是移除
绝对值较小的关键点的结果：
2.4消除边界响应（关键点的进一步筛选）
由于DoG对图像中的边缘有比较强的响应值，而一旦特征落在图形的边缘上，这些点就是不稳定的点。根据Harris[4]角点可以知道，一个角点在任何方向上平移都应该保证局部窗口内的像素值的剧烈变化，而边缘上的点沿着边缘方向移动时局部窗口内的像素值基本没有什么变化。如下图所示：
同样，一个平坦的DoG响应峰值往往在横跨边缘的地方有较大的主曲率，而在垂直的方向有较小的主曲率。而主曲率可以通过2×2的Hessain矩阵H求出：
H的特征值与D的主曲率成正比例。可以避免求取具体的特征值，因为我们只关心特征值的比例。令α为较大的特征值，β为较小的特征值。
那么
如果γ为较大特征值与较小特征值之间的比值，α=γβ，这样便有
上式的结果只与两个特征的比例有关，而与具体特征值无关。两个特征值相等时，
最小，随着γ增大，其值也增加。所以，要想检查主曲率的比值小于某一阈值γ，只要检查下式是否成立：
D.Lowe在论文中给出γ=10，对于主曲率比值大于10的特征点将被删除，否则，这些特征点将被保留。下面是移除边界响应的结果：
到此为止，特征点的位置已经被确定，假设在D(x,y,σ(o,s))上位置为(xk,yk)处确定了一个特征点，在插值后的偏移Xo=(xo,yo,so)(注：尺度σ的导数在DoG金字塔中是根据同一组的图像前后差分求得，也就是s+1层图像减去s-1层图像再除以2，因此这里用so替代σo)，那么先前给出的描述SIFT特征的结构体将如下被填充：
feat.x = (xk+xo) × 2^o
feat.y = (yk+yo) × 2^o
feat.r = yk;
feat.c = xk;
feat.octv = o;
feat.intvl = s;
feat.subintvl = so;
（注：如果最先开始图像的尺寸被加倍了，那么feat.x和feat.y要除以2）
2.5关键点的尺度计算
这一步计算出关键点的尺度，根据2.1节给出的公式(6)，如下：
（注：如果最先开始图像的尺寸被加倍了，那么feat.scl要除以2）
2.6特征点的方向分配
为了实现图像旋转的不变性，需要根据检测到的特征点的局部图像结构求得一个方向基准。我们使用图像梯度的方法求取局部结构的稳定方向。对于已经检测到的特征点，我们知道该特征点在DoG金字塔中的确切位置和尺度（由feat.r, feat.c, feat.octv, feat.intvl, feat.scl_octv确定），在与尺度相应的高斯图像L(x,y,σ(feat.octv, feat.intvl))上使用有限差分，计算以特征点为中心，以3×1.5×feat.scl_octv为半径的区域内图像梯度的幅角和幅值（模值），如下图所示，幅角和幅值的计算公式如下：
在完成特征点邻域的高斯图像的梯度计算后，使用直方图统计邻域内像素的梯度方向和幅值。梯度方向直方图的横轴是梯度方向角，纵轴是梯度方向角对应的梯度幅值累加。梯度方向直方图将0-360度的范围分为36个柱(bins)，每10度一个柱。直方图的峰值代表该特征点处邻域内图像梯度的主方向，也即该特征点的主方向，如下图所示：
每个累加到梯度方向直方图的采样点的梯度幅值都要进行权重处理，加权采用圆形高斯加权函数，其标准偏差σ为1.5×feat.scl_octv，（上面的局部邻域半径通过3σ定理得到）。如下图所示，由于SIFT算法只考虑了尺度和旋转不变性，并没有考虑仿射不变性。通过高斯加权，使特征点附近的梯度幅值有较大的权重，这样可以部分弥补因没有仿射不变性而产生的特征点不稳定的问题。
另外，当存在一个相当于主峰值能量80%能量的峰值时，则将这个方向认为是该特征点的辅方向。一个特征点可能会被指定具有多个方向（一个主方向，一个以上辅方向），这可以增强匹配的鲁棒性，如下图所示。实际编程实现中，就是把该特征点复制成多份特征点，并将方向值分别赋给这些复制后的特征点，并且，离散的梯度方向直方图要进行插值拟合处理，来求得更精确的方向角度值。
下面是特征方向和尺度分配的结果（箭头指向代表方向，长度代表尺度）：
2.7特征描述子表示（特征点的特征矢量生成）
由SIFT描述子h(x,y,θ)是对特征点附近邻域高斯图像梯度统计结果的一种表示，它是一个三维的阵列，但通常在编程实现时只表示成一个矢量（1维向量）。特征描述子与特征点所在的尺度有关，因此，对梯度的求取应在特征点对应的高斯图像上进行。将特征点附近的邻域划分成Bp*Bp个子区域，每个子区域的尺寸为m×feat.scl_octv个子像素，其中，m=3，Bp=4。考虑到实际计算时，需要采用双线性插值，计算的图像区域为m × feat.scl_octv × (Bp+1)。如果再考虑旋转因素，那么实际计算的图像区域应为m × feat.scl_octv × (Bp+1) × √2。如下图所示：
为了保证特征矢量具有旋转不变形，需要以特征点为中心，将特征点附近邻域内(m × feat.scl_octv × (Bp+1) × √2 × m × feat.scl_octv × (Bp+1) × √2)图像的梯度的位置和方向旋转一个方向角feat.ori。如下图所示：
在特征点附近邻域图像梯度的位置和角度旋转后，再以特征点为中心，取一个m × feat.scl_octv × Bp × m × feat.scl_octv × Bp大小的图像区域，并将它等间隔划分成Bp×Bp个子区域，每个间隔为m × feat.scl_octv像素，在每个子区域内计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，形成一个种子点。求与特征点主方向时有所不同，此时，每个子区域的梯度方向直方图将0-360度的范围划分为8个方向范围，每个范围为45度，这样，每个种子点共有8个方向的梯度强度信息。由于存在4×4（Bp × Bp）个子区域，所以，共有4×4×8=128个数据，即特征向量的长度feat.d=128。如下图所示：
同样，要对特征矢量需要信息高斯加权处理，加权采用标准偏差为m × feat.scl_octv × Bp / 2的标准高斯函数。特征矢量形成后，为了去除光照变化的影响，需要对它们进行归一化处理。在归一化处理后，对于特征矢量中值大于0.2的要进行截断，即大于0.2的值只取0.2，然后，再进行一次归一化处理，其目的是提高特征的鉴别性。
3疑问
(1)在每次需要确定特征点邻域半径大小时，使用的尺度并不是公式(6)确定的尺度，而是去掉所在组序号o计算出的尺度，这是为什么？
(2)前面说到
的绝对值小于0.3的点要被舍弃，在Rob Hess的SIFT代码中，实际上是
的绝对值小于0.3/S的点被舍弃，为何这个对比度阈值要除以S？
4参考文献
[1]D.Lowe, "Distinctive Image Features from Scale-Invariant Keypoints", January 5, 2004
[2]T.Lindeberg, "Scale-space theory:A basic tool for analysing structures at different scales"
[3]王永明，王贵锦，“图像局部不变性特征与描述”，国防工业出版社
[4]C. Harris and M. Stephens (1988). "A combined corner and edge detector". Proceedings of the 4th Alvey Vision Conference. pp. 147–151.
本文出自 “UnderTheHood” 博客，请务必保留此出处http://underthehood.blog.51cto.com/2531780/658350