高光谱拼接算法（三）SIFT 特征点检测

本篇内容关于SIFT 特征点检测，由于内容较多，本篇只包含 SIFT 的检测逻辑，特征匹配逻辑和代码会放在下一篇。

1. 从 Harris 到 SIFT

在上一篇中，我们详细介绍了 Harris 角点探测，其首次系统性地用结构张量来描述局部特征。
但其终究只是一个“原始基座”，一个很明显的问题是：

无人机从不同高度拍摄同一片农田，或者用不同分辨率传感器对同一区域成像时，同一个角点为什么就检测不到了？

这就是 Harris 最核心的问题：没有尺度不变性。
Harris 的窗口大小是固定的，比如 7×7。当我们把图像放大，原来 7×7 窗口里的细节就变了：一个角点在原图里看是角点，放大后可能变成了平滑的边缘，或者变成了更复杂的纹理。
尺度变了，特征就变了。

除此之外，Harris 还有另一个问题：

Harris 能告诉你"这里有个角点"，但说不出来"这个角点长什么样"，因此也就没法在另一张图像中找到同一个角点。

换句话说，它只有检测器，没有描述子，而拼接算法的核心步骤之一就是跨图像匹配同一特征。

因此，在 Harris 诞生后的十几年里，计算机视觉领域一直在寻找一种更好的局部特征，于是，而 1999 年 Lowe 提出了 SIFT（Scale-Invariant Feature Transform） 。其核心贡献在于：

它同时解决了"尺度不变性"和"特征描述"两个问题，并且对旋转、光照、仿射变化都有很强的鲁棒性。

2004 年，其正式发表定稿版论文 Distinctive Image Features from Scale-Invariant Keypoints，即使在今天 DL 方法盛行的背景下，SIFT 在图像拼接、SLAM、三维重建等领域依然被广泛使用。

2. 尺度空间理论（Scale Space Theory）

展开 SIFT 具体步骤前的一个基础问题是：

"尺度不变"这件事到底怎么做到的？

直观上，一个物体在不同尺度下看起来是不同的。远处的树是一个轮廓，近处的树能看到每一片叶子。
以此出发，问题就是：

能不能构造一种表示方法，让同一个特征在不同尺度下都能被稳定识别？

这就引出了尺度空间理论，其核心想法很简单：

引入一个连续参数 \(\sigma\)（尺度参数），通过逐渐模糊图像，模拟出"从近到远、从细节到轮廓"的视觉变化过程。

模糊的程度由 \(\sigma\) 控制：\(\sigma\) 越小，保留的细节越多；\(\sigma\) 越大，图像越模糊，只剩大尺度结构。

而实现这一模糊过程最常用的工具还是上一篇里见过的高斯卷积：

\[L(x,y,\sigma) = G(x,y,\sigma) * I(x,y) \]

其中：

\[G(x,y,\sigma)=\frac{1}{2\pi\sigma^2}\exp\left(-\frac{x^2+y^2}{2\sigma^2}\right) \]

概括来说：高斯核是唯一能保证"尺度空间生成"的线性核。 列举几个它的重要性质：

1. 半群性质：先用 \(\sigma_1\) 模糊，再用 \(\sigma_2\) 模糊，等价于直接用 \(\sigma=\sqrt{\sigma_1^2+\sigma_2^2}\) 模糊。这保证了不同尺度的变换可以叠加。
2. 非增局部极值：高斯模糊不会产生新的局部极值，只会消除已有的极值。这保证了"特征随尺度变化而消失"是单调的，不会凭空冒出新特征。
3. 各向同性：高斯核是圆对称的，不会引入方向偏好。

但其实本质上还是加权平均，来看一个例子：对于当前中心点 \((0,0)\) 和偏移点 \((2,0)\) 的权重比例：

\[\frac{G(2,0,\sigma)}{G(0,0,\sigma)} = e^{-\frac{4}{2\sigma^2}}= e^{-\frac{2}{\sigma^2}} \]

代几个值如下：

\(\sigma\)	\(\frac{G(2,0,\sigma)}{G(0,0,\sigma)}\)	含义
\(0.5\)	\(0.03\%\)	几乎只有中心像素参与计算
\(1\)	\(13.5\%\)	邻近区域开始产生影响
\(2\)	\(60.7\%\)	更大范围像素参与平均
\(4\)	\(88.2\%\)	远处像素与中心像素影响接近

显然， 随着 \(\sigma\) 增大，中心权重不断下降，说明像素值不再主要由自身决定，而是越来越多地受到周围更大范围像素的影响，因此图像会变得越来越模糊。

总结来说就是：

随着 \(\sigma\) 增大，高斯核的权重分布会逐渐变得平坦。更大范围内的像素都会参与平均，从而逐步“抹平”图像细节，只保留较大的结构轮廓。

这样，我们就能模拟出一个特征在不同分辨率，不同距离下的视觉效果。
尺度空间理论还有更多细节，之后遇到再展开。

3. 从 Harris 的角点到 SIFT 的尺度特征

现在我们知道了尺度空间的概念，那一个很容易得到的想法就是：

Harris 已经能检测角点了，SIFT 是不是直接把 Harris 搬到尺度空间里，对同一张图的多个尺度空间内检测到的角点进行分析处理来获得尺度不变性？

实际上，这种思路确实存在，例如后来的 Harris-Laplace ，就是在多尺度下运行 Harris 检测器，但其更像是一种“遍历”，并没有较根本地解决尺度不变性问题。

展开来说，由于不进行尺度空间的相关操作，Harris 输出实际上只是特征点坐标：\((x,y)\)。
而一个角点无论尺度如何变化，Harris 都可能给出较强响应。

那么问题来了：

哪一个尺度才是真正对应这个结构的尺度？或者说哪个尺度才是描述这个角点的最优尺度？

就像拍照时调整镜头焦距一样。焦距太近时，只能看到树叶的纹理；焦距太远时，整棵树又缩成一个轮廓。只有在某个合适的观察尺度下，目标结构才会呈现出最稳定、最容易识别的形态。

因此，为了真正获得尺度不变性，SIFT 的思路是这样的：

找到某个局部结构最自然、最稳定的观察尺度。只有确定了这个尺度，后续才能在与尺度相匹配的邻域内构造特征描述，从而实现不同分辨率图像之间的稳定匹配。

最终，SIFT 在这一环节的输出变成了：\((x,y,\sigma)\)。
这里的 \(\sigma\) 不是人为指定的参数，而是需要由算法自动发现的特征属性。是特征最稳定的尺度。

现在就到了下一个问题：

如何在尺度空间中，自动找到每个局部结构对应的最佳尺度 \(\sigma\)？

4. 高斯差分 DoG

我们已经知道：尺度不变性的关键在于找到局部结构最稳定的观察尺度。

现在同样对比来看：

Harris 本来就能检测特征点，为什么不直接用 Harris 在尺度空间里的最大响应值作为最佳尺度？

问题出在 Harris 所依赖的本质信息上，简单来说就是：最大相应值并不代表最优尺度，二者的语义不是对应的。

回顾结构张量：

\[M= \begin{bmatrix} I_x^2&I_xI_y\\ I_xI_y&I_y^2 \end{bmatrix} \]

其中：

\[I_x=\frac{\partial I}{\partial x}, \qquad I_y=\frac{\partial I}{\partial y} \]

全部来自图像的一阶导数，代表亮度变化有多剧烈。检测的是梯度能量聚集的位置。
如果尺度变化，那么 \(I_x\) 和 \(I_y\) 都会变化。最终响应值 \(R\) 也会变化。
问题是，这个变化没有一个统一规律。 即使在多个尺度下运行 Harris 并选取响应最大的结果，也不能保证找到结构真正对应的最佳尺度。

因此，我们需要一个具有明确尺度意义响应的检测器。

4.1 高斯拉普拉斯 LoG

理论证明，高斯拉普拉斯（Laplacian of Gaussian，LoG） 恰好具备这样的性质，其过程如下：
首先对图像进行我们上面提到的高斯平滑：

\[L(x,y,\sigma) = G(x,y,\sigma)*I(x,y) \]

然后计算拉普拉斯算子，即在各坐标方向的二阶偏导和：

\[\nabla^2L = \frac{\partial^2L}{\partial x^2} + \frac{\partial^2L}{\partial y^2} \]

合起来就是：

\[LoG = \nabla^2(G*I) \]

这是比较理论的步骤，实际计算中，会先对高斯函数求二阶导得到 LoG 核后截断，只对原图像做一次卷积得到响应值：

\[\nabla^2G = \frac{1}{\pi\sigma^4} \left( \frac{x^2+y^2}{2\sigma^2}-1 \right) e^{-\frac{x^2+y^2}{2\sigma^2}} \]

\[LoG = (\nabla^2G)*I \]

两者等价，而在真正应用时，还会加入一步尺度归一化：

\[LoG_{norm} = \sigma^2(\nabla^2G)*I \]

这是因为二阶导天然会随着尺度增大而衰减，例如比如同一个位置：

尺度	LoG
\(\sigma=1\)	50
\(\sigma=2\)	20
\(\sigma=4\)	5

即使结构完全一样，响应也会因为模糊越来越强而变小，这样我们根本没法比较不同尺度。
尺度空间理论证明：对于 \(n\) 阶导数，应乘 \(\sigma^n\) 进行归一化。

这便是 LoG 在理论上的完整步骤，现在再来看看原理。

4.2 为什么 LoG 响应峰值对应最优尺度？

LoG 本质上是对高斯平滑后的图像计算二阶导数，而二阶导数描述的并不是亮度大小，而是图像灰度变化本身的变化速度，也就是曲率。
因此 LoG 本质上是在找灰度变化最剧烈的局部特征，现在再联系到尺度：

假设图像中存在一个亮圆斑，我们现在在其中心使用一个很小的 LoG：

\[\sigma=1 \]

按照惯例，其对应检测器尺寸为：

\[W \approx 6\sigma + 1=7 \]

如果这个亮斑的尺寸很大，那么检测器只能看到局部区域，看不到亮斑边缘，此时检测器内部几乎全是亮区域，中心与周围差异不大。因此响应很弱：

\[LoG \approx 0 \]

继续增大尺度：

\[\sigma=3 \]

LoG 检测器同样变大，开始覆盖整个圆斑结构。
而 LoG 核具有典型的“墨西哥帽”结构：中心权重为正，外围环形区域权重为负：

此时中心正权重区域恰好覆盖亮斑主体。得到大量正贡献。
而外围负权重区域落在边缘背景，尺度合适的情况下，背景接近 0，负贡献就会很小。此时 \(LoG\) 达到极大值。

而如果尺度继续增大，外围负环就会开始覆盖大量亮区域和背景，正负贡献互相抵消，响应就会开始下降。
因此同一个位置就会在不同尺度下出现不同响应：

尺度	LoG响应
1	8
2	25
3	65
4	120
5	95
6	60
7	20

因此，LoG 值其实反应的是在当前尺度下，目标与背景之间的对比结构是否最明显。

这里很容易有这一个问题：

LoG 好像只是不断换一个更大的检测器去扫描图像，那为什么不直接遍历检测器尺寸，反而要依靠 \(\sigma\) 呢？

但这其实是搞错了主角，截断 LoG 的尺寸只是因为高斯核在 \(|3\sigma|\) 外无限接近 0，我们真正是在遍历 \(\sigma\) ，在不同的模糊程度下，去找到这个位置的局部结构和背景对比最明显的尺度。

4.3 实际方案：DoG

LoG 的理论很好，但就和 Harris 为什么要对结构张量进行工程优化一样，其计算代价相当高。
因此，我们在实际工程中使用的是其近似版本：高斯差分（Difference of Gaussian，DoG）
对于尺度空间：

\[L(x,y,\sigma) = G(x,y,\sigma)*I(x,y) \]

如果取两个相邻尺度：

\[L(x,y,k\sigma) \]

和

\[L(x,y,\sigma) \]

直接做差：

\[D(x,y,\sigma) = L(x,y,k\sigma)-L(x,y,\sigma) \]

看起来只是简单相减，但根据尺度空间理论可以证明：

\[D(x,y,\sigma) \approx (k-1)\sigma^2\nabla^2L \]

也就是说：

\[DoG \propto Normalized\ LoG \]

它实际上近似于尺度归一化后的 LoG，这也是 SIFT 中真正使用的方法，在下面的具体步骤中，我们就可以看到它的位置。

5. SIFT 特征点检测

下面就开始 SIFT 的实际构建了，这部分逻辑十分复杂，步骤也较多，但可以说是最标准化的流程，后续的几个出名算法其实本质上就是针对不同目标对 SIFT 的简化改进。

5.1 DoG 金字塔

一个很明显的道理是尺度空间不可能无限大，如果一直增大 \(\sigma\)，图像最终会被模糊成一团灰色，继续计算没有意义。

因此 SIFT 采用了一个方法：

当尺度增大到一定程度后，不再继续增大 \(\sigma\)，而是直接缩小图像。

一个需要注意的是原论文中，原始图像在开始处理前会先通过插值放大一倍来进行更密集的采样，而随着现代成像技术的提升，这点并不强制了。

回到主线，假如尺寸为 1024×1024 的初始图，在处理中会逐渐减小到 512×512，再到 256×256、128×128 等。而这样做的效果其实与继续增大尺度非常接近。

假设图像中有一个半径 \(r=64\) 的圆斑，在初始图中它占据很大面积，但如果直接缩小图像一半，那么圆斑半径也变成 \(32\)，相当于目标自动变小了。

其核心逻辑是：

目标变小与检测器变大，本质上都是在改变二者的相对尺度。

这便形成了 SIFT 的金字塔结构，其中存在两个核心设计：Octave 和 Interval。

SIFT 将尺度空间划分成多个组，每个组称为一个 Octave（八度）。
每经过一个 Octave，图像尺寸缩小一半：

\[1024 \rightarrow 512 \rightarrow 256 \rightarrow 128 \rightarrow \cdots \]

这样整个尺度空间被分解成多个层级，小尺度特征主要在前几个 Octave 中检测，而大尺度特征则在后几个 Octave 中检测。

但仅有 Octave 还不够，因为同一个 Octave 内仍然需要搜索最佳尺度，所以每个 Octave 内还会划分多个 Interval（尺度间隔）。

原论文中设：

\[s=3 \]

注意：这里的 \(s\) 并不是后面马上要提到的 Gaussian 层数或者 DoG 层数，而是和下一部分检测有关，它是指：

一个 Octave 中真正参与尺度搜索的 Interval 数量。

展开来说，为了让一个 Octave 恰好覆盖尺度翻倍：

\[\sigma \rightarrow 2\sigma \]

我们定义一个尺度倍率：

\[k=2^{1/s} \]

当 \(s=3\) 时：

\[k=2^{1/3} \approx1.26 \]

因此尺度会按照如下方式增长：

\[\sigma, \quad k\sigma, \quad k^2\sigma, \quad k^3\sigma=2\sigma \]

例如初始尺度取 \(\sigma=1.6\) ，则对应：

层号	\(\sigma\)
0	1.60
1	2.02
2	2.54
3	3.20

这样便定义出了一个 Octave 中的尺度划分方式。

接下来需要构建对应的高斯图像：

\[L(x,y,\sigma) = G(x,y,\sigma)*I(x,y) \]

如果只考虑这 3 个 Interval，那么理论上只需要生成对应尺度的 Gaussian 图像即可。

但后面还需要构建 DoG 图像，并进行三维尺度空间极值检测，因此 SIFT 会额外生成若干高斯层作为边界缓冲。

最终每个 Octave 实际会生成：

\[s+3 \]

个 Gaussian 图像。

当 \(s=3\) 时，生成结果就如下：

Gaussian层	\(\sigma\)	作用
\(L_0\)	1.60	有效尺度层
\(L_1\)	2.02	有效尺度层
\(L_2\)	2.54	有效尺度层
\(L_3\)	3.20	有效尺度层，达到 \(2\sigma\)，也是下一 Octave 的起始来源
\(L_4\)	4.03	边界缓冲层，用于构建 DoG
\(L_5\)	5.08	边界缓冲层，用于构建 DoG 与极值检测

最终，多个 Octave 的 Gaussian 图像共同组成高斯金字塔（Gaussian Pyramid）。

得到高斯金字塔后，在每个 Octave 内对相邻 Gaussian 图像做差：

\[D_i = L_{i+1} - L_i \]

例如：

Gaussian层	Gaussian尺度	DoG层	DoG对应尺度
\(L_0\)	\(\sigma_0\)	-	-
\(L_1\)	\(k\sigma_0\)	\(D_0=L_1-L_0\)	\(k\sigma_0\)
\(L_2\)	\(k^2\sigma_0\)	\(D_1=L_2-L_1\)	\(k^2\sigma_0\)
\(L_3\)	\(k^3\sigma_0\)	\(D_2=L_3-L_2\)	\(k^3\sigma_0\)
\(L_4\)	\(k^4\sigma_0\)	\(D_3=L_4-L_3\)	\(k^4\sigma_0\)
\(L_5\)	\(k^5\sigma_0\)	\(D_4=L_5-L_4\)	\(k^5\sigma_0\)

一般情况下：

\[(s+3) \text{ 个 Gaussian} \rightarrow (s+2) \text{ 个 DoG} \]

而这些 DoG 图像会共同组成 DoG 金字塔（Difference of Gaussian Pyramid）。

之所以需要额外生成两层 DoG，是因为后面寻找极值点时，需要同时比较当前层、上一层和下一层总计\(3\times3\times3\) 的三维邻域。
最顶层和最底层 DoG 无法参与极值检测，只能作为边界层存在，这样以来，最终真正参与尺度空间搜索的 DoG 层数恰好为：

\[(s+2)-2=s \]

这也正好对应最开始设定的 Interval 数量，这样，我们就完成了极值检测的所有准备。

5.2 尺度空间极值检测

这一步相对简单，得到 DoG 金字塔后，其最顶层和最底层 DoG 不参与检测，假设当前考察点位于：

\[D(x,y,\sigma) \]

现在，这个点周围有当前层 8 个邻居、上一层 9 个邻居以及下一层 9 个邻居，总计26个邻居。

SIFT 要求：

当前点必须同时大于（或小于）全部 26 个邻居。

遍历整个金字塔，这样的点才能被保留。这一步就是：尺度空间极值检测（Scale Space Extrema Detection）

也就是说：只有中心点同时成为空间和尺度两个方向上的极值，才能成为候选关键点。

此时得到的结果已经包含了 \((x,y,\sigma)\) ，前两维说明位置显著，第三维说明尺度显著，并且，我们还可以通过尺度系数和检测器大小的关系，估计出特征尺寸。

到这里，我们已经成功提取出了候选关键点，但 SIFT 仍然还有相当一部分内容，我们继续。

5.3 关键点精确定位

现在，我们已经获得了一批候选关键点，但这里有一个问题：

DoG 金字塔本质上是离散采样的尺度空间，而真实极值点的位置并不一定刚好落在采样网格上。

比如某个位置附近的 DoG 响应如下：

位置	DoG
\(x-1\)	80
\(x\)	100
\(x+1\)	95

按照规则： \(D(x)=100\) 已经是局部最大值。
但真正的峰值有可能位于：

\[x+0.7 \]

这时，更符合关键点要求的应该是 \(x+1\)，离散坐标导致了这一误差。

同样的问题也会出现在尺度方向：

\[\sigma_i \rightarrow \sigma_{i+1} \rightarrow \sigma_{i+2} \]

由于离散取值，真实最佳尺度可能位于两层之间。

所以，如果直接把离散极值点作为最终结果，会导致定位误差较大，且容易保留大量不稳定噪声点。

最终，从候选点到最终特征点间， SIFT 会再进一步进行亚像素级优化，具体可划分为三步：

5.3.1 亚像素定位

这步涉及到的数学原理和公式较多，这里先概述一下其核心逻辑：

使用泰勒展开和差分近似候选点附近的特征值，来估计真正极值点相对于当前采样点的偏移量，由此进行后续的修正和筛选。

展开来，设候选点为

\[\mathbf x_0=(x_0,y_0,\sigma_0)^T \]

则 DoG 函数在该点附近可以近似写成

\[D(\mathbf x) \approx D(\mathbf x_0) + \frac{\partial D^T}{\partial \mathbf x} (\mathbf x-\mathbf x_0) + \frac12 (\mathbf x-\mathbf x_0)^T \frac{\partial^2D}{\partial\mathbf x^2} (\mathbf x-\mathbf x_0) \]

记当前位置到真实极值点的偏移量为

\[\hat{\mathbf x} = \mathbf x-\mathbf x_0 \]

现在，这个式子可以写成：

\[D(\mathbf x) \approx D(\mathbf x_0) + \frac{\partial D^T}{\partial \mathbf x} \hat{\mathbf x} + \frac12 \hat{\mathbf x}^T \frac{\partial^2D}{\partial\mathbf x^2} \hat{\mathbf x} \]

关于 \(\hat{\mathbf x}\) 求导使结果为 0 ：

\[\frac{\partial D}{\partial\mathbf x} + \frac{\partial^2D}{\partial\mathbf x^2} \hat{\mathbf x} = 0 \]

整理即可得到真实极值相对于当前离散采样点的偏移量：

\[\hat{\mathbf x} =- \left( \frac{\partial^2D}{\partial\mathbf x^2} \right)^{-1} \frac{\partial D}{\partial\mathbf x} = (\hat x,\hat y,\hat\sigma)^T \]

我们以求解这个式子为目标继续，用离散采样点间的差分进行近似计算一阶导数和二阶导数，比如一阶偏导可由中心差分近似为：

\[\frac{\partial D}{\partial x} \approx \frac{D(x+1,y,\sigma)-D(x-1,y,\sigma)}2 \]

同理，\(\frac{\partial D}{\partial y},\frac{\partial D}{\partial\sigma}\) 也采用完全相同的方法计算。

二阶偏导则为：

\[\frac{\partial^2D}{\partial x^2} \approx D(x+1,y,\sigma) - 2D(x,y,\sigma) + D(x-1,y,\sigma) \]

\(\frac{\partial^2D}{\partial y^2},\frac{\partial^2D}{\partial\sigma^2}\) 也同理。
而混合偏导则采用二维中心差分：

\[\frac{\partial^2D}{\partial x\partial y} \approx \frac{ D(x+1,y+1,\sigma) - D(x+1,y-1,\sigma) - D(x-1,y+1,\sigma) + D(x-1,y-1,\sigma) }{4} \]

\(D_{x\sigma}, D_{y\sigma}\) 也同理。

最终，这些二阶偏导组成 DoG 函数的 Hessian 矩阵（注意不是 Harris）：

\[\frac{\partial^2D}{\partial\mathbf x^2} = H = \begin{bmatrix} D_{xx} & D_{xy} & D_{x\sigma}\\ D_{yx} & D_{yy} & D_{y\sigma}\\ D_{\sigma x} & D_{\sigma y} & D_{\sigma\sigma} \end{bmatrix} \]

结合一阶梯度向量：

\[\frac{\partial D}{\partial\mathbf x} = \mathbf g = \begin{bmatrix} D_x \\ D_y \\ D_\sigma \end{bmatrix} \]

即可计算亚像素偏移量：

\[\hat{\mathbf x} = - \left( \frac{\partial^2D}{\partial\mathbf x^2} \right)^{-1} \frac{\partial D}{\partial\mathbf x} = -H^{-1}\mathbf g \]

求得偏移量后，下一步就是判断当前离散采样点是否确实是距离真实极值最近的采样点。

这部分不复杂，如果：

\[|\hat x|, |\hat y|, |\hat\sigma| \le0.5 \]

这说明在三维采样空间里，真实极值仍距离当前采样点最近，此时就可以接受该定位结果。

而如果某一维度满足：

\[|\hat x|>0.5, \quad |\hat y|>0.5, \quad \text{或} \quad |\hat\sigma|>0.5 \]

则说明真实极值实际上已经更接近相邻的离散采样点，因此继续以当前点作为展开中心已经不再合适。

此时，SIFT 会移动到对应方向的邻居点作为新的展开中心，重新计算梯度、Hessian 矩阵以及偏移量，并重复上述过程，这时会有如下两种情况：

1. 在新展开点上的所有维度的偏移量都小于 \(0.5\)，使用该点进行后续操作。
2. 达到最大迭代次数仍未符合要求，删除当前候选点。

总结来说，亚像素定位这步的主要目的不是筛选候选点，而是以更高精度修正每个候选点的位置。
而更明显的筛选逻辑体现在下面两步：

5.3.2 剔除低对比度点

经过亚像素定位后，我们已经得到了更加精确的关键点位置，但这并不意味着这些关键点都是稳定可靠的。
比如，某个候选点的峰值十分平缓，很容易受到噪声、光照变化等因素影响，在另一张图像中甚至可能完全消失。

因此，SIFT 会进一步判断：

该关键点的 DoG 响应值是否足够大。

这步比较简单粗暴，直接利用上面计算的偏移量代入泰勒展开估计真实极值处的 DoG 响应值：

\[D(\hat{\mathbf x}) = D(\mathbf x_0) + \frac12 \left( \frac{\partial D}{\partial\mathbf x} \right)^T \hat{\mathbf x} \]

其中

\[\hat{\mathbf x} = H^{-1}\mathbf g \]

即可得到连续空间中的 DoG 极值响应。

随后就是阈值判断：

\[|D(\hat{\mathbf x})| < T_c \]

其中 \(T_c\) 表示对比度阈值，在原论文中，其默认值为 \(0.03\) 。
若响应值低于该阈值，则认为该极值仅由微弱纹理或噪声产生，稳定性较差。
因此直接删除该候选点。

经过这一步后，大量由噪声产生的弱响应关键点都会被过滤掉，仅保留具有明显局部结构的特征点。

5.3.3 剔除边缘响应点

完成低对比度筛选后，还有另一类熟悉的不稳定关键点：边缘响应点。

这部分和 上一篇 Harris 角点探测的形式有些相似，但目的不同：对于角点来说，无论图像沿哪个方向发生微小移动，都会导致响应迅速下降，因此定位十分稳定。
而对于边缘来说，只有垂直于边缘方向变化明显，沿着边缘方向移动时几乎没有变化，因此关键点的位置很容易发生漂移。

为了区分二者，SIFT 同样利用亚像素定位时计算的 Hessian 矩阵：

\[H = \begin{bmatrix} D_{xx} & D_{xy}\\ D_{yx} & D_{yy} \end{bmatrix} \]

这里仅保留空间方向的二阶导数，因为边缘判断只发生在当前尺度层内，不涉及尺度方向。
Hessian 矩阵的两个特征值 \(\alpha\) 和 \(\beta\)（令 \(\alpha > \beta\)）分别代表两个主曲率方向上的变化强度。对于边缘点，\(\alpha \gg \beta\)；对于非边缘点，两者大小相当。

同样为了避免直接计算特征值，我们利用迹和行列式来间接判断：

\[\operatorname{Tr}(H) = D_{xx}+D_{yy} = \alpha+\beta \]

\[\det(H) = D_{xx}D_{yy} - D_{xy}^2 = \alpha\beta \]

在这两个式子上，我们设计 \(r = \alpha/\beta\)，则：

\[\frac{\operatorname{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\) 有关，且 \(r\) 越大比值越大：

\(r\)	\(\frac{(r+1)^2}{r}\)	含义
1	4.0	两个方向曲率相等（理想角点）
5	7.2	一个方向曲率远大于另一个
10	12.1	非常明显的边缘响应
20	22.05	极强边缘响应

原论文取阈值 \(r=10\)，即：

\[\frac{\operatorname{Tr}(H)^2}{\det(H)} < \frac{(r+1)^2}{r} \]

不满足该条件的候选点就会被识别为边缘响应被去除。

经过这步保留下来的关键点就是最终特征点。
它们不仅具有较高的对比度，而且在两个方向上都具有良好的定位稳定性，更适合作为后续特征描述与匹配的基础。

本篇的内容先到这里，下一篇的内容，也是 SIFT 的后半部分关于在得到特征点的基础上，如何构造具有旋转不变性的描述子并在不同图像间进行匹配。