高光谱拼接算法（四）SIFT 特征匹配

本篇代码仓位：SIFT

1. 特征描述

在上一篇中，我们完成了 SIFT 特征点检测的全部流程,得到了一批位置精确、尺度明确、对比度足够且非边缘的稳定特征点，每个特征点携带的信息是：

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

现在便正式来到下一部分：

拿到了这些关键点后，怎么在另一张图像中找到它？也就是构造描述子。

1.1 获取特征主方向

在 SIFT 的流程里，我们需要为每个特征点赋予一个主方向，这是为了让描述子具有旋转不变性。

我们直接展开具体过程来理顺其原理，首先：

以单个特征点为例，下一步是在该特征点对应的 Gaussian 图像 \(L(x,y,\sigma)\) 上，以特征点为中心取一个邻域窗口，对窗口内每个像素计算梯度幅值和方向。

说明两个细节：

1. 操作对象不是 DoG 金字塔，因为 DoG 是差分响应图，Gaussian 才具有原始图像信息。
2. 邻域窗口的大小一般以特征点尺度为基准，取约覆盖 \(16\sigma\) 大小的二维邻域。

现在再来看出现的两个新词：梯度幅值和梯度方向，我们同样通过中心差分近似计算，公式如下：

1. 梯度幅值：

\[m(x,y) = \sqrt{ (L(x+1,y)-L(x-1,y))^2 + (L(x,y+1)-L(x,y-1))^2 } \]

2. 梯度方向：

\[\theta(x,y) = \arctan\left( \frac{L(x,y+1)-L(x,y-1)} {L(x+1,y)-L(x-1,y)} \right) \]

虽然二者出现的有些突兀，但是从梯度本身的概念出现其实不难理解，这是为了得到特征点所代表的局部特征的变化强度和方向，用于后续的描述子构建。

1.1.1 方向直方图

现在，对于窗口内的每个像素，我们都有了一个响应值和一个代表方向的角度，SIFT 用直方图来统计这些方向，具体如下：

1. 设计直方图共 36 个 bin，每个 bin 覆盖 \(10^\circ\)（\(360^\circ/36=10^\circ\)）。
2. 每个像素根据其梯度方向投票到对应的 bin 中。
3. 投票权重不是简单的 +1，而是加权投票：

\[\text{贡献} = m(x,y) \times w(x,y) \]

其中， \(m(x,y)\) 是刚刚计算的梯度幅值，\(w(x,y)\) 则是以关键点为中心的高斯权重（取 \(\sigma=1.5\sigma_{scale}\)）。

需要补充的是，在实际实现中，SIFT 并不会将一个像素的全部贡献直接投给某一个方向 bin，而是采用线性插值进行分配。

这是因为梯度方向通常不是某个 bin 的中心值。如果全部投给其中一个 bin，当方向仅发生极小变化时（如 \(69^\circ\rightarrow71^\circ\)），投票对象便会突然切换，导致方向直方图出现不连续的跳变。
因此，SIFT 会根据梯度方向距离相邻两个 bin 中心的远近，将梯度贡献按比例分配给这两个 bin。例如，对于方向为 \(73^\circ\) 的像素，则总投票量会按照和区间中心的距离比例分别累加到 相邻两个方向 bin（如覆盖 \(60^\circ\sim70^\circ\) 与 \(70^\circ\sim80^\circ\) 的两个 bin） 中，而不是全部投给其中一个。

这样，即使梯度方向发生轻微变化，方向直方图也会平滑变化，从而减少量化误差，提高主方向估计的稳定性。

再解释两个设计细节：

1. 为什么用梯度幅值加权：梯度幅值代表"这个方向的变化有多强"。如果某个像素的梯度幅值很小，说明它附近的灰度变化很微弱，其方向信息受噪声影响较大，不应赋予过高的投票权重。
2. 为什么还要高斯加权：离关键点中心越远的像素，越容易受到图像变形、遮挡等因素的影响。高斯加权是为了确保关键点中心附近的信息最可靠，边缘区域仅供参考。

1.1.2 主方向与辅助方向

现在，刚刚得到的直方图的最高峰值对应的方向就是该特征点的主方向。

然后，为了获得更精确的角度值，SIFT 会对峰值及其相邻两个 bin 做抛物线插值：

\[\hat{\theta} = \theta_{peak} + \frac{h_{-1}-h_{+1}}{2(2h_{peak}-h_{-1}-h_{+1})} \times\Delta\theta \]

其中 \(h_{-1}\)、\(h_{peak}\)、\(h_{+1}\) 分别是峰值 bin 及其左右相邻 bin 的投票数，\(\Delta\theta=10^\circ\)。

举个例子：

步骤	计算公式 / 说明	代入数值	计算结果
1. 确定峰值与相邻值	提取峰值 bin 及其左右相邻 bin 的投票数	\(h_{-1} = 80\), \(h_{peak} = 100\), \(h_{+1} = 90\)	\(\theta_{peak} = 50^\circ\)
2. 计算分子	左侧投票数减去右侧投票数	\(h_{-1} - h_{+1}\)	\(80 - 90 = -10\)
3. 计算分母	峰值投票数的两倍减去两侧投票数之和，再乘 2	\(2 \times (2h_{peak} - h_{-1} - h_{+1})\)	\(2 \times (200 - 170) = 60\)
4. 计算修正偏移量	(分子 / 分母) × bin 宽度 \(\Delta\theta\)	\(\frac{-10}{60} \times 10^\circ\)	\(\approx -1.67^\circ\)
5. 得出精确角度	峰值角度 + 修正偏移量	\(50^\circ + (-1.67^\circ)\)	\(48.33^\circ\)

至此，每个特征点携带的信息扩展为四维：

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

但这还不够，如果一个关键点位于"十字交叉"这类双向结构上，可能有两个同样显著的方向。因此，SIFT 为这张情况进行了专门处理：

如果存在某个非主峰 bin，其值达到主峰值的 80% 以上，则为该关键点创建一个新的辅助方向。

也就是说：一个特征点可能会对应多组 \((x,\ y,\ \sigma,\ \theta)\) ， 这会在后续匹配中增加成功的概率。

到这里，我们为每个特征点扩展了一维方向信息，下一步，它就可以帮助我们生成具有旋转不变性的描述子。

1.2 特征描述子

现在我们有了位置、尺度、方向这些关键信息，最后一步就是：如何用一个向量描述这个关键点周围的信息？

SIFT 描述子的设计思路概括如下：

在特征点坐标周围取邻域，进行旋转归一获取旋转不变性，再次使用梯度方向直方图投票组成描述子，并进行合理的处理尽可能使其对光照影响不敏感。

1.2.1 生成描述子

现在，再次取约 \(16\sigma\) 大小邻域窗口，将其划分为 \(4\times4\) 个子区域（此区域数一般固定）。

下一步，对整个邻域中的像素旋转矩阵，以特征值为中心旋转到主方向 \(\theta\)：

\[\Delta x=x-x_0,\qquad \Delta y=y-y_0 \]

\[\begin{bmatrix} x' \\ y' \end{bmatrix} = \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix} \]

这样无论原始图像怎么旋转，描述子看到的都是"正面朝上"的局部区域。
这一步就是让匹配具有旋转不变性的关键，如果同一个局部区域在另一幅图中被旋转，其会被旋转归一化后再生成对应描述子，这样就抹平了旋转带来的影响。

继续，在每个子区域内，计算各像素的梯度方向并减去主方向，保证描述子始终以关键点主方向作为参考方向。继续投票到一个 8 个 bin 的方向直方图中，每个 bin 覆盖 \(45^\circ\)。

\[\theta' = \theta-\theta_{keypoint} \]

投票权重依然是梯度幅值乘以高斯权重，但这次高斯窗口覆盖整个区域，\(\sigma\) 取区域宽度的一半。

而和上面相同的，这里实际上也使用了插值。我们从一个问题出发：

如果某个像素刚好位于两个子区域的边界，或者梯度方向恰好落在两个方向 bin 的交界处，该如何投票？

如果直接把它全部投给其中一个子区域或一个方向 bin，那么像素只发生极小的位移，就可能导致投票对象突然改变，从而使描述子发生明显跳变，同样的问题也存在于方向统计中。

因此，SIFT 并不会将一个像素的贡献全部投给某一个子区域和某一个方向 bin，而是采用三线性插值（Trilinear Interpolation）进行加权分配。

一个像素的梯度贡献会按照距离大小，分别分配给邻近的 4 个子区域和每个子区域中相邻的 2 个方向 bin，最多影响 8 个统计单元。

现在，每个 \(4\times4\) 子区域投票完的直方图中每个 bin 都有对应投票结果，共 8 个值，组成一个 8 维向量，总共：

\[4\times4\times8 = 128 \]

这个128 维的特征向量就是 SIFT 描述子。

1.2.1 优化描述子

到这里还没有结束：原始 128 维向量还不能直接用，因为光照差异会导致梯度幅值整体放大或缩小。

所以首先要对结果做 L2 归一化：

\[\mathbf h_{norm} = \frac{\mathbf h}{\|\mathbf h\|_2} = \frac{\mathbf h}{\sqrt{\sum_{i=1}^{128}h_i^2}} \]

这样就能抵消光照带来的整体线性缩放问题。
但问题还没完全解决：如果图像局部出现过曝或过暗，某些方向上的梯度会异常大，归一化后依然会扭曲整体分布。

因此 SIFT 加了一步截断：

将 \(\mathbf h_{norm}\) 中大于 0.2 的元素截断为 0.2，然后再次 L2 归一化。

注意这是第二次 L2 归一化，因为截断后向量长度再次发生变化，需要重新归一化。
0.2 的作用就是限制单个方向上的最大贡献，避免光照饱和导致的梯度异常主导整个描述子。

最终，SIFT 特征点的完整信息如下：

\[(x,\ y,\ \sigma,\ \theta,\ \mathbf d_{128}) \]

这便是 SIFT 特征描述的完整逻辑，现在，我们就可以应用其进行匹配。

2. 特征匹配

有了描述子，特征匹配就转化为了一个向量匹配问题：

对于第一张图像中的每个特征点，在第二张图像中找到"最相似"的特征点。

由于 SIFT 描述子经过了两次 L2 归一化，各维度已经处于统一的尺度，因此可以直接使用欧氏距离衡量两个描述子的相似程度，距离越小，说明两者的局部结构越接近：

\[d(\mathbf d_i,\mathbf d_j) = \sqrt{\sum_{k=1}^{128}(d_i^{(k)}-d_j^{(k)})^2} \]

但就像给大模型搭 Agent，怎么使用这个公式同样关键：

2.1 暴力匹配

最简单的匹配方式就是暴力匹配：

对于图像 \(A\) 中的每个特征点，计算它与图像 \(B\) 中所有特征点的欧氏距离，取距离最小的那个作为匹配结果。

整个过程十分直接，因此也称为最近邻匹配（Nearest Neighbor Matching）。
其时间复杂度为：\(\mathcal O(N_A\times N_B)\)， 如果两幅图各有几千个特征点，总计算量通常在百万到千万次距离计算量级，现代 CPU 一般可以在几十毫秒内完成。

当特征点数量较少时，暴力匹配已经足够高效。但随着特征点数量不断增加，逐一计算所有距离的代价会迅速增长。因此，工程中通常会采用 FLANN 等近似最近邻搜索算法，来提高匹配速度。

不过，暴力匹配有一个明显的问题：

如果某个特征点在另一张图像中根本没有真正的对应点怎么办？

暴力匹配仍然会返回一个距离最小的结果，即使这个"最近邻"实际上只是所有候选中最不像那么不像的那个。

因此，SIFT 原论文还提出了一个的筛选策略：比率测试（Ratio Test）。

2.2 Ratio Test 比率测试

比率测试的逻辑并不复杂：

对于每个查询特征点，找到描述子空间中距离最近的两个候选点，其距离分别记为 \(d_1\) 和 \(d_2\)，仅当 \(\frac{d_1}{d_2}<T\) 时，才认为该匹配有效。

其中，\(d_1\) 为最近邻距离，\(d_2\) 为次近邻距离，论文中推荐 \(T=0.7\) 或 \(T=0.8\) 。

如图所示，从理解上来说，比率测试的逻辑十分巧妙：

如果一个特征点在另一张图像中有清晰的对应点，它与最近邻之间的距离一定远小于与其它候选点的距离；反之，如果根本不存在真正对应点，那么最近邻和次近邻之间的差别通常不会太大，因为它们都只是"矮子里拔将军"。

可以看到，Ratio Test 实际衡量的是最近邻是否具有明显的唯一性，而不是最近邻本身距离有多小。

因此，比率测试是 SIFT 匹配流程中关键的筛选环节，仅一次测试便可以过滤掉大量误匹配，这里简单测试一组结果如下：

可以发现，即使经过了比率测试，仍然可能存在误匹配。例如，两幅图像中不同位置恰好具有十分相似的局部纹理，它们的描述子距离依然可能非常接近，从而通过比率测试。

因此，后续通常还需要结合几何一致性约束进一步剔除误匹配，这也是后续算法要解决的问题了。