你对SIFT算法了解多少—原理详解与演示

1. 简介

  SIFT算法，全称为尺度不变性特征匹配算法（Scale Invariant Feture Transform），它是计算机视觉中最常应用的特征匹配算法之一，用来检测与描述图像中的局部性特征，这种描述具有尺度不变性，对旋转、缩放、亮度变化等保持不变性。该算法由David Lowe1999年发表，在2004年进行完善，后续在其基础上有很多改进的算法，目前用到的特征匹配算法可能不是SIFT算法本身而是改进后的算法，但是无论用什么改进的算法，其基本思想是不变的，只不过是在计算性能和效率上有所提升。

2. SIFT的特点及性能

  SIFT局部特征描述有助于辨识物体，这些特征描述子基于物体上的局部特征点生成而与物体的大小及旋转无关。对噪声、旋转、亮度变化的鲁棒性较好，而且SIFT特征描述在一定程度上也能对目标物体部分遮蔽进行检测，在少数的目标物体上也能够产生很多的SIFT特征向量，特征信息量大，经过筛选优化后可进行快速匹配，SIFT特点总结如下：

• 图像的局部特征描述子具有较好的稳定性，能够适应旋转、尺度缩放、亮度的变化，能在一定程度上不受视角变化、仿射变换、噪声的干扰。
• 多量性，即使少数的几个物体也可以产生大量的SIFT特征向量。
• 独特性好，特征信息量丰富，适用于海量特征库进行快速、准确的匹配。
• 高速性，经优化的SIFT匹配算法能快速进行匹配。
• 可扩展性，能够与其他形式的特征向量进行联合。

  目标物体的自身状态、所处环境和成像设备等因素影响图像配准、目标识别跟踪的性能。而SIFT算法在一定程度上可解决：

• 目标的旋转、缩放、平移（RST）
• 图像仿射/投影变换（视点viewpoint）
• 光照影响（illumination）
• 目标遮挡（occlusion）
• 复杂场景（clutter）
• 噪声干扰

3. SIFT特征检测步骤

• 尺度空间的极值检测： 搜索所有尺度空间上的图像，通过高斯微分函数来识别潜在的对尺度和旋转不变性的兴趣点。
• 关键点定位：在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。关键点的选择依据于它们的稳定程度。
• 特征方向确定：基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向。所有后面的对图像数据的操作都相对于关键点的方向、尺度和位置进行变换，从而提供对于这些变换的不变性。
• 关键点描述：在每个特征点周围的邻域内，在选定的尺度上测量图像的局部梯度，这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变换。

  下面根据SIFT特征检测步骤对该算法进行通俗地解读，SIFT算法涉及到的数学概念及实现流程较多，下面分模块逐步进行解读。

4. 尺度空间的极值检测

4.1 图像尺度空间的概念

  在一定的范围内，无论是大目标物体还是小目标物体，人眼都可以分辨出来，然而计算机要实现不仅在清晰的层面上辨别出物体，还要在模糊的层面上辨别出物体有什么特点却有一定的难度，竟然有难度该如何解决呢？

  解决办法并不复杂，人类在什么层面上观察到的目标，让计算机也在相同条件下辨别目标，所以要让机器能够对物体在不同尺度上有一个统一的认知，就需要考虑图像在不同尺度下都存在的特点。Lindeberg等人已证明高斯卷积核是实现尺度变换的唯一变换核，并且是唯一的线性核。SIFT算法在构建尺度空间时候也是采取高斯核函数进行滤波，使原始图像保存最多的细节特征，经过高斯滤波后细节特征逐渐减少来模拟大尺度情况下的特征表示。如下图1所示，有6副Lena图像，第1幅图像为原图，对图像做一些改变后逐渐变得模糊看不清图像中的具体内容。

图1 高斯模糊图像

4.2 尺度空间的表示

  早期图像的多尺度通常使用图像金字塔来表示形式，图像金字塔是同一图像在不同的分辨率下得到的多组结果，其生成过程首先是对原始图像进行平滑，然后对处理后的图像进行降采样，降采样后得到一系列尺寸不断缩小的图像。

  一个图像的尺度空间$L(x,y,\sigma )$定义为一个变化尺度的高斯函数$G(x,y,\sigma )$与原图像$I(x,y)$的卷积。
$$L(x,y,\sigma )=G(x,y,\sigma )\ast I(x,y)\text{\hspace{1em}(1)}$$
其中，$\ast$表示卷积，高斯核函数$G(x,y,\sigma )=\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{(x-m/2{)}^2+(y-n/2{)}^2}{2{\sigma }^2}}$，$m,n$为高斯模板的维度，$(x,y)$表示图像像素点的位置，$\sigma$表示尺度空间因子，对应高斯正态分布中的标准差，该值越大图像越模糊，对应的尺度也就越大，反之，值越小表示图像被平滑的越少，相应的尺度也就越小。大尺度对应于图像的概貌特征，小尺度对应于图像的细节特征。

4.3 高斯金字塔的构建

  高斯金字塔的构建包含两部分：

• 对图像做不同尺度的高斯模糊。
• 对图像做降采样，通常是采样前图像的水平、垂直方向的1/2。

尺度空间使用高斯金字塔实现如下图2所示。

图2 尺度空间的高斯金字塔实现

图像的金字塔模型将原始图像不断降阶采样得到一系列大小不一的图像，从下到上由大变小构成金字塔模型。原图像为金字塔的第一层，每次降采样得到的新图像为金字塔的某一层（每层一张图像）。

  为体现尺度空间的连续性（下文有解释），高斯金字塔在简单降采样的同时还进行了高斯模糊。如上图所示，对金字塔的每层图像中的每张图像使用不同参数做高斯模糊，使得金字塔的每层含有多张高斯模糊后的图像，在这里将金字塔每层多张图像合称为一组（Octave），金字塔的每层含有一组图像，不难发现组数与金字塔的层数相等，每组中的每张图像称为该组中的层（Interval）。

说明： 由于每个Octave中的多张图像是按层次叠放的，为避免与高斯图像金字塔中的层数混淆，本文对“层”未特别说明的地方一般都指每个Octave中的各层图像。

  不难发现，高斯金字塔有多组，每组又有多层。一组中的多个层之间的尺度是不一样的（也就是使用的高斯参数$\sigma$是不同的），相邻两层之间的尺度相差一个比例因子$k$。如果每组有$S$层，则$k=2^{\frac{1}{S}}$。上一组图像的最底层图像是由前一组中尺度为$2\sigma$的图像进行因子为2的降采样得到的（高斯金字塔先从底层建立）。

  假设金字塔的组数共有$o$组，其组数根据图像的大小和塔顶图像的大小共同决定，其计算公式如下：
$$o=lo{g}_{2}min(M,N)-a\text{\hspace{1em}(2)}$$
其中，$M,N$分别是原图像的行和列，$a$可以在$0-lo{g}_{2}min(M,N)$之间任意取值（即塔顶图像的最小维数的对数值）。

  以一个$512\times 512$的图像$I$为例，构建高斯金字塔步骤：(从0开始计数，倒立的金字塔）

1. 金字塔的组数，$lo{g}_{2}512=9$，减去因子3，构建的金字塔的组数为6。取每组的层数为3（一般为3~5层）。
2. 构建第0组，将图像的宽和高都增加一倍，变成$1024\times 1024（I_0）$。第0层$I_0\ast G(x,y,{\sigma }_0)$，第1层$I_0\ast G(x,y,k{\sigma }_0)$，第2层$I_0\ast G(x,y,k^2{\sigma }_0)$。
3. 构建第1组，对$I_0$降采样变成$512\times 512（I_1）$。第0层$I_1\ast G(x,y,2{\sigma }_0)$，第1层$I_1\ast G(x,y,2k{\sigma }_0)，第2层I_1\ast G(x,y,2k^2{\sigma }_0)$。
4. …
5. 构建第$o$组，第$s$层$I_o\ast G(x,y,2^o{k}^s{\sigma }_0)$。

  根据上述分析，高斯模糊参数$\sigma$（尺度空间），可由下面关系式得到
$$\sigma (o,s)={\sigma }_0\cdot 2^{o+\frac{s}{S}}\text{\hspace{1em}(3)}$$
其中$o$为所在的组，$s$为所在的层，${\sigma }_0$为初始的尺度，$S$为每组的层数。

  在Lowe论文的算法实现中${\sigma }_0=1.6$，$o_{min}=-1$，$S=3$，$o_{min}=-1$就是首先将原图像的长和宽各扩展一倍。

  从上面分析可以得知相邻组之间的尺度关系
$${\sigma }_{o+1}=2{\sigma }_o\text{\hspace{1em}(4)}$$
  同一组内相邻层的图像尺度关系
$${\sigma }_{s+1}=k\cdot {\sigma }_s=2^{\frac{1}{S}}\cdot {\sigma }_s\text{\hspace{1em}(5)}$$

  为了有效的在尺度空间检测到稳定的关键点，Lowe提出了高斯差分尺度空间（DoG Scale-Space），将每一组相邻的两层相减就可以得到DoG金字塔，利用不同尺度的高斯差分核与图像卷积生成。
$$D(x,y,\sigma )=(G(x,y,k\sigma )-G(x,y,\sigma ))\ast I(x,y)=L(x,y,k\sigma )-L(x,y,\sigma )\text{\hspace{1em}(6)}$$

注： 虽然使用Laplacian of Gaussian（LoG）能够很好地找到找到图像中的兴趣点，但是需要大量的计算量，而Difference of Gaussian（DoG）算子计算简单，是尺度归一化的LoG算子的近似，所以使用DoG图像的极大极小值近似寻找特征点。有关DoG寻找特征点的介绍及方法详见DoG (Difference of Gaussian)角点检测，极值点检测用的非极大值抑制。

  在实际计算时，使用高斯金字塔每组中相邻上下两层图像相减，得到高斯差分图像，如下图3所示，进行极值检测。

图3 高斯差分过程

4.4 检测DoG尺度空间极值点

  图像中的特征点是由DoG尺度空间的局部极值点组成，为寻找DoG空间中的极值点，每一个像素点要和它所有的相邻点进行比较，看其是否比它的图像域和尺度域的相邻点大或者小。如下图4所示，中间的检测点和它同尺度的8个相邻点及上下相邻尺度对应的9×2=18个点共26个点比较，以确保在二维图像空间和尺度空间都检测到极值点。 一个点如果在DoG尺度空间在本层以及上下两层的26个点中是最大或最小值时，就认为该点是图像在该尺度下的一个特征点。

图4 极值点检测

  在极值比较的过程中，每一组图像的首末两层是无法进行极值比较的，为了满足尺度变化的连续性，我们在每一组图像的顶层继续用高斯模糊生成了3幅图像，高斯金字塔有每组S+3层图像，DoG金字塔每组有S+2层图像。解释如下：

图5 高斯差分金字塔图像

  同样假设S=3，也就是每组有3层，则$k=2^{\frac{1}{s}}=2^{\frac{1}{3}}$，按照上图可得Gauss Space和DoG space分别有3个（S个）和2个（S-1个）分量，在DoG space中，第一层两项分别是$\sigma ,k\sigma$; 第二层两项分别是$2\sigma ,2k\sigma$，由于无法比较极值，我们必须在高斯空间继续添加高斯模糊项，使得形成$\sigma ,k\sigma ,k^2\sigma ,k^3\sigma ,k^4\sigma$，这样就可以选择DoG space中的中间三项$k\sigma ,k^2\sigma ,k^3\sigma$（只有左右都有才能有极值），那么下一组中（由上一层降采样获得）所得三项即为$2k\sigma ,2k^2\sigma ,2k^3\sigma$，其首项$2k\sigma =2\cdot 2^{\frac{1}{3}}\sigma =2^{\frac{4}{3}}\sigma$。刚好与上一组末项$k^3\sigma =2\cdot (2^{\frac{1}{3}}{)}^3\sigma =2\cdot 2^{\frac{3}{3}}\sigma$尺度变化连续起来，所以每次要在Gaussian space添加3项，每组共S+3层图像，相应的DoG金字塔有S+2层图像。

5. 关键点定位

  上述检测到的极值点是离散空间的极值点，以下通过拟合三维二次函数来精确确定关键点的位置和尺度，同时去除低对比度的关键点和不稳定的边缘响应点(因为DoG算子会产生较强的边缘响应)，以增强匹配稳定性、提高抗噪声能力。

5.1 关键点的精确定位

  离散空间的极值点并不是真正的极值点，下图6显示了二维函数离散空间得到的极值点与连续空间极值点的差别。利用已知的离散空间点插值得到的连续空间极值点的方法叫做子像素插值。

图6 离散空间点插值

  为了提高关键点的稳定性，需要对尺度空间DoG函数进行曲线插值。利用DoG函数在尺度空间的Taylor展开式(插值函数)，任意一极值点在$(x_0,y_0,{\sigma }_0)$处的泰勒展开如下（去掉2阶后面的高次项）：

$$f\left(\left[\begin{array}{c}x\\ y\\ \sigma \end{array}\right]\right)\approx f\left(\left[\begin{array}{c}{x}_0\\ y_0\\ {\sigma }_0\end{array}\right]\right)+\left[\frac{\partial f}{\partial x},\frac{\partial f}{\partial y},\frac{\partial f}{\partial \sigma }\right]\left(\left[\begin{array}{c}x\\ y\\ \sigma \end{array}\right]-\left[\begin{array}{c}{x}_0\\ y_0\\ \sigma \end{array}\right]\right)\text{\hspace{1em}(7)}$$

上式的矩阵表示如下：
$$D(X)=D+\frac{\partial D^T}{\partial X}X+\frac{1}{2}{X}^T\frac{{\partial }^{2}D}{\partial X^2}X\text{\hspace{1em}(8)}$$

其中，对$X$求导并令方程等于0，可以得到极值点的偏移量为
$$\hat{X}=-\frac{{\partial }^2{D}^{-1}}{\partial X^2}\frac{\partial D}{\partial X}\text{\hspace{1em}(9)}$$

代入上述矩阵方程，对应极值点的值为
$$D(\hat{X})=D+\frac{1}{2}\frac{\partial D^T}{\partial X}\hat{X}\text{\hspace{1em}(10)}$$

其中, $\hat{X}$代表相对插值中心的偏移量，当它在任一维度上的偏移量大于0.5时（即$x$或$y$或 $\sigma$），意味着插值中心已经偏移到它的邻近点上，所以必须改变当前关键点的位置。同时在新的位置上反复插值直到收敛；也有可能超出所设定的迭代次数或者超出图像边界的范围，此时这样的点应该删除，在Lowe中进行了5次迭代。另外，过小的点易受噪声的干扰而变得不稳定，所以将小于某个经验值(Lowe论文中使用0.03，Rob Hess等人实现时使用0.04/S)的极值点删除。同时，在此过程中获取特征点的精确位置(原位置加上拟合的偏移量)以及尺度($\sigma$)。

5.2 消除边缘响应

  一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的主曲率，而在垂直边缘的方向有较小的主曲率。DoG算子会产生较强的边缘响应，需要剔除不稳定的边缘响应点。获取特征点处的Hessian矩阵，主曲率通过一个2x2的Hessian矩阵$H$求出（$D$的主曲率和$H$的特征值成正比）：
$$H=\left[\begin{array}{cc}{D}_{xx}& D_{xy}\\ D_{xy}& D_{yy}\end{array}\right]\text{\hspace{1em}(11)}$$

假设$H$的特征值为$\alpha$和$\beta$（$\alpha$和$\beta$代表$x$和$y$方向的梯度）且$\alpha >\beta$，则

$$T_r(H)=D_{xx}+D_{yy}=\alpha +\beta \text{\hspace{1em}(12)}$$
$$Det(H)=D_{xx}{D}_{yy}-(D_{xy}{)}^2=\alpha \beta \text{\hspace{1em}(13)}$$
令$\alpha =\gamma \beta$，则
$$\frac{Tr(H{)}^2}{Det(H)}=\frac{(\alpha +\beta {)}^2}{\alpha \beta }=\frac{(\gamma \beta +\beta {)}^2}{\gamma {\beta }^2}=\frac{(\gamma +1{)}^2}{\gamma }\text{\hspace{1em}(14)}$$

$(\gamma +1{)}^2/\gamma$的值在两个特征值相等时最小，随着$\gamma$的增大而增大。值越大，说明两个特征值的比值越大，即在某一个方向的梯度值越大，而在另一个方向的梯度值越小，而边缘恰恰就是这种情况。所以为了剔除边缘响应点，需要让该比值小于一定的阈值，因此，为了检测主曲率是否在某域值$\gamma$下，只需检测:
$$\frac{Tr(H{)}^2}{Det(H)}<\frac{(\gamma +1{)}^2}{\gamma }\text{\hspace{1em}(15)}$$

如果$(\alpha +\beta {)}^2/\alpha \beta >(\gamma +1{)}^2/\gamma$则舍弃，在Lowe文章中，取$\gamma =10$。

6. 特征方向确定

  上一步中确定了每幅图中的特征点，为每个特征点计算一个方向，依照这个方向做进一步的计算， 利用关键点邻域像素的梯度方向分布特性为每个关键点指定方向参数，使算子具备旋转不变性。

6.1 梯度计算

  对于在DoG金字塔中检测出的关键点，采集其所在高斯金字塔图像$3\sigma$邻域窗口内像素的梯度和方向分布特征。梯度的模值和方向如下：
$$m(x,y)=\sqrt{(L(x+1,y)-L(x-1,y){)}^2+(L(x,y+1)-L(x,y-1){)}^2}\text{\hspace{1em}(16)}$$
$$\theta (x,y)=ta{n}^{-1}((L(x,y+1)-L(x,y-1))/L(x+1,y)-L(x-1,y))\text{\hspace{1em}(17)}$$

$L$为关键点所在的尺度空间值，按Lowe的建议，梯度的模值$m(x,y)$按 $\sigma =1.5{\sigma }_{oct}$的高斯分布加成，按尺度采样的$3\sigma$原则，领域窗口半径为$3\times 1.5{\sigma }_{oct}$。

6.2 梯度直方图

  梯度直方图的范围是0~360°，其中每45°一个$bin$，总共8个$bins$，或每10°为一个$bin$，总共有36个$bins$。随着距中心点越远的邻域对梯度直方图的贡献越小。在实际计算时，在以关键点为中心的邻域窗口内采样，并用直方图统计邻域像素的梯度方向。Lowe论文中还提到要使用高斯函数对直方图进行平滑，减少突变的影响。直方图的峰值则代表了该关键点处邻域梯度的主方向，即作为该关键点的方向。如下图7所示，为简化，图中只画了八个方向的直方图。

图7 关键点方向直方图

方向直方图的峰值则代表了该特征点处邻域梯度的方向，以直方图中最大值作为该关键点的主方向。为了增强匹配的鲁棒性，只保留峰值大于主方向峰值80％的方向作为该关键点的辅方向。因此，对于同一梯度值的多个峰值的关键点位置，在相同位置和尺度将会有多个关键点被创建但方向不同。仅有15％的关键点被赋予多个方向，但可以明显的提高关键点匹配的稳定性。实际编程实现中，就是把该关键点复制成多份关键点，并将方向值分别赋给这些复制后的关键点，并且，离散的梯度方向直方图要进行插值拟合处理，来求得更精确的方向角度值。

图8 梯度直方图插值处理

如果在第$i$个$bin$寻找一个确定的方向，由以上分析做出以下假设：
假设插值抛物线方程为$H(u)=a{u}^2+bu+c$，其中$a,b,c$为抛物线的系数，$u\in [-1,1]$，对此抛物线求导并令其等于0。

由$H^{\prime }(u)=0$得$u_{max}=-\frac{b}{2a}$，即抛物线的顶点位置横坐标。现分别将$H(-1),H(0),H(1)$代入抛物线插值方程

$$\{\begin{array}{l}H(-1)=a-b+c\\ H(0)=c\\ H(1)=a+b+c\end{array}\text{\hspace{1em}(18)}$$

解得

$$\{\begin{array}{l}a=[H(1)+H(-1)]/2-H(0)\\ b=[H(1)-H(-1)]/2\\ c=H(0)\end{array}\text{\hspace{1em}(19)}$$

进一步可得

$u_{max}=-\frac{b}{2a}=\frac{H(-1)-H(1)}{2[H(-1)+H(1)-2H(0)]}$，即局部坐标系中的值。$bin$在直方图中的索引号为$i^{\prime }=i+u_{max}=i+\frac{H(-1)-H(1)}{2[H(-1)+H(1)-2H(0)]}$。

至此，图像的关键点已检测完毕，每个关键点有三个信息：位置、所处尺度、方向。由此可以确定一个SIFT特征区域。

7. 关键点描述

  通过以上步骤，对于每一个关键点，拥有三个信息：位置、尺度以及方向。接下来就是为每个关键点建立一个描述符，使其不随各种变化而改变，比如光照变化、视角变化等等。并且描述符应该有较高的独特性，以便于提高特征点正确匹配的概率。

  SIFT描述子是关键点邻域高斯图像梯度统计结果的一种表示。通过对关键点周围图像区域分块，计算块内梯度直方图，生成具有独特性的向量，这个向量是该区域图像信息的一种抽象，具有唯一性。Lowe建议描述子使用在关键点尺度空间内4×4的窗口中计算的8个方向的梯度信息，共4*4*8=128维向量表征。

7.1 确定描述子所需的区域

  特征描述子与特征点所在的尺度有关，所以对梯度的求取应在特征点对应的高斯图像上进行。将关键点附近的邻域划分为$d\times d$(Lowe建议$d=4$)个子区域，每个子区域做为一个种子点，每个种子点有8个方向。每个子区域的大小与关键点方向分配时相同，即每个区域有$3{\sigma }_{oct}$个子像素，为每个子区域分配边长为$3{\sigma }_{oct}$的矩形区域进行采样($3{\sigma }_{oct}$个子像素实际用边长为$\sqrt{3{\sigma }_{oct}}$的矩形区域即可包含，为了简化计算取其边长为$3{\sigma }_{oct}$，并且采样点宜多不宜少)。考虑到实际计算时，需要采用三线性插值，所需图像窗口边长为$3{\sigma }_{oct}\times (d+1)$。在考虑到旋转因素(方便下一步将坐标轴旋转到关键点的方向)，如下图9所示，实际计算所需的图像区域半径为：

图9 邻域半径变化

7.2 坐标轴旋转至主方向及分配采样点

  将坐标轴旋转为关键点的方向，以确保旋转不变性。

图10 坐标轴旋转

旋转后邻域内采样点的新坐标为
$$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)=\left(\begin{array}{cc}cos\theta & -sin\theta \\ sin\theta & cos\theta \end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)\text{\hspace{1em}(20)}$$
其中$(x,y)\in [-radius,radius]$

  将邻域内的采样点分配到对应的子区域内，将子区域内的梯度值分配到8个方向上，计算其权值。旋转后的采样点坐标在半径为$radius$的圆内被分配到的子区域，计算影响子区域的采样点的梯度和方向，分配到8个方向上，旋转后的采样点落在子区域的坐标为

$$\left(\begin{array}{c}{x}^{\prime \prime }\\ y^{\prime \prime }\end{array}\right)=\frac{1}{3{\sigma }_{oct}}\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)+\frac{d}{2}\text{\hspace{1em}(21)}$$

该坐标为归一化后的坐标，在特征点局部坐标系中把$3\sigma$长度作为1个单元长度，$d/2$是把坐标系由特征点处平移至左上角的边界点上。Lowe建议子区域的像素的梯度大小按$\sigma =0.5d$的高斯加权计算，即

$$w=m(a+x,b+y)\ast e^{-\frac{(x^{\prime }{)}^2+(y^{\prime }{)}^2}{2\times (0.5d{)}^2}}\text{\hspace{1em}(22)}$$

其中$a，b$为关键点在高斯金字塔图像中的位置坐标。

7.3 插值计算种子点八个方向的梯度

  采样点在子区域中的坐标$(x^{\prime \prime },y^{\prime \prime })$(图中蓝色窗口内红色点)线性插值，计算其对每个种子点的贡献。如图11中的红色点，落在第0行和第1行之间，对这两行都有贡献。对第0行第3列种子点的贡献因子为$dr$，对第1行第3列的贡献因子为$1-dr$，同理，对邻近两列的贡献因子为$dc$和$1-dc$，对邻近两个方向的贡献因子为$do$和$1-do$。则最终累加在每个方向上的梯度大小为：

$$weight=w\ast d{r}^k\ast (1-dr{)}^{1-k}\ast d{c}^m\ast (1-dc{)}^{1-m}\ast d{o}^n\ast (1-do{)}^{1-n}\text{\hspace{1em}(23)}$$

其中$k，m，n$为0（像素点超出了对要插值区间的四个邻近子区间所在范围）或为1（像素点处在对要插值区间的四个邻近子区间之一所在范围）。

图11 每个种子点八个方向的梯度描述

7.4 特征描述子及门限化

  如上统计的4*4*8=128个梯度信息即为该关键点的特征向量。特征向量形成后，为了去除光照变化的影响，需要对它们进行归一化处理，对于图像灰度值整体漂移，图像各点的梯度是邻域像素相减得到，所以也能去除。得到的描述子向量为$H=(h_1,h_2,...,h_{128})$，归一化后的特征向量为$L=(L_1,L_2,....,L_{128})$，则

$$L_i=\frac{h_i}{\sqrt{{\sum }_{j=1}^{128}{h}_j}}j=1,2,3,...\text{\hspace{1em}(24)}$$

描述子向量门限。非线性光照，相机饱和度变化对造成某些方向的梯度值过大，而对方向的影响微弱。因此设置门限值(向量归一化后，一般取0.2)截断较大的梯度值。然后，再进行一次归一化处理，提高特征的鉴别性。

7.5 SIFT特征提取算法流程总结

经过上述对SIFT特征提取的详细分析与推导，SIFT特征提取算法流程总结如下：

图12 SIFT特征提取算法流程

8. SIFT特征匹配

8.1 特征匹配策略

  生成了两幅图的描述子$k1\ast 128$维和$k2\ast 128$维后，将两图中各个scale（所有scale）的描述子进行匹配，匹配上128维即可表示两个特征点匹配成功。

  实际计算过程中，为了增强匹配的稳健性，Lowe建议对每个关键点使用4×4共16个种子点来描述，这样对于一个关键点就可以产生128个数据，即最终形成128维的SIFT特征向量。此时SIFT特征向量已经去除了尺度变化、旋转等几何变形因素的影响，再继续将特征向量的长度归一化，则可以进一步去除光照变化的影响。当两幅图像的SIFT特征向量生成后，下一步我们采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量。取图像1中的某个关键点，并找出其与图像2中欧式距离最近的前两个关键点，在这两个关键点中，如果最近的距离除以次近的距离少于某个比例阈值，则接受这一对匹配点。降低这个比例阈值，SIFT匹配点数目会减少，但更加稳定。为了排除因为图像遮挡和背景混乱而产生的无匹配关系的关键点，Lowe提出了比较最近邻距离与次近邻距离的方法，距离比率ratio小于某个阈值认为是正确匹配。因为对于错误匹配，由于特征空间的高维性，相似的距离可能有大量其他的错误匹配，从而它的ratio值比较高。Lowe推荐ratio的阈值为0.8。但作者对大量任意存在尺度、旋转和亮度变化的两幅图片进行匹配，结果表明ratio取值在0. 4~0. 6之间最佳，小于0. 4的很少有匹配点，大于0. 6的则存在大量错误匹配点。作者建议ratio的取值原则如下:

• ratio=0. 4　对于准确度要求高的匹配；
• ratio=0. 6　对于匹配点数目要求比较多的匹配；
• ratio=0. 5　一般情况下。

也可按如下原则：当最近邻距离<200时，ratio=0. 6，反之ratio=0. 4。ratio的取值策略能排分错误匹配点。

8.2 实验演示

8.2.1 特征向量的生成

  如下图所示，图13（加有噪声）与图15为模板图像（参考图像），图17为待匹配图像。在正常尺度大小情况下，它们的特征点检测如下图14、图16与图18所示。

图13 参考图像1

图14 参考图像1的特征向量(859个)

图15 参考图像2

图16 参考图像2的特征向量(384个)

图17 待匹配图像

图18 待匹配图像的特征向量(2302个)

  当两幅图像的SIFT特征向量生成后，下一步我们采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量。取图像1中的某个关键点，并找出其与图像2中欧式距离最近的前两个关键点，在这两个关键点中，如果最近的距离除以次近的距离少于某个比例阈值（本文为0.6），则接受这一对匹配点。

8.2.2 尺度变换特征匹配

  为体现SIFT现特征描述子基于物体上的局部特征点生成而与物体的大小无关的特点，现将待匹配图像以不同大小比例尺度（0.1-0.5）进行缩放后再与参考图像进行匹配，遍历参考图像，找到匹配数比较合适的缩放级别，观察待匹配图像与哪一幅参考图像更加匹配，即哪两幅图像相似信息量更多。其匹配效果如下图18所示。

图19 实验演示

待匹配图像（scale=0.1）与参考图像1

待匹配图像（scale=0.2）与参考图像1

待匹配图像（scale=0.3）与参考图像1

待匹配图像（scale=0.4）与参考图像1

待匹配图像（scale=0.1）与参考图像2

待匹配图像（scale=0.2）与参考图像2

待匹配图像（scale=0.3）与参考图像2

待匹配图像（scale=0.4）与参考图像2

  从图18实验结果演示可以发现，通过缩放遍历，待匹配图像检测到的特征向量数逐渐增多，同时与参考图像的匹配对数也逐渐增多，缩放尺度scale=0.4时，参考图像共有384个特征向量，待匹配图像检测到472个特征特征向量，匹配数达到62对，此时已找到匹配数比较适合的缩放级别，完成特征匹配过程。

9. SIFT的缺点

  SIFT在图像的不变特征提取方面拥有无与伦比的优势，但并不完美，仍然存在：

• 实时性不高；
• 对边缘光滑的目标无法准确提取特征点；

近来不断有人改进，其中最著名的有SURF和CSIFT。

10. 结尾

  本篇文章详细介绍了SIFT算法的来源，特点，检测性能及具体检测步骤，并且详细讲解了SIFT算法理论原理。还根据检测到的SIFT特征向量进行了实验演示，验证了图像局部SIFT特征描述子具有较好的稳定性，独特性好，特征信息量丰富，匹配准确。最后简单地介绍了该算法存在的缺点。

参考资料

1. David G.Lowe Distinctive Image Features from Scale-Invariant Keypoints. January 5, 2004.
2. David G.Lowe Object Recognition from Local Scale-Invariant Features. 1999
3. Tony Lindeberg Scale-space theory: A basic tool for analysing structures at different scales.1994
4. SIFT特征点提取
5. SIFT算法详解
6. SIFT特征提取分析
7. 宋丹 10905056 尺度不变特征变换匹配算法Scale Invariant Feature Transform （SIFT）(PPT)
8. 线性插值与抛物线插值
9. Matlab图像处理学习笔记（六）：基于sift特征点的人民币识别