本文已参与「新人创作礼」活动，一起开启掘金创作之路

引言

SIFT算法是为了解决图片的匹配问题，想要从图像中提取一种对图像的大小和旋转变化保持鲁棒的特征，从而实现匹配。这一算法的灵感也十分的直观：人眼观测两张图片是否匹配时会注意到其中的典型区域（特征点部分），如果我们能够实现这一特征点区域提取过程，再对所提取到的区域进行描述就可以实现特征匹配了。于是问题就演变成了以下几个子问题：

应该选取什么样的点作为特征点呢？：人眼对图像中的高频区域更加的敏感，由此我们应该选择变化剧烈的边缘或者角点进行检测，这里我们选择检测角点。（这里的intuition不是很明白，但感觉来说可能是我们提取的是特征点进行匹配，而边缘往往是由多个点组成。）
怎样使得选取的特征点有尺度不变性呢？：使用高斯金字塔获取不同尺寸下的图像变体，由这些变体获得尺度不变特征。
怎样使得选取的特征点有缩放不变性呢？：使用高斯金字塔获取不同尺寸下的图像变体，每个变体里都提取特征点再缩放回原大小以获得图像不同尺寸下的特征点。
怎样使得选取的特征点有旋转不变性呢？：后续采取将特征点区域旋转到主方向的设定以获得旋转不变性。
如何描述特征点区域呢？：使用特征点区域内每个方向的梯度赋值，类似HOG算子。

以上几个问题在SIFT算法里都用了很有意思的trick，后续会一一介绍。

一、高斯金字塔

引入高斯金字塔的目的在引言中已经介绍过了--解决图片缩放及尺度变化下特征提取的问题，高斯金字塔的Intuition有两个：1. 人看物体时近大远小，可以对图片下采样实现（金字塔->组）；2. 人看物体时近处清晰，远处模糊，可以对图像高斯平滑实现（高斯->层）；具体的推导可以参见我的另一篇博客：

Opencv学习笔记（六）图像金字塔 - 掘金 (juejin.cn)

在SIFT里，高斯金字塔的层数和组数有着如下设定：
组数： $O=[log_2min(M,N)]-3$
层数： $S=n+3$
组数的设定是来自于提出SIFT算法的原始论文给出的经验值，理论上来说知要 $O\leq[log_2min(M,N)]$ 即可，层数的设定则是有着理论依据的，这里的 $n$ 是我们想要提取特征点的图片层数，由于提取出高斯金字塔后需要计算层间差分以获得高斯差分金字塔(DOG, Difference of Gausssian)，高斯金字塔层数需要比DOG层数多1，而计算特征值时要求在尺度层面，即上下相邻层间计算，则DOG层数要比特征层数多2，则要求 $S=n+3$ 。

附：在SIFT中对于高斯金字塔有几个需要补充的知识点。

SIFT算法中的S指的是"scale"即尺度的概念，这一概念有别于传统的图像尺寸，而是一种连续变化的参数，在高斯金字塔里的 $\sigma$ 就是这样一个尺度空间的参数，更特殊的，高斯核是唯一可以产生连续多尺度变化的线性核，这就是为什么我们用高斯金字塔，这一部分内容涉及到尺度空间的概念，可以自行了解。
图像金子塔一文中提到了两个基础参数，层间尺度变化系数 $k$ 和初始尺度 $\sigma_0$ ，在原始算法中规定 $k=2^{o+\frac{r}{n}}，r={0,1...,n+2}$ ，由于每组的起始层是由前一组倒数第三层降采样得到，实际上保证了每组开始的尺度为 $\sigma_0,2\sigma_0,3\sigma_0...$ 。此外我们提到 $\sigma_0=1.6$ ，而实际拍摄图片时相机已经原始景物进行了一次尺度变换，此变换的系数原文设定经验值为 $\sigma'=0.5$ ，由此我们得到 $\sigma_0=\sqrt{1.6^2-\sigma'^2}=1.52$ 。这一式子的来历为：高斯滤波体现为高斯核与原图进行卷积 $f(x)\bigotimes G'(x)$ ，两次高斯滤波的级联相当于 $f(x)\bigotimes (G'(x) \bigotimes G(x))$ ，而两个高斯卷积的结果仍为高斯卷积，且方差满足平方和关系，可用时域卷积或者傅里叶变换证明，详见参考文献。

二、高斯差分金字塔

通过高斯金字塔，我们获取了不同尺度的图片，接下来的问题是如何获取高频区域呢，一个很简单的思路就是按照边缘检测的算法使用差分滤波器如拉普拉斯滤波器、sobel滤波器在图片上滑动找到灰度值变化剧烈的区域。而经前人研究，归一化的高斯拉普拉斯算子的极大值极小值相较于其他特征提取函数可以获得最稳定的图像特征，因此我们打算使用归一化的高斯拉普拉斯算子在多尺度图片上提取特征，但是使用这种方式提取复杂度会很高，又尺度归一化高斯拉普拉斯算子和DOG函数有着如下关系: $G(x,y,k\sigma)-G(x,y,\sigma)\approx(k-1)\sigma^2 \nabla^2 G$ 证明如下： 忽略高斯函数系数:

\begin{aligned} &G(x,y,\sigma)=\frac{1}{\sigma^2}exp(-\frac{x^2+y^2}{2\sigma^2})\\ &\frac{\partial{G}}{\partial x}=-\frac{x}{\sigma^4}exp(-\frac{x^2+y^2}{2\sigma^2})\\ &\frac{\partial^2{G}}{\partial x^2}=\frac{-\sigma^2+x^2}{\sigma^6}exp(-\frac{x^2+y^2}{2\sigma^2})\\ &\nabla^2 G(x,y)=\frac{\partial^2{G}}{\partial x^2}+\frac{\partial^2{G}}{\partial y^2}\\ &\qquad \qquad \ =\frac{-2\sigma^2+x^2+y^2}{\sigma^6}exp(-\frac{x^2+y^2}{2\sigma^2})\\ &\frac{\partial{G}}{\partial \sigma}=\frac{-2\sigma^2+x^2+y^2}{\sigma^5}exp(-\frac{x^2+y^2}{2\sigma^2})\\ &\Rightarrow\\ &\sigma\nabla^2G=\frac{\partial G}{\partial \sigma} \end{aligned}

又对于差分高斯金字塔有：

\begin{aligned} &DOG=\frac{G(x,y,k\sigma)-G(x,y,\sigma)}{(k-1)\sigma}\approx \frac{\partial G}{\partial \sigma}\\ & DOG \approx(k-1)\sigma^2\nabla^2G \end {aligned}

由此我们不再需要在高斯金子塔上进行卷积操作，只需要计算在尺度层面 $\sigma$ 进行层间差分得到高斯差分金字塔(DOG)，即完成了特征提取的操作。在这里插入图片描述

三、特征点处理

得到DOG后，我们理论上来说就已经获得了特征值，但我们需要对特征值进行一些处理，就像我们在cannny边缘检测中使用的非极大值抑制等思路，去掉没那么特征的特征值。

1.阈值化

简单的阈值化去除掉变换没有那么剧烈的点，这些点就有可能是噪声引起的，算是在高斯拉普拉斯之外又加了一层去噪措施。 $val=\begin{cases} val & abs(val)>0.5\frac{T}{n}\\ 0 & otherwise \end{cases}$ 这里的 $T$ 为经验值0.04， $n$ 为之前提过的提取特征点数目。

2.非极大值抑制

非极大值一致的思路就和在其他算法中的一样，我们要求选取出的特征值应该是其领域范围内的极值，不同之处在于其他算法要求的只是图像这一二维平面内的极值，而这里要求还要在尺度 $\sigma$ 这一层面上也是极值，这也承接了前文DOG层数要比提取特征所用层数多2的设定。在这里插入图片描述

3. 二阶泰勒修正

由于我们的图片在 $x,y,\sigma$ 方向上都只能取到离散值，即使通过前两个步骤取到了一些特征点，但这些特征点都不够精确，我们需要引入二阶泰勒函数对齐进行修正，使得特征点可以出现在亚像素（亚尺度）位置。在这里插入图片描述

$f(X)=f(X_0)+\frac{\partial f^T}{\partial X}\hat{X}+\frac{1}{2}\hat{X^T}\frac{\partial^2 f}{\partial X^2}\hat{X}$

上式给出了特征点 $X_0$ 附近函数的近似值 $f(x)$ ，其中 $\hat{X}=X-X_0$ ，对该式求取一阶导数零点，即可得函数实际极值点位置与 $X_0$ 的距离，从而对离散特征点修正到亚尺度处。

$f'(X)=\frac{\partial f^T}{\partial X}+\frac{\partial f^2}{\partial X^2}\hat{X}$

$\hat{X}_{ex}=-\frac{\partial^2 f^{-1}}{\partial X^2}\frac{\partial f}{\partial X}$

带入前式可求取新的特征点的灰度值： $f(X')=f(X_0)+\frac{1}{2}\frac{\partial f^T}{\partial X'}\hat{X}_{ex}$

需要注意的是，上式是一个不断迭代的过程，也即根据当前特征点二阶泰勒求取新的特征点这一过程会不断重复直到满足终止条件，如 $\hat{X}$ 过小。另需注意当所得解超出离散极值点一定范围时需要舍去，因为二阶泰勒拟合只在其附近有效。

附：一个挺有意思的点在于这里的迭代目的和梯度下降法里的迭代并不是同一个目的，梯度下降里我们求取的是函数在某点的梯度，因为我们很难对复杂非线性的神经网络求除其一阶导数零点，即使求取到了也可能落入局部极值；而这里我们是直接求到了函数的一阶零点，迭代的目的是因为函数本身是一种近似，泰勒展开只取到了第二项，需要不断对原函数进行近似。这两种迭代的目的并不相同。

4.低对比度去除

目的和之间的阈值化类似，同样是去除掉没那么剧烈变化的特征点。要求： $|f(x)|\geq\frac{T}{n}$

5.边缘效应去除

引言中提过我们想要提取的特征点为角点而非边缘，而前述一系列措施只能保证取到灰度值变换剧烈的点，而边缘点同样符合这一特征，因此我们将通过以下方式去除边缘点。

计算黑森矩阵 $H(x,y)=\begin{bmatrix} D_{xx}(x,y)&D_{xy}(x,y)\\ D_{yx}(x,y)&D_{yy}(x,y) \end{bmatrix}$
若矩阵行列式 $Det(H)<0$ ，舍去该特征点
若矩阵行列式和迹不满足： $\frac{Tr(H)}{Det(H)}<\frac{(\gamma_0+1)^2}{\gamma_0}$ ，舍去该特征点， $\gamma_0$ 为有实际意义的经验值，通常设定为10。

接下来解释以下这几个步骤的意义，角点和边缘点的区别在于边缘在图像中表现为一条线，垂直于线的方向频率高，沿着线的方向频率比较低；而角点则在多（大于2）个方向方向出现强高频分量。黑森矩阵实际上是函数的二阶偏导构成的矩阵，可以反应函数的曲率变化状况。又对于二次型矩阵有如下性质：

假定二次型矩阵 $H$ 两个特征值为 $\alpha,\beta$ ，则 $Det(H)=\alpha\beta$ ， $Tr(H)=\alpha+\beta$ ;
实二次型矩阵的特征值异号时，该矩阵为不定矩阵，黑森矩阵为不定矩阵时，该临界点为非极值点；
黑森矩阵的特征值标定了函数在相应特征向量方向上变化的快慢。

由性质1，2我们可以推导处当 $Det(H)<0$ 时，特征点为非极值点，舍去；由性质1，3我们可以推导出当 $\frac{Tr(H)}{Det(H)}=\frac{(\alpha+\beta)^2}{\alpha^2\beta^2}$ 过小时，由两特征特征向量的比值 $\gamma$ 构成的式子 $\frac{(\gamma+1)^2}{\gamma}$ 同样较小，且对勾函数在 $\gamma>1$ 时单增，我们可以根据 $\frac{Tr(H)}{Det(H)}$ 大小判断特征向量的相对大小，该值过小，说明函数在该点不同方向上的变化非常不均匀，类似于边缘，舍去。

四、特征点描述子

通过前述步骤，我们已经获得了不同尺寸层面上的稳定特征点，接下来需要对其进行描述。

1. 确定特征点区域方向

引言中提到为了使特征点拥有旋转不变性，我们会将特征点区域统一旋转到特定方向，这一方向即为特征点区域的主方向，因此我们需要先确定主方向。

计算方式为：统计在离该特征点尺度 $\sigma^*$ 最近的尺度层 $\sigma_{oct}上，$ 以特征点为中心 $3\sigma'=3*1.5\sigma_{oct}$ 范围内的像素的梯度幅值及方向，对于范围内的梯度幅值用 $1.5\sigma^*$ 高斯核滤波以实现距离加权。得到一系列梯度幅值和方向对 $pairs=\{amp,ang\}$ 后，将 $360°$ 方向划分为多个bins，将包含相应 $ang$ 的 $pairs$ 中的 $amp$ 值累加到对应的bins上，如果 $ang$ 在两个bins划分间则根据距离分配 $amp$ 。（此处和HOG算法类似，可自行查阅）。最终可以获得特征点区域内幅值-方向直方图，选取其中幅值最大的方向作为主方向。在这里插入图片描述

附：这里的几个参数像是统计区域的半径，高斯加权的方差，有着不同的设置方式，但内在思路是一致的。此外如果在主方向之外出现了幅值达到了主方向幅值80%的其他方向，我们会将其作为辅方向，后续匹配时会出现两个位置、尺度相同但方向不同的两个特征点区域。

2. 特征点区域描述子

获得了特征点区域的主方向之后，我们就可以利用该值计算出有旋转不变性的描述子了。首先还是和前一步类似，统计该特征点所在尺度层面一定区域内的梯度幅值和方向，实施起来略有不同。

首先在特征点所处尺度层面内划定特征区域，半径为： $r=\frac{3\sqrt{2}\sigma_{oct}(d+1)}{2}$ ，其中 $d$ 为我们在一个维度上划分的子块数目，通常为4；
将该区域划分为共 $d*d$ 个子块，每个子块内包含多个像素点；
将划定出来的区域旋转到该区域的主方向（前一步已求出）；
在每个子块内构筑幅值-方向直方图，统计8个方向的梯度幅值。由此该区域可用 $8*d*d$ 维向量表示，由此完成了特征点区域描述子的构建；附：需要注意的是以上的计算都是发生在相应的高斯金字塔尺度层面上，而非原图或者DOG上，此外由于旋转产生的像素值丢失通过插值算法解决，可以自行了解。如果想在原图上可视化SIFT特征，需要将我们获得的稳定特征点坐标变换回原始的图像尺寸上，简单的乘以原下采样倍数即可。