SIFT 尺度不变的特征变化匹配算法

参考：
博客
ppt
百度百科

一、 Introduction

事实上，SIFT特征就是对图像的一种描述，这种描述具有尺度不变性，通过在图像中检测出关键点，得到局部特征的描述子，得到图像的SIFT特征向量，作为该图像的描述。

SIFT算法具有这样的特性：

• SIFT特征是图像的局部特征，对尺度（scale），旋转（rotation）和亮度（Intensity）保持不变性，对视角变化，仿射变换和噪声也保持一定的不变性
• 区分性（Distinctiveness）好，信息量丰富，适用于在海量特征数据库中进行快速，准确的匹配
• 少数几个物体也可以产生大量的SIFT特征向量
• 优化的SIFT算法可以达到实时的要求
• 方便与其他形式的特征向量联合

基本步骤如下：

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

二、Step

更加细致的，我们可以具体分成下面这几个步骤：

• 构建尺度空间
• 获得高斯差分金字塔
• 确定DoG极值点
• 子像素定位潜在特征点
• 去除对比度低的点，降低边缘响应
• 建立关键点描述
• 去吧，皮卡丘！
  

1. 尺度空间极值检测

尺度空间理论是这样的：在图像空间中引入一个尺度的概念，使得原来的单尺度

1.1 尺度空间

一个图像的尺度空间 $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}(0)}$$
其中， $\ast$ 表示卷积运算，高斯函数表示如下：
$$G(x,y,\sigma )=\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{(x-m/2{)}^2+(y-n/2{)}^2}{2{\sigma }^2}}\text{\hspace{1em}(1)}$$
其中，$m,n$ 都表示为高斯模糊模板的维度，$\sigma$ 为尺度空间因子。也就是说，通过尺度空间因子的变化，我们可以得到模糊程度从大到小的一个不同尺度的尺度空间序列，称为高斯金字塔，塔的高度由图像的原始尺寸和塔顶模糊程度最大的图像的尺寸决定：
n = log ⁡ 2 { min ⁡ ( M , N ) } − t , t ∈ [ 0 , min ⁡ ( M , N ) ] (2) n=\log_{2}\{\min{(M,N)}\}-t,\quad t\in[0,\min{(M,N)}]\tag{2} n=log2​{min(M,N)}−t,t∈[0,min(M,N)](2)

1.2 高斯金字塔

由前面所讲得到的高斯金字塔如下：

Figure 1. 高斯金字塔

从上图可以看到，每层金子塔含有若干张图片，一层中的一组图片称为Octave，这是为了保证尺度的连续性，因此在一层之中使用了不同的尺度进行高斯滤波，使得在进行下采样得到下一层的尺度图片的时候跨度较小，通常是对倒数第三张图片进行下采样来得到下一层的第一张图片。

1.3 高斯差分金字塔

事实上，高斯差分算子能够得到更加稳定的图像特征，高斯差分算子和尺度归一化的拉普拉斯算子非常近似：

$$G(x,y,k\sigma )-G(x,y,\sigma )=\sigma {\nabla }^{2}G\text{\hspace{1em}(3)}$$
$$Laplacian=\sigma {\nabla }^{2}G\text{\hspace{1em}(4)}$$
如果我们用差分来近似表示微分的话，则有：
$$\frac{\partial G}{\partial \sigma }=\frac{G(x,y,k\sigma )-G(x,y,\sigma )}{k\sigma -\sigma }\approx \sigma {\nabla }^{2}G\text{\hspace{1em}(5)}$$
即我们可以得到拉普拉斯算子和高斯差分算子之间的关系：
$$G(x,y,k\sigma )-G(x,y,\sigma )\approx (k-1){\sigma }^2{\nabla }^{2}G\text{\hspace{1em}(6)}$$
其中 $k$ 是一个常数，并不会影响我们确定极值点，而求取极值点是我们找到关键点的前提步骤，而通过这样的替代我们得到了更好的效果。

最终我们就可以得到：
$$\begin{array}{rl}D(x,y,k\sigma )& =(G(x,y,k\sigma )-G(x,y,\sigma ))\ast I(x,y)\\ & =L(x,y,k\sigma )-L(x,y,\sigma )\end{array}\text{\hspace{1em}(7)}$$

1.4 空间极值点检测

构造好了尺度空间之后，接下来就是在尺度空间中进行极值点的检测了。所谓极值点指的是局部极值点，SIFT的关键点就是由尺度空间中的局部极值点构成的。

通常，寻找局部极值点很简单，只需要对每个像素和其8领域中的像素进行比较，判断是否为局部最大值即可。然而，在尺度空间中不是这样的，尺度空间中的局部极值点还需要和其相邻的两个尺度连续的图片对应的18个（上一层的图片中 9 个，加上下一层的图片中的 9 个）像素点的值进行比较，如下图所示：

即，上图中绿色的18个像素点。显然，这样的判断方式只能在一个Octave的中间层中进行，也就是在组内进行，因此通常需要 $S+2$ 个层，即，比所需层数多两层，来使得我们需要的每一层都能进行局部极值点的检测

1.5 尺度空间确定的参数

首先，我们需要这个图像是在哪个Octave中，在这个Octave中的哪一层，以及初始的尺度是多少，因此吧需要确定的参数有：

• $\sigma$：尺度空间坐标 ，决定模板的大小，即，$\sigma$ 的值越小，模板越小，从而平滑程度也就越小，模糊程度也就越小，相对的细节会比较多，反之亦然。通常模板的大小为 $(6\sigma +1)\times (6\sigma +1)$
• $O$：Octave的数量（多少组）
• $S$：组内的层数

其中，参数的关系如下：
$$\sigma (o,s)={\sigma }_0{2}^{o+\frac{s}{S}},o\in [0,\dots ,O-1],s\in [0.\dots ,S-1]\text{\hspace{1em}(8)}$$

其中：

• ${\sigma }_0$：是基准层尺度
• $o$：Octave的索引
• $s$：一个Octave内部的索引

通过确定初始的尺度，在哪一个Octave中的哪一层，根据上面的公式，我们就可以得到对应的每一层的尺度坐标。

而如果已知Octave，希望得到这个Octave的所有图像的时候，我们可以用下面的式子来实现同一个Octave中的图像的连续求取：
$$\sigma (s)=\sqrt{(k^s{\sigma }_0{)}^2-(k^{s-1}{\sigma }_0{)}^2}\text{\hspace{1em}(9)}$$

其中，$s-1$ 就是上一层的索引，然后上一层的图像通过上式得到的 $\sigma (s)$ 进行高斯模糊得到下一层的图像，因此通过这条式子，我们可以从这个Octave的第一层开始，连续求出Octave中的全部图像。

当然，你还可以通过下式直接求出某个Octave中具体一层的图像的尺度：
$$\sigma \_oct(s)={\sigma }_0{2}^{\frac{s}{S}},s\in [0,\dots ,S+2]\text{\hspace{1em}(10)}$$
但是，正如我们前面所说的，我们是对每个Octave中间的图像进行极值点的确定，而如果直接这样显然丢失了第一层，即 $\sigma ={\sigma }_0$ 的分析，因此我们考虑：

• 先对初始图像作一个 $\sigma =0.5$ 的高斯模糊
• 然后对图像进行上采样，得到尺度为初始图象两倍的图像，记尺度为 $\sigma ={\sigma }_{-1}$则此时 ${\sigma }_{-1}=1.0$

我们将此图像作为第一个Ocrtave的第一层图像，这样我们就能保证我们对所需要的每一层图像都进行极值点的确定。

2. 关键点定位

我们需要找到在这个尺度中的极值点，而显然这是在一个离散的空间中进行求取的，而离散空间中的极值点不一定就是精确的极值点，假如我们把这些离散的像素点看成对于一个三维二次曲线的采样的话，我们并不能保证值最大的点就一定对应这个曲线的一个极值点上，而只是在离散空间中是局部最大值。因此我们需要通过拟合一个三维二次曲线来确定真正的价值点的位置和尺度，并去除对比度低的关键点和不稳定的边缘响应点（DoG算子会产生较强的边缘响应）。

利用DoG函数在尺度空间中的泰勒展开式：
$$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}(11)}$$

令式子值为 0，得到：
$$\hat{x}=\frac{{\partial }^2{D}^{-1}}{\partial x^2}\cdot \frac{\partial D}{\partial x}\text{\hspace{1em}(12)}$$

即，根据 $(x,y,\sigma )$ 来获得极值点的位置。

2.1 去除对比度低的点

从前面的式子中我们可以得到对应极值点方程的值为：
$$D(\hat{x})=D+\frac{1}{2}\frac{\partial D^T}{\partial x}\hat{x}\text{\hspace{1em}(13)}$$

通过设定某个经验值，去除 $D$ 小于这个经验值的极值点。

2.2 降低边缘响应

另一方面，DoG算子会产生较强的边缘响应，如上图所示，因此需要去除不稳定的边缘相应点，通过一个 Hessian矩阵：
$$H=\left[\begin{array}{cc}{D}_{xx}& D_{xy}\\ D_{xy}& D_{yy}\end{array}\right]\text{\hspace{1em}(14)}$$
其中：

• $D_{xx}$：水平方向的梯度
• $D_{yy}$：垂直方向的梯度

假如我们令 $D_{xx}=\alpha ,D_{yy}=\beta$，则我们可以得到Hessian矩阵的对角线之和 $Tr$ 以及行列式 $Det$ 为：
$$\begin{gathered} Tr(H)=D_{xx}+D_{yy}=\alpha +\beta \\ Det(H)=D_{xx}{D}_{yy}-D_{xy}^2=\alpha \beta \end{gathered}$$

然后，我们假设 $\alpha =r\beta$，其中，$\alpha >\beta$，则有：
$$\frac{Tr(H{)}^2}{Det(H)}=\frac{(\alpha +\beta {)}^2}{\alpha \beta }=\frac{(r\beta +\beta {)}^2}{r{\beta }^2}=\frac{(r+1{)}^2}{r}$$
显然，当且仅当两个特征值相等，即 $r=1$ 的时候这个值最小，而值越大的时候，说明两个特征值之间的比值越大，也就是说，此时在某个方向上的梯度越大，而另外一个方向的梯度越小，显然这就是边缘的情况，也就是我们需要去除边缘响应的目标像素。因此，通过设定一个阈值 $r_0$，若：
$$\frac{Tr(H{)}^2}{Det(H)}>\frac{(r_0+1{)}^2}{r_0}$$

则剔除该点以降低边缘响应。

3. 方向确定

为每个模糊图像计算其梯度：

$$\{\begin{array}{l}m(x,y)=\sqrt{(L(x+1,y)-L(x-1,y){)}^2+(L(x,y+1)-L(x,y-1){)}^2}\\ \theta (x,y)={\tan }^{-1}\left[\frac{L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)}\right]\end{array}$$
其中：

• $L(x,y)$：关键点所在位置的尺度空间值
• $m(x,y)$：为梯度的大小
• $\theta (x,y)$：为梯度的方向

此外，需要使用 $\sigma =1.5$ 的高斯窗口对每个点的邻域范围内进行高斯模糊加成。
然后，在完成梯度的计算之后，对领域范围内的像素的梯度方向进行直方图统计，将 $0°\backsim 360°$ 的方向等分为 $36$ 份，然后为每个统计值大于最大值 $0.8$，即比直方图中值最高的统计量的 $80\%$ 的方向建立关键点（当然，包括了最大值的那个方向），通过赋予一些关键点（根据经验统计，大概 $15\%$ 左右的关键点会出现多个方向）多个方向可以明显的提高关键点匹配的稳定性（即在一个尺度位置上有两个不同方向的关键点）。

4. 关键点描述

现在，每个关键点包含的信息有：

• 尺度空间中的位置
• 尺度
• 方向

接下来需要为每个关键点建立一个描述符，用一组向量将这个关键点描述出来，使其具有：

• 旋转不变性
• （一定的）光照不变性
• 视角不变性

并且，这个描述独特性较高，可以提高我们匹配的正确率。
**SIFT关键点描述子是在关键点邻域中高斯图像梯度的一种统计结果的表示。**通过对关键点周围像素进行分快，计算块内的梯度直方图，生成具有独特性的向量，作为这个图像区域的一个抽象的表示，具有一定的独特性。

通常，在关键点的尺度空间内 $4\times 4$ 的窗口中计算 $8$ 个方向的梯度，即，最终得到一个 $4\times 4\times 8=128$ 维的向量：

• 首先，确定描述子需要的图像区域
• 将坐标轴旋转为关键点的方向
• 将邻域内的采样点分配到对应的子区域内，将子区域内的梯度值分配到8个方向上，计算其权值
• 插值计算每个种子点八个方向的梯度
• 对它们进行归一化处理
• 设置描述子向量阈值
• 按特征点的尺度对特征描述向量进行排序

最终，生成了我们的 SIFT 描述向量。