文章目录
【SIFT介绍】Scale-Invariant Feature Transform——尺度不变特征变换(三)
简介
许多实际应用需要在一张或多张图像中定位参考位置,例如图像对齐、去除畸变、物体跟踪、3D重建等。我们已经看到,角点可以相当可靠地定位,并且不依赖于方向。然而,典型的角点检测器仅提供每个候选点的位置和强度,它们不提供任何有关其特征或“身份”的信息,这些信息可以用于匹配。另一个限制是大多数角点检测器仅在特定尺度或分辨率下工作,因为它们基于一组固定的滤波器。
本章介绍了局部特征检测的尺度不变特征变换(SIFT)技术,该技术最初由D. Lowe [152]提出,并自此成为成像行业的“主力”方法。其目标是定位能够鲁棒识别的图像特征,以便在多幅图像和图像序列中进行匹配,以及在不同视角条件下进行物体识别。SIFT采用了“尺度空间”的概念[151],在多个尺度层次或图像分辨率下捕捉特征,这不仅增加了可用特征的数量,还使得该方法对尺度变化具有高度的容忍度。这使得在例如物体向相机移动并因此连续改变其尺度的情况下跟踪特征成为可能,或将使用不同变焦设置拍摄的图像拼接在一起。
通过简化尺度空间计算和特征检测或使用GPU硬件[20, 90, 218],已经实现了SIFT算法的加速变种。
原则上,SIFT的工作方式类似于多尺度角点检测器,具有亚像素定位精度,并为每个候选点附加了旋转不变的特征描述符。这个(通常是128维的)特征描述符总结了对应特征点周围空间邻域中的梯度方向分布,因此可以像“指纹”一样使用。SIFT特征计算涉及的主要步骤如下:
- 在拉普拉斯-高斯(LoG)尺度空间中检测极值点,以定位潜在的兴趣点。
- 通过拟合连续模型来精确确定位置和尺度,从而对关键点进行优化。
- 通过周围图像梯度方向的主导方向为特征点分配方向。
- 通过归一化局部梯度直方图来形成特征描述符。
这些步骤都将在本章的其余部分详细描述。
我们在这里如此详细地解释SIFT技术有几个原因。首先,这是迄今为止我们所讨论的最复杂的算法,其各个步骤精心设计并相互依赖,涉及众多需要考虑的参数。因此,深入理解其内部工作原理和局限性对于成功使用以及在结果不如预期时分析问题非常重要。
3. 创建局部描述符
对于在分层DoG尺度空间中检测到的每个局部最大值,都会创建一个候选关键点,并按照前述步骤将其细化到连续位置(见公式(25.56)-(25.64))。然后,对于每个细化后的关键点 k = ( p , q , x , y ) k = (p, q, x, y) k=(p,q,x,y),计算一个或多个(最多四个)局部描述符。如果局部方向不唯一,则可能为一个位置创建多个(最多四个)描述符。此过程涉及以下步骤:
- 从相应的高斯尺度空间级别的梯度分布中找到关键点 k k k的主要方向。
- 对于每个主要方向,在关键点 k k k处创建一个单独的SIFT描述符。
3.1 寻找主要方向
从高斯尺度空间获取局部方向
通过对分层高斯尺度空间 G p , q ( u , v ) \mathbf G_{p,q}(u, v) Gp,q(u,v)(见公式(25.32))的梯度值进行采样,得到方向向量。对于八度 p p p和尺度级别 q q q下任意晶格位置 ( u , v ) (u, v) (u,v),局部梯度计算如下:
∇ p , q ( u , v ) = ( d x d y ) = 0.5 ⋅ ( G p , q ( u + 1 , v ) − G p , q ( u − 1 , v ) G p , q ( u , v + 1 ) − G p , q ( u , v − 1 ) ) (25.73) \mathbf{\nabla}_{p,q}(u, v) = \begin{pmatrix} d_x \\ d_y \end{pmatrix} = 0.5 \cdot \begin{pmatrix} \mathbf G_{p,q}(u+1, v) - \mathbf G_{p,q}(u-1, v) \\ \mathbf G_{p,q}(u, v+1) - \mathbf G_{p,q}(u, v-1) \end{pmatrix} \tag{25.73} ∇p,q(u,v)=(dxdy)=0.5⋅(Gp,q(u+1,v)−Gp,q(u−1,v)Gp,q(u,v+1)−Gp,q(u,v−1))(25.73)
从这些梯度向量中,得到梯度幅度和方向(即极坐标):
E p , q ( u , v ) = ∥ ∇ p , q ( u , v ) ∥ = d x 2 + d y 2 , (25.74) E_{p,q}(u, v) = \left\|\mathbf{\nabla}_{p,q}(u, v)\right\| = \sqrt{d_x^2 + d_y^2}, \tag{25.74} Ep,q(u,v)=∥∇p,q(u,v)∥=dx2+dy2,(25.74)
ϕ p , q ( u , v ) = ∠ ∇ p , q ( u , v ) = tan − 1 ( d y d x ) . (25.75) \phi_{p,q}(u, v) = \angle \mathbf{\nabla}_{p,q}(u, v) = \tan^{-1}\left(\frac{d_y}{d_x}\right). \tag{25.75} ϕp,q(u,v)=∠∇p,q(u,v)=tan−1(dxdy).(25.75)
这些标量场 E p , q E_{p,q} Ep,q和 ϕ p , q \phi_{p,q} ϕp,q通常预先计算,用于高斯尺度空间 G \mathbf G G的所有相关八度 p p p和级别 q q q。
方向直方图
为了找到给定关键点的主要方向,计算梯度向量从关键点中心周围的正方形窗口中收集的方向角直方图 h ϕ h_{\phi} hϕ。通常,直方图有 n orient = 36 n_{\text{orient}} = 36 norient=36个箱,即角度分辨率为10°。方向直方图从使用各向同性高斯加权函数的正方形区域中收集,其宽度 σ w \sigma_w σw与关键点的尺度级别 q q q的减小尺度 σ ˙ q \dot{\sigma}_q σ˙q成比例(见公式(25.37))。通常使用“ σ \sigma σ是关键点尺度的1.5倍”的高斯加权函数[153],即:
σ w = 1.5 ⋅ σ ˙ q = 1.5 ⋅ σ 0 ⋅ 2 q / Q . (25.76) \sigma_w = 1.5 \cdot \dot{\sigma}_q = 1.5 \cdot \sigma_0 \cdot 2^{q/Q}. \tag{25.76} σw=1.5⋅σ˙q=1.5⋅σ0⋅2q/Q.(25.76)
注意, σ w \sigma_w σw与八度指数 p p p无关,因此每个八度使用相同的加权函数。为了计算方向直方图,从大小为 2 r w × 2 r w 2r_w \times 2r_w 2rw×2rw的正方形区域收集关键点周围的高斯梯度,其中
r w = ⌈ 2.5 ⋅ σ w ⌉ (25.77) r_w = \left\lceil 2.5 \cdot \sigma_w \right\rceil \tag{25.77} rw=⌈2.5⋅σw⌉(25.77)
充分维度以避免数值截断效应。对于表25.3中列出的参数( σ 0 = 1.6 , Q = 3 \sigma_0 = 1.6, Q = 3 σ0=1.6,Q=3),高斯加权函数的 σ w \sigma_w σw值(以八度的坐标单位表示)如下:
q 0 1 2 3 σ w 1.6000 2.0159 2.5398 3.2000 r w 4 5 6 7 (25.78) \begin{array}{cccc} q & 0 & 1 & 2 & 3 \\ \sigma_w & 1.6000 & 2.0159 & 2.5398 & 3.2000 \\ r_w & 4 & 5 & 6 & 7 \\ \end{array} \tag{25.78} qσwrw01.6000412.0159522.5398633.20007(25.78)
在算法25.7中,高斯加权函数的 σ w \sigma_w σw和 r w r_w rw分别在第7行和第8行计算。在八度 p p p和级别 q q q的高斯尺度空间 G \mathbf G G中,每个晶格点 ( u , v ) (u, v) (u,v)的梯度向量 ∇ p , q ( u , v ) \mathbf{\nabla}_{p,q}(u, v) ∇p,q(u,v)在第16行计算。从中得到梯度幅度 E p , q ( u , v ) E_{p,q}(u, v) Ep,q(u,v)和方向 ϕ p , q ( u , v ) \phi_{p,q}(u, v) ϕp,q(u,v)(第29-30行)。相应的高斯权重从网格点 ( u , v ) (u, v) (u,v)与兴趣点 ( x , y ) (x, y) (x,y)之间的空间距离计算(第18行):
w G ( u , v ) = exp ( − ( u − x ) 2 + ( v − y ) 2 2 σ w 2 ) . (25.79) w_G(u, v) = \exp \left( -\frac{(u-x)^2+(v-y)^2}{2\sigma_w^2} \right). \tag{25.79} wG(u,v)=exp(−2σw2(u−x)2+(v−y)2).(25.79)
对于网格点 ( u , v ) (u, v) (u,v),累积到方向直方图中的量为
z = E p , q ( u , v ) ⋅ w G ( u , v ) , (25.80) z = E_{p,q}(u, v) \cdot w_G(u, v), \tag{25.80} z=Ep,q(u,v)⋅wG(u,v),(25.80)
即局部梯度幅度乘以高斯窗口函数(算法25.7第19行)。
方向直方图 h ϕ h_{\phi} hϕ由 n orient n_{\text{orient}} norient个箱组成,因此角度 ϕ ( u , v ) \phi(u, v) ϕ(u,v)的连续箱号为
κ ϕ = n orient 2 π ⋅ ϕ ( u , v ) (25.81) \kappa_{\phi} = \frac{n_{\text{orient}}}{2\pi} \cdot \phi(u, v) \tag{25.81} κϕ=2πnorient⋅ϕ(u,v)(25.81)
(见算法25.7第20行)。为了将连续方向收集到具有离散箱的直方图中,需要进行量化。最简单的方法是选择“最近”箱(通过四舍五入),并将相关量(记作 z z z)全部添加到选定的箱中。另一种常用技术是将量 z z z分配到两个最接近的箱中。给定连续的箱值 κ ϕ \kappa_{\phi} κϕ,两个最接近的离散箱的索引分别为
k 0 = ⌊ κ ϕ ⌋ m o d n orient , k 1 = ( ⌊ κ ϕ ⌋ + 1 ) m o d n orient , (25.82) k_0 = \lfloor \kappa_{\phi} \rfloor \mod n_{\text{orient}}, \quad k_1 = (\lfloor \kappa_{\phi} \rfloor + 1) \mod n_{\text{orient}}, \tag{25.82} k0=⌊κϕ⌋modnorient,k1=(⌊κϕ⌋+1)modnorient,(25.82)
然后量 z z z(公式25.80)被分配并累积到方向直方图 h ϕ h_{\phi} hϕ的相邻箱 k 0 k_0 k0和 k 1 k_1 k1中,形式为
h ϕ ( k 0 ) ← h ϕ ( k 0 ) + ( 1 − α ) ⋅ z , h_{\phi}(k_0) \leftarrow h_{\phi}(k_0) + (1 - \alpha) \cdot z, hϕ(k0)←hϕ(k0)+(1−α)⋅z,
h ϕ ( k 1 ) ← h ϕ ( k 1 ) + α ⋅ z , (25.83) h_{\phi}(k_1) \leftarrow h_{\phi}(k_1) + \alpha \cdot z, \tag{25.83} hϕ(k1)←hϕ(k1)+α⋅z,(25.83)
其中 α = κ ϕ − ⌊ κ ϕ ⌋ \alpha = \kappa_{\phi} - \lfloor \kappa_{\phi} \rfloor α=κϕ−⌊κϕ⌋。这个过程在图25.18中有示例说明(见算法25.7第21-25行)。
图 25.18
通过线性插值累积到多个直方图箱中。假设某个量 z z z(蓝色条)要添加到离散直方图 h ϕ h_{\phi} hϕ 的连续位置 κ ϕ \kappa_{\phi} κϕ。 κ ϕ \kappa_{\phi} κϕ 相邻的直方图箱是 k 0 = ⌊ κ ϕ ⌋ k_0 = \lfloor \kappa_{\phi} \rfloor k0=⌊κϕ⌋ 和 k 1 = ⌊ κ ϕ ⌋ + 1 k_1 = \lfloor \kappa_{\phi} \rfloor + 1 k1=⌊κϕ⌋+1。累积到箱 k 1 k_1 k1 的部分 z 1 = z ⋅ α z_1 = z \cdot \alpha z1=z⋅α(红色条),其中 α = κ ϕ − k 0 \alpha = \kappa_{\phi} - k0 α=κϕ−k0。同样地,添加到箱 k 0 k_0 k0 的量是 z 0 = z ⋅ ( 1 − α ) z_0 = z \cdot (1 - \alpha) z0=z⋅(1−α)(绿色条)。图 25.19
方向直方图示例。36 个径向条中的每一个对应于方向直方图 h ϕ h_{\phi} hϕ 中的一个条目。具有索引 k k k 的每个径向条的长度(半径)与相应箱 h ϕ ( k ) h_{\phi}(k) hϕ(k) 中累积的值成正比,其方向为 ϕ k \phi_k ϕk。图 25.20
通过反复应用具有 1D 核 H = 1 4 ⋅ ( 1 , 2 , 1 ) H = \frac{1}{4} \cdot (1, 2, 1) H=41⋅(1,2,1) 的圆形低通滤波器来平滑方向直方图(从图 25.19 中)。
平滑方向直方图
图 25.19 展示了方向直方图的几何渲染,解释了单元索引(离散角度 ϕ k \phi_k ϕk)和累积量( z z z)的相关性。在计算主要方向之前,通常通过应用(圆形)低通滤波器来平滑原始方向直方图 h ϕ h_{\phi} hϕ,通常使用简单的3抽头高斯或盒型滤波器(见算法 25.7 中的 SmoothCircular() 程序,第6-16行)。通过多次应用滤波器可以实现更强的平滑效果,如图 25.20 所示。实际操作中,两到三次平滑迭代就足够了。
定位和插值方向峰值
在平滑方向直方图之后,下一步是检测 h ϕ h_{\phi} hϕ 中的峰值条目。如果直方图箱 k k k 是局部最大值且其值不低于最大直方图条目的一定比例,则认为它是显著的方向峰值,即,仅当满足以下条件时:
h ϕ ( k ) > h ϕ ( ( k − 1 ) m o d n orient ) ∧ h ϕ ( k ) > h ϕ ( ( k + 1 ) m o d n orient ) ∧ h ϕ ( k ) > t domor ⋅ max i h ϕ ( i ) , (25.84) h_{\phi}(k) > h_{\phi}((k - 1) \mod n_{\text{orient}}) \wedge \\ h_{\phi}(k) > h_{\phi}((k + 1) \mod n_{\text{orient}}) \wedge \\ h_{\phi}(k) > t_{\text{domor}} \cdot \max_i h_{\phi}(i), \tag{25.84} hϕ(k)>hϕ((k−1)modnorient)∧hϕ(k)>hϕ((k+1)modnorient)∧hϕ(k)>tdomor⋅imaxhϕ(i),(25.84)
其中 t domor = 0.8 t_{\text{domor}} = 0.8 tdomor=0.8 是一个典型的限值。
为了获得比方向直方图箱(通常间隔10°)更高的角度分辨率,通过对相邻直方图值进行二次插值来计算连续的峰值方向。给定离散峰值索引 k k k,通过拟合二次函数到连续三个直方图值 h ϕ ( k − 1 ) , h ϕ ( k ) , h ϕ ( k + 1 ) h_{\phi}(k-1), h_{\phi}(k), h_{\phi}(k+1) hϕ(k−1),hϕ(k),hϕ(k+1) 得到插值(连续)峰值位置 k ˘ \breve{k} k˘:
k ˘ = k + h ϕ ( k − 1 ) − h ϕ ( k + 1 ) 2 ⋅ ( h ϕ ( k − 1 ) − 2 h ϕ ( k ) + h ϕ ( k + 1 ) ) , (25.85) \breve{k} = k + \frac{h_{\phi}(k-1) - h_{\phi}(k+1)}{2 \cdot \left( h_{\phi}(k-1) - 2h_{\phi}(k) + h_{\phi}(k+1) \right)}, \tag{25.85} k˘=k+2⋅(hϕ(k−1)−2hϕ(k)+hϕ(k+1))hϕ(k−1)−hϕ(k+1),(25.85)
所有索引都取模 n orient n_{\text{orient}} norient。从公式 (25.81) 得到的(连续)主要方向角 θ ∈ [ 0 , 2 π ) \theta \in [0, 2\pi) θ∈[0,2π) 为:
θ = ( k ˘ m o d n orient ) ⋅ 2 π n orient , (25.86) \theta = (\breve{k} \mod n_{\text{orient}}) \cdot \frac{2\pi}{n_{\text{orient}}}, \tag{25.86} θ=(k˘modnorient)⋅norient2π,(25.86)
其中 θ ∈ [ 0 , 2 π ) \theta \in [0, 2\pi) θ∈[0,2π)。这样可以以远超方向直方图粗分辨率的精度估计主要方向。注意,在某些情况下,一个给定的关键点可能会获得多个直方图峰值(见算法 25.6 中的 FindPeakOrientations() 程序,第18-31行)。在这种情况下,对于每个主要方向,在相同的关键点位置创建单独的 SIFT 描述符(见算法 25.3,第8行)。
图 25.21 显示了在对不同图像中的一组检测到的关键点应用不同数量的平滑步骤后,方向直方图以及从方向直方图(公式 (25.86))计算出的插值主要方向
θ
\theta
θ 通过相应的向量表示的结果。
图 25.21
方向直方图和主要方向(示例)。对方向直方图进行了 n = 0 , … , 3 n = 0,\ldots, 3 n=0,…,3 次平滑迭代。插值后的主要方向显示为从每个特征中心点发出的径向线。直方图图表的大小与检测到相应关键点的绝对尺度( σ p , q \sigma_{p,q} σp,q,见表25.3)成比例。颜色表示包含尺度空间八度 p p p 的索引(红色 = 0,绿色 = 1,蓝色 = 2,品红色 = 3)。
3.2 SIFT 描述符构建
对于每个关键点 k ′ = ( p , q , x , y ) \mathbf k' = (p, q, x, y) k′=(p,q,x,y) 和每个主导方向 θ \theta θ,通过在高斯尺度空间 G \mathbf G G 的八度 p p p 和层次 q q q 采样周围的梯度,获得对应的 SIFT 描述符。
描述符几何
描述符计算的几何结构如图 25.22 所示。描述符结合了大小为 w d × w d w_d \times w_d wd×wd 的正方区域内的梯度方向和幅度,该区域以相关特征点的(连续)位置 ( x , y ) (x, y) (x,y) 为中心,并与其主导方向 θ \theta θ 对齐。描述符的边长设定为 w d = 10 ⋅ σ ˙ q w_d = 10 \cdot \dot{\sigma}_q wd=10⋅σ˙q,其中 σ ˙ q \dot{\sigma}_q σ˙q 表示关键点的缩减尺度(内圈的半径)。它取决于关键点的尺度层次 q q q(见表 25.4)。
该区域被划分为 n s p a t × n s p a t n_{spat} \times n_{spat} nspat×nspat 个大小相同的子正方形;通常 n s p a t = 4 n_{spat} = 4 nspat=4(见表 25.5)。每个梯度样本的贡献通过一个宽度为 σ d = 0.25 ⋅ w d \sigma_d = 0.25 \cdot w_d σd=0.25⋅wd 的圆形高斯函数(蓝色圆)衰减。权重径向下降,在 r d = 2.5 ⋅ σ d r_d = 2.5 \cdot \sigma_d rd=2.5⋅σd 处几乎为零(图 25.22 中的绿色圆)。因此,仅需包括该区域外的样本来计算描述符统计数据。
为了实现旋转不变性,描述符区域与在前面步骤中确定的关键点主导方向对齐。为了使描述符对尺度变化不敏感,其大小
w
d
w_d
wd(以八度
p
p
p 的网格坐标单位表示)设置为与关键点的缩减尺度
σ
˙
q
\dot{\sigma}_q
σ˙q 成比例(见公式 (25.37)),即:
w
d
=
s
d
⋅
σ
˙
q
=
s
d
⋅
σ
0
⋅
2
q
/
Q
,
(25.87)
w_d = s_d \cdot \dot{\sigma}_q = s_d \cdot \sigma_0 \cdot 2^{q/Q}, \tag{25.87}
wd=sd⋅σ˙q=sd⋅σ0⋅2q/Q,(25.87)
其中
s
d
s_d
sd 是一个常数大小因子。对于
s
d
=
10
s_d = 10
sd=10(见表 25.5),描述符大小
w
d
w_d
wd 从 16.0(在层次 0)到 25.4(在层次 2)不等,如表 25.4 所示。注意,描述符大小
w
d
w_d
wd 仅取决于尺度层次索引
q
q
q,与八度索引
p
p
p 无关。因此,同样的描述符几何适用于尺度空间的所有八度。
描述符的空间分辨率由参数
n
s
p
a
t
n_{spat}
nspat 指定。通常
n
s
p
a
t
=
4
n_{spat} = 4
nspat=4(如图 25.22 所示),因此总空间桶数为
n
s
p
a
t
×
n
s
p
a
t
=
16
n_{spat} \times n_{spat} = 16
nspat×nspat=16(在这种情况下)。每个空间描述符桶对应一个大小为
(
w
d
/
n
s
p
a
t
)
×
(
w
d
/
n
s
p
a
t
)
(w_d/n_{spat}) \times (w_d/n_{spat})
(wd/nspat)×(wd/nspat) 的区域。例如,在任何八度的尺度层次
q
=
0
q = 0
q=0,
σ
˙
0
=
1.6
\dot{\sigma}_0 = 1.6
σ˙0=1.6,相应的描述符大小为
w
d
=
s
d
⋅
σ
˙
0
=
10
⋅
1.6
=
16.0
w_d = s_d \cdot \dot{\sigma}_0 = 10 \cdot 1.6 = 16.0
wd=sd⋅σ˙0=10⋅1.6=16.0(见表 25.4)。在这种情况下(如图 25.23 所示),描述符覆盖 16 × 16 个梯度样本,如 [153] 所建议。图 25.24 显示了一个示例,其中 M 形特征点标记对齐到主导方向并缩放到相关尺度层次的描述符区域宽度
w
d
w_d
wd。
图 25.22
SIFT 描述符的几何结构。描述符从一个以关键点位置 ( x , y ) (x, y) (x,y) 为中心、与关键点的主导方向 θ \theta θ 对齐的正方支持区域计算而来,并被划分为 n spat × n spat n_{\text{spat}} \times n_{\text{spat}} nspat×nspat(4 × 4)个子正方形。内圈(灰色)的半径对应特征点的缩减尺度值 σ q \sigma_q σq。蓝色圆表示应用于梯度的高斯权重函数的宽度 σ d \sigma_d σd;在绿色圆( r d r_d rd)外,其值实际上为零。
表 25.4
不同尺度级别 q q q 下的 SIFT 描述符尺寸(对于大小因子 s d = 10 s_d = 10 sd=10 和每八度 Q = 3 Q = 3 Q=3 个级别)。 σ ˙ q \dot{\sigma}_q σ˙q 是关键点的缩减尺度, w d w_d wd 是描述符的大小, σ d \sigma_d σd 是高斯权重函数的宽度, r d r_d rd 是描述符支持区域的半径。对于 Q = 3 Q = 3 Q=3,仅与尺度级别 q = 0 , 1 , 2 q = 0, 1, 2 q=0,1,2 相关。所有长度都以八度的(即缩减的)坐标单位表示。
图 25.23
SIFT 描述符与相关八度(级别 q = 0 q = 0 q=0,参数 s d = 10 s_d = 10 sd=10)的离散样本网格的几何关系。在这种情况下,缩减尺度是 σ ˙ 0 = 1.6 \dot{\sigma}_0 = 1.6 σ˙0=1.6,描述符的宽度是 w d = s d ⋅ σ ˙ 0 = 10 ⋅ 1.6 = 16.0 wd = s_d \cdot \dot{\sigma}_0 = 10 \cdot 1.6 = 16.0 wd=sd⋅σ˙0=10⋅1.6=16.0。图 25.24
与主导方向对齐的标记关键点。注意,在具有多个主导方向的关键点位置插入了多个特征实例。标记的大小与检测到相应关键点的绝对尺度( σ p , q \sigma_{p,q} σp,q,见表 25.3)成正比。颜色表示包含八度 p p p 的尺度空间的索引(红色 = 0,绿色 = 1,蓝色 = 2,品红色 = 3)。
梯度特征
实际的SIFT描述符是通过对高斯尺度级别内描述符空间支持区域的梯度方向进行直方图统计得到的特征向量。这需要一个3D直方图 h ∇ ( i , j , k ) h_\nabla(i, j, k) h∇(i,j,k),其中两个空间维度 ( i , j ) (i, j) (i,j) 对应于 n spat × n spat n_{\text{spat}} \times n_{\text{spat}} nspat×nspat 子区域,另外一个维度 ( k ) (k) (k) 对应于 n angl n_{\text{angl}} nangl 个梯度方向。因此,该直方图包含 n spat × n spat × n angl n_{\text{spat}} \times n_{\text{spat}} \times n_{\text{angl}} nspat×nspat×nangl 个箱。图25.25展示了这一结构的典型设置, n spat = 4 n_{\text{spat}} = 4 nspat=4 和 n angl = 8 n_{\text{angl}} = 8 nangl=8(见表25.5)。在这种安排中,八个方向箱 k = 0 , … , 7 k = 0, \ldots, 7 k=0,…,7 附加到16个空间位置箱 (A1, … , D4) 上,总共有128个直方图箱。
对于给定的关键点 k ′ = ( p , q , x , y ) k' = (p, q, x, y) k′=(p,q,x,y),直方图 h ∇ h_\nabla h∇ 积累了高斯尺度空间级别 G p , q \mathbf G_{p,q} Gp,q 在支持区域中心坐标 ( x , y ) (x, y) (x,y) 周围的梯度方向(角度)。在该区域内的每个网格点 ( u , v ) (u, v) (u,v),估计梯度向量 ∇ G \nabla _G ∇G(如方程25.73所述),并计算梯度幅值 E ( u , v ) E(u, v) E(u,v) 和方向 φ ( u , v ) \varphi(u, v) φ(u,v)(见方程25.74-25.75和算法25.7的27-31行)。出于效率考虑,通常会为所有相关的尺度级别预先计算 E ( u , v ) E(u, v) E(u,v) 和 φ ( u , v ) \varphi(u, v) φ(u,v)。
每个梯度样本根据梯度幅值
E
E
E 和样本点
(
u
,
v
)
(u, v)
(u,v) 与关键点中心
(
x
,
y
)
(x, y)
(x,y) 的距离对梯度直方图
h
∇
h_\nabla
h∇ 作出特定的贡献量
z
z
z。同样使用高斯加权函数(宽度
σ
d
\sigma_d
σd)来衰减随空间距离增加的样本量;因此,累积的量为
z
(
u
,
v
)
=
R
(
u
,
v
)
⋅
w
G
=
R
(
u
,
v
)
⋅
exp
(
−
(
u
−
x
)
2
+
(
v
−
y
)
2
2
σ
d
2
)
.
(25.88)
z(u, v) = R(u, v) \cdot w_G = R(u, v) \cdot \exp \left(-\frac{(u-x)^2+(v-y)^2}{2\sigma^2_d}\right). \tag{25.88}
z(u,v)=R(u,v)⋅wG=R(u,v)⋅exp(−2σd2(u−x)2+(v−y)2).(25.88)
高斯函数
w
G
(
)
w_G()
wG() 的宽度
σ
d
\sigma_d
σd 与描述符区域的边长成正比,
σ
d
=
0.25
⋅
w
d
=
0.25
⋅
s
d
⋅
σ
˙
q
.
(25.89)
\sigma_d = 0.25 \cdot w_d = 0.25 \cdot s_d \cdot \dot{\sigma}_q. \tag{25.89}
σd=0.25⋅wd=0.25⋅sd⋅σ˙q.(25.89)
加权函数从中心开始径向衰减,在距离
r
d
=
2.5
⋅
σ
d
r_d = 2.5 \cdot \sigma_d
rd=2.5⋅σd 处几乎为零。因此,在计算梯度直方图时,仅需考虑距离关键点中心较近的梯度样本(图25.22中的绿色圆圈,见算法25.8的第7和17行)。对于给定的关键点
k
′
=
(
p
,
q
,
x
,
y
)
k' = (p, q, x, y)
k′=(p,q,x,y),高斯梯度的采样可限制在由
x
±
r
d
x \pm r_d
x±rd 和
y
±
r
d
y \pm r_d
y±rd 界定的方形区域内的网格点
(
u
,
v
)
(u, v)
(u,v)(见算法25.8的第8-10行和15-16行)。然后将每个样本点
(
u
,
v
)
(u, v)
(u,v) 进行仿射变换
(
u
′
v
′
)
=
1
w
d
(
cos
(
−
θ
)
−
sin
(
−
θ
)
sin
(
−
θ
)
cos
(
−
θ
)
)
(
u
−
x
v
−
y
)
,
(25.90)
\begin{pmatrix} u' \\ v' \end{pmatrix} = \frac{1}{w_d} \begin{pmatrix} \cos(-\theta) & -\sin(-\theta) \\ \sin(-\theta) & \cos(-\theta) \end{pmatrix} \begin{pmatrix} u-x \\ v-y \end{pmatrix}, \tag{25.90}
(u′v′)=wd1(cos(−θ)sin(−θ)−sin(−θ)cos(−θ))(u−xv−y),(25.90)
该变换执行一个主方向
θ
\theta
θ 的旋转,并将原始(旋转后)的
w
d
×
w
d
w_d \times w_d
wd×wd 大小的方形映射到单位方形,坐标为
u
′
,
v
′
∈
[
−
0.5
,
+
0.5
]
u', v' \in [-0.5, +0.5]
u′,v′∈[−0.5,+0.5](见图25.23)。
为了使特征向量具有旋转不变性,各个梯度方向
ϕ
(
u
,
v
)
\phi(u, v)
ϕ(u,v) 旋转主方向,即
ϕ
′
(
u
,
v
)
=
(
ϕ
(
u
,
v
)
−
θ
)
m
o
d
2
π
,
(25.91)
\phi'(u, v) = (\phi(u, v) - \theta) \mod 2\pi, \tag{25.91}
ϕ′(u,v)=(ϕ(u,v)−θ)mod2π,(25.91)
其中
ϕ
′
(
u
,
v
)
∈
[
0
,
2
π
)
\phi'(u, v) \in [0, 2\pi)
ϕ′(u,v)∈[0,2π),以保持相对方向不变。对于每个梯度样本,具有连续坐标
(
u
′
,
v
′
,
ϕ
′
)
(u', v', \phi')
(u′,v′,ϕ′) 的相应量
z
(
u
,
v
)
z(u, v)
z(u,v)(方程25.88)被累积到3D梯度直方图
h
∇
h_\nabla
h∇ 中。有关该步骤的完整描述,请参见算法25.9中的过程UpdateGradientHistogram()。首先将坐标
(
u
′
,
v
′
,
ϕ
′
)
(u', v', \phi')
(u′,v′,ϕ′)(见方程25.90)映射到连续直方图位置
(
i
′
,
j
′
,
k
′
)
(i', j', k')
(i′,j′,k′)
i
′
=
n
spat
⋅
u
′
+
0.5
⋅
(
n
spat
−
1
)
,
j
′
=
n
spat
⋅
v
′
+
0.5
⋅
(
n
spat
−
1
)
,
k
′
=
ϕ
′
⋅
n
angl
2
π
,
(25.92)
i' = n_{\text{spat}} \cdot u' + 0.5 \cdot (n_{\text{spat}}-1), \\ j' = n_{\text{spat}} \cdot v' + 0.5 \cdot (n_{\text{spat}}-1), \\ k' = \phi' \cdot \frac{n_{\text{angl}}}{2\pi}, \tag{25.92}
i′=nspat⋅u′+0.5⋅(nspat−1),j′=nspat⋅v′+0.5⋅(nspat−1),k′=ϕ′⋅2πnangl,(25.92)
使得
i
′
,
j
′
∈
[
−
0.5
,
n
spat
−
0.5
]
i', j' \in [-0.5, n_{\text{spat}}-0.5]
i′,j′∈[−0.5,nspat−0.5] 和
k
′
∈
[
0
,
n
angl
)
k' \in [0, n_{\text{angl}})
k′∈[0,nangl)。
类似于通过线性插值在两个箱之间插入连续位置的1D直方图(见图25.18),量
z
z
z 通过三线性插值分布到八个相邻的直方图箱中。
z
z
z 贡献给各个直方图箱的分位数由坐标
(
i
′
,
j
′
,
k
′
)
(i', j', k')
(i′,j′,k′) 与受影响的直方图箱的离散索引
(
i
,
j
,
k
)
(i, j, k)
(i,j,k) 的距离决定。索引
(
i
,
j
,
k
)
(i, j, k)
(i,j,k) 是可能组合的集合
{
i
0
,
i
1
}
×
{
j
0
,
j
1
}
×
{
k
0
,
k
1
}
\{i_0, i_1\} \times \{j_0, j_1\} \times \{k_0, k_1\}
{i0,i1}×{j0,j1}×{k0,k1},
i
0
=
⌊
i
′
⌋
,
i
1
=
(
i
0
+
1
)
,
j
0
=
⌊
j
′
⌋
,
j
1
=
(
j
0
+
1
)
,
k
0
=
⌊
k
′
⌋
m
o
d
n
angl
,
k
1
=
(
k
0
+
1
)
m
o
d
n
angl
,
(25.93)
i_0 = \lfloor i' \rfloor,\qquad i_1 = (i_0 + 1), \\ j_0 = \lfloor j' \rfloor, \qquad j1 = (j_0 + 1), \tag{25.93} \\ k_0 = \lfloor k' \rfloor \mod n_{\text{angl}},\qquad k_1 = (k_0 + 1) \mod n_{\text{angl}},
i0=⌊i′⌋,i1=(i0+1),j0=⌊j′⌋,j1=(j0+1),k0=⌊k′⌋modnangl,k1=(k0+1)modnangl,(25.93)
相应的分位数(权重)为
α
0
=
⌊
i
′
⌋
+
1
−
i
′
=
i
1
−
i
′
,
α
1
=
1
−
α
0
,
β
0
=
⌊
j
′
⌋
+
1
−
j
′
=
j
1
−
j
′
,
β
1
=
1
−
β
0
,
γ
0
=
⌊
k
′
⌋
+
1
−
k
′
,
γ
1
=
1
−
γ
0
,
(25.94)
\alpha_0 = \lfloor i' \rfloor + 1 - i' = i_1 - i', \qquad \alpha_1 = 1 - \alpha_0, \\ \beta_0 = \lfloor j' \rfloor + 1 - j' = j_1 - j',\qquad \beta_1 = 1 - \beta_0, \tag{25.94} \\ \gamma_0 = \lfloor k' \rfloor + 1 - k', \qquad \gamma_1 = 1 - \gamma_0,
α0=⌊i′⌋+1−i′=i1−i′,α1=1−α0,β0=⌊j′⌋+1−j′=j1−j′,β1=1−β0,γ0=⌊k′⌋+1−k′,γ1=1−γ0,(25.94)
最后更新梯度直方图的(八个)受影响箱为
h
∇
(
i
0
,
j
0
,
k
0
)
←
+
z
⋅
α
0
⋅
β
0
⋅
γ
0
,
h
∇
(
i
1
,
j
0
,
k
0
)
←
+
z
⋅
α
1
⋅
β
0
⋅
γ
0
,
h
∇
(
i
0
,
j
1
,
k
0
)
←
+
z
⋅
α
0
⋅
β
1
⋅
γ
0
,
⋮
⋮
h
∇
(
i
1
,
j
1
,
k
1
)
←
+
z
⋅
α
1
⋅
β
1
⋅
γ
1
.
(25.95)
h_\nabla(i_0, j_0, k_0) \stackrel{+}{\gets} z \cdot \alpha_0 \cdot \beta_0 \cdot \gamma_0, \\ h_\nabla(i_1, j_0, k_0) \stackrel{+}{\gets} z \cdot \alpha_1 \cdot \beta_0 \cdot \gamma_0, \\ h_\nabla(i_0, j_1, k_0) \stackrel{+}{\gets} z \cdot \alpha_0 \cdot \beta_1 \cdot \gamma_0, \\ \vdots \qquad \qquad \qquad\quad \vdots \\ h_\nabla(i_1, j_1, k_1) \stackrel{+}{\gets} z \cdot \alpha_1 \cdot \beta_1 \cdot \gamma_1. \tag{25.95}
h∇(i0,j0,k0)←+z⋅α0⋅β0⋅γ0,h∇(i1,j0,k0)←+z⋅α1⋅β0⋅γ0,h∇(i0,j1,k0)←+z⋅α0⋅β1⋅γ0,⋮⋮h∇(i1,j1,k1)←+z⋅α1⋅β1⋅γ1.(25.95)
需要注意的是,坐标
k
k
k 代表一个方向,因此必须以循环方式处理,如图25.26所示(也见算法25.9的第11-12行)。
对于每个直方图箱,贡献梯度样本的范围覆盖相邻箱的一半,即相邻箱的支持区域重叠,如图25.27所示。
图25.25
对于 n spat = 4 n_{\text{spat}} = 4 nspat=4 和 n angl = 8 n_{\text{angl}} = 8 nangl=8 的SIFT描述符结构。每个16个空间箱 i j = A 1 , … , D 4 ij = A1, \ldots, D4 ij=A1,…,D4 提供8个方向箱 k = 0 , … , 7 k = 0, \ldots, 7 k=0,…,7。因此,梯度直方图 h ∇ h_\nabla h∇ 包含128个单元,这些单元排列成一个1D特征向量(如(b)所示的 A 10 , A 12 , … , D 46 , D 47 A10, A12, \ldots, D46, D47 A10,A12,…,D46,D47)。
图25.26
梯度直方图的3D结构,其中空间维度 ( i , j ) (i, j) (i,j) 为 n spat × n spat = 4 × 4 n_{\text{spat}} \times n_{\text{spat}} = 4 \times 4 nspat×nspat=4×4 个箱,方向轴 ( k ) (k) (k) 上有 n angl = 8 n_{\text{angl}} = 8 nangl=8 个箱。为了将量 z z z 累积到某个连续位置 ( i ′ , j ′ , k ′ ) (i', j', k') (i′,j′,k′),八个相邻的箱接收 z z z 的不同分位数,这些分位数由三线性插值确定 (a)。注意,方向轴 ϕ \phi ϕ 上的箱以循环方式处理;例如, k = 0 k = 0 k=0 处的箱也被视为与 k = 7 k = 7 k=7 处的箱相邻 (b)。图25.27
梯度场中的重叠支持区域。由于在直方图计算中使用了三线性插值,方向直方图 h ∇ h_\nabla h∇ 的单元相关的空间区域重叠。圆圈的阴影表示由高斯加权函数分配给每个样本的权重 w G w_G wG,其值取决于每个样本与关键点中心的距离(见方程25.88)。
SIFT描述符的归一化
梯度直方图
h
∇
h_\nabla
h∇ 的元素是SIFT特征向量
f
sift
f_{\text{sift}}
fsift 的原材料。从梯度直方图计算特征向量的过程在算法25.10中描述。首先,将大小为
n
spat
×
n
spat
×
n
angl
n_{\text{spat}} \times n_{\text{spat}} \times n_{\text{angl}}
nspat×nspat×nangl 的3D梯度直方图
h
∇
h_\nabla
h∇(包含连续值)展平为长度为
n
spat
2
⋅
n
angl
n_{\text{spat}}^2 \cdot n_{\text{angl}}
nspat2⋅nangl(通常为128)的1D向量
f
f
f,其中
f
[
(
i
⋅
n
spat
+
j
)
⋅
n
angl
+
k
]
←
h
∇
(
i
,
j
,
k
)
,
(25.96)
f[(i \cdot n_{\text{spat}} + j) \cdot n_{\text{angl}} + k] \leftarrow h_\nabla(i, j, k), \tag{25.96}
f[(i⋅nspat+j)⋅nangl+k]←h∇(i,j,k),(25.96)
其中
i
,
j
=
0
,
…
,
n
spat
−
1
i, j = 0, \ldots, n_{\text{spat}}-1
i,j=0,…,nspat−1,
k
=
0
,
…
,
n
angl
−
1
k = 0, \ldots, n_{\text{angl}}-1
k=0,…,nangl−1。因此,
f
f
f 中的元素按图25.25所示的顺序排列,方向索引
k
k
k 变化最快,空间索引
i
i
i 变化最慢(见算法25.10,第3-8行)。
图像对比度的变化对梯度幅度有线性影响,因此也对特征向量
f
f
f 的值有影响。为了消除这些影响,随后对向量
f
f
f 进行归一化:
f
(
m
)
←
1
∥
f
∥
⋅
f
(
m
)
,
(25.97)
f(m) \leftarrow \frac{1}{\|f\|} \cdot f(m), \tag{25.97}
f(m)←∥f∥1⋅f(m),(25.97)
对于所有
m
m
m,使得
f
f
f 具有单位范数(见算法25.10,第9行)。由于梯度是从局部像素差异中计算得出的,绝对亮度的变化不会影响梯度幅度,除非发生饱和。这样的非线性照明变化往往会产生峰值梯度值,通过将
f
f
f 的值裁剪到预定义的最大值
t
f
c
l
i
p
t_{fclip}
tfclip 来补偿,即
f
(
m
)
←
min
(
f
(
m
)
,
t
f
c
l
i
p
)
,
(25.98)
f(m) \leftarrow \min(f(m), t_{fclip}), \tag{25.98}
f(m)←min(f(m),tfclip),(25.98)
通常
t
f
c
l
i
p
=
0.2
t_{fclip} = 0.2
tfclip=0.2,如文献[153]中所建议(见算法25.10,第10行)。此步骤后,
f
f
f 再次按照方程 (25.97) 进行归一化。
最后,将实值特征向量
f
f
f 转换为整数向量:
f
sift
(
m
)
←
min
(
round
(
s
f
s
c
a
l
e
⋅
f
(
m
)
)
,
255
)
,
(25.99)
f_{\text{sift}}(m) \leftarrow \min(\text{round}(s_{fscale} \cdot f(m)), 255), \tag{25.99}
fsift(m)←min(round(sfscale⋅f(m)),255),(25.99)
其中
s
f
s
c
a
l
e
s_{fscale}
sfscale 是预定义常数(通常
s
f
s
c
a
l
e
=
512
s_{fscale} = 512
sfscale=512)。
f
sift
f_{\text{sift}}
fsift 的元素范围为 [0, 255],以方便编码和存储为字节序列(见算法25.10,第12行)。
给定关键点
k
′
=
(
p
,
q
,
x
,
y
)
k' = (p, q, x, y)
k′=(p,q,x,y) 的最终SIFT描述符是一个元组
s
=
(
x
′
,
y
′
,
σ
,
θ
,
f
sift
)
,
(25.100)
s = (x', y', \sigma, \theta, f_{\text{sift}}), \tag{25.100}
s=(x′,y′,σ,θ,fsift),(25.100)
其中包含关键点的插值位置
x
′
,
y
′
x', y'
x′,y′(在原始图像坐标中)、绝对尺度
σ
\sigma
σ、其主方向
θ
\theta
θ 以及相应的整数值梯度特征向量
f
sift
f_{\text{sift}}
fsift(见算法25.8,第27行)。请记住,对于位于同一关键点位置的不同主方向,可能会生成多个SIFT描述符。这些描述符将具有相同的位置和尺度值,但不同的
θ
\theta
θ 和
f
sift
f_{\text{sift}}
fsift 数据。