SIFT 特征点提取

SIFT 是一种从图像中提取独特不变特征的方法，其特点为基于图像的一些局部特征，而与图像整体的大小和旋转无关。并且该方法对于光照、噪声、仿射变换具有一定鲁棒性，同时能生成大量的特征点。

SIFT 的具体步骤

尺度空间极值检测：使用差分高斯函数识别潜在的兴趣点
特征点定位：剔除对比度不高和处于边界位置的特征点
分配方向：计算特征点的方向用于下一步构建描述
特征点描述：

尺度空间极值检测

尺度空间极值检测的作用就是发掘同一图像在不同尺度下都存在的特征点。通过对原始图像进行不断地降采样和平滑操作，生成一串不同尺度的图片，然后再找出每张图片都存在的特征点。

高斯尺度空间

作者在论文中提出，高斯核是唯一可以产生多尺度的核。对图像使用不同尺度空间因子（标准差）的高斯核进行卷积，能够得到不同模糊程度的图像，也就是将图像转换到高斯尺度空间。
L(x,y,σ)=G(x,y,σ)∗I(x,y)L(x, y, \sigma) = G(x, y,\sigma) * I(x, y) L(x,y,σ)=G(x,y,σ)∗I(x,y)
其中，G（x,y,σ）G（x,y,\sigma）G（x,y,σ）是高斯核函数:
G(x,y,σ)=12πσ2e−(x2+y2)/2σ2G(x, y, \sigma) = \frac{1}{2 \pi \sigma^2} e^{-(x^2 + y^2)/2\sigma^2} G(x,y,σ)=2πσ21e−(x2+y2)/2σ2

而 $ * $ 代表卷积， L(x,y,σ)L(x, y, \sigma)L(x,y,σ) 代表图像的高斯尺度空间。检测特征点的算子中较好的算子是高斯拉普拉斯算子 σ2∇2G\sigma^2\nabla^2Gσ2∇2G ，但是其运算代价较大，故使用差分高斯（DoG, Difference of Gaussina）来近似替代。
KaTeX parse error: No such environment: align at position 8: \begin{̲a̲l̲i̲g̲n̲}̲ D(x, y, \sigma…
在上述公式中， kkk 一个比例因子，通过将其设置为不同的值从而控制高斯核的大小。

高斯金字塔与高斯差分金字塔

为了从不同的空间尺度获取图片的特征点信息，首先要在同一图片大小上使用不同的高斯核进行卷积运算，再进行差分计算。然后再进行将图片降采样，到更小的尺度上再次进行高斯卷积运算，并同样进行差分计算。重复上述过程，即可得到高斯金字塔与高斯差分金字塔。

DoG 空间极值检测

即找出差分金字塔中，每一层中每一个点的 3×3×33\times 3 \times 33×3×3 邻域中的极值（包括其本身）。

参数取值

为了使每一层之间都只差一个比例因子 kkk ，同时高斯金字塔下方一组的最上层和上方一组的最下层也只差一个比例因子 kkk ，则设置
k=21sk = 2^{\frac{1}{s}} k=2s1
则每一组的最上方一层的比例因子即为 222，直接让对这层进行降采样，即可得到下一组最下方一层的原始图片，如图所示：

在 DoG 空间极值检测中可以发现，若最后指定需要寻找 sss 层的特征点，那么最上层与最下层由于在边缘，其邻域将会少一层。那么为了能够在这两层也能有完整的邻域，则高斯差分金字塔中每组都需要有 s+2 层，那么高斯金字塔则需要更多一层，即 s+3s+3s+3 层。那么在刚才的基础上，高斯金字塔需要再加上三层，如图：

需要注意的是，由于在最后又添加了两层，所以向下一组进行下采样的不是最后一层，而是倒数第三层了。

那么一共需要多少组呢，由于第一组的图片是原始尺寸，但每组都会使图片的尺寸缩小一倍，且在极值检测时还需要考虑 3 × 3 的邻域，则最后一组的图片尺寸需要大于 3，则有：
min⁡(w,h)2n−1≥3\frac{\min(w, h)}{2^{n-1}} \ge 3 2n−1min(w,h)≥3
其中 nnn 表示组数，则对两侧求对数，即可得到 nnn 的值：

int(floor(log2(min(w, h) * 2 / 3)))

则对于一张尺寸为 (w,h)(w, h)(w,h) 的图片，构建的高斯金字塔参数如下：

o = int(floor(log2(min(w, h) * 2 / 3)))         # 组数
s = 3                                          # 极值检测层数，作者设定为 3
ipo = S + 3                                       # 每组的层数
sigma_0 = 1.6                                  # 高斯核方差基数，作者设定为1.6

特征点定位

特征点定位中包含了三个操作

特征点真实坐标与真实值的计算
低对比度特征点消除
边缘效应消除

有限差分求导法

∂f∂x=f3−f12h∂f∂y=f4−f22h∂2f∂x2=f1+f3−2f0h2∂2f∂y2=f2+f4−2f0h2∂2f∂x∂y=(f5+f7)−(f6+f8)4h2\frac{\partial f}{\partial x} = \frac{f_3 - f_1}{2h}\\[1em] \frac{\partial f}{\partial y} = \frac{f_4 - f_2}{2h}\\[1em] \frac{\partial ^2 f}{\partial x^2} = \frac{f_1 +f_3 -2f_0}{h^2}\\[1em] \frac{\partial ^2 f}{\partial y^2} = \frac{f_2 +f_4 -2f_0}{h^2}\\[1em] \frac{\partial ^2 f}{\partial x \partial y} = \frac{(f_5 + f_7)-(f_6 + f_8)}{4h^2}\\[1em] ∂x∂f=2hf3−f1∂y∂f=2hf4−f2∂x2∂2f=h2f1+f3−2f0∂y2∂2f=h2f2+f4−2f0∂x∂y∂2f=4h2(f5+f7)−(f6+f8)

x−σx-\sigmax−σ，y−σy-\sigmay−σ 轴同理。

特征点定位与低对比度特征点消除

使用DoG空间极值检测所得到的特征点坐标是离散的，并不是极值点真正的坐标。这里使用了泰勒展开公式，对极值点进行拟合。已知二阶展开的泰勒展开式为：
f(x^)=f(x0)+f′(x0)(x^−x0)+f′′(x0)2(x^−x0)2f(\hat x) = f(x_0) + f'(x_0)(\hat x-x_0) + \frac{f''(x_0)}{2}(\hat x-x_0)^2 f(x^)=f(x0)+f′(x0)(x^−x0)+2f′′(x0)(x^−x0)2
现在使已得到的每个点为向量 x0={x,y,σ}x_0 = \{x, y, \sigma\}x0={x,y,σ} ，$\hat x $ 则为想要得到的极值的真实坐标，将泰勒展开改写为矩阵形式：
D(x^)=D(x0)+∂D(x0)T∂x^(x^−x0)+12(x^−x0)T∂2D(x0)∂x^2(x^−x0)D(\hat x) = D(x_0) + \frac{\partial D(x_0)^T}{\partial \hat x}(\hat x-x_0) + \frac{1}{2}(\hat x - x_0)^T\frac{\partial ^2 D(x_0)}{\partial \hat x^2}(\hat x-x_0) D(x^)=D(x0)+∂x^∂D(x0)T(x^−x0)+21(x^−x0)T∂x^2∂2D(x0)(x^−x0)
由于 x^\hat xx^ 处为极值，D(x^)D(\hat x)D(x^) 在此处的导数为 0，则

0=∂D(x0)∂x^+∂2D(x0)∂x^2(x^−x0)(x^−x0)=−∂2D(x0)∂x^2−1∂D(x0)∂x^0 = \frac{\partial D(x_0)}{\partial \hat x} + \frac{\partial ^2 D(x_0)}{\partial \hat x^2}(\hat x-x_0)\\[1em] (\hat x - x_0) = -\frac{\partial ^2 D(x_0)}{\partial \hat x^2}^{-1} \frac{\partial D(x_0)}{\partial \hat x} 0=∂x^∂D(x0)+∂x^2∂2D(x0)(x^−x0)(x^−x0)=−∂x^2∂2D(x0)−1∂x^∂D(x0)

则得到了极值点的坐标，接着将 $ \hat x - x_0$ 带入泰勒展开式即可得到：

D(x^)=D(x0)+12∂D(x0)T∂x^(x^−x0)D(\hat x) = D(x_0) + \frac{1}{2}\frac{\partial D(x_0)^T}{\partial \hat x}(\hat x-x_0) D(x^)=D(x0)+21∂x^∂D(x0)T(x^−x0)

显然 (x^−x0)(\hat x - x_0)(x^−x0) 表示偏移，若偏移的任何一项的绝对值超过了 0.5，说明该点已经偏移到了其他的点上，需要对该点的坐标进行更新。更新位置后再次重复以上操作，直到收敛（若偏移到图片外则舍弃该点）。

通过以上方法，就可以得到该极值点的真实值，该点含义就是通过差分操作得到的邻域中最具有特征的点，该点的值可以理解为对比度。而不是每个极值都具有大的对比度，如白墙上的一个不那么白的点，虽然它是极值点，但是很容易受到噪声的影响，所以设定一个阈值 TTT，当 ∣D(x^)∣<T|D(\hat x)| < T∣D(x^)∣<T 时，则认为该点的对比度不够大，则舍弃该点。

消除边缘效应

单纯地去除低对比度的点是不够的，考虑一条垂直的边界线，在横向上的变化很大，但是在纵向的变化却很小。这种点不能很好地反映特征，所以要把这些点给去除。去除的方法是使用黑塞矩阵，黑塞矩阵即为多元二阶导数矩阵，处于边缘的点特征值的比值较大，则特征值可以反应该点是否处于边缘 (此处不考虑 σ\sigmaσ轴 )。
H=[DxxDxyDxyDyy]H = \begin{bmatrix} D_{xx} &D_{xy}\\ D_{xy} &D_{yy} \end{bmatrix} H=[DxxDxyDxyDyy]
该黑塞矩阵有两个特征值 α,β\alpha,\betaα,β ，则有：
Tr(H)=Dxx+Dyy=α+βDet(H)=DxxDyy−Dxy2Tr(H) = D_{xx} +D_{yy} = \alpha + \beta \\ Det(H) = D_{xx} D_{yy} - D_{xy}^2 Tr(H)=Dxx+Dyy=α+βDet(H)=DxxDyy−Dxy2
假设 α>β\alpha > \betaα>β ，同时 α=γβ\alpha = \gamma \betaα=γβ ，此处需要的就是 γ\gammaγ 小于某个阈值，让两个特征值的比值不要太大，有
Tr(H)2Det(H)=(α+β)2αβ=(γβ+β)2γβ2=(1+γ)2γ\frac{Tr(H)^2}{Det(H)} = \frac{(\alpha + \beta)^2}{\alpha\beta} = \frac{(\gamma \beta + \beta)^2}{\gamma \beta^2} = \frac{(1+\gamma)^2}{\gamma} Det(H)Tr(H)2=αβ(α+β)2=γβ2(γβ+β)2=γ(1+γ)2
当比值最小时，γ=1\gamma = 1γ=1 ，则设置一个阈值，只要 γ\gammaγ 小于这个阈值即可说明该点非边缘位置。由于 $\gamma $ 增大时，Tr(H)2Det(H)\frac{Tr(H)^2}{Det(H)}Det(H)Tr(H)2 也增大，则只需要观测当前黑塞矩阵的 Tr(H)2Det(H)\frac{Tr(H)^2}{Det(H)}Det(H)Tr(H)2 是否小于 γ=\gamma =γ= 阈值时，(1+γ)2γ\frac{(1+\gamma)^2}{\gamma}γ(1+γ)2 的值即可，作者设置该阈值 γ=10\gamma = 10γ=10 ，即有效的特征点需要满足：
Tr(H)2Det(H)<(1+γ)2γ\frac{Tr(H)^2}{Det(H)}< \frac{(1+ \gamma)^2}{\gamma} Det(H)Tr(H)2<γ(1+γ)2

分配方向

获得特征点后，需要根据特征点附近的数据，来给特征点分配方向。

如图，在高斯金字塔中，以特征点像素左上角坐标为中心，2×3×σ2\times 3 \times \sigma2×3×σ 为边长，画出一个正方形，计算该正方形中每一个像素位置的梯度的方向角与幅值，计算公式如下：
m(x,y)=(L(x+1,y)−L(x−1,y))2+(L(x,y+1)−L(x,y−1))2θ(x,y)=arctan⁡((L(x,y+1)−L(x,y−1))/(L(x+1,y)−L(x−1,y)))m(x, y) = \sqrt{(L(x+1, y ) - L(x-1,y) )^2 + (L(x, y+1) -L(x, y-1))^2 } \\[1em] \theta(x, y) = \arctan((L(x, y+1) - L(x, y-1) ) / (L(x+1, y) -L(x-1,y))) m(x,y)=(L(x+1,y)−L(x−1,y))2+(L(x,y+1)−L(x,y−1))2θ(x,y)=arctan((L(x,y+1)−L(x,y−1))/(L(x+1,y)−L(x−1,y)))

将每个点的角度和幅值计算出来以后，再将幅值乘以一个权值作为新的幅值，从而反应与特征点的远近关系，如图中蓝色圈。离圆心越远权值越小，离圆心越近，权值越大。该权值的计算公式为：
Wx,y=−(x2+y2)2×1.5×σW_{x, y} = \frac{-(x^2 + y^2)}{2 \times 1.5 \times \sigma} Wx,y=2×1.5×σ−(x2+y2)
将每个点的幅值与该点对应的权值相乘得到新的幅值以后，就可以对这些点进行统计。统计的方式为，以这些点的角度作为横坐标，幅值与权值之积作为纵坐标。由于统计时角度是离散的。所以将角度分为 36 组，每 10° 为一组，如图：

图中仅分为了8组以供参考，实际操作过程中常分为 36 组。完成统计工作后，需要对统计的数据进行平滑操作，平滑操作的方法如公式：
H(i)=(6H(i)+4(H(i−1)+H(i+1))+H(i−2)+H(i+2))/12H(i) = (6H(i) + 4(H(i-1)+H(i+1)) + H(i-2)+H(i+2)) / 12 H(i)=(6H(i)+4(H(i−1)+H(i+1))+H(i−2)+H(i+2))/12
接着就可以在统计的结果中，找到幅值之和最大的那个方向，这个方向就是该特征点的主方向。而其他方向中，所有大于主方向幅值之和 80% 的方向，被称为辅方向。具有辅方向的特征点约占 15% 左右，但这些点对于特征匹配的稳定性起到了很大的作用。对于这种具有辅方向的特征点，我们需要额外创建一个位置相同的特征点，而其梯度数据则是辅方向的角度与幅值。

最后还需要做的是，将离散的角度进行连续化，让最接近每个峰值的3个统计值拟合成抛物线，以插值出峰值位置以获得更好的精度。拟合公式为：
B=i+H(i−1)−H(i+1)2(H(i−1)+H(i+1)−2H(i))B = i +\frac{H(i-1) - H(i+1)}{2(H(i-1) + H(i+1) -2H(i))} B=i+2(H(i−1)+H(i+1)−2H(i))H(i−1)−H(i+1)
将得出的值还原到 360° 的尺度上，并更新特征点的角度信息，这一部分就完成了。

特征点描述