SIFT（Scale-invariant feature transform，尺度不变性特征变换）是一种检测局部特征的算法，该算法通过求一幅图中的特征点（interest points,or corner points）及其有关scale和orientation的描述子得到特征并进行图像特征点匹配，获得了良好效果，详细解析如下：

---

算法描述

整个算法分为以下几个部分：

1.构建尺度空间

尺度空间理论目的是模拟图像数据的多尺度特性，高斯卷积核是实现尺度变换的唯一卷积核，于是一副二维图像的尺度空间定义为：

$$L(x,y,\sigma) = G(x,y,\sigma)*I(x,y)$$ 其中，$G(x,y,\sigma)$是尺度可变的高斯函数，$G(x,y,\sigma)=\frac{1}{2\pi \sigma^2}e^{-(x^2+y^2)}/2\sigma^2$，$(x,y)$是空间坐标，$\sigma$是尺度坐标。$\sigma$大小决定图像的平滑程度，大尺度对应图像的概貌特征，小尺度对应图像的细节特征。大的$\sigma$值对应粗糙尺度(低分辨率)，反之，对应精细尺度(高分辨率)。

为了有效的在尺度空间检测到稳定的关键点，提出了高斯差分尺度空间（DOG scale-space）。利用不同尺度的高斯差分核与图像卷积生成。构造高斯差分尺度空间(DOG scale-space):

$$D(x,y,\sigma) = (G(x,y,k\sigma)-G(x,y,\sigma))*I(x,y) = L(x,y,k\sigma)-L(x,y,\sigma)$$

下图所示不同σ下图像尺度空间：

关于尺度空间的理解说明：$2kσ$中的2是必须的，尺度空间是连续的。在 Lowe的论文中 ，将第0层的初始尺度定为1.6（最模糊），图片的初始尺度定为0.5（最清晰）. 在检测极值点前对原始图像的高斯平滑以致图像丢失高频信息，所以 Lowe 建议在建立尺度空间前首先对原始图像长宽扩展一倍，以保留原始图像信息，增加特征点数量。尺度越大图像越模糊。

图像金字塔的建立：对于一幅图像I，建立其在不同尺度(scale)的图像，也成为子八度（octave），这是为了scale-invariant，也就是在任何尺度都能够有对应的特征点，第一个子八度的scale为原图大小，后面每个octave为上一个octave降采样的结果，即原图的1/4（长宽分别减半），构成下一个子八度（高一层金字塔）。

由图片size决定建几个塔，每塔几层图像(S一般为3-5层)。0塔的第0层是原始图像(或你double后的图像)，往上每一层是对其下一层进行Laplacian变换（高斯卷积，其中$σ$值渐大，例如可以是$σ, k*σ, k*k*σ$…），直观上看来越往上图片越模糊。塔间的图片是降采样关系，例如1塔的第0层可以由0塔的第3层downsample得到，然后进行与0塔类似的高斯卷积操作。

---

2.检测尺度空间极值点

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

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

---

3.精确定位极值点

通过拟和三维二次函数以精确确定关键点的位置和尺度（达到亚像素精度），同时去除低对比度的关键点和不稳定的边缘响应点(因为DoG算子会产生较强的边缘响应)，以增强匹配稳定性、提高抗噪声能力，在这里使用近似Harris Corner检测器。

(1) 空间尺度函数泰勒展开式如下：

$$D(x,y,\sigma) = D(x,y,\sigma) + \frac{\partial D^T}{\partial x} x + \frac{1}{2}x^T \frac{\partial^2 D^T}{\partial x^2} x$$ 求导，并令其为0，得到精确解$\hat x$： $$\hat x = - \frac{\partial^{2} D^{-1}}{\partial x^2} \frac{\partial D^T}{\partial x}$$

(2) 在已经检测到的特征点中,要去掉低对比度的特征点和不稳定的边缘响应点。去除低对比度的点：把上面公式代入上上面公式，即在DoG Space的极值点处$D(x)$取值，只取前两项可得：

$$D(\hat x) = D + \frac{1}{2} \frac{\partial D^T}{\partial x} \hat x$$ 若$|D(\hat x)| \geq 0.03$，该特征点就保留下来，否则丢弃。

(3) 边缘响应的去除。
一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的主曲率，而在垂直边缘的方向有较小的主曲率。主曲率通过一个2×2 的Hessian矩阵$H$求出：

$$H = \begin{bmatrix} D_{xx} & D_{xy} \\ D_{xy} & D_{yy} \end{bmatrix}$$ 导数由采样点相邻差估计得到。$D$的主曲率和H的特征值成正比，令$\alpha$为较大特征值，$\beta$为较小的特征值，则 $$Tr(H) = D_{xx}+D_{yy} = \alpha+\beta\\ Det(H) = D_{xx}D_{yy}-(D_{xy})^2 = \alpha\beta$$ 令$\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}$$ $\frac{(\gamma+1)^2}{\gamma}$的值在两个特征值相等的时候最小，随着$\gamma$的增大而增大，因此，为了检测主曲率是否在某域值$\gamma$下，只需检测 $$\frac{Tr(H)^2}{Det(H)} = \frac{(\gamma+1)^2}{\gamma}$$ if $\frac{(\alpha+\beta)}{\alpha\beta}> \frac{(\gamma+1)^2}{\gamma}$, throw it out. 在Lowe的文章中，取$r＝10$。

---

4.为每个关键点指定方向参数

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

$$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)=a\tan2\frac{(L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)}$$ 式中为$(x,y)$处梯度的模值和方向公式。其中$L$所用的尺度为每个关键点各自所在的尺度。
至此，图像的关键点已检测完毕，每个关键点有三个信息：位置、所处尺度、方向。由此可以确定一个SIFT特征区域。

在实际计算时，我们在以关键点为中心的邻域窗口内采样，并用直方图统计邻域像素的梯度方向。梯度直方图的范围是0～360度，其中每45度一个柱，总共8个柱, 或者每10度一个柱，总共36个柱。Lowe论文中还提到要使用高斯函数对直方图进行平滑，减少突变的影响。直方图的峰值则代表了该关键点处邻域梯度的主方向，即作为该关键点的方向。

---

5.关键点描述子的生成

首先将坐标轴旋转为关键点的方向，以确保旋转不变性。以关键点为中心取8×8的领域作为采样窗口。将采样点和特征点的相对方向通过高斯加权后归入包含8个方向直方图，最后获得2*2*8的32维特征描述子。示意图如下。

图左部分的中央为当前关键点的位置，每个小格代表关键点邻域所在尺度空间的一个像素，利用公式求得每个像素的梯度幅值与梯度方向，箭头方向代表该像素的梯度方向，箭头长度代表梯度模值，然后用高斯窗口对其进行加权运算。

图中蓝色的圈代表高斯加权的范围（越靠近关键点的像素梯度方向信息贡献越大）。然后在每4×4的小块上计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，即可形成一个种子点，如图右部分示。此图中一个关键点由2×2共4个种子点组成，每个种子点有8个方向向量信息。这种邻域方向性信息联合的思想增强了算法抗噪声的能力，同时对于含有定位误差的特征匹配也提供了较好的容错性。

这样就可以对每个feature形成一个4*4*8=128维的描述子，每一维都可以表示4*4个格子中一个的scale/orientation。将这个向量归一化之后，就进一步去除了光照的影响。

---

6.128维关键点描述子的生成

1.确定计算描述子所需的图像区域

描述子梯度方向直方图由关键点所在尺度的模糊图像计算产生。图像区域的半径通过下式计算：

$$radius=\sqrt{\frac{3\sigma_{oct} \times \sqrt2 \times (d+1)+1}{2}}$$ $\sigma_{oct}$是关键点所在组(octave)的组内尺度，$d=4$。

2.将坐标移至关键点主方向

那么旋转角度后新坐标为

$$\begin{pmatrix} \hat x \\ \hat y \end{pmatrix} = \begin{pmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{pmatrix} \times \begin{pmatrix} x \\ y \end{pmatrix}$$

3.生成方向直方图

在图像半径区域内对每个像素点求其梯度幅值和方向，然后对每个梯度幅值乘以高斯权重参数，生成方向直方图。

$$weight = |grad(I_{\sigma}(x,y))| \times exp(-\frac{x_k^2 + y_k^2}{2\sigma_w}) \times (1-d_r) \times (1-d_c) \times (1-d_o)$$ $x_k$：该点与关键点的列距离；
$y_k$：该点与关键点的行距离；
$\sigma_w$：等于描述子窗口宽度$3\sigma \times$ 直方图列数(取4)的一半。

4.生成种子点

在窗口宽度为$2 \times 2$的区域内计算8个方向的梯度直方图，绘制每个梯度方向的累加值，即可形成一个种子点，然后再在下一个$2 \times 2$的区域内进行直方图统计，形成下一个种子点，共生成16个种子点。

5.描述子向量元素门限化

方向直方图每个方向上梯度幅度限制在一定门限值以下（门限一般取0.2）。

6.描述子向量元素规范化

$W = (w_1,w_2,...,w_{128})$ 为得到的128描述子向量
$L = (l_1,l_2,...,l_{128})$ 为规范化后的向量
$l_j = w_j / \sqrt{\sum_{i=1}^{128}w_i}$，其中$j=1,2,...,128$

---