opencv2.4.9 源码分析——SIFT
赵春江
blog.csdn.net/zhaocj
一、SIFT 算法
SIFT(尺度不变特征变换,Scale‐Invariant Feature Transform)是在计算机视觉领域中检测
和描述图像中局部特征的算法,该算法于 1999 年被 David Lowe 提出,并于 2004 年进行了
补充和完善。该算法应用很广,如目标识别,自动导航,图像拼接,三维建模,手势识别,
视频跟踪等。不幸的是,该算法已经在美国申请了专利,专利拥有者为 Lowe 所在的加拿大
不列颠哥伦比亚大学,因此我们不能随意使用它。
用 SIFT 算法所检测到的特征是局部的,而且该特征对于图像的尺度和旋转能够保持不
变性。同时,这些特征对于亮度变化具有很强的鲁棒性,对于噪声和视角的微小变化也能保
持一定的稳定性。SIFT 特征还具有很强的可区分性,它们很容易被提取出来,并且即使在低
概率的不匹配情况下也能够正确的识别出目标来。因此鲁棒性和可区分性是 SIFT 算法最主
要的特点。
SIFT 算法分为 4 个阶段:
1、尺度空间极值检测:该阶段是在图像的全部尺度和全部位置上进行搜索,并通过应
用高斯差分函数可以有效地识别出尺度不变性和旋转不变性的潜在特征点来;
2、特征点的定位:在每个候选特征点上,一个精细的模型被拟合出来用于确定特性点
的位置和尺度。而特征点的最后选取依赖的是它们的稳定程度;
3、方向角度的确定:基于图像的局部梯度方向,为每个特性点分配一个或多个方向角
度。所有后续的操作都是相对于所确定下来的特征点的角度、尺度和位置的基础上进行的,
因此特征点具有这些角度、尺度和位置的不变性;
4、特征点的描述符:在所选定的尺度空间内,测量特征点邻域区域的局部图像梯度,
将这些梯度转换成一种允许局部较大程度的形状变形和亮度变化的描述符形式。
下面就详细阐述 SIFT 算法的这 4 个阶段:
1、尺度空间极值检测
特征点检测的第一步是能够识别出目标的位置和尺度,对于同一个目标在不同的视角下
这些位置和尺度可以被重复的分配。并且这些检测到的位置是不随图像尺度的变化而改变
的,因为它们是通过搜索所有尺度上的稳定特征得到的,所应用的工具就是被称为尺度空间
的连续尺度函数。
1
真实世界的物体只有在一定尺度上才有意义,例如我们能够看到放在桌子上的水杯,但
对于整个银河系,这个水杯是不存在的。物体的这种多尺度的本质在自然界中是普遍存在的。
尺度空间就是试图在数字图像领域复制这个概念。又比如,对于某幅图像,我们是想看到叶
子还是想看到整棵树,如果是树,那么我们就应该有意识的去除图像的细节部分(如叶子、
细枝等)。在去除细节部分的过程中,我们一定要确保不能引进新的错误的细节。因此在创
建尺度空间的过程中,我们应该对原始图像逐渐的做模糊平滑处理。进行该操作的唯一方法
是高斯模糊处理,因为已经被证实,高斯函数是唯一可能的尺度空间核。
图像的尺度空间用 L(x, y, σ)函数表示,它是由一个变尺度的高斯函数 G(x, y, σ)与图像
I(x, y)通过卷积产生,即
,,,,,
,, 12
其中,表示在 x 和 y 两个方向上进行卷积操作,而 G(x, y, σ)为
(1)
(2)
(3)
(4)
σ 是尺度空间因子,它决定着图像模糊平滑处理的程度。在大尺度下(σ 值大)表现的是图
像的概貌信息,在小尺度下(σ 值小)表现的是图像的细节信息。因此大尺度对应着低分辨
率,小尺度对应着高分辨率。(x, y)则表示在 σ 尺度下的图像像素坐标。
需要说明的是,公式 1 中的图像 I(x, y)是具有无限分辨率的图像,也就是说它的尺度
σ=0,即 I(x, y) = L(x, y, 0)。因此公式 1 得到的尺度空间图像 L(x, y, σ)是由尺度为 0 的图像
L(x, y, 0)生成的。很显然,尺度为 0,即无限分辨率的图像是无法获得的,Lowe 就是把初始
图像的尺度设定为 0.5。那么由 L(x, y, σ1)得到 L(x, y, σ2),即由尺度为 σ1 的图像生成尺度为
σ2 的图像的公式为:
,,,, ,,,
,,
1
2
其中,
由于尺度为 0 的图像无法得到,因此在实际应用中要想得到任意尺度下的图像,一定是
利用公式 3 生成的,即由一个已知尺度(该尺度不为 0)的图像生成另一个尺度的图像,并
且一定是小尺度的图像生成大尺度的图像。
利用 LoG(高斯拉普拉斯方法,Laplacian of Gaussian),即图像的二阶导数,能够在不同
的尺度下检测到图像的斑点特征,从而可以确定图像的特征点。但 LoG 的效率不高。因此
2
SIFT 算法进行了改进,通过对两个相邻高斯尺度空间的图像相减,得到一个 DoG(高斯差分,
Difference of Gaussians)的响应值图像 D(x, y, σ)来近似 LoG:
,,,,,,,,,,,
其中,k 为两个相邻尺度空间倍数的常数。
(5)
可以证明 DoG 是对 LoG 的近似表示,并且用 DoG 代替 LoG 并不影响对图像斑点位置的
检测。而且用 DoG 近似 LoG 可以实现下列好处:第一是 LoG 需要使用两个方向的高斯二阶
微分卷积核,而 DoG 直接使用高斯卷积核,省去了卷积核生成的运算量;第二是 DoG 保留
了个高斯尺度空间的图像,因此在生成某一空间尺度的特征时,可以直接使用公式 1(或公
式 3)产生的尺度空间图像,而无需重新再次生成该尺度的图像;第三是 DoG 具有与 LoG 相
同的性质,即稳定性好、抗干扰能力强。
为了在连续的尺度下检测图像的特征点,需要建立 DoG 金字塔,而 DoG 金字塔的建立
又离不开高斯金字塔的建立,如下图所示,左侧为高斯金字塔,右侧为 DoG 金字塔:
图 1 高斯金字塔和 DoG 金字塔
高斯金字塔共分 O 组(Octave),每组又分 S 层(Layer)。组内各层图像的分辨率是相
同的,即长和宽相同,但尺度逐渐增加,即越往塔顶图像越模糊。而下一组的图像是由上一
组图像按照隔点降采样得到的,即图像的长和宽分别减半。高斯金字塔的组数 O 是由输入
图像的分辨率得到的,因为要进行隔点降采样,所以在执行降采样生成高斯金字塔时,一直
到不能降采样为止,但图像太小又毫无意义,因此具体的公式为:
logmin,2
(6)
3
其中,X 和 Y 分别为输入图像的长和宽,⌊ ⌋表示向下取整。
金字塔的层数 S 为:
S = s + 3
(7)
Lowe 建议 s 为 3。需要注意的是,除了公式 7 中的第一个字母是大写的 S 外,后面出现的都
是小写的 s。
高斯金字塔的创建是这样的:设输入图像的尺度为 0.5,由该图像得到高斯金字塔的第
0 组的第 0 层图像,它的尺度为 σ0,我们称 σ0 为基准层尺度,再由第 0 层得到第 1 层,它
的尺度为 kσ0,第 2 层的尺度为 k2σ0,以此类推。这里的 k 为:
2
我们以 s=3 为例,第 0 组的 6(s+3=6)幅图像的尺度分别为:
σ0,kσ0,k2σ0,k3σ0,k4σ0,k5σ0
写成更一般的公式为:
σ = krσ0 r∈[0,…,s+2]
(8)
(9)
(10)
第 0 组构建完成后,再构建第 1 组。第 1 组的第 0 层图像是由第 0 组的倒数第 3 层图
像经过隔点采样得到的。由公式 10 可以得到,第 0 组的倒数第 3 层图像的尺度为 ksσ0,k 的
值代入公式 8,得到了该层图像的尺度正好为 2σ0,因此第 1 组的第 0 层图像的尺度仍然是
2σ0。但由于第 1 组图像是由第 0 组图像经隔点降采样得到的,因此相对于第 1 组图像的分
辨率来说,第 0 层图像的尺度为 σ0,即尺度为 2σ0 是相对于输入图像的分辨率来说的,而尺
度为 σ0 是相对于该组图像的分辨率来说的。这也就是为什么我们称 σ0 为基准层尺度的原因
(它是每组图像的基准层尺度)。第 1 组其他层图像的生成与第 0 组的相同。因此可以看出,
第 1 组各层图像的尺度相对于该组分辨率来说仍然满足公式 10。这样做的好处就是编程的
效率会提高,并且也保证了高斯金字塔尺度空间的连续性。而之所以会出现这样的结果,是
因为在参数选择上同时满足公式 7、公式 8 以及对上一组倒数第 3 层图像降采样这三个条件
的原因。
那么第 1 组各层图像相对于输入图像来说,它们的尺度为:
σ = 2krσ0 r∈[0,…,s+2]
(11)
该公式与公式 10 相比较可以看出,第 1 组各层图像的尺度比第 0 组相对应层图像的尺
度大了一倍。高斯金字塔的其他组的构建以此类推,不再赘述。下面给出相对于输入图像的
4
各层图像的尺度公式:
σ(o,r) = 2okrσ0 o∈[0,…,O-1], r∈[0,…,s+2]
(12)
其中,o 表示组的坐标,r 表示层的坐标,σ0 为基准层尺度。k 用公式 8 代入,得:
,2 ∈0,…,1, ∈0,…,2
(13)
在高斯金字塔中,第 0 组第 0 层的图像是输入图像经高斯模糊后的结果,模糊后的图像
的高频部分必然会减少,因此为了最大程度的保留原图的信息量,Lowe 建议在创建尺度空
间前首先对输入图像的长宽扩展一倍,这样就形成了高斯金字塔的第‐1 组。设输入图像的尺
度为 0.5,那么相对于输入图像,分辨率扩大一倍后的尺度应为 1,由该图像依次进行高斯
平滑处理得到第‐1 组的各个层的尺度图像,方法与其他组的一样。由于增加了第‐1 组,因此
公式 13 重新写为:
,2 ∈1,0,…,1, ∈0,…,2
(14)
DoG 金字塔是由高斯金字塔得到的,即高斯金字塔组内相邻两层图像相减得到 DoG 金
字塔。如高斯金字塔的第 0 组的第 0 层和第 1 层相减得到 DoG 金字塔的第 0 组的第 0 层图
像,高斯金字塔的第 0 组的第 1 层和第 2 层相减得到 DoG 金字塔的第 0 组的第 1 层图像,
以此类推。需要注意的是,高斯金字塔的组内相邻两层相减,而两组间的各层是不能相减的。
因此高斯金字塔每组有 s+3 层图像,而 DoG 金字塔每组则有 s+2 层图像。
极值点的搜索是在 DoG 金字塔内进行的,这些极值点就是候选的特征点。
在搜索之前,我们需要在 DoG 金字塔内剔除那些像素值过小的点,因为这些像素具有
较低的对比度,它们肯定不是稳定的特征点。
极值点的搜索不仅需要在它所在尺度空间图像的邻域内进行,还需要在它的相邻尺度空
间图像内进行,如图 2 所示。
每个像素在它的尺度图像中一共有 8 个相邻点,而在它的下一个相邻尺度图像和上一
个相邻尺度图像还各有 9 个相邻点(图 2 中绿色标注的像素),也就是说,该点是在 3×3×3
的立方体内被包围着,因此该点在 DoG 金字塔内一共有 26 个相邻点需要比较,来判断其是
否为极大值或极小值。这里所说的相邻尺度图像指的是在同一个组内,因此在 DoG 金字塔
内,每一个组的第 0 层和最后一层各只有一个相邻尺度图像,所以在搜索极值点时无需在这
两层尺度图像内进行,从而使极值点的搜索就只在每组的中间 s 层尺度图像内进行。
搜索的过程是这样的:从每组的第 1 层开始,以第 1 层为当前层,对第 1 层的 DoG 图
5
像中的每个点取一个 3×3×3 的立方体,立方体上下层分别为第 0 层和第 2 层。这样,搜索得
到的极值点既有位置坐标(该点所在图像的空间坐标),又有尺度空间坐标(该点所在层的
尺度)。当第 1 层搜索完成后,再以第 2 层为当前层,其过程与第 1 层的搜索类似,以此类
推。
图 2 DoG 中极值点的搜索
2、特征点的定位
通过上一步,我们得到了极值点,但这些极值点还仅仅是候选的特征点,因为它们还存
在一些不确定的因素。首先是极值点的搜索是在离散空间内进行的,并且这些离散空间还是
经过不断的降采样得到的。如果把采样点拟合成曲面后我们会发现,原先的极值点并不是真
正的极值点,也就是离散空间的极值点并不是连续空间的极值点。在这里,我们是需要精确
定位特征点的位置和尺度的,也就是要达到亚像素精度,因此必须进行拟合处理。
我们使用泰勒级数展开式作为拟合函数。如上所述,极值点是一个三维矢量,即它包括
极值点所在的尺度,以及它的尺度图像坐标,即 X = (x, y, σ)T,因此我们需要三维函数的泰
勒级数展开式,设我们在 X0 = (x0, y0, σ0)T 处进行泰勒级数展开,则它的矩阵形式为:
6
12
∂∂12
(15)
公式 15 为舍去高阶项的形式,而它的矢量表示形式为:
合后连续空间下的插值点坐标,设,则表示相对于插值中心,插值后的偏移量。
在这里 X0 表示离散空间下的插值中心(在离散空间内也就是采样点)坐标,X 表示拟
(16)
因此公式 16 经过变量变换后,又可写成:
∂∂12
∂∂ 12
对上式求导,得
(17)
(18)
(19)
让公式 17 的导数为 0,即公式 18 为 0,就可得到极值点下的相对于插值中心 X0 的偏移量:
把公式 19 得到的极值点带入公式 17 中,就得到了该极值点下的极值:
∂
∂
∂∂12
∂∂12
∂∂12
∂∂12
12
(20)
对于公式 19 所求得的偏移量如果大于 0.5(只要 x、y 和 σ 任意一个量大于 0.5),则表
明插值点已偏移到了它的临近的插值中心,所以必须改变当前的位置,使其为它所偏移到的
7
插值中心处,然后在新的位置上重新进行泰勒级数插值拟合,直到偏移量小于 0.5 为止(x、
y 和 σ 都小于 0.5),这是一个迭代的工程。当然,为了避免无限次的迭代,我们还需要设置
一个最大迭代次数,在达到了迭代次数但仍然没有满足偏移量小于 0.5 的情况下,该极值点
就要被剔除掉。另外,如果由公式 20 所得到的极值过小,即0.03时(假设图
像的灰度值在 0~1.0 之间),则这样的点易受到噪声的干扰而变得不稳定,所以这些点也应
该剔除。而在 opencv 中,使用的是下列公式来判断其是否为不稳定的极值:
(21)
其中 T 为经验阈值,系统默认初始化为 0.04。
极值点的求取是在 DoG 尺度图像内进行的,DoG 图像的一个特点就是对图像边缘有很
强的响应。一旦特征点落在图像的边缘上,这些点就是不稳定的点。这是因为一方面图像边
缘上的点是很难定位的,具有定位的歧义性;另一方面这样的点很容易受到噪声的干扰而变
得不稳定。因此我们一定要把这些点找到并剔除掉。它的方法与 Harris 角点检测算法相似,
即一个平坦的 DoG 响应峰值往往在横跨边缘的地方有较大的主曲率,而在垂直边缘的方向
上有较小的主曲率,主曲率可以通过 2×2 的 Hessian 矩阵 H 求出:
,, ,
, ,
(22)
其中 Dxx(x, y)、Dyy(x, y)和 Dxy(x, y)分别表示对 DoG 图像中的像素在 x 轴方向和 y 轴方向上求
二阶偏导和二阶混合偏导。在这里,我们不需要求具体的矩阵 H 的两个特征值——α 和 β,
而只要知道两个特征值的比例就可以知道该像素点的主曲率。
矩阵 H 的直迹和行列式分别为:
Tr(H) = Dxx + Dyy = α + β
Det(H) = DxxDyy – (Dxy)2 = α β
(23)
(24)
我们首先剔除掉那些行列式为负数的点,即 Det(H) < 0,因为如果像素的曲率有不同的
符号,则该点肯定不是特征点。设 α > β,并且 α = γ β,其中 γ > 1,则
Tr
Det
1
(25)
上式的结果只与两个特征值的比例有关,而与具体的特征值无关。我们知道,当某个像
素的 H 矩阵的两个特征值相差越大,即 γ 很大,则该像素越有可能是边缘。对于公式 25,
当两个特征值相等时,等式的值最小,随着 γ 的增加,等式的值也增加。所以,要想检查主
曲率的比值是否小于某一阈值 γ,只要检查下式是否成立即可:
8