SIFT算法详解
大纲
- 引言
- 一、高斯金字塔
- 二、高斯差分金字塔
- 三、特征点处理
- 1.阈值化
- 2.非极大值抑制
- 3. 二阶泰勒修正
- 4.低对比度去除
- 5.边缘效应去除
- 四、特征点描述子
- 1. 确定特征点区域方向
- 2. 特征点区域描述子
- 总结
- 参考:
引言
SIFT算法是为了解决图片的匹配问题,想要从图像中提取一种对图像的大小和旋转变化保持鲁棒的特征,从而实现匹配。这一算法的灵感也十分的直观:人眼观测两张图片是否匹配时会注意到其中的典型区域(特征点部分),如果我们能够实现这一特征点区域提取过程,再对所提取到的区域进行描述就可以实现特征匹配了。于是问题就演变成了以下几个子问题:
以上几个问题在SIFT算法里都用了很有意思的trick,后续会一一介绍。
一、高斯金字塔
引入高斯金字塔的目的在引言中已经介绍过了–解决图片缩放及尺度变化下特征提取的问题,高斯金字塔的Intuition有两个:1. 人看物体时近大远小,可以对图片下采样实现(金字塔->组);2. 人看物体时近处清晰,远处模糊,可以对图像高斯平滑实现(高斯->层);具体的推导可以参见我的另一篇博客:
Opencv学习笔记(六)图像金字塔
在SIFT里,高斯金字塔的层数和组数有着如下设定:
组数:O=[log2min(M,N)]−3O=[log_2min(M,N)]-3O=[log2min(M,N)]−3
层数:S=n+3S=n+3S=n+3
组数的设定是来自于提出SIFT算法的原始论文给出的经验值,理论上来说知要O≤[log2min(M,N)]O\leq[log_2min(M,N)]O≤[log2min(M,N)]即可,层数的设定则是有着理论依据的,这里的nnn是我们想要提取特征点的图片层数,由于提取出高斯金字塔后需要计算层间差分以获得高斯差分金字塔(DOG, Difference of Gausssian),高斯金字塔层数需要比DOG层数多1,而计算特征值时要求在尺度层面,即上下相邻层间计算,则DOG层数要比特征层数多2,则要求S=n+3S=n+3S=n+3。
附:在SIFT中对于高斯金字塔有几个需要补充的知识点。
二、高斯差分金字塔
通过高斯金字塔,我们获取了不同尺度的图片,接下来的问题是如何获取高频区域呢,一个很简单的思路就是按照边缘检测的算法使用差分滤波器如拉普拉斯滤波器、sobel滤波器在图片上滑动找到灰度值变化剧烈的区域。而经前人研究,归一化的高斯拉普拉斯算子的极大值极小值相较于其他特征提取函数可以获得最稳定的图像特征,因此我们打算使用归一化的高斯拉普拉斯算子在多尺度图片上提取特征,但是使用这种方式提取复杂度会很高,又尺度归一化高斯拉普拉斯算子和DOG函数有着如下关系:
G(x,y,kσ)−G(x,y,σ)≈(k−1)σ2∇2GG(x,y,k\sigma)-G(x,y,\sigma)\approx(k-1)\sigma^2 \nabla^2 GG(x,y,kσ)−G(x,y,σ)≈(k−1)σ2∇2G
证明如下:
忽略高斯函数系数:
G(x,y,σ)=1σ2exp(−x2+y22σ2)∂G∂x=−xσ4exp(−x2+y22σ2)∂2G∂x2=−σ2+x2σ6exp(−x2+y22σ2)∇2G(x,y)=∂2G∂x2+∂2G∂y2=−2σ2+x2+y2σ6exp(−x2+y22σ2)∂G∂σ=−2σ2+x2+y2σ5exp(−x2+y22σ2)⇒σ∇2G=∂G∂σ\begin{aligned} &G(x,y,\sigma)=\frac{1}{\sigma^2}exp(-\frac{x^2+y^2}{2\sigma^2})\\ &\frac{\partial{G}}{\partial x}=-\frac{x}{\sigma^4}exp(-\frac{x^2+y^2}{2\sigma^2})\\ &\frac{\partial^2{G}}{\partial x^2}=\frac{-\sigma^2+x^2}{\sigma^6}exp(-\frac{x^2+y^2}{2\sigma^2})\\ &\nabla^2 G(x,y)=\frac{\partial^2{G}}{\partial x^2}+\frac{\partial^2{G}}{\partial y^2}\\ &\qquad \qquad \ =\frac{-2\sigma^2+x^2+y^2}{\sigma^6}exp(-\frac{x^2+y^2}{2\sigma^2})\\ &\frac{\partial{G}}{\partial \sigma}=\frac{-2\sigma^2+x^2+y^2}{\sigma^5}exp(-\frac{x^2+y^2}{2\sigma^2})\\ &\Rightarrow\\ &\sigma\nabla^2G=\frac{\partial G}{\partial \sigma} \end{aligned}G(x,y,σ)=σ21exp(−2σ2x2+y2)∂x∂G=−σ4xexp(−2σ2x2+y2)∂x2∂2G=σ6−σ2+x2exp(−2σ2x2+y2)∇2G(x,y)=∂x2∂2G+∂y2∂2G =σ6−2σ2+x2+y2exp(−2σ2x2+y2)∂σ∂G=σ5−2σ2+x2+y2exp(−2σ2x2+y2)⇒σ∇2G=∂σ∂G
又对于差分高斯金字塔有:
DOG=G(x,y,kσ)−G(x,y,σ)(k−1)σ≈∂G∂σDOG≈(k−1)σ2∇2G\begin{aligned} &DOG=\frac{G(x,y,k\sigma)-G(x,y,\sigma)}{(k-1)\sigma}\approx \frac{\partial G}{\partial \sigma}\\ & DOG \approx(k-1)\sigma^2\nabla^2G \end {aligned}DOG=(k−1)σG(x,y,kσ)−G(x,y,σ)≈∂σ∂GDOG≈(k−1)σ2∇2G
由此我们不再需要在高斯金子塔上进行卷积操作,只需要计算在尺度层面σ\sigmaσ进行层间差分得到高斯差分金字塔(DOG),即完成了特征提取的操作。
三、特征点处理
得到DOG后,我们理论上来说就已经获得了特征值,但我们需要对特征值进行一些处理,就像我们在cannny边缘检测中使用的非极大值抑制等思路,去掉没那么特征的特征值。
1.阈值化
简单的阈值化去除掉变换没有那么剧烈的点,这些点就有可能是噪声引起的,算是在高斯拉普拉斯之外又加了一层去噪措施。
val={valabs(val)>0.5Tn0otherwiseval=\begin{cases} val & abs(val)>0.5\frac{T}{n}\\ 0 & otherwise \end{cases}val={val0abs(val)>0.5nTotherwise
这里的TTT为经验值0.04,nnn为之前提过的提取特征点数目。
2.非极大值抑制
非极大值一致的思路就和在其他算法中的一样,我们要求选取出的特征值应该是其领域范围内的极值,不同之处在于其他算法要求的只是图像这一二维平面内的极值,而这里要求还要在尺度σ\sigmaσ这一层面上也是极值,这也承接了前文DOG层数要比提取特征所用层数多2的设定。
3. 二阶泰勒修正
由于我们的图片在x,y,σx,y,\sigmax,y,σ方向上都只能取到离散值,即使通过前两个步骤取到了一些特征点,但这些特征点都不够精确,我们需要引入二阶泰勒函数对齐进行修正,使得特征点可以出现在亚像素(亚尺度)位置。
f(X)=f(X0)+∂fT∂XX^+12XT^∂2f∂X2X^f(X)=f(X_0)+\frac{\partial f^T}{\partial X}\hat{X}+\frac{1}{2}\hat{X^T}\frac{\partial^2 f}{\partial X^2}\hat{X}f(X)=f(X0)+∂X∂fTX^+21XT^∂X2∂2fX^
上式给出了特征点X0X_0X0附近函数的近似值f(x)f(x)f(x),其中X^=X−X0\hat{X}=X-X_0X^=X−X0,对该式求取一阶导数零点,即可得函数实际极值点位置与X0X_0X0的距离,从而对离散特征点修正到亚尺度处。
f′(X)=∂fT∂X+∂f2∂X2X^f'(X)=\frac{\partial f^T}{\partial X}+\frac{\partial f^2}{\partial X^2}\hat{X}f′(X)=∂X∂fT+∂X2∂f2X^
X^ex=−∂2f−1∂X2∂f∂X\hat{X}_{ex}=-\frac{\partial^2 f^{-1}}{\partial X^2}\frac{\partial f}{\partial X}X^ex=−∂X2∂2f−1∂X∂f
带入前式可求取新的特征点的灰度值:
f(X′)=f(X0)+12∂fT∂X′X^exf(X')=f(X_0)+\frac{1}{2}\frac{\partial f^T}{\partial X'}\hat{X}_{ex}f(X′)=f(X0)+21∂X′∂fTX^ex
需要注意的是,上式是一个不断迭代的过程,也即根据当前特征点二阶泰勒求取新的特征点这一过程会不断重复直到满足终止条件,如X^\hat{X}X^过小。另需注意当所得解超出离散极值点一定范围时需要舍去,因为二阶泰勒拟合只在其附近有效。
附:一个挺有意思的点在于这里的迭代目的和梯度下降法里的迭代并不是同一个目的,梯度下降里我们求取的是函数在某点的梯度,因为我们很难对复杂非线性的神经网络求除其一阶导数零点,即使求取到了也可能落入局部极值;而这里我们是直接求到了函数的一阶零点,迭代的目的是因为函数本身是一种近似,泰勒展开只取到了第二项,需要不断对原函数进行近似。这两种迭代的目的并不相同。
4.低对比度去除
目的和之间的阈值化类似,同样是去除掉没那么剧烈变化的特征点。要求:
∣f(x)∣≥Tn|f(x)|\geq\frac{T}{n}∣f(x)∣≥nT
5.边缘效应去除
引言中提过我们想要提取的特征点为角点而非边缘,而前述一系列措施只能保证取到灰度值变换剧烈的点,而边缘点同样符合这一特征,因此我们将通过以下方式去除边缘点。
- 计算黑森矩阵H(x,y)=[Dxx(x,y)Dxy(x,y)Dyx(x,y)Dyy(x,y)]H(x,y)=\begin{bmatrix} D_{xx}(x,y)&D_{xy}(x,y)\\ D_{yx}(x,y)&D_{yy}(x,y) \end{bmatrix}H(x,y)=[Dxx(x,y)Dyx(x,y)Dxy(x,y)Dyy(x,y)]
- 若矩阵行列式Det(H)<0Det(H)<0Det(H)<0,舍去该特征点
- 若矩阵行列式和迹不满足:Tr(H)Det(H)<(γ0+1)2γ0\frac{Tr(H)}{Det(H)}<\frac{(\gamma_0+1)^2}{\gamma_0}Det(H)Tr(H)<γ0(γ0+1)2,舍去该特征点,γ0\gamma_0γ0为有实际意义的经验值,通常设定为10。
接下来解释以下这几个步骤的意义,角点和边缘点的区别在于边缘在图像中表现为一条线,垂直于线的方向频率高,沿着线的方向频率比较低;而角点则在多(大于2)个方向方向出现强高频分量。黑森矩阵实际上是函数的二阶偏导构成的矩阵,可以反应函数的曲率变化状况。又对于二次型矩阵有如下性质:
由性质1,2我们可以推导处当Det(H)<0Det(H)<0Det(H)<0时,特征点为非极值点,舍去;
由性质1,3我们可以推导出当Tr(H)Det(H)=(α+β)2α2β2\frac{Tr(H)}{Det(H)}=\frac{(\alpha+\beta)^2}{\alpha^2\beta^2}Det(H)Tr(H)=α2β2(α+β)2过小时,由两特征特征向量的比值γ\gammaγ构成的式子(γ+1)2γ\frac{(\gamma+1)^2}{\gamma}γ(γ+1)2同样较小,且对勾函数在γ>1\gamma>1γ>1时单增,我们可以根据Tr(H)Det(H)\frac{Tr(H)}{Det(H)}Det(H)Tr(H)大小判断特征向量的相对大小,该值过小,说明函数在该点不同方向上的变化非常不均匀,类似于边缘,舍去。
四、特征点描述子
通过前述步骤,我们已经获得了不同尺寸层面上的稳定特征点,接下来需要对其进行描述。
1. 确定特征点区域方向
引言中提到为了使特征点拥有旋转不变性,我们会将特征点区域统一旋转到特定方向,这一方向即为特征点区域的主方向,因此我们需要先确定主方向。
计算方式为:统计在离该特征点尺度σ∗\sigma^*σ∗最近的尺度层σoct上,\sigma_{oct}上,σoct上,以特征点为中心3σ′=3∗1.5σoct3\sigma'=3*1.5\sigma_{oct}3σ′=3∗1.5σoct范围内的像素的梯度幅值及方向,对于范围内的梯度幅值用1.5σ∗1.5\sigma^*1.5σ∗高斯核滤波以实现距离加权。得到一系列梯度幅值和方向对pairs={amp,ang}pairs=\{amp,ang\}pairs={amp,ang}后,将360°360°360°方向划分为多个bins,将包含相应angangang的pairspairspairs中的ampampamp值累加到对应的bins上,如果angangang在两个bins划分间则根据距离分配ampampamp。(此处和HOG算法类似,可自行查阅)。最终可以获得特征点区域内幅值-方向直方图,选取其中幅值最大的方向作为主方向。
附:这里的几个参数像是统计区域的半径,高斯加权的方差,有着不同的设置方式,但内在思路是一致的。此外如果在主方向之外出现了幅值达到了主方向幅值80%的其他方向,我们会将其作为辅方向,后续匹配时会出现两个位置、尺度相同但方向不同的两个特征点区域。
2. 特征点区域描述子
获得了特征点区域的主方向之后,我们就可以利用该值计算出有旋转不变性的描述子了。首先还是和前一步类似,统计该特征点所在尺度层面一定区域内的梯度幅值和方向,实施起来略有不同。
r=32σoct(d+1)2r=\frac{3\sqrt{2}\sigma_{oct}(d+1)}{2}r=232σoct(d+1),其中ddd为我们在一个维度上划分的子块数目,通常为4;
附:需要注意的是以上的计算都是发生在相应的高斯金字塔尺度层面上,而非原图或者DOG上,此外由于旋转产生的像素值丢失通过插值算法解决,可以自行了解。如果想在原图上可视化SIFT特征,需要将我们获得的稳定特征点坐标变换回原始的图像尺寸上,简单的乘以原下采样倍数即可。
总结
至此SIFT算法就讲解完毕了,匹配的部分根据提的特征采用其他的聚类算法即可,总的来说这个算法还是有一定难度,本文也只是针对其他博客没有提到的细节引入了数学推导,使得整个思考过程更加连贯,更加详细的数学证明如有限差分法还请移步参考。
参考:
SIFT原理部分:
SIFT算法详解
SIFT算法原理详解
SIFT算法
SIFT(尺度不变特征变换)
LOG和DOG关系部分:
LOG算子
黑森矩阵意义部分:
Hessian矩阵以及在图像中的应用
Hessian 矩阵的特征值有什么含义?
Hessian矩阵与多元函数极值
Hessian矩阵(黑塞矩阵)
为何矩阵特征值乘积等于矩阵行列式值?
高斯相乘部分:
两个高斯函数的卷积仍为一高斯函数
总结
- 上一篇: python查询斐波那契数列通项公式_分
- 下一篇: mxm智能教育机器人无法智能对话_关于智