SIFT(尺度不变特征变换,Scale-lnvariantFeatureTransform)是计算机视觉领域,检测和描述图像局部特征的算法,该算法于1999年被DavidLowe提出,并于2004年进行了补充和完善。之所以称它为变换,是因为它会将图像数据变换为相对于局部图像特征的坐标。
当两幅图像性质相似时,角检测核MSER适合作为整体图像特征。然而,当出现变化(旋转、光照变化、视点变化)时,就应使用SIFT。SIFT特征(称为关键点)对图像尺度和旋转是不变的,并且对仿射失真、三维视点变化、噪声和光照变化具有很强的鲁棒性。SIFT的输入是一副图像,输出是一个n维特征向量。
该算法应用很广,如目标识别,自动导航,图像拼接,三维建模,手势识别,视频跟踪等。
SIFT算法分为4个阶段:
- 尺度空间极值检测: 在所有图像上进行搜索,应用高斯差分函数识别出特征点,这些特征点具有尺度不变性和旋转不变性;
- 改进特征点位置精度:一个连续函数被取样时,它真正的最大值或最小值实际上可能位于样本点之间。为了提高特征点位置精度,我们可以在每个极值点处拟合一个内插函数,利用内插函数寻找精度更高的极值点;
- 确定特征点方向:根据图像局部性质,为每个特征点分配一个方向。SIFT使用了一种简单的方法,使用特征点的尺度来确定高斯平滑图像L,如此一来,所有方向的计算就都以尺度不变的方式执行。对该尺度的每个图像样本 L ( x , y ) L(x,y) L(x,y),使用像素差计算梯度幅度和方向角 M ( x , y ) M(x,y) M(x,y)和 θ ( x , y ) \theta(x,y) θ(x,y)。
- 特征点描述符:计算每个特征点的特征向量。
下面详细阐述这4个阶段。
尺度空间极值检测
特征点检测的第一步是识别目标的位置和尺度,同一个目标在不同的视角下,位置和尺度是不同的。理论上,位置不随图像尺度变化而变化。
现实世界中,物体只有在一定尺度上才有意义。例如,对于整个银河系来说,桌子上的水杯是不存在的;但对于人类来说,这个水杯是存在的。尺度空间即试图在数字图像领域描述该意义。
在一副图像中,我们更想看到叶子还是更想看到整棵树?如果答案是后者,那么我们应该去除图像的细节部分(如叶子、细枝)。在去除细节部分的过程中,我们一定要确保不能引进错误的细节。因此在创建尺度空间的过程中,我们应该对原始图像逐渐的做模糊平滑处理。进行该操作的唯一方法是高斯模糊,因为高斯函数是唯一的尺度空间核。
图像的尺度空间用 L ( x , y , σ ) L(x,y, σ) L(x,y,σ)表示,公式如下:
L ( x , y , σ ) = G ( x , y , σ ) ⊕ I ( x , y ) L(x,y,σ)=G(x,y, σ)⊕I(x,y) L(x,y,σ)=G(x,y,σ)⊕I(x,y)
(1)
其中,⊕表示卷积,G表示二维高斯函数:
G ( x , y , σ ) = 1 2 π σ 2 e − x 2 + y 2 2 σ 2 G(x,y,σ)=\frac{1}{2\pi\sigma^2} e^{-\frac{x^2 + y^2}{2\sigma^2}} G(x,y,σ)=2πσ21e−2σ2x2+y2
(2)
σ是尺度空间因子,它决定了图像模糊的程度。大尺度(σ值大)表现的是图像的概貌信息,小尺度(σ值小)表现的是图像的细节信息。因此大尺度对应着低分辨率,小尺度对应着高分辨率。
需要说明的是,公式中的图像I具有无限分辨率,也就是说它的尺度空间因子为0。很显然,无限分辨率的图像是无法获得的,Lowe把初始图像的尺度空间因子设定为0.5。由 L ( x , y , σ 1 ) L(x,y,σ_1) L(x,y,σ1)得到 L ( x , y , σ 2 ) L(x,y,σ_2) L(x,y,σ2)的公式为:
L ( x , y , σ 2 ) = G ( x , y , σ 2 2 − σ 1 2 ) ⊕ L ( x , y , σ 1 ) L(x,y,σ_2)=G(x,y, \sqrt{σ^2_2-σ^2_1})⊕L(x,y,σ_1) L(x,y,σ2)=G(x,y,σ22−σ12)⊕L(x,y,σ1)
(3)
其中,
G ( x , y , σ 2 2 − σ 1 2 ) = 1 2 π ( σ 2 2 − σ 1 2 ) e − x 2 + y 2 2 ( σ 2 2 − σ 1 2 ) G(x,y, \sqrt{σ^2_2-σ^2_1})=\frac{1}{2\pi(σ^2_2-σ^2_1)} e^{-\frac{x^2 + y^2}{2(σ^2_2-σ^2_1)}} G(x,y,σ22−σ12)=2π(σ22−σ12)1e−2(σ22−σ12)x2+y2
(4)
因为我们无法得到尺度为0 的图像,所以利用公式3生成,用小尺度的图像生成大尺度的图像。
利用LoG(高斯拉普拉斯方法,Laplacian of Gaussian),即图像的二阶导数,能够在不同的尺度下检测到图像特征,从而确定特征点。但是LoG的效率不高,因此SIFT算法进行了改进,采用另外一种方式确定特征点。将尺度空间中相邻图像相减,得到DoG(高斯差分,Difference of Gaussians)图像。
D ( x , y , σ ) = ( G ( x , y , k σ ) − G ( x , y , σ ) ) ⊕ I ( x , y ) = L ( x , y , k σ ) − L ( x , y , σ ) D(x,y, σ)=(G(x,y, kσ)-G(x,y, σ))⊕I(x,y)=L(x,y, kσ)-L(x,y, σ) D(x,y,σ)=(G(x,y,kσ)−G(x,y,σ))⊕I(x,y)=L(x,y,kσ)−L(x,y,σ)
(5)
其中,k为两个相邻尺度空间倍乘的常数。
可以证明DoG是LoG的近似,并且用DoG 代替LoG 并不影响特征点检测。用DoG 近似LoG还有以下好处:
- LoG 需要使用两个方向的高斯二阶微分卷积核,而DoG 直接使用高斯卷积核,省去了卷积核生成时间;
- DoG 保留了高斯尺度空间的图像,因此在生成某一空间尺度的特征时,可以直接使用公式1(或公式3)产生的尺度空间图像,而无需重新生成该尺度的图像;
- DoG 稳定性高、抗干扰能力强。
为了在连续的尺度下检测图像的特征点,需要建立 DoG 金字塔,而DoG金字塔的建立又离不开高斯金字塔,如下图所示,左侧为高斯金字塔,右侧为 DoG金字塔。
高斯金字塔共分O组(Octave),每组又分S层(Layer)。组内各层图像的分辨率是相同的,尺度逐渐增加,即越往塔顶图像越模糊。而下一组的图像是由上一组图像按照隔点降采样得到的,即图像的长和宽分别减半。高斯金字塔的组数O是由输入图像的分辨率得到的。O应足够大,以得到数量充足的图像,考虑到O过大(组数过多),某个组内图像分辨率会过低,所以使用如下公式确定O的取值:
O = ⌊ l o g 2 m i n ( X , Y ) − 2 ⌋ O = \lfloor log_2min(X,Y)-2\rfloor O=⌊log2min(X,Y)−2⌋(6)
其中,X 和Y 分别为输入图像的长和宽,⌊ ⌋表示向下取整。 S = s + 3 S = s + 3 S=s+3(7)
建立高斯金字塔的过程如下:输入图像的尺度为0.5,由该图像得到高斯金字塔第0组的第0层图像,它的尺度为
σ
0
σ_0
σ0,我们称
σ
0
σ_0
σ0为基准层尺度,再由第0层得到第1层,它的尺度为
k
σ
0
kσ_0
kσ0,第2层的尺度为
k
2
σ
0
k^2σ_0
k2σ0,以此类推。这里的k为:
k
=
2
1
s
k=2^{\frac1s}
k=2s1
(8)
我们以s=3为例,第0组的6幅图像的尺度分别为:
σ 0 , k σ 0 , k 2 σ 0 , k 3 σ 0 , k 4 σ 0 , k 5 σ 0 σ_0,kσ_0,k^2σ_0,k^3σ_0,k^4σ_0,k^5σ_0 σ0,kσ0,k2σ0,k3σ0,k4σ0,k5σ0
(9)
写成更一般的公式:
σ = k r σ 0 σ=k^rσ_0 σ=krσ0
(10)
各层图像的尺度公式为:
σ ( o , r ) = 2 o k r σ 0 σ(o,r)=2^ok^rσ_0 σ(o,r)=2okrσ0
(11)
将公式8代入公式11得:
σ ( o , r ) = 2 o + r s σ 0 σ(o,r)=2^{o+\frac rs}σ_0 σ(o,r)=2o+srσ0
(12)
为了最大程度的保留原图的信息量,论文中建议,在创建尺度空间前首先对输入图像的长宽扩展一倍。如果输入图像的尺度为0.5,那么长宽扩展一倍后的图像尺度为1。
DoG金字塔是由高斯金字塔得到的,高斯金字塔组内相邻两层图像相减得到DoG金字塔,两组间的各层不能相减。高斯金字塔每组有s+3层图像,DoG 金字塔每组有s+2层图像。
搜索极值点是在DoG金字塔内进行的,在搜索之前应剔除像素值过小的点,因为这些点对比度较低,必定不是稳定的特征点。搜索范围为尺度空间图像的邻域。
改进特征点位置精度
第一步我们得到极值点,但这些极值点还需要经过筛选。原因是:极值点的搜索是在离散空间内进行的,如果把采样点拟合成曲面,我们会发现,第一步得到的极值点并不是真正的极值点(也即精度不够高)。离散空间的极值点并不是连续空间的极值点。如果我们想获得更精确的极值点,必须进行拟合。
使用泰勒级数展开式作为拟合函数。极值点是一个三维矢量
X
=
(
x
,
y
,
σ
)
T
X = (x, y, σ)^T
X=(x,y,σ)T,因此需要三维函数的泰勒级数展开式。在
X
0
=
(
x
0
,
y
0
,
σ
0
)
T
X_0 = (x_0, y_0, σ_0)^T
X0=(x0,y0,σ0)T 处进行泰勒级数展开,矩阵形式为:
x
y
σ
=
x
0
y
0
σ
0
+
[
∂
f
∂
x
∂
f
∂
y
∂
f
∂
σ
]
(
x
y
σ
−
x
0
y
0
σ
0
)
+
1
2
(
[
x
y
σ
]
−
[
x
0
y
0
σ
0
]
)
∂
2
f
∂
x
∂
x
∂
2
f
∂
x
∂
y
∂
2
f
∂
x
∂
σ
∂
2
f
∂
x
∂
y
∂
2
f
∂
y
∂
y
∂
2
f
∂
y
∂
σ
∂
2
f
∂
x
∂
σ
∂
2
f
∂
y
∂
σ
∂
2
f
∂
σ
∂
σ
(
x
y
σ
−
x
0
y
0
σ
0
)
\begin{matrix} x \\y \\σ \end{matrix} =\begin{matrix} x_0 \\y_0 \\σ_0 \end{matrix}+ [\frac{\partial f} {\partial x} \frac{\partial f} {\partial y}\frac{\partial f} {\partial σ}] (\begin{matrix} x \\y \\σ \end{matrix} -\begin{matrix} x_0 \\y_0 \\σ_0 \end{matrix}) + \frac12([x\:y\:σ]-[x_0\:y_0\:σ_0]) \begin{matrix} \frac{\partial ^2f}{\partial x\partial x}\ \frac{\partial ^2f} {\partial x\partial y}\ \frac{\partial ^2f} {\partial x\partial σ}\\ \frac{\partial ^2f} {\partial x\partial y}\ \frac{\partial ^2f} {\partial y\partial y}\ \frac{\partial ^2f} {\partial y\partial σ}\\ \frac{\partial ^2f} {\partial x\partial σ}\ \frac{\partial ^2f} {\partial y\partial σ}\ \frac{\partial ^2f} {\partial σ\partial σ} \end{matrix} (\begin{matrix} x \\y \\σ \end{matrix} -\begin{matrix} x_0 \\y_0 \\σ_0 \end{matrix})
xyσ=x0y0σ0+[∂x∂f∂y∂f∂σ∂f](xyσ−x0y0σ0)+21([xyσ]−[x0y0σ0])∂x∂x∂2f ∂x∂y∂2f ∂x∂σ∂2f∂x∂y∂2f ∂y∂y∂2f ∂y∂σ∂2f∂x∂σ∂2f ∂y∂σ∂2f ∂σ∂σ∂2f(xyσ−x0y0σ0)
(13)
式13舍去了高阶项,其矢量表示形式为:
f ( X ) = f ( X 0 ) + ∂ f T ∂ X ( X − X 0 ) + 1 2 ( X − X 0 ) T ∂ 2 f ∂ X 2 ( X − X 0 ) f(X)=f(X_0)+\frac {\partial f^T}{\partial X}(X-X_0) + \frac12(X-X_0)^T \frac{\partial^2f}{\partial X^2}(X-X_0) f(X)=f(X0)+∂X∂fT(X−X0)+21(X−X0)T∂X2∂2f(X−X0)
(14)
对上式求导得:
f
(
X
)
∂
X
=
f
(
X
0
)
+
∂
f
T
∂
X
+
1
2
∂
2
f
∂
X
2
(
X
−
X
0
)
+
1
2
(
X
−
X
0
)
T
∂
2
f
∂
X
2
\frac {f(X)}{\partial X}=f(X_0)+\frac {\partial f^T}{\partial X} + \frac12 \frac{\partial^2f}{\partial X^2}(X-X_0) +\frac12(X-X_0)^T \frac{\partial^2f}{\partial X^2}
∂Xf(X)=f(X0)+∂X∂fT+21∂X2∂2f(X−X0)+21(X−X0)T∂X2∂2f
(15)
令上式为0,可求得偏移量为:
X
−
X
0
=
−
∂
2
f
∂
X
2
∂
f
∂
X
X-X_0=-\frac {\partial^2f}{\partial X^2}\frac {\partial f}{\partial X}
X−X0=−∂X2∂2f∂X∂f
(16)
将公式16代入公式14得到:
f
(
X
)
=
f
(
X
0
)
+
1
2
∂
f
T
∂
X
(
X
−
X
0
)
f(X)=f(X_0)+\frac 12 \frac {\partial f^T}{\partial X}(X-X_0)
f(X)=f(X0)+21∂X∂fT(X−X0)
(17)
在DoG图像内搜索极值点(特征点),该图像的边缘有很强的响应,边缘上的特征点是不稳定的。
确定特征点方向
经过上述操作,我们已找出图像的所有特征点,这些特征点具有尺度不变性。为了实现旋转不变性,还需要为特征点分配一个方向。
每个特征点除了分配一个主方向外,有可能需要分配一个或多个辅方向,增加辅方向的目的是为了增强图像匹配的鲁棒性。
特征点描述符
经过上述步骤,每个特征点被分配了坐标位置、尺度和方向。在图像局部区域内,下面将计算局部区域的描述符,描述符既具有可区分性,又具有不变性。最终得到128维的特征向量。