本文对 SIFT 算法进行了详细梳理。SIFT即尺度不变特征变换(Scale-Invariant Feature Transform),是一种用于检测和描述图像局部特征的算法。该算法对图像的尺度和旋转具有不变性,并且在一定程度上能够抵御亮度变化和视角变化,具备较强的鲁棒性。此外,SIFT检测流程还提出了一些技巧,以加速特征提取过程。作为传统特征提取领域中极为重要的算法之一,SIFT 算法值得深入研究。
本文首先介绍高斯二阶偏导在尺度检测中的作用,进而引入二维拉普拉斯高斯算子用于blob特征提取。接着,阐述了两种提升卷积效率的方法。最后,对尺度、旋转和仿射变换进行统一处理,构建出最终的 SIFT 描述子。
参考资料:
计算机视觉 鲁鹏 清晰完整合集
opencv学习笔记(二十二):图像特征检测之SIFT算法
尺度不变特征转换-SIFT
理想特征曲线
在之前角点检测的文章中,我们发现Harris角点对于尺度是不敏感的,在不同尺度下角点特征可能检测不到。那么如何让我们的检测特征在不同尺度上都有效呢?
如下图,我们希望对房子这一特征在不同的缩放尺度上都能检测,一个理想的特征曲线(w.r.t尺度)应该是在刚好能够框住房子的尺度下,特征检测响应结果达到最大。这样,在不同尺度下都能够准确检测出特征。
高斯二阶偏导的尺度检测作用
我们知道,用高斯一阶偏导可以进行边缘的检测,如图所示,边缘也就是信号发生剧烈变化的地方,导数会出现极值。用高斯偏导核对信号做卷积,并找出极值点即可得到边缘信息。
而高斯的二阶偏导同样也能用于边缘的检测,如图,其与信号卷积的结果在边缘处的值为0。
这个高斯二阶偏导称为高斯Laplacian算子。
key idea
关键点来了,当我们使用一个确定σ的高斯二阶偏导和不同尺度的信号进行卷积,其响应会在某一特定尺度下达到最大。如图:
图中σ为1的Laplacian和不同尺度信号卷积,我们发现在最右边那个尺度下其中心响应达到最大。
也就是说,在σ确定的时候,这个高斯Laplacian能够对不同的尺度的信号进行一个筛选作用。其对应的那个中心响应最大的尺度,就是筛选出来的对应尺度。
确定σ,筛选不同尺度方波
我们可以从波形叠加的角度进一步分析这个过程,可以看到上图中的方波信号其实是由两个边缘信号组成的,我们可以把高斯拉普拉斯算子作用于方波信号看作是两个步骤:
第一步,将拉普拉斯分别作用于两个边缘信号;
第二步,将响应结果的波形进行叠加。
这方便我们理解,也有利于后续的分析过程。
从前面的分析中我们可以知道Laplacian算子作用于左边那个边缘信号上时会产生这样的一个波形:
作用于右边那个边缘信号上会产生一个水平翻转的波形:
随着方波信号的尺度不断减小,两个边缘不断靠近,这两个波形也逐步靠近,当足够靠近的时候,两个波形产生了叠加,中心的响应会增大;当两个响应的波谷峰值重叠到一起时,整体响应中心达到了一个最大值;那如果尺度进一步缩小呢?两个波形再往里叠加,又会相互抵消,所以响应又会减少;直到两个边缘缩小到重叠在一起,两个响应波形也重叠在一起,从波形图上可知,两个波形完全抵消掉了,此时响应为0。
上述过程中响应中心最大的那个方波的尺度,就是匹配当前Laplacian算子(由sigma决定)的方波尺度。
我们用一个确定的sigma,从不同尺度的方波中找出了一个特定尺度的方波信号。两者match了。
用这个理念,就能从信号中筛选出来特定尺度的信号。
同理,如果我们用很多不同σ的Laplacian对同一个方波信号做卷积,看哪个sigma能产生最大响应,那么这个sigma就能表示该方波的尺度。
确定方波,用多个σ进行匹配
然而,当我们用各种不同sigma的Laplacian去对同一个方波做卷积时,会产生一个问题:用Laplacian去做卷积,随着尺度的增加是有一个衰减的效应的,如图:
从这张图可以看到,随着sigma的增加,两个边缘的响应波形位置不变(因为方波确定,其两个边缘就确定),但波形的重叠度越来越高,中心的响应应该越来越大才是,但这里由于Laplacian随尺度增大而衰减的效应,导致即使两波型的波谷重叠,也达不到最大值,反而近乎消失了。
也就是说Laplacian的衰减作用太大,抵消了波形叠加的作用。
我们期望的情况是,无论sigma是多少,对于单个边缘的响应波形,峰值是一定的。单个边缘的响应峰值一定,考虑两个边缘的响应叠加才有意义。
如何让σ增加而响应峰值保持一定?我们先给出结论:
对于高斯偏导的响应,我们对其乘以σ进行补偿。
对于高斯二阶偏导的响应(Laplacian),我们对其乘以σ的平方进行补偿。
下图是补偿了σ的平方后,响应结果展示,可以看到,随σ变大,两个边缘的laplacian响应结果峰值不再衰减。此时影响响应结果的只有波形叠加的多少,可以看到,在σ=8时,叠加达到了最大,那么此方波尺度所对应的σ就是8。完成了该方波的尺度特性提取。
接下来我们从数学上对响应结果的σ补偿其进行解释。
响应结果归一化的数学分析
通过上述的讨论,我们知道一个Laplacian对方波信号的响应实际是其对两个阶跃信号(边缘)的响应波形叠加,只要这个阶跃信号的响应波形的峰值不随σ变化,我们的目的就达到了。
所以我们把目标定位到研究阶跃信号和高斯偏导、高斯二阶偏导的响应结果。
对高斯偏导的响应归一化
我们先分析高斯偏导的响应归一化。如图
高斯偏导和阶跃函数做卷积,其最高峰值是在信号突变处,也就是x = 0这里,在这里卷积得到的最大值是什么呢?从图像中我们可以看到在x < 0 时阶跃信号为0,所以卷积结果为0,在x>0时,阶跃信号为1,其卷积结果是高斯偏导函数和x轴的面积,也就是(0, ∞)的积分结果。
我们把这个积分写出来,设高斯函数为f(x),则高斯偏导函数:
(0, ∞)的积分结果:
这个积分结果就是高斯偏导在阶跃信号边缘处的卷积结果,也是整体阶跃信号响应的峰值处。
现在要让这个最大值在不同的σ下保持一致,就需要乘以一个σ,使其变为常数。
不同的σ高斯偏导卷积响应结果都乘以σ,就实现了各个响应结果的归一化。
对高斯二阶偏导的响应归一化
同样的,我们也希望在不同的σ下,高斯二阶偏导响应结果的峰值是一样的。和高斯偏导不同,这个响应结果峰值不在x=0处,图中可知x=0时响应结果为0,因为此时高斯二阶偏导在0上0下都有值,且面积相等,相互抵消。而当高斯二阶偏导和阶跃信号如下图所示时,卷积响应是最大的:
我们用面积来理解是比较方便的,此时的面积没有0之下的”负面积“抵消,所以这时候是最大的。由高斯二阶偏导的函数:
令其为0,可知x = ±σ时,卷积响应最大。
我们再把这个积分计算一下,这时我们要找的就是(σ, ∞)的积分表达式了。
这个积分是响应最大处的值,要保持其与σ无关,则需要将响应乘以σ的平方。
由此,我们就得到了前文所述的补偿结论。
2D Blob检测
拉普拉斯高斯算子
和一维空间一样,我们在二维空间中也可以用高斯二阶偏导来进行尺度的提取。这个高斯二阶偏导也叫做Laplacian of Gaussian,拉普拉斯高斯算子。它长这样:
公式为:
其中二维高斯函数:
代入求导得到:
当然,我们也得用σ平方做补偿:
尺度的选择
我们用不同的σ的算子和图像进行卷积,通过寻找响应最大的那个σ,就能确定当前的信号尺度。
这个是二维空间中一个信号:
那么什么情况下响应能达到最大值呢?σ和具体的图像信号尺度有什么关系呢?
我们看下图三种不同的信号尺度与同一σ进行卷积,这里假定信号最高值为1。黑色的信号线对应上图中的圆的截面,红色线对应Laplacian高斯算子截面。
我们还是用面积角度来思考(实际应该是通过体积,但这里由于其是圆形,中心对称,用面积考虑也成立)。信号为0时的卷积结果为0,信号为1的卷积结果由laplacian算子与x轴之间的面积决定,从图中各部分的面积可以很容易地看出,只有中间这种情况下,也就是拉普拉斯高斯算子的0值处和信号圆的边缘处刚好重合时,其卷积结果最大。也就是当前σ对应的信号尺度。
我们令\(\nabla^{2} G=0\),从上文中的表达式可以得到,其等效于\(\frac{x^{2}+y^{2}}{\sigma^{2}} - 2=0\),也就是
这是一个半径为\(\sqrt{2σ}\)的圆。
这个圆和上图中半径为r的信号圆刚好重叠时,卷积响应最大。因此我们可以得出σ和信号尺度的对应关系:
也就是用σ去和图像做卷积,检测尺度为 \(r = \sqrt{2} \sigma\)
如上右图,横轴为不同的σ,纵轴为和图片中心的卷积响应。可以看到σ在4-5之间时,响应最大,由此得到该处的信号半径r,绘制左图的圆。
尺度空间中提取特征
用不同的σ的拉普拉斯高斯算子和原图像做卷积
然后在图像每个点上都形成一个尺度空间:
在这个尺度空间中我们去找极大值,就可以知道这个当前点存在多大尺度的blob特征。
这里我们在实践时需要注意两点:
- 这里并不是找该点尺度空间中的最大值,而是每个尺度都和相邻的两个尺度比对,如果该尺度最大,则判定该点处有一个这个尺度的特征。这是因为在某个点上可能有同心圆的情况,那么这个点就可能在多个尺度下都有特征,在多个尺度下都有极值点存在。
- 在同一尺度下邻近的像素点可能检测到同样的特征信息,因此,在一个尺度下我们要对每个点在其局部范围内判断是否为最大值,如果该点是其领域范围内最大值,则将该点输出,否则不输出。
整体来讲,我们对一个像素点要考虑其上下两个尺度和其当前尺度周围的值,共27个(3X3X3局部区域内)进行比对,这里采用了非最大化抑制的思想。
我们把找到的特征输出到图像上,其中圆的大小代表不同的σ对应的不同尺度:
效率问题
上述的检测流程是没问题的,最大的问题是效率。试想,每个尺度下都需要和整个图像进行一次卷积操作;而且随着σ的增大,对应的卷积核也在增大。核越大,计算开销就越大。有没有什么方法可以减小所需的卷积核大小呢。这里采用了两种方法改进,我们下面介绍。
高斯差分(DoG)替代拉普拉斯高斯算子(LoG)
上文中我们一直采用的是拉普拉斯高斯算子对图像进行卷积操作,而高斯差分和它具有着非常相似的函数曲线,如图:
拉普拉斯高斯算子(Laplacian of Gaussian, LoG)的表达式:
高斯差分(Differece of Gaussian, DoG)的表达式:
他们之间有一个(k - 1)倍的近似关系:
也就是说原先的LoG构造的算子,可以用DoG构造的算子替换。
这种替换有什么好处呢?我们构造如图所示的尺度空间:
上图中左侧我们构造了一系列的Gaussian算子,得到了黄色的卷积响应,它们的σ以k为比例等比增加。通过相邻两项相减,就能得到对应的DoG,也就近似的得到了LoG算子,得到了蓝色的卷积响应,对应的算子的σ也同样是以k等比增加。
这里可以忽略 k - 1这个系数,是因为我们在不同尺度下只是比大小,k确定后这是一个常量,不会影响我们响应的相对大小。
这样,我们原先用Laplacian的不同σ去做卷积,现在我们使用Gaussian的不同σ去做卷积。而且好像还多了一个尺度,这种能提升效率吗?
这里的关键是:高斯卷积是可以分解的! 一个大的高斯核可以分解成两个小的高斯核依次进行卷积。
如果一个高斯核σ = a,另一个高斯核σ = b,那么他们两个对图像依次卷积,就等效于用一个\(σ = \sqrt{a^2+b^2}\)的高斯核对图像进行一次卷积。同理,一个大的高斯核,可以等效于两个小的高斯核依次卷积,它们的σ满足勾股定理。
这样的话,上图中左边那一堆高斯核就不需要对原图进行卷积了,而是可以通过对下一层的卷积结果再卷积一次得到。
比如,最下层以σ对原图做了一次卷积,而第二层不需要通过kσ和原图卷积得到,而是可以通过一个\(\sqrt{(kσ)^2-\sigma^2}\) 的高斯核对第一层的卷积结果再卷积一次得到;同样的第三层是通过对第二层的卷积结果,再来一个\(\sqrt{(k^2\sigma)^2-(k\sigma)^2}\)的卷积得到;以此类推。
这样的做法,虽然整体上多了一次卷积,但每次卷积的核大大减小了,其效率的提升是显著的。
Octave分组
通过高斯差分的替代,我们减小了每次卷积所需要的卷积核大小,然而当尺度进一步增大,层数再往上走,高斯核还是会不可避免地增加,有没有再减少高斯核大小的方法呢?
这里需要引入一个idea,如果当前图像上有一个半径为r的特征,用\(\sigma = r/\sqrt{2}\)检测出来了。那么如果我们把图像缩小一倍,特征的半径变成了r/2,那么我们检测这个特征需要的σ也就缩小到原来的一半σ/2,也就是所需要的高斯核也减小了一半(高斯核的半窗宽一般设置为3σ)。而最终检测出来的特征尺度我们再放大一倍,就能得到原图的特征了。
图像尺寸---特征尺寸 r --- σ --- 窗宽,这几项之间都是固定比例,它们的放大缩小都是同步的。
原图缩小了,所需要的高斯核就缩小了。这里面缩小原图会带来一些精度损失,但对于特征检测来说这点损失不算什么。
当尺度再往上走,我们虽然没法缩小高斯核了,但我们可以缩小图像,间接达到缩小高斯核的目的。
如图左侧的高斯卷积所示,下面5层作为一组,我们称为一个octave,这个里面我们是以原图为基础进行卷积;上面5层是另一个octave,在这里我们把原图像缩小一倍,如果我们沿用下方的5个σ对缩小的图分别进行卷积,这个操作就相当于用2倍的下方的σ对原图进行卷积,实际上检测的尺度为下方octave的2倍。
注意:这里用高斯核卷积的时候,σ在octave之间是不连续的,不是等比数列关系。
这样处理,我们继续往上扩展尺度的话,就再也不用增加高斯核的大小了,我们只需要每往上一个octave,把图像缩小一倍,然后沿用第一个octave中的高斯核,最后再把尺度放大回去即可。
我们再看上图右侧的实际特征输出,由于尺度空间中我们只对比当前尺度和上下两个尺度的响应大小,决定是否输出特征。在第一个octave中,我们实际输出的尺度空间为\(kσ\)和\(k^2σ\) ;同理,在上面这个octave中,输出的是\(2kσ\)和\(2k^2σ\)。
我们期望最终输出的每个点上的尺度空间是连续的,也就是这两个octave输出的总共四个尺度最好是连续的。那么如何设置octave的大小,k值设定为多少才能达到这种效果呢?
这是一个纯数学问题,我们直接使用数学结论:如果每个octave的输出尺度数量为s,则令
就能保证输出尺度空间的连续性。
我们把这个结论带到上图中验证一下,上图中每个octave输出2个尺度,也就是s = 2,那么\(k=\sqrt{2}\) 。我们带入可知这两个octave的输出尺度为\(\sqrt{2}\sigma\),\(2\sigma\),\(2\sqrt{2}\sigma\),\(4\sigma\) ,满足等比连续要求。
关键点定位
从上面的讨论中,我们已经得到了许多极值点,然而这些极值点却不是真正的极值点,这是因为我们取的尺度都是离散的尺度,是以k的倍数增长的一个离散尺度空间;并且,我们的像素也都是离散的像素。如果把上文中得到的响应函数表示的话,应该是 \(f(x, y, \sigma)\),一个关于像素位置和尺度的函数,我们前面的操作,只是得到了这个函数离散的一些数据点,这些点更有可能在真正的极值点附近而已。如下图所示:
泰勒展开
如何通过这些“伪极值点”得到真极值点呢?
我们的问题是:已知f的离散点,且知道哪些点可能在极值点附近,求f的真极值点。
我们可以通过泰勒展开来近似函数在某个点附近的信息,在harris角点的文章中就运用到泰勒展开去了解E(u, v)在(0, 0)附近的函数情况。这里,我们还是通过函数在伪极值点附近的泰勒展开去了解其附近的情况,从而求得真极值点:
在harris角点文章中我们用的是二元二阶泰勒展开,这里是三元二阶泰勒展开。
有限差分求导
但是这个里面由很多求导数的,函数我们都不知道如何求导数呢?
答案是用差分,在图像中我们求导数通常都是用差分来做,叫做有限差分求导。
如果一个图像的像素编号如下:
那么它的有限差分求导为:
多个尺度层间的编号如图:
差分求导:
差分求导:
极值点的确定
我们对上文的泰勒展开式进行简化:
令其导数为0:
得到:
这个就是当前的伪极值点到真极值点的一个偏移量。
用这个偏移量就能到真极值点了吗?还不行,因为我们用的泰勒展开毕竟是一个近似,有限差分也是真实导数的一个近似,所以只能说我们应该离真实极值点近了一些。实际算法实现过程中这是一个迭代的过程。不断地逼近真实极值点。在迭代时有如下几种情况:
- 当位移量三个分量都小于0.5时,说明位移量已经很小了,可以认为已经收敛。得到的就是精确极值点。
- 当迭代次数超限,位移量的三个分量还不是全都小于0.5,则认为不能收敛,该点不在极值点附近,把该点丢弃掉。
- 如果收敛了,我们把极值点带回原函数得到的值超出原先极值点的值一定范围,则表示泰勒展开没有很好的拟合原函数,则该点也被丢弃。
低对比度点舍弃
在得到真极值点后,我们代入算出极值,如果|f(X)|<T/s,则认为该点对比度太低,需要舍弃。在SIFT作者的原文中,这里T给了一个经验值0.04,s是上文中提到的单个octave输出尺度个数。
去除边缘效应
由于高斯差分函数对于边缘也有很强的响应,从而产生噪音。因此上面的得到的真极值点中,有一些是边缘点,而边缘点我们是不期望得到的。我们更希望得到角点。因此我们又要利用Harris角点检测的一些特性。我们把该点处的M矩阵求解出来,这个矩阵也叫海森矩阵(Hessian Matrix),用这个矩阵,我们就能判定哪些是边缘点了。我们把这些点也去掉。
这里放一张Harris角点检测文章中的判断图,R是由M矩阵得到的响应值:
经过层层筛选,我们终于把不同尺度的blob特征检测出来了,接下来的问题是,我们如何去描述这些特征,使得同一特征在不同的图像中具有相同的描述,继而完成后续的匹配任务。
我们需要首先解决一些变换问题。
尺度缩放问题
在这两张图中,我们通过前文的特征检测已经成功地将不同尺度下的特征检测出来了,然而在匹配时,毕竟它们的尺度不一致,我们没办法直接进行匹配,所以通常的做法是将它们缩放到大小一致的圆形区域内,也就是针对缩放进行归一化处理。这个是最好处理的变换问题。
仿射变换问题
如下图所示,当相机视角发生变化时,图像中的物体可能会出现拉伸变形,这种变换可以用仿射变换表示。
当图像发生仿射变换时,特征也发生了仿射变换,也被拉伸了,然而我们之前的算法中输出的特征区域是一个圆形的区域,如图:
这就导致在之后的匹配中,这两个特征区域无法匹配到一起。
而我们期望通过一种自适应的算法,形成这样的一个椭圆特征区域,也就是能刚好把特征框起来的区域:
这里我们用区域内信号在各方向上变化的快慢来解决这个问题。在原来的原型区域内,各个方向的强度变化是不一致的。而在新的椭圆形区域内,各个方向上强度变化大体是一致的。(两者相对而言)
如何描述一个区域内信号强度在各个方向上变化的快慢呢?
在Harris角点检测一文中我们已经对此进行了讨论。使用区域图像的二阶矩矩阵M来描述各个方向信号变化快慢问题。
当我们得到M矩阵后,对其进行对角化,R矩阵描述了最大强度变化和最小强度变化的方向。λ1和λ2分别描述了这两个方向上强度变化的大小。λ越大,等值面椭圆的轴越短,表示这个方向强度变化越快。反之在这个方向变化越慢。
具体的说采用一种迭代的方法:
从这张图中,我们对圆中的像素求取它的M矩阵,并通过对角化得到它的个方向强度变化情况。直观判断一下,图中沿着大概50°这个方向强度变化大,140°这个方向强度变化小。M矩阵对应的椭圆应该是在50°方向上为短轴,140°方向上为长轴。此时,我们就把这个框住特征区域的圆沿着强度变化大的方向(50°)缩小一些。再重新计算缩小后的区域的M矩阵,新的M矩阵长短轴之间的差距理论上来讲应该是缩小了一些,由此继续迭代下去。
随着迭代的进行,框住特征区域的圆逐渐“椭圆化”,而M矩阵对应的椭圆逐渐趋近于圆。直到M矩阵对应的两个λ相等,表示强度变化在各个方向上都一样了。此时,框住特征区域的椭圆就是我们要找的。
同时,我们仍然需要像处理尺度缩放问题一样,将该椭圆再通过仿射变换转化为单位圆。这也是归一化的操作。
旋转变换问题
在完成尺度归一化后,我们会发现本应匹配的两个特征,它们的旋转角度不一致,这个角度也得进行统一。
这里我们的做法是计算特征区域内各个像素的梯度方向,统计一个梯度直方图,这里梯度的方向和长度计算可以参见边缘检测文章:
选择这个直方图中最大的那个方向,作为该区域的”主方向“,按照这个主方向,将特征区域旋转到统一的方向,比如x正向。
把所有的特征区域都按照它们各自的主方向旋转到统一的x正向,即可完成角度的统一。
SIFT描述子
上述的尺度缩放和仿射变换问题是为了把特征归一化到一个统一大小的圆形区域,旋转问题是处理的角度统一。现在我们可以对区域的特征进行真正的描述(向量化)了。
我们把特征区域分成16个小区域,在每一个小区域内建立梯度直方图,如果把360°分成8分,直方图可以向量化为8维向量,16个区域总计\(16*8=128\),也就是128维向量。这个向量就是我们的SIFT描述子。
可以看到,SIFT描述子本质上是对特征点邻域内的梯度信息进行统计和编码 ,最终组合成描述子向量,用于表征图像局部特征的方向和相对强度分布。并不是以信号的强度进行编码,所以图像的强度(如光线的明暗)并不会对其产生影响。