SIFT特征提取分析
读‘D. G. Lowe. Distinctive Image Features from Scale-Invariant Keypoints[J],IJCV,2004’ 笔记
关键点是指图像中或者视觉领域中明显区别于其周围区域的地方，这些关键点对于光照，视角相对鲁棒，所以对图像关键点提取特征的好坏直接影响后续分类、识别的精度。
特征描述子就是对关键点提取特征的过程，应该具备可重复性、可区分性、准确性、有效性和鲁棒性。
SIFT(Scale-Invariant Feature Transform)是一种特征描述子。该描述子具有尺度不变性和光照不变性。
检测尺度空间极值
构建尺度空间
这里的尺度可以理解为图像的模糊程度，就是眼睛近视的度数。尺度越大细节越少，SIFT特征希望提取所有尺度上的信息，所以对图像构建尺度空间, 也就是实用不同的平滑核对图像进行平滑。
It has been shown by Koenderink (1984) and Lindeberg (1994) that under a variety of reasonable assumptions the only possible scale-space kernel is the Gaussian function.
所以，这里的平滑核选用高斯核，空间尺度有高斯核尺度决定
其中
是原图像，*是卷积符号，
对应尺度下的尺度图像。
是高斯核。
其中左上角为原图，右上角为
尺度下的图像，下行左右分别是
尺度的图像。
尺度空间描述了不同模糊程度上的特征，而并没有刻画图像尺寸大小的情形，所以文中将尺度空间和图像金字塔结合起来。图像的金字塔是通过对图像逐步平滑下采样得到的。
令每一层金字塔(octave)中尺度变化范围为
,每层金字塔测量s个尺度，那么第t层金字塔的尺度范围就是
,第
层金字塔的第一幅图像由上层塔中
尺度下的下采样的到，下采样比例为0.5.
最终的尺度金字塔图像如下:
每一层金字塔中多层尺度变化因子是
,第n层尺度是
这样能够保证从金字塔底部到顶部的尺度是连续的
,
DoG极值
现在已经将图像变到尺度空间了，希望在这个空间找到比较显著的点，即变化明显的点，也就是说在图像逐渐模糊的过程中引起你注意的点，这个点突不突出可以通过尺度轴上变化趋势来体现，即尺度轴上的梯度极值。
在尺度轴上计算梯度就是图像沿着尺度轴差分呗
其示意图如下:
这里DoG近似等于尺度归一化的高斯拉普拉斯算子
,而尺度归一化的高斯拉普拉斯算子相较于其他角点检测算子，如梯度，Hessian或Harris焦点特征能够更稳定的图像特征[1],这篇参考文献我没看。
现在就要找极值点了，在DoG空间中寻找极值，也就是说比周围点都大或者都小的点认为是关键点。既然DoG是差分图像叠在一起的，这里的周围当然是3D的周围，共包括相邻的26个邻点。如下图
标注的点要和周围26个绿色的点比较。
那现在就是沿着尺度轴，逐步查找对应尺度上的极值。我们现在在图像的一个像素点沿着尺度轴看过去，就成了一维图像，让我们看看这个问题。
每层金字塔中有s层尺度图像，就是说有s个尺度，我们希望在s个尺度上都寻找极值点，那么应该有s+2层DoG图像才可以，因为这个过程是三层数据比较嘛。 而DoG图像是两层图像差分得来的，所以s+2层DoG图像应该需要s+3层尺度图像。
所以在每一层金字塔上都应该计算s+3层尺度图像。
注：这里是沿着尺度轴类比到一维上观测的示意图。图中单位尺度使用a表示，每一层金字塔包含3个尺度，即
。发现第一层金字塔最后寻找尺度为到
,而第二层金字塔中最先寻找到的尺度为
,所以尺度是连续的。
这样我们就从DoG空间中粗略的找到了关键点。但为了进一步精化结果，需要对检
的结果进行筛选，去掉噪声点。
去除干扰
无论是LoG还是DoG在做关键点检测时都会受到边缘的影响，如下图
图中最高点实际上是在一个一条脊线上相对较高的地方，那么认为这个点并不是blob，而是边缘的一部分，但是这种点再DoG空间中是会被找极值找到的。
对于这种边缘，虽然是一个极值点，但是会发现这个极值垂直于边缘方向变换剧烈而平行于边缘变化相对平缓，那么就可以通过这两个方向变化幅度的比值判定找到的最大值点是山脊上的一处制高点还是一处独立的山峰即blob的中心点。下图是山脊线在不同方向的剖面，其中上一张子图是平行于脊线，下图垂直于脊线。
当然了也会存在平面上的一个很小的土丘被检测出来，这被认为是噪声点，这种噪声在DoG空间的幅值都比较小，可以通过阈值化处理。
要去除干扰就要分别处理这两种情况。
去除干扰1: 较小的极值
当然可以直接对找到的极值点的值进行阈值化过滤，但是为了更加精确的定位关键点的位置，对关键点进行3D二次函数的拟合(可以类比最优化里牛顿法的做法)，然后对函数求极值。
将每个关键点的一段函数
进行泰勒二次展开:
其中x是相对于关键点的距离。
和
分别表示三元函数`!$D(x,y,\sigma)'在关键点的取值和一阶导数，二阶导数。这个展开就是在关键点周围拟合出的3D二次曲线，找极值就令其导数为零。于是：
此时的极值为
文中剔除
的点。
去除干扰2:边缘噪声
一维情形下曲线的变化趋势可以通过曲率反应。
那么在检测最大值点是不是边缘噪声时可以选择曲面变化最大和最小的两个剖面，判断其剖面曲线在该点的曲率，曲率越大表示越陡峭，曲率越小表示越平缓，而对于blob区域，即山峰而言，理想状态是任意方向的曲率都相同且比较大。而山脊线即边缘则会在沿线方向的变化较缓慢，即曲率较小，在垂直脊线的方向变化较剧烈，曲率较大。可以通过两个曲率的比值判定是山峰还是山脊线。
对于二维曲面，曲面点处的二阶矩矩阵，即Hessian矩阵刻画了该点周围的变化趋势。说白了二阶矩矩阵就是二元函数在该点的协方差矩阵，协方差矩阵的特征值对应着在特征向量方向上的投影，该值越大就反应函数在该方向变化趋势越缓慢，也就是曲率越大。所以说Hessian矩阵的特征值和特征值所在特征向量方向上的曲率成正比。我们希望使用曲率的比值剔除边缘点，就可以使用Hessian的特征值的比值替代曲率的比值。
Hessian矩阵定义如下：
可以由二阶差分计算得到，那么接下来就可以计算特征值的比值了。
先等一等，假设两个特征值分别为
,看下面两个式子：
其中
分别是矩阵的迹和矩阵的行列式，这样发现特征值的比值可以通过这两个两计算得到，而不需要特征值分解，这样简单多了。
假设
为较大特征值，且
,则
当
是
最小，所以当
越大时，对应的
越大。所以我们要将
的点剔除就相当于将
的点剔除。论文中
。
关键点的表示
ok,在DoG空间中已经找到并细化出若干关键点，那么这些关键点该怎么表示呢？
可以使用三元组
表示关键点，其中三个值分别表示关键点的位置、尺度和方向。为什么要计算特征点的方向呢？这是为了使特征描述子具有旋转不变特性(下一节有讲到，为什么就旋转不变了)。
样本点的方向定义(其实就是对应尺度图像上对应位置的梯度)：
是对应尺度的尺度图像。
而关键点的方向是由他周围区域内样本点决定的，比如
区域内的所有样本点计算梯度值和方向，将方向分为若干个bins，然后使用高斯函数加权将周围样本点统计出方向直方图，比如论文中将一周360°分成36个bins，那么我们统计方向落在每个bins中样本，将该样本梯度值乘以高斯权重后加入这个bins中，就得到了长度为36的方向直方图。在这个方向直方图中，最大峰值对应的bins就是关键点的方向。如果存在多个峰值或者存在大于0.8倍最大峰值的bins，那么就在该点创建多个关键点，这些关键点的位置、尺度一致，但是方向不同(参考[2]中此处也存在错误)。
论文中采用的高斯权重函数的尺度
,其中
是该关键点所在的图像在尺度空间中的尺度。而由高斯函数
原则，关键点的邻域半径应该设置为
.
局部图象描述子
找到了图像在不同尺度的关键点，我们还希望刻画关键点周围的特征，以利于后续的分类或者匹配等操作。
局部特征肯定离不开关键点周围的区域了，将关键点附近半径为
的邻域(论文中选取边长为16的正方形邻域)划分为
的子区域(论文中推荐
)，在每个子区域统计长度为
的方向直方图(论文中直方图长度为8)，每个直方图被称为一个种子点，这样一个关键点的特征描述子就是一个长度为
的向量(论文中是
维的特征向量)。
如下图所示
其中中间的点是检测到的关键点，蓝色的点表示该尺度图像中的像素点，红色方格表示划分出的子区域
，在每个子区域统计方向直方图得到种子点。这里我们只是示意，具体的区域大小接下来详细讨论。我们还注意到关键点发出的一条橙黄色箭头，这表示关键点的方向。
那么该描述子怎么体现旋转不变性呢？让我们来看下面这个过程
这两个区域应该是一个图像，只不过旋转了90°，而如果按照‘从最上边的点顺时针统计E点的邻域’，分别是
计算出来的差别就会好大，但是如果我们将E点的方向固定，都指向右侧，那么它周围的邻域的相对位置也就固定了，这时候再按照‘从最上边的点顺时针统计E点的邻域’那么结果就是是相同的了，我们可以认为图中左右两个区域是同一区域的不同视角。
现在再来看图像，图像是规则的矩形区域，我们在统计方向直方图是也是有固定的规则的，矩形区域的选择一般都是平行于图像边缘的，所以为了去除旋转的影响，可以，将每个关键的方向固定到相同方向，那么相同区域经过固定规则检测出来的结果就非常相近，这就避免了旋转角度的影响。一般可以将图像旋转，使关键点的方向统一对齐到图像x轴的方向[3]，然后再对旋转后的图像划分子区域统计方向直方图。有些资料说将坐标轴方向旋转到关键点方向上，可能意思和我一样，但我觉得没有“旋转图像使将关键点方向和x轴方向重合”这种表述清晰。我刚开始看到那种表述一直搞不明白在干嘛。。。
如下图，我们将图像旋转使关键点的方向对齐到x轴方向，然后在使用红色的网格对旋转后的图像统计方向直方图。
坐标旋转后的值
这里的
是关键点方向和x轴方向的夹角，顺时针旋转该角度为负值，逆时针为正值。
ok，现在保证描述子旋转不变了，让我们按照红色网格区域划分子区域来统计方向直方图吧。可是这时候又发现了个问题，这怎么统计啊？对于每旋转的图像正好像素点落在红色网格内，可现在旋转后落得乱七八糟，怎么办呢？这时候就需要采用双线性插值的方法插出网格点的值了~
我们从上图中取出一小块来分析：
其中蓝色的点就是旋转后的点，而绿色的点就是网格的中心，我们现在要插值插出绿色点的值，就要看蓝色点在这些点的影响了。图中
分别表示蓝色的点对左上角绿色点的在两个方向的影响率，最终的蓝色点分布到左上角点的影响是
,这里
是绿色点的值。关于双线性插值见[5][4].
假设绿色的点的梯度值为
,方向为60°，那么其在左上角的影响为
,方向为60°，
为高斯权重。那么对于每一个网格点都要计算他周围的点对他的影响，然后累加起来，这时候会产生不同的点在同一个网格点产生影响的方向不同，此时应该将方向插分到与之最近的统计方向上，使得方向直方图更加准确。比如直方图长度为8，那么每45°一个方向，所以60°的方向影响应该插分到45°和90°方向上，所以影响力有
落在了45°方向上，
落在90°方向上。
综上，方向直方图统计网格点
处的值时，在方向o上增量为
这里，
是p点周围的旋转后的样本点，因为是逐像素的所以限制其距离在单位1内。
是
的坐标，
是高斯权重，
分别是
对该网格点在两个方向的影响率和对所求方向上的影响率。
现在旋转不变性还剩下一个一直没有说的问题就是区域大小的选择和高斯权重尺度的选择。下面就和说一说~
在统计方向直方图串联成局部特征描述子时，每个子区域的选择应该和计算关键点方向时的区域大小一致，即
,其中
是图像在尺度空间的尺度。考虑到旋转的问题，我们不希望旋转后半径为
的区域内有部分是空的，所以我们选的区域要能够在旋转后仍然后
的区域可以检测到。
另外由于需要采用双线性插值的方法获得网格点的值，所以每一个子区域的半径为
总体区域的半径应该为
去除光照影响
为了去除光照影响，一般将关键点生成的特征向量
归一化。
是方差。
描述子向量门限。 非线性光照，相机饱和度变化会造成某些方向的梯度值过大，而对方向的鉴别性微弱。因此一般归一化后，取0.2截断较大的梯度值，然后再次进行归一化，提高特征的鉴别性。[3]
总结
将所有尺度上关键点提取的局部描述子取出来，就完成了图像sift特征的提取。在图像匹配时，对两幅图像分别提取sift特征，然后在每个尺度上进行描述子的匹配。在匹配中如果一个关键点的最小匹配距离和第二小匹配距离小于某个阈值，则认为是正确匹配。这里参见参考文献3.
SIFT特征特性：
尺度不变性
旋转不变性
光照不变性
鉴别性强，信息量丰富
SIFT特征的缺点
实时性不高，因为不断的下采样和插值等操作
对于边缘光滑的目标无法准确提取特征点
SIFT特征的用途
主要解决图像配准和目标识别跟踪中下述问题
目标的旋转、缩放、平移
图像的仿射/射影变换
部分减轻光照影响
目标的部分遮挡
杂物场景(crowd scene)
噪声
后序
关于LoG和DoG的关系以及什么是尺度归一化的问题，还有我想了好久的问题：LoG和DoG都是通过zero-acrossing检测边缘的，这里为什么使用极值来检测blob？这些将在下一篇文章斑点检测分析。
由于最近在挖目标跟踪的东西，就挖到了sift特征，所以花了近一周时间看了Lowe的论文，加上网上的资料，总的来说应该还算清楚的捋顺了SIFT特征的原理。
关于具体的实验及其结果，我暂时没有实现，因为我了解SIFT只是为了更好的理解现在看的东西，所以暂时就不实现了，等忙完了这一段，抽空把代码补上。
Mikolajczyk,K. 2002. Detection of local features invariant to affine transformation, Ph.D. thesis, Institut National Polytechnique de Grenoble, France. ↩