【SIFT介绍】Scale-Invariant Feature Transform——尺度不变特征变换（二）

---

【SIFT介绍】Scale-Invariant Feature Transform——尺度不变特征变换（二）

---

简介

许多实际应用需要在一张或多张图像中定位参考位置，例如图像对齐、去除畸变、物体跟踪、3D重建等。我们已经看到，角点可以相当可靠地定位，并且不依赖于方向。然而，典型的角点检测器仅提供每个候选点的位置和强度，它们不提供任何有关其特征或“身份”的信息，这些信息可以用于匹配。另一个限制是大多数角点检测器仅在特定尺度或分辨率下工作，因为它们基于一组固定的滤波器。

本章介绍了局部特征检测的尺度不变特征变换（SIFT）技术，该技术最初由D. Lowe [152]提出，并自此成为成像行业的“主力”方法。其目标是定位能够鲁棒识别的图像特征，以便在多幅图像和图像序列中进行匹配，以及在不同视角条件下进行物体识别。SIFT采用了“尺度空间”的概念[151]，在多个尺度层次或图像分辨率下捕捉特征，这不仅增加了可用特征的数量，还使得该方法对尺度变化具有高度的容忍度。这使得在例如物体向相机移动并因此连续改变其尺度的情况下跟踪特征成为可能，或将使用不同变焦设置拍摄的图像拼接在一起。

通过简化尺度空间计算和特征检测或使用GPU硬件[20, 90, 218]，已经实现了SIFT算法的加速变种。

原则上，SIFT的工作方式类似于多尺度角点检测器，具有亚像素定位精度，并为每个候选点附加了旋转不变的特征描述符。这个（通常是128维的）特征描述符总结了对应特征点周围空间邻域中的梯度方向分布，因此可以像“指纹”一样使用。SIFT特征计算涉及的主要步骤如下：

1. 在拉普拉斯-高斯（LoG）尺度空间中检测极值点，以定位潜在的兴趣点。
2. 通过拟合连续模型来精确确定位置和尺度，从而对关键点进行优化。
3. 通过周围图像梯度方向的主导方向为特征点分配方向。
4. 通过归一化局部梯度直方图来形成特征描述符。

这些步骤都将在本章的其余部分详细描述。

我们在这里如此详细地解释SIFT技术有几个原因。首先，这是迄今为止我们所讨论的最复杂的算法，其各个步骤精心设计并相互依赖，涉及众多需要考虑的参数。因此，深入理解其内部工作原理和局限性对于成功使用以及在结果不如预期时分析问题非常重要。

---

2 关键点的选择和优化

关键点通过三个步骤识别：（1）在DoG尺度空间中检测极值点，（2）通过局部插值进行位置优化，以及（3）消除边缘响应。以下详细介绍这些步骤，并在算法25.3-25.6中总结。

2.1 局部极值检测

在第一步中，候选兴趣点被检测为我们在上一节描述的3D DoG尺度空间中的局部极值。极值检测在每个八度音阶$p$中独立进行。为了方便，我们定义3D尺度空间坐标$\mathbf{c}=(u,v,q)$，由空间位置$(u,v)$和层次索引$q$组成，以及函数
$$D(\mathbf{c}):={\mathbf{D}}_{p,q+k}(u,v)\text{\hspace{1em}(25.51)}$$
作为从给定八度音阶$p$中选择DoG值的简写表示法。此外，为了收集尺度空间位置$\mathbf{c}$周围3D邻域的DoG值，我们定义映射
$${\mathsf{N}}_{\mathbf{c}}(i,j,k):=D\left(\mathbf{c}+i\cdot {\mathbf{e}}_i+j\cdot {\mathbf{e}}_j+k\cdot {\mathbf{e}}_k\right)\text{\hspace{1em}(25.52)}$$
其中 i , j , k ∈ { − 1 , 0 , 1 } i, j, k \in \{-1, 0, 1\} i,j,k∈{−1,0,1}以及3D单位向量
$${\mathbf{e}}_i=(1,0,0{)}^{\top },{\mathbf{e}}_j=(0,1,0{)}^{\top },{\mathbf{e}}_k=(0,0,1{)}^{\top }\text{\hspace{1em}(25.53)}$$
邻域${\mathsf{N}}_{\mathbf{c}}$包括中心值$D(\mathbf{c})$及其26个紧邻的值（见图25.15(a)）。这些值将用于估计尺度空间位置$c$的3D梯度向量和Hessian矩阵，如下所述。

如果DoG尺度空间位置$\mathbf{c}$的关联值$D(\mathbf{c})={\mathsf{N}}_{\mathbf{c}}(0,0,0)$为负且小于所有相邻值，或者为正且大于所有相邻值，则$\mathbf{c}$被接受为局部极值（最小值或最大值）。此外，可以指定最小差异$t_{extrm}\ge 0$，指示中心值必须至少偏离周围值的程度。因此，判断给定邻域${\mathsf{N}}_{\mathbf{c}}$是否包含局部最小值或最大值可以表示为
$$\mathsf{IsLocalMin}({\mathsf{N}}_{\mathbf{c}}):={\mathsf{N}}_{\mathbf{c}}(0,0,0)<0\wedge {\mathsf{N}}_{\mathbf{c}}(0,0,0)+t_{extrm}<\underset{(i,j,k)\ne (0,0,0)}{\min }{\mathsf{N}}_{\mathbf{c}}(i,j,k),\text{\hspace{1em}(25.54)}$$
$$\mathsf{IsLocalMax}({\mathsf{N}}_{\mathbf{c}}):={\mathsf{N}}_{\mathbf{c}}(0,0,0)>0\wedge {\mathsf{N}}_{\mathbf{c}}(0,0,0)-t_{extrm}<\underset{(i,j,k)\ne (0,0,0)}{\max }{\mathsf{N}}_{\mathbf{c}}(i,j,k)\text{\hspace{1em}(25.55)}$$
（见算法25.5中的过程$\mathsf{IsExtremum}({\mathsf{N}}_{\mathbf{c}})$）。如图25.15(b-c)所示，可以为极值检测指定具有18或10个单元的替代3D邻域。

图 25.15
用于检测 DoG 尺度空间中局部极值的不同 3D 邻域。红色立方体表示参考坐标 $\mathbf{c}=(u,v,q)$ 处的 DoG 值，该坐标位于空间位置 $(u,v)$ 和尺度层次 $q$ （在某个八度音阶 $p$ 内）。完整的 3 × 3 × 3 邻域包含 26 个元素（a）；其他类型的邻域分别包含 18 个（b）或 10 个（c）元素，也常被使用。如果中心的 DoG 值大于/小于所有相邻值（绿色立方体），则检测到局部最大值/最小值。

2.2 位置优化

在DoG尺度空间中检测到局部极值后，只知道其离散的3D坐标$\mathbf{c}=(u,v,q)$，由空间网格位置$(u,v)$和相关尺度层次的索引$q$组成。在第二步中，通过拟合二次函数到局部邻域来估计每个候选关键点的更精确的连续位置，如[37]中所提议。这在尺度空间的较高八度音阶尤为重要，因为连续抽取导致空间分辨率越来越粗糙。位置优化基于离散DoG函数的局部二阶泰勒展开式，产生一个连续近似函数，其最大值或最小值可以解析地找到。附录C.3.2部分提供了更多细节和说明性示例。

在分层DoG尺度空间$\mathbf{D}$的八度音阶$p$中任一极值位置$\mathbf{c}=(u,v,q)$处，使用相应的3×3×3邻域${\mathcal{N}}_D(c)$来估计连续3D梯度的元素，即
$${\nabla }_D(\mathbf{c})=\left(\begin{array}{c}dx\\ dy\\ d\sigma \end{array}\right)\approx \frac{1}{2}\cdot \left(\begin{array}{c}D(\mathbf{c}+{\mathbf{e}}_i)-D(\mathbf{c}-{\mathbf{e}}_i)\\ D(\mathbf{c}+{\mathbf{e}}_j)-D(\mathbf{c}-{\mathbf{e}}_j)\\ D(\mathbf{c}+{\mathbf{e}}_k)-D(\mathbf{c}-{\mathbf{e}}_k)\end{array}\right),\text{\hspace{1em}(25.56)}$$
其中$D()$如公式(25.51)所定义。同样，位置$\mathbf{c}$的3×3 Hessian矩阵获得如下
$${\mathbf{H}}_D(\mathbf{c})=\left(\begin{array}{ccc}{d}_{xx}& d_{xy}& d_{x\sigma }\\ d_{xy}& d_{yy}& d_{y\sigma }\\ d_{x\sigma }& d_{y\sigma }& d_{\sigma \sigma }\end{array}\right),\text{\hspace{1em}(25.57)}$$
所需的二阶导数估计为
$$\begin{array}{rl}{d}_{xx}& =D(\mathbf{c}-{\mathbf{e}}_i)-2\cdot D(\mathbf{c})+D(\mathbf{c}+{\mathbf{e}}_i),\\ d_{yy}& =D(\mathbf{c}-{\mathbf{e}}_j)-2\cdot D(\mathbf{c})+D(\mathbf{c}+{\mathbf{e}}_j),\\ d_{\sigma \sigma }& =D(\mathbf{c}-{\mathbf{e}}_k)-2\cdot D(\mathbf{c})+D(\mathbf{c}+{\mathbf{e}}_k),\\ d_{xy}& =\frac{D(\mathbf{c}+{\mathbf{e}}_i+{\mathbf{e}}_j)-D(\mathbf{c}-{\mathbf{e}}_i+{\mathbf{e}}_j)-D(\mathbf{c}+{\mathbf{e}}_i-{\mathbf{e}}_j)+D(\mathbf{c}-{\mathbf{e}}_i-{\mathbf{e}}_j)}{4},\\ d_{x\sigma }& =\frac{D(\mathbf{c}+{\mathbf{e}}_i+{\mathbf{e}}_k)-D(\mathbf{c}-{\mathbf{e}}_i+{\mathbf{e}}_k)-D(\mathbf{c}+{\mathbf{e}}_i-{\mathbf{e}}_k)+D(\mathbf{c}-{\mathbf{e}}_i-{\mathbf{e}}_k)}{4},\\ d_{y\sigma }& =\frac{D(\mathbf{c}+{\mathbf{e}}_j+{\mathbf{e}}_k)-D(\mathbf{c}-{\mathbf{e}}_j+{\mathbf{e}}_k)-D(\mathbf{c}+{\mathbf{e}}_j-{\mathbf{e}}_k)+D(\mathbf{c}-{\mathbf{e}}_j-{\mathbf{e}}_k)}{4}.\end{array}\text{\hspace{1em}(25.58)}$$

参见算法25.5中的程序Gradient(${\mathsf{N}}_{\mathbf{c}}$)和Hessian(${\mathsf{N}}_{\mathbf{c}}$)以获取更多详细信息。根据梯度向量${\nabla }_D(\mathbf{c})$和Hessian矩阵${\mathbf{H}}_D(\mathbf{c})$，点$\mathbf{c}$周围的二阶泰勒展开式为
$${\tilde{D}}_{\mathbf{c}}(\mathbf{x})=D(\mathbf{c})+{\nabla }_D^{\top }(\mathbf{c})\cdot (\mathbf{x}-\mathbf{c})+\frac{1}{2}(\mathbf{x}-\mathbf{c}{)}^{\top }\cdot {\mathbf{H}}_D(\mathbf{c})\cdot (\mathbf{x}-\mathbf{c}),\text{\hspace{1em}(25.59)}$$
其中连续位置$\mathbf{x}=(x,y,\sigma {)}^{\top }$。标量值函数${\tilde{D}}_{\mathbf{c}}(\mathbf{x})\in \mathbb{R}$，其中$\mathbf{c}=(u,v,q{)}^{\top }$和$\mathbf{x}=(x,y,\sigma {)}^{\top }$，是八度音阶$p$、尺度层次$q$和空间位置$u,v$处离散DoG函数$D_{p,q}(u,v)$的局部连续近似函数。这是一个二次函数，在位置
$$\breve{\mathbf{x}}=\left(\begin{array}{c}\breve{x}\\ \breve{y}\\ \breve{\sigma }\end{array}\right)=\mathbf{c}+\mathbf{d}=\mathbf{c}\underset{\mathbf{d}=\breve{\mathbf{x}}-\mathbf{c}}{\underbrace{-{\mathbf{H}}_D^{-1}(\mathbf{c})\cdot {\nabla }_D(\mathbf{c})}}\text{\hspace{1em}(25.60)}$$
处具有极值（最大值或最小值），假设Hessian矩阵${\mathbf{H}}_D$的逆存在。通过将极值位置$\breve{\mathbf{x}}$代入公式(25.59)，可以找到连续近似函数$\tilde{D}$的峰值（最小值或最大值）
$$\begin{gathered} D_{\text{peak}}(\mathbf{c})={\tilde{D}}_c(\breve{\mathbf{x}})=D(\mathbf{c})+\frac{1}{2}{\nabla }_D^{\top }(\mathbf{c})\cdot (\breve{\mathbf{x}}-\mathbf{c}) \\ =D(\mathbf{c})+\frac{1}{2}{\nabla }_D^{\top }(\mathbf{c})\cdot \mathbf{d},\text{\hspace{1em}(25.61)} \end{gathered}$$
其中$\mathbf{d}=\breve{\mathbf{x}}-\mathbf{c}$（参见公式(25.60)）表示邻域的离散中心位置$\mathbf{c}$和连续极值位置$\breve{\mathbf{x}}$之间的3D向量。

只有当DoG的估计幅度超过给定阈值$t_{\text{peak}}$时，尺度空间位置$\mathbf{c}$才保留为候选兴趣点，即
$$|D_{\text{peak}}(\mathbf{c})|>t_{\text{peak}}.\text{\hspace{1em}(25.62)}$$
如果公式(25.60)中的$\mathbf{c}$到估计的（连续）峰值位置$\hat{x}$的距离$\mathbf{d}=(x^{\prime },y^{\prime },{\sigma }^{\prime }{)}^{\top }$在任何空间方向上大于预定义的限制（通常为0.5），则中心点$\mathbf{c}=(u,v,q{)}^{\top }$将通过最大±1单位步长沿$u,v$轴移动到相邻的一个DoG单元，即
$$\mathbf{c}\leftarrow \mathbf{c}+\left(\begin{array}{c}\min (1,\max (-1,\text{round}(x^{\prime })))\\ \min (1,\max (-1,\text{round}(y^{\prime })))\\ 0\end{array}\right).\text{\hspace{1em}(25.63)}$$
$\mathbf{c}$的$q$分量在此版本中不作修改，即搜索继续在原来的尺度层次进行。基于新点周围的3D邻域，再次执行泰勒展开（公式(25.60)）以估计新的峰值位置。这一过程重复进行，直到峰值位置位于当前DoG单元内部或达到允许的重新定位步数$n_{\text{refine}}$（通常设置为4或5）。如果成功，这一步的结果是候选特征点
$$\breve{\mathbf{c}}=(\breve{x},\breve{y},\breve{q}{)}^{\top }=\mathbf{c}+(x^{\prime },y^{\prime },0{)}^{\top }.\text{\hspace{1em}(25.64)}$$
请注意（在此实现中）即使3D泰勒展开表明估计的峰值位于另一个尺度层次，尺度层次$q$仍保持不变。参见算法25.4中的过程RefineKeyPosition()以获取这些步骤的简明总结。

需要提到的是，原始文献[153]对上述位置优化过程没有特别详细的说明，因此各种开源SIFT实现中使用的方式略有不同。例如，只要$|x^{\prime }|$或$|y^{\prime }|$大于0.6，VLFeat[241]的实现就移动到同一尺度层次的直接邻居，如前所述。S. Nowozin的AutoPano-SIFT[18]计算空间位移$d=||(x^{\prime },y^{\prime })||$的长度，如果$d>2$则丢弃当前点。否则，它以${\Delta }_u=\text{round}(x^{\prime })$，${\Delta }_v=\text{round}(y^{\prime })$移动，而不限制位移到±1。OpenCV使用的开源SIFT库[106]也在空间方向上进行全移动，并且在每次迭代中还可能通过${\Delta }_q=\text{round}({\sigma }^{\prime })$改变尺度层次。

2.3 抑制对类似边缘结构的响应

在前一步中，候选兴趣点被选择为DoG尺度空间中泰勒近似具有局部最大值且外推的DoG值高于给定阈值（$t_{peak}$）的位置。然而，DoG滤波器也会对类似边缘的结构产生强烈响应。在这些位置上，兴趣点无法以足够的稳定性和重复性进行定位。为了消除边缘附近的响应，Lowe建议使用2D DoG结果沿空间$x$、$y$轴的主曲率，利用函数的主曲率与函数在给定点的Hessian矩阵的特征值成比例这一事实。

对于DoG尺度空间中的特定点$\mathbf{c}=(u,v,q)$，其邻域为${\mathsf{N}}_D$（见公式(25.52)），其空间坐标的$2\times 2$Hessian矩阵为

$${\mathbf{H}}_{xy}(\mathbf{c})=\left(\begin{array}{cc}{d}_{xx}& d_{xy}\\ d_{xy}& d_{yy}\end{array}\right),\text{\hspace{1em}(25.65)}$$

其中$d_{xx}$、$d_{xy}$、$d_{yy}$的定义见公式(25.58)，即这些值可以从相应的$3\times 3$Hessian矩阵${\mathbf{H}}_D(\mathbf{c})$（见公式(25.57)）中提取。

矩阵${\mathbf{H}}_{xy}(\mathbf{c})$有两个特征值${\lambda }_1,{\lambda }_2$，我们定义为按大小排序，使得${\lambda }_1$具有更大的绝对值（$|{\lambda }_1|\ge |{\lambda }_2|$）。如果某点$\mathbf{c}$的两个特征值的大小相似，函数沿两个正交方向具有高曲率，在这种情况下，$\mathbf{c}$很可能是一个可以可靠定位的良好参考点。在最佳情况下（例如在角点附近），特征值比率$\rho ={\lambda }_1/{\lambda }_2$接近1。相反，如果比率$\rho$很高，则可以得出结论，在此位置上单一方向占主导地位，这通常发生在边缘附近。

要估计比率$\rho$，不需要实际计算特征值本身。根据[153]中的描述，特征值${\lambda }_1,{\lambda }_2$的和与积可以表示为

$${\lambda }_1+{\lambda }_2=\mathrm{trace}({\mathbf{H}}_{xy}(\mathbf{c}))=d_{xx}+d_{yy},\text{\hspace{1em}(25.66)}$$

$${\lambda }_1\cdot {\lambda }_2=\det ({\mathbf{H}}_{xy}(\mathbf{c}))=d_{xx}\cdot d_{yy}-d_{xy}^2.\text{\hspace{1em}(25.67)}$$

如果行列式$\det ({\mathbf{H}}_{xy})$为负，则基础2D函数的主曲率符号相反，因此可以将点$\mathbf{c}$舍弃，因为它不是极值点。否则，如果两个特征值${\lambda }_1,{\lambda }_2$的符号相同，则比率

$${\rho }_{1,2}=\frac{{\lambda }_1}{{\lambda }_2}\text{\hspace{1em}(25.68)}$$

为正（其中${\lambda }_1={\rho }_{1,2}\cdot {\lambda }_2$），因此表达式

$$a=\frac{[\mathrm{trace}({\mathbf{H}}_{xy}(\mathbf{c})){]}^2}{\det ({\mathbf{H}}_{xy}(\mathbf{c}))}=\frac{({\lambda }_1+{\lambda }_2{)}^2}{{\lambda }_1\cdot {\lambda }_2}\text{\hspace{1em}(25.69)}$$

$$=\frac{({\rho }_{1,2}\cdot {\lambda }_2+{\lambda }_2{)}^2}{{\rho }_{1,2}\cdot {\lambda }_2^2}=\frac{{\lambda }_2^2\cdot ({\rho }_{1,2}+1{)}^2}{{\rho }_{1,2}\cdot {\lambda }_2^2}=\frac{({\rho }_{1,2}+1{)}^2}{{\rho }_{1,2}}\text{\hspace{1em}(25.70)}$$

仅依赖于比率${\rho }_{1,2}$。如果Hessian矩阵${\mathbf{H}}_{xy}$的行列式为正，量$a$在${\rho }_{1,2}=1$时达到最小值（4.0），即两个特征值相等（见图25.16）。注意比率$a$对于${\rho }_{1,2}={\lambda }_1/{\lambda }_2$或${\rho }_{1,2}={\lambda }_2/{\lambda }_1$是相同的，因为

$$a=\frac{({\rho }_{1,2}+1{)}^2}{{\rho }_{1,2}}=\frac{{\left(\frac{1}{{\rho }_{1,2}}+1\right)}^2}{\frac{1}{{\rho }_{1,2}}}.\text{\hspace{1em}(25.71)}$$

要验证给定位置$\mathbf{c}$处的特征值比率${\rho }_{1,2}$是否低于指定限值${\rho }_{\max }$（使$\mathbf{c}$成为良好候选点），只需检查条件

$$a\le a_{\max },\text{其中}{a}_{\max }=\frac{({\rho }_{\max }+1{)}^2}{{\rho }_{\max }},\text{\hspace{1em}(25.72)}$$

而无需实际计算个别特征值${\lambda }_1$和${\lambda }_2$。

${\rho }_{\max }$应大于1，通常选择在3到10的范围内（[153]中建议${\rho }_{\max }=10$）。在公式(25.72)中的结果值$a_{\max }$是常数，只需计算一次（见算法25.3，第2行）。对于不同${\rho }_{\max }$值的检测示例如图25.17所示。注意，随着${\rho }_{\max }$从3提高到40，边缘附近出现的候选点明显增多。

图 25.16
通过指定 $a_{\max }$ 限制主曲率比率（边缘比率）${\rho }_{1,2}$。当特征值比率 ${\rho }_{1,2}=\frac{{\lambda }_1}{{\lambda }_2}$ 为1时，即当两个特征值 ${\lambda }_1,{\lambda }_2$ 相等时，量$a$（蓝线）达到最小值，这表示一个类似角点的事件。通常，在图像线附近，只有一个特征值占主导地位，因此 ${\rho }_{1,2}$ 和 $a$ 的值显著增加。在此示例中，通过设置 $a_{\max }=\frac{(5+1{)}^2}{5}=7.2$（红线），将主曲率比率 ${\rho }_{1,2}$ 限制为 ${\rho }_{\max }=5.0$。

图 25.17
通过控制最大曲率比率 ${\rho }_{\max }$ 来拒绝类似边缘的特征。圆的大小与检测到相应关键点的尺度水平成比例，颜色表示包含的八度（0 = 红色，1 = 绿色，2 = 蓝色，3 = 洋红色）。