SIFT, Scale Invariant Feature Transform,尺度不变特征变换,由加拿大教授David G.Lowe提出的. SIFT特征对旋转、尺度缩放、亮度变化等保持不变性,是一种非常稳定的局部特征.
注:该文是学习记录,比较杂七杂八,比较多的复制粘贴,这里提前把相关博文备注如下:
[1] - SIFT特征详解-2015.10.03
[2] - SIFT - 百度百科
[3] - SIFT算法详解-2012.04.28
[4] - SIFT特征提取分析 - 2012.06.06
[5] - SIFT特征点提取 - 2018.01.15
[6] - sift特征提取算法 - 2015.05.23
[7] - 【动手学计算机视觉】第七讲:传统目标检测之SIFT特征 - 2019.10.30
[8] - python利用opencv实现SIFT特征提取与匹配 - 2020.11.05
[9] - 基于SIFT特征的图像检索 vs CNN - 2018.06.12
[10] - Implementing sift in python a complete guide-2019
[11] - 超强大的SIFT图像匹配技术详细指南(附Python代码)- 2019.11.20
[12] - SIFT(尺度不变特征变换) - bilibili
1. SIFT 特点
- 图像的局部特征,对旋转、尺度缩放、亮度变化保持不变,对视角变化、仿射变换、噪声也保持一定程度的稳定性.
- 独特性好,信息量丰富,适用于海量特征库进行快速、准确的匹配.
- 多量性,即使是很少几个物体也可以产生大量的SIFT特征
- 高速性,经优化的SIFT匹配算法甚至可以达到实时性
- 扩招性,可以很方便的与其他的特征向量进行联合.
SIFT算法可以解决的问题:
目标的自身状态、场景所处的环境和成像器材的成像特性等因素影响图像配准、目标识别跟踪的性能. 而SIFT算法在一定程度上可解决:
- 目标的旋转、缩放、平移(RST)
- 图像仿射/投影变换(视点viewpoint)
- 光照影响(illumination)
- 目标遮挡(occlusion)
- 杂物场景(clutter)
- 噪声
SIFT算法的实质是在不同的尺度空间上查找关键点(特征点),并计算出关键点的方向. SIFT所查找到的关键点是一些十分突出,不会因光照,仿射变换和噪音等因素而变化的点,如角点、边缘点、暗区的亮点及亮区的暗点等.
2. SIFT 特征检测步骤
主要分四步:
尺度空间的极值检测
搜索所有尺度空间上的图像,通过高斯微分函数来识别潜在的对尺度和选择不变的兴趣点.
特征点定位
在每个候选的位置上,通过一个拟合精细模型来确定位置尺度,关键点的选取依据他们的稳定程度.
特征方向赋值
基于图像局部的梯度方向,分配给每个关键点位置一个或多个方向,后续的所有操作都是对于关键点的方向、尺度和位置进行变换,从而提供这些特征的不变性.
特征点描述
在每个特征点周围的邻域内,在选定的尺度上测量图像的局部梯度,这些梯度被变换成一种表示,这种表示允许比较大的局部形状的变形和光照变换.
3. 尺度空间
在一定的范围内,无论物体是大还是小,人眼都可以分辨出来. 然而计算机要有相同的能力却不是那么的容易,在未知的场景中,计算机视觉并不能提供物体的尺度大小,其中的一种方法是把物体不同尺度下的图像都提供给机器,让机器能够对物体在不同的尺度下有一个统一的认知. 在建立统一认知的过程中,要考虑的就是在图像在不同的尺度下都存在的特征点.
3.1. 多分辨率图像金字塔
在早期图像的多尺度通常使用图像金字塔表示形式. 图像金字塔是同一图像在不同的分辨率下得到的一组结果,其生成过程一般包括两个步骤:
- 对原始图像进行平滑
对处理后的图像进行降采样(通常是水平、垂直方向的1/2)
降采样后得到一系列不断尺寸缩小的图像. 显然,一个传统的金字塔中,每一层的图像是其上一层图像长、高的各一半. 多分辨率的图像金字塔虽然生成简单,但其本质是降采样,图像的局部特征则难以保持,也就是无法保持特征的尺度不变形.
3.2. 高斯模糊
使用高斯模糊技术(Gaussian Blur) 可以用来降低图像中的噪点.
对于图像中的每个像素,高斯模糊技术会基于其相邻像素计算一个值. 以下是应用高斯模糊之前和之后的图像示例.
如图所示,纹理和次要细节将从图像中删除,并且仅保留诸如形状和边缘之类的相关信息.
高斯模糊成功地去除了图像中的噪点,强调了图像的重要特征.
4. 高斯尺度空间
可以通过图像的模糊程度来模拟人在距离物体由远到近时物体在视网膜上成像过程,距离物体越近其尺寸越大图像也越模糊,这就是高斯尺度空间,使用不同的参数模糊图像(分辨率不变),是尺度空间的另一种表现形式.
图像和高斯函数进行卷积运算能够对图像进行模糊,使用不同的“高斯核”可得到不同模糊程度的图像. 一副图像其高斯尺度空间可由其和不同的高斯卷积得到:
$$ L(x, y, \sigma) = G(x, y, \sigma) * I(x, y) $$
其中,$x$ 和 $y$ 是图像坐标,$G(x, y, \sigma)$ 是高斯函数:
$$ G(x, y, \sigma) = \frac{1}{2 \pi \sigma^2} e ^{\frac{x^2 + y^2}{2 \sigma^2}} $$
$\sigma$ 称为尺度空间因子,它是高斯正态分布的标准差,反映了图像被模糊的程度,其值越大图像越模糊,对应的尺度也就越大.
$L(x, y, \sigma)$ 代表着图像的高斯尺度空间.
构建尺度空间的目的是为了检测出在不同的尺度下都存在的特征点,而检测特征点较好的算子是 $\Delta ^2 G$(高斯拉普拉斯, LoG),
$$ \Delta ^2 = \frac{\partial ^2}{\partial x^2} + \frac{\partial ^2}{\partial y^2} $$
使用LoG虽然能较好的检测到图像中的特征点,但是其运算量过大,通常可使用DoG(差分高斯,Difference of Gaussina) 来近似计算 LoG.
设 k 为相邻两个高斯尺度空间的比例因子,则 DoG 定义为:
$$ D(x, y, \sigma) = [G(x, y, k \sigma) - G(x, y, \sigma)] * I(x, y) $$
$$ D(x, y, \sigma) = L(x, y, k \sigma) - L(x, y, \sigma) $$
从上式可以知道,将相邻的两个高斯空间的图像相减就得到了DoG的响应图像.
为了得到DoG图像,先要构建高斯尺度空间,而高斯的尺度空间可以在图像金字塔降采样的基础上加上高斯滤波得到,也就是对图像金字塔的每层图像使用不同的参数 $\sigma$ 进行高斯模糊,使每层金字塔有多张高斯模糊过的图像,这样一组图像就是一个 Octave.
降采样时,金字塔上边一组图像的第一张是由其下面一组图像倒数第三张降采样得到.
尺度空间是从单个图像生成的具有不同比例的图像的集合.
易知,高斯金字塔有多组,每组又有多层. 一组中的多个层之间的尺度是不一样的(也就是使用的高斯参数 $\sigma$ 是不同的),相邻的两层之间的尺度相差一个比例因子 k. 如果每组有 S 层,则$k = 2^{\frac{1}{S}}$. 上一组图像的最底层图像是由下一组中尺度为 $2 \sigma$ 的图像进行因子为 2 的降采样得到的(高斯金字塔)
高斯金字塔的组数一般是:
$$ o = [log_2 min(m, n)] - a $$
其中,o 表示金字塔的层数,m、n 分别为图像的行数和列数. 减去的系数 $a$ 可以在 $[0, log_2(m, n)]$ 区间内任意取值,和具体需要的金字塔的顶层图像的大小有关.
高斯模糊参数 $\sigma$ (尺度空间) 可由下面关系式得到:
$$ \sigma(o, s) = \sigma_0 \cdot 2 ^{\frac{o+s}{S}} $$
其中,o 为所在的组,s 为所在的层,$\sigma_0$ 为初始的尺度,S 为每组的层数.
在 Lowe 的算法实现中,$\sigma_0 = 1.6, o_{min} = -1, S=3$,$o_{min}=-1$ 就是首先讲原图像的长和宽各扩展一倍.
从上面可知,同一组内相邻层的图像尺度关系为:
$$ \sigma _{s+1} = k \cdot \sigma_s = 2^{\frac{1}{S}} \cdot \sigma_s $$
相邻组之间的尺度关系为:
$$ \sigma_{o+1} = 2 \sigma_{o} $$
示例:假设原始图像尺寸为 275x183,缩放图像尺寸为 138x93,对于该两个图像,将创建两个模糊图像,
然而,要缩放多少次呢,而且每个缩放图像需要创建多少个后续的模糊图像呢?理想的数量为 4 个,并且对于每个 Octave,模糊图像的数目应该为 5 个.
4.1. 高斯金字塔构建示例
假设图像 I 为 512x512 图像,构建高斯金字塔步骤为:(从 0 开始计数,倒立的金字塔)
[1] - 金字塔的组数,$log_2(512)=9$,减去因子 3,则构建的金字塔的组数为 6. 取每层的层数为 3.
[2] - 构建第 0 组,将图像宽和高增加一倍,变成 1024x1024($I_0$)
第 0 层 $I_0 * G(x, y, \sigma_0)$
第 1 层 $I_1 * G(x, y, k \sigma_0)$
第 2 层 $I_0 * G(x, y, k^2 \sigma_0)$
[3] - 构建第 1 组,对 $I_0$ 降采样,变成 512x512 ($I_1$)
第 0 层 $I_1 * G(x, y, 2 \sigma_0)$
第 1 层 $I_1 * G(x, y, 2k \sigma_0)$
第 2 层 $I_1 * G(x, y, 2k^2 \sigma_0)$
[4] - 类似地,依次进行
[5] - 构建第 o 组,第 s 层 $I_o * G(x, y, 2^o k^s \sigma_0)$
高斯金字塔构建成功后,将每一组相邻的两层想减,即可得到 DoG 金字塔.
相关源码:
From: SIFT特征点提取 - 2018.01.15
for (o = 0; o < octvs; o++)//金字塔组数为octvs,
for (i = 0; i < intvls + 3; i++)//每一组有intvls + 3 层,intvls一般为3
{
if (o == 0 && i == 0)//如果是第一组第1层
gauss_pyr[o][i] = cvCloneImage(base);//base 为原始灰度图像经过升采样或降采样得到的图像
/* base of new octvave is halved image from end of previous octave */
else if (i == 0)//建立非第一组的第1层
gauss_pyr[o][i] = downsample(gauss_pyr[o - 1][intvls]);//降采样图像
/* blur the current octave's last image to create the next one */
else//建立非第一组的非第1层
{
gauss_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i - 1]),IPL_DEPTH_32F, 1);
cvSmooth(gauss_pyr[o][i - 1], gauss_pyr[o][i],CV_GAUSSIAN, 0, 0, sig[i], sig[i]);// sig[i]为模糊系数
}//cvSmooth 为平滑处理函数,也即模糊处理. CV_GAUSSIAN 为选用高斯函数对图像模糊
return gauss_pyr;//返回建好的金字塔
4.2. 高斯差分金字塔构建
差分金字塔的是在高斯金字塔的基础上操作的,其建立过程是:在高斯金子塔中的每组中相邻两层相减(下一层减上一层)就生成高斯差分金字塔.
如图:
示例: 以第一个 Octave 为例,图像都具有相同比例,通过在前一个图像上应用高斯模糊,来创建每个后续图像,如图:
其源码实现如:
for (o = 0; o < octvs; o++)//octvs为高斯金字塔组数
for (i = 0; i < intvls + 2; i++)//因为相减,故高斯金字塔中每组有(intvls + 2)层图像
{
dog_pyr[o][i] = cvCreateImage(cvGetSize(gauss_pyr[o][i]),IPL_DEPTH_32F, 1);
cvSub(gauss_pyr[o][i + 1], gauss_pyr[o][i], dog_pyr[o][i], NULL);//cvSub为opencv内置相减函数
}
return dog_pyr;//返回高斯差分金字塔
5. DoG 空间极值点检测
为了寻找尺度空间的极值点,每个像素点要和其图像域(同一尺度空间)和尺度域(相邻的尺度空间)的所有相邻点进行比较,当其大于(或者小于)所有相邻点时,改点就是极值点.
极值点(关键点)是由DOG空间的局部极值点组成的,关键点的初步探查是通过同一组内各DoG相邻两层图像之间比较完成的. 为了寻找DoG函数的极值点,每一个像素点要和它所有的相邻点比较,看其是否比它的图像域和尺度域的相邻点大或者小. 如图下图所示,中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的 9×2 个点共26个点比较,以确保在尺度空间和二维图像空间都检测到极值点.
一个点如果在DOG尺度空间本层以及上下两层的26个邻域中是最大或最小值时,就认为该点是图像在该尺度下的一个特征点. 下图中将叉号点要比较的26个点都标为了绿色.
因为首尾两张图都只有一个相邻图,无法进行极值比较(因为极值定义中,极值点是连续的). 为了满足尺度变化连续性,要在Octave每一组的顶层使用高斯模糊再生成3幅图像,则高斯金字塔有S+3幅图像,而DOG金字塔有S+2幅图像.
5.1. 尺度变化的连续性
假设 S=3,也就是每组有 3 层,则 $k = 2 ^{\frac{1}{S}} = 2 ^{\frac{1}{3}}$,也就是高斯金字塔每组有 (S-1)3 层图像,DoG金字塔每组有 (S-2)2 层图像.
在 DoG 金字塔的第一组有两层尺度分别为 $\sigma, k\sigma$,第二组有两层的尺度分别为 $2\sigma, 2k\sigma$,由于只有两项是无法比较取得极值的(只有左右两边均有值才能有极值).
由于无法取得极值,那么就需要继续对每组的图像进行高斯模糊,使得尺度行形成 $\sigma, k\sigma, k^2\sigma, k^3 \sigma, k^4 \sigma$,这样就可以选择中间的三项 $k\sigma, k^2\sigma, k^3 \sigma$.
对应的下一组由上一组降采样得到的三项是 $2k\sigma, 2k^2\sigma, 2k^3 \sigma$,其首项 $2k \sigma = 2 \cdot 2^{\frac{1}{3}}\sigma=2^{\frac{4}{3}}\sigma$,刚好与上一组的最后一项 $k^3 \sigma = 2^{\frac{3}{3}} \sigma$ 的尺度连续起来.
5.2. DoG极值点检测实现
参考实现如:
if (val > 0)//极大值检测
{
for (i = -1; i <= 1; i++)
for (j = -1; j <= 1; j++)
for (k = -1; k <= 1; k++)
if (val < pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))//pixval32f为提取图像像素位置上的灰度值
return 0;}
else /* check for minimum */
{
for (i = -1; i <= 1; i++)
for (j = -1; j <= 1; j++)
for (k = -1; k <= 1; k++)
if (val > pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))//r c为图像的行数和列数,dog_pyr为高斯差分图
return 0;}
6. 删除不好的极值点(特征点)
通过比较检测得到的DoG的局部极值点是在离散的空间搜索得到的,由于离散空间是对连续空间采样得到的结果,因此在离散空间找到的极值点不一定是真正意义上的极值点,如图,
因此要设法将不满足条件的点剔除掉. 可以通过尺度空间DoG函数进行曲线拟合寻找极值点,这一步的本质是去掉DoG局部曲率非常不对称的点.
要剔除掉的不符合要求的点主要有两种:
- 低对比度的特征点
- 不稳定的边缘响应点
6.1. 剔除低对比度的特征点
候选特征点$x$,其偏移量定义为 $\Delta x$,其对比度为 $D(x)$的绝对值$|D(x)|$,对$D(x)$ 应用二阶泰勒展开式:
$$ D(x) = D + \frac{\partial D^T}{\partial x} \Delta x + \frac{1}{2} \Delta x^T \frac{\partial ^2 D }{\partial x^2} \Delta x $$
由于 $x$ 是 $D(x)$ 的极值点,所以对上式求导,并令其为 0,得到:
$$ \Delta = - \frac{\partial ^2 D^{-1}}{\partial x^2} \frac{\partial D(x)}{\partial x} $$
然后,再把求得的 $\Delta x$ 代入到 $D(x)$ 的泰勒展开式中,有:
$$ D(\hat{x}) = D + \frac{1}{2} \frac{\partial D^T}{\partial x} \hat{x} $$
设对比度的阈值为 T,若 $|D(\hat{x})| \ge T$,则该特征点保留,否则剔除.
注:
任意极值点在其 $(x_0, y_0, \sigma_0)$ 处的二阶泰勒展开式如:
6.2. 剔除不稳定的边缘响应点
在边缘梯度的方向上主曲率值比较大,而沿着边缘方向则主曲率值较小. 候选特征点的DoG函数D(x) 的主曲率与 2x2 Hessian 矩阵 H 的特征值成正比.
$$ H= \left[ \begin{matrix} D_{xx} & D_{yx} \\ D_{xy} & D_{yy} \end{matrix} \right] $$
其中, $D_{xx}, D_{xy}, D_{yx}$ 是候选点邻域对应位置的差分求得的.
为了避免求具体的值,可以使用 H 特征值的比例. 设 $\alpha = \lambda _{max}$ 为 H 的最大特征值,$\beta = \lambda _{min}$ 为 H 的最小特征值,则:
$$ Tr(H) = D_{xx} + D_{yy} = \alpha + \beta $$
$$ Det(H) = D_{xx} + D_{yy} - D_{xy}^2 = \alpha \cdot \beta $$
其中,$Tr(H)$ 为矩阵 H 的迹,$Det(H)$ 为矩阵 H 的行列式.
设 $\gamma = \frac{\alpha}{\beta}$ 表示最大特征值和最小特征值的比值,则:
$$ \frac{Tr(H)^2}{Det(H)} = \frac{(\alpha + \beta)^2}{\alpha \beta} = \frac{(\gamma \beta + \beta)^2}{\gamma \beta^2} = \frac{(\gamma + 1)^2}{\gamma} $$
上式的结果与两个特征值的比例有关,和具体的大小无关,当两个特征值想等时其值最小,并且随着 $\gamma$ 的增大而增大. 因此为了检测主曲率是否在某个阈值下,只需检测:
$$ \frac{Tr(H)^2}{Det(H)} > \frac{(T_{\gamma} + 1)^2}{T_{\gamma}} $$
如果上式成立,则剔除该特征点,否则保留. (Lowe 论文中取 $T_{\gamma} =10$).
6.3. 泰勒插值的源码实现
while (i < SIFT_MAX_INTERP_STEPS)//SIFT_MAX_INTERP_STEPS=5为最大迭代次数,避免长时迭代
{
interp_step(dog_pyr, octv, intvl, r, c, &xi, &xr, &xc);// 泰勒展开拟合,xi,xr,xc依次为x、y、σ方向偏移量,
if (ABS(xi) < 0.5 && ABS(xr) < 0.5 && ABS(xc) < 0.5)//如果当前偏移量绝对值中的每个值均小于0.5,退出迭代
break;
c += cvRound(xc);//计算行坐标,cvRound 为四舍五入.
r += cvRound(xr);
intvl += cvRound(xi);
if (intvl < 1 ||//不在计算的图像层中
intvl > intvls ||//高斯差分每组的层数为intvls
c < SIFT_IMG_BORDER ||//靠近图像边缘5个像素的区域不做检测,SIFT_IMG_BORDER=5,
r < SIFT_IMG_BORDER ||
c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER ||//靠近图像边缘5个像素的区域不做检测
r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER)
{
return NULL;
}
i++;//迭代计数
}
#
static void interp_step(IplImage*** dog_pyr, int octv, int intvl, int r, int c,double* xi, double* xr, double* xc)
{
CvMat* dD, *H, *H_inv, X;
double x[3] = { 0 };
dD = deriv_3D(dog_pyr, octv, intvl, r, c);//一阶偏导数
H = hessian_3D(dog_pyr, octv, intvl, r, c);//Hessian 矩阵即二阶导数组成的矩阵
H_inv = cvCreateMat(3, 3, CV_64FC1);
cvInvert(H, H_inv, CV_SVD);//求Hessian矩阵的逆矩阵
cvInitMatHeader(&X, 3, 1, CV_64FC1, x, CV_AUTOSTEP);
cvGEMM(H_inv, dD, -1, NULL, 0, &X, 0); //cvGEMM为矩阵乘法,//第一个矩阵的系数;//H_inv、dD第一二个矩阵//-1矩阵前的常数//X结果矩阵
cvReleaseMat(&dD);
cvReleaseMat(&H);
cvReleaseMat(&H_inv);
*xi = x[2];
*xr = x[1];
*xc = x[0];
}
7. 求取特征点的主方向
经过上面的步骤已经找到了在不同尺度下都存在的特征点,为了实现图像旋转不变性,需要给特征点的方向进行赋值. 利用特征点邻域像素的梯度分布特性来确定其方向参数,再利用图像的梯度直方图求取关键点局部结构的稳定方向.
找到了特征点,也就可以得到该特征点的尺度$\sigma$,也就可以得到特征点所在的尺度图像
$$ L(x, y) = G(x, y, \sigma) * I(x, y) $$
计算以特征点为中心,以 $3 \times 1.5\sigma$ 为半径的区域图像的幅角和幅值,每个点 $L(x, y)$ 的梯度的模 $m(x, y)$ 以及方向 $\theta (x, y)$ 可通过如下公式求得:
$$ m(x, y) = \sqrt{[L(x+1, y) - L(x-1, y)]^2 + [L(x, y+1) - L(x, y-1)]^2} $$
$$ \theta(x, y) = \text{arctan} \frac{L(x, y+1) - L(x, y-1)}{L(x+1, y) - L(x-1, y)} $$
计算得到梯度方向后,就要使用直方图统计特征点邻域内像素对应的梯度方向和幅值. 梯度方向的直方图的横轴是梯度方向的角度(梯度方向的范围是0到360度,直方图每36度一个柱共10个柱,或者没45度一个柱共8个柱),纵轴是梯度方向对应梯度幅值的累加,在直方图的峰值就是特征点的主方向.
在完成关键点的梯度计算后,使用直方图统计领域内像素的梯度和方向. 梯度直方图将0~360度的方向范围分为36个柱(bins),其中每柱10度. 如图所示,直方图的峰值方向代表了关键点的主方向,(为简化,图中只画了八个方向的直方图).
在Lowe的论文还提到了使用高斯函数对直方图进行平滑以增强特征点近的邻域点对关键点方向的作用,并减少突变的影响. 为了得到更精确的方向,通常还可以对离散的梯度直方图进行插值拟合.
对梯度方向直方图的平滑处理:Opencv 采用的平滑公式为:
$$ H(i) = \frac{h(i-2) + h(i+2)}{16} + \frac{4 \times (h(i-1) + h(i+1))}{16} + \frac{6 \times h(i)}{16} $$
其中,$i \in [0, 35]$,$h$ 和 $H$ 分别表示平滑前和平滑后的直方图. 由于角度是循环的,即 $0^o = 360^o$,如果出现 $h(i)$,$i$ 超出了 [0, 35] 的范围,则可以通过圆周循环的方法找到它所对应的、在 $0^o = 360^o$ 之间的值,如 $h(-1) = h(35)$.
具体而言,关键点的方向可以由和主峰值最近的三个柱值通过抛物线插值得到.
梯度直方图抛物线插值:
在梯度直方图中,当存在一个相当于主峰值80%能量的柱值时,则可以将这个方向认为是该特征点辅助方向. 所以,一个特征点可能检测到多个方向(也可以理解为,一个特征点可能产生多个坐标、尺度相同,但是方向不同的特征点). Lowe在论文中指出: 15%的关键点具有多方向,而且这些点对匹配的稳定性很关键.
直方图峰值代表该关键点邻域内图像梯度的主方向,当存在另一个相当于主峰值 80% 能量的峰值时,则认为这个方向是该关键点的辅方向. 所以一个关键点可能检测得到多个方向,这可以增强匹配的鲁棒性.
得到特征点的主方向后,对于每个特征点可以得到三个信息$(x, y, \sigma, \theta)$,即位置、尺度和方向. 由此可以确定一个SIFT特征区域,一个SIFT特征区域由三个值表示,中心表示特征点位置,半径表示关键点的尺度,箭头表示主方向. 具有多个方向的关键点可以被复制成多份,然后将方向值分别赋给复制后的特征点,一个特征点就产生了多个坐标、尺度相等,但是方向不同的特征点.
8. 生成特征描述
通过以上的步骤已经找到了SIFT特征点位置、尺度和方向信息,下面就需要使用一组向量来描述关键点也就是生成特征点描述子,这个描述符不只包含特征点,也含有特征点周围对其有贡献的像素点. 描述子应具有较高的独立性,以保证匹配率.
特征描述符的生成大致有三个步骤:
- 校正旋转主方向,确保旋转不变性.
- 生成描述子,最终形成一个128维的特征向量
- 归一化处理,将特征向量长度进行归一化处理,进一步去除光照的影响.
为了保证特征矢量的旋转不变性,要以特征点为中心,在附近邻域内将坐标轴旋转 $\theta$ (特征点的主方向)角度,即将坐标轴旋转为特征点的主方向.
旋转后邻域内像素的新坐标为:
$$ \left[ \begin{matrix} x' \\ y' \end{matrix} \right] = \left[ \begin{matrix} cos \theta & -sin \theta \\ sin \theta & cos \theta \end{matrix} \right] \left[ \begin{matrix} x \\ y \end{matrix} \right] $$
旋转后以主方向为中心取 8x8 的窗口. 下图所示,左图的中央为当前关键点的位置,每个小格代表为关键点邻域所在尺度空间的一个像素,求取每个像素的梯度幅值与梯度方向,箭头方向代表该像素的梯度方向,长度代表梯度幅值,然后利用高斯窗口对其进行加权运算. 最后在每个 4×4 的小块上绘制8个方向的梯度直方图,计算每个梯度方向的累加值,即可形成一个种子点,如图所示. 每个特征点由4个种子点组成,每个种子点有8个方向的向量信息. 这种邻域方向性信息联合增强了算法的抗噪声能力,同时对于含有定位误差的特征匹配也提供了比较理性的容错性.
与求主方向不同,此时每个种子区域的梯度直方图在0-360之间划分为8个方向区间,每个区间为45度,即每个种子点有8个方向的梯度强度信息.
在实际的计算过程中,为了增强匹配的稳健性,Lowe建议:对每个关键点使用4×4共16个种子点来描述,这样一个关键点就可以产生128维(4x4x8)的SIFT特征向量.
通过对特征点周围的像素进行分块,计算块内梯度直方图,生成具有独特性的向量,这个向量是该区域图像信息的一种抽象,具有唯一性.
特征向量形成后,为了去除光照变化的影响,需要对它们进行归一化处理.
主要处理概括如图:
9. SIFT 特征匹配
在 SIFT 特征提取完成以后,即可生成图像的描述子.
假设两幅图像 A、B 的 SIFT 描述子分别为 k1*128 维、k2*128 ,可以将两图的各个尺度的描述子进行匹配,如果两个特征点匹配上 128 维,则表示其匹配上了.
当两幅图像的SIFT特征向量生成后,下一步我们采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量. 取图像 A 中的某个关键点,并找出其与图像 B 中欧式距离最近的前两个关键点,在这两个关键点中,如果最近的距离除以次近的距离少于某个比例阈值,则接受这一对匹配点. 降低这个比例阈值,SIFT匹配点数目会减少,但更加稳定.
为了排除因为图像遮挡和背景混乱而产生的无匹配关系的关键点, Lowe提出了比较最近邻距离与次近邻距离的方法, 距离比率ratio小于某个阈值的认为是正确匹配. 因为对于错误匹配,由于特征空间的高维性, 相似的距离可能有大量其他的错误匹配,从而它的ratio值比较高. Lowe推荐ratio的阈值为0.8. 但对大量任意存在尺度、旋转和亮度变化的两幅图片进行匹配,结果表明ratio取值在0. 4~0. 6之间最佳,小于0. 4的很少有匹配点,大于0. 6的则存在大量错误匹配点. 故,建议ratio的取值原则如下:
ratio=0. 4 对于准确度要求高的匹配;
ratio=0. 6 对于匹配点数目要求比较多的匹配;
ratio=0. 5 一般情况下.
也可按如下原则:当最近邻距离<200时 ratio=0. 6,反之ratio=0. 4. ratio的取值策略能排分错误匹配点.
10. OpenCV 实现
#!/usr/bin/python3
#!--*-- coding: utf-8 --*--
import cv2
#
imgfile1 = 'test1.jpg'
imgfile2 = 'test2.jpg'
img1 = cv2.imread(imgfile1)
img2 = cv2.imread(imgfile2)
#
sift = cv2.xfeatures2d.SIFT_create()
# 计算关键点和描述子
# kp为关键点keypoints, des为描述子descriptors
kp1, des1 = sift.detectAndCompute(img1, None)
kp2, des2 = sift.detectAndCompute(img2, None)
# 关键点可视化
img3 = cv2.drawKeypoints(img1, kp1, img1, color=(0, 255, 255))
img4 = cv2.drawKeypoints(img2, kp2, img2, color=(0, 255, 255))
# M1:匹配方法一
# 描述子匹配
index_params = dict(algorithm=1, trees=6)
search_params = dict(checks=50)
flann = cv2.FlannBasedMatcher(index_params, search_params)
# knn 计算描述子的匹配
matche = flann.knnMatch(des1, des2, k=2)
matchesMask = [[0, 0] for i in range(len(matche))]
# 匹配结果可视化
result = []
for m, n in matche:
if m.distance < 0.6 * n.distance:
result.append([m])
img5 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, matche, None, flags=2)
cv2.imshow("SIFTMatchResult", img5)
cv2.waitKey(0)
## M2: 匹配方法二
## K近邻算法求取在空间中距离最近的K个数据点,并将这些数据点归为一类
#matcher = cv2.BFMatcher()
#raw_matches = matcher.knnMatch(des1, des2, k = 2)
#good_matches = []
#for m1, m2 in raw_matches:
# # 如果最接近和次接近的比值大于一个既定的值,那么保留这个最接近的值,认为它和其匹配的点为good_match
# if m1.distance < ratio * m2.distance:
# good_matches.append([m1])
#img5 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, matche, None, flags=2)
11. SIFT 的缺点
SIFT在图像的不变特征提取方面拥有无与伦比的优势,但并不完美,仍然存在:
- 实时性不高.
- 有时特征点较少.
- 对边缘光滑的目标无法准确提取特征点,等.
如图所示,对模糊的图像和边缘平滑的图像,检测出的特征点过少,对圆更是无能为力. 近来不断有人改进,其中最著名的有SURF和CSIFT.
.1. 尺度变化的连续性
假设 S=3,也就是每组有 3 层,则 $k = 2 ^{\frac{1}{S}} = 2 ^{\frac{1}{3}}$,也就是高斯金字塔每组有 (S-1)3 层图像,DoG金字塔每组有 (S-2)2 层图像.
在 DoG 金字塔的第一组有两层尺度分别为 $\sigma, k\sigma$,第二组有两层的尺度分别为 $2\sigma, 2k\sigma$,由于只有两项是无法比较取得极值的(只有左右两边均有值才能有极值).
由于无法取得极值,那么就需要继续对每组的图像进行高斯模糊,使得尺度行形成 $\sigma, k\sigma, k^2\sigma, k^3 \sigma, k^4 \sigma$,这样就可以选择中间的三项 $k\sigma, k^2\sigma, k^3 \sigma$.
对应的下一组由上一组降采样得到的三项是 $2k\sigma, 2k^2\sigma, 2k^3 \sigma$,其首项 $2k \sigma = 2 \cdot 2^{\frac{1}{3}}\sigma=2^{\frac{4}{3}}\sigma$,刚好与上一组的最后一项 $k^3 \sigma = 2^{\frac{3}{3}} \sigma$ 的尺度连续起来.
5.2. DoG极值点检测实现
参考实现如:
if (val > 0)//极大值检测
{
for (i = -1; i <= 1; i++)
for (j = -1; j <= 1; j++)
for (k = -1; k <= 1; k++)
if (val < pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))//pixval32f为提取图像像素位置上的灰度值
return 0;}
else /* check for minimum */
{
for (i = -1; i <= 1; i++)
for (j = -1; j <= 1; j++)
for (k = -1; k <= 1; k++)
if (val > pixval32f(dog_pyr[octv][intvl + i], r + j, c + k))//r c为图像的行数和列数,dog_pyr为高斯差分图
return 0;}
6. 删除不好的极值点(特征点)
通过比较检测得到的DoG的局部极值点是在离散的空间搜索得到的,由于离散空间是对连续空间采样得到的结果,因此在离散空间找到的极值点不一定是真正意义上的极值点,如图,
因此要设法将不满足条件的点剔除掉. 可以通过尺度空间DoG函数进行曲线拟合寻找极值点,这一步的本质是去掉DoG局部曲率非常不对称的点.
要剔除掉的不符合要求的点主要有两种:
- 低对比度的特征点
- 不稳定的边缘响应点
6.1. 剔除低对比度的特征点
候选特征点$x$,其偏移量定义为 $\Delta x$,其对比度为 $D(x)$的绝对值$|D(x)|$,对$D(x)$ 应用二阶泰勒展开式:
$$ D(x) = D + \frac{\partial D^T}{\partial x} \Delta x + \frac{1}{2} \Delta x^T \frac{\partial ^2 D }{\partial x^2} \Delta x $$
由于 $x$ 是 $D(x)$ 的极值点,所以对上式求导,并令其为 0,得到:
$$ \Delta = - \frac{\partial ^2 D^{-1}}{\partial x^2} \frac{\partial D(x)}{\partial x} $$
然后,再把求得的 $\Delta x$ 代入到 $D(x)$ 的泰勒展开式中,有:
$$ D(\hat{x}) = D + \frac{1}{2} \frac{\partial D^T}{\partial x} \hat{x} $$
设对比度的阈值为 T,若 $|D(\hat{x})| \ge T$,则该特征点保留,否则剔除.
注:
任意极值点在其 $(x_0, y_0, \sigma_0)$ 处的二阶泰勒展开式如:
6.2. 剔除不稳定的边缘响应点
在边缘梯度的方向上主曲率值比较大,而沿着边缘方向则主曲率值较小. 候选特征点的DoG函数D(x) 的主曲率与 2x2 Hessian 矩阵 H 的特征值成正比.
$$ H= \left[ \begin{matrix} D_{xx} & D_{yx} \\ D_{xy} & D_{yy} \end{matrix} \right] $$
其中, $D_{xx}, D_{xy}, D_{yx}$ 是候选点邻域对应位置的差分求得的.
为了避免求具体的值,可以使用 H 特征值的比例. 设 $\alpha = \lambda _{max}$ 为 H 的最大特征值,$\beta = \lambda _{min}$ 为 H 的最小特征值,则:
$$ Tr(H) = D_{xx} + D_{yy} = \alpha + \beta $$
$$ Det(H) = D_{xx} + D_{yy} - D_{xy}^2 = \alpha \cdot \beta $$
其中,$Tr(H)$ 为矩阵 H 的迹,$Det(H)$ 为矩阵 H 的行列式.
设 $\gamma = \frac{\alpha}{\beta}$ 表示最大特征值和最小特征值的比值,则:
$$ \frac{Tr(H)^2}{Det(H)} = \frac{(\alpha + \beta)^2}{\alpha \beta} = \frac{(\gamma \beta + \beta)^2}{\gamma \beta^2} = \frac{(\gamma + 1)^2}{\gamma} $$
上式的结果与两个特征值的比例有关,和具体的大小无关,当两个特征值想等时其值最小,并且随着 $\gamma$ 的增大而增大. 因此为了检测主曲率是否在某个阈值下,只需检测:
$$ \frac{Tr(H)^2}{Det(H)} > \frac{(T_{\gamma} + 1)^2}{T_{\gamma}} $$
如果上式成立,则剔除该特征点,否则保留. (Lowe 论文中取 $T_{\gamma} =10$).
6.3. 泰勒插值的源码实现
while (i < SIFT_MAX_INTERP_STEPS)//SIFT_MAX_INTERP_STEPS=5为最大迭代次数,避免长时迭代
{
interp_step(dog_pyr, octv, intvl, r, c, &xi, &xr, &xc);// 泰勒展开拟合,xi,xr,xc依次为x、y、σ方向偏移量,
if (ABS(xi) < 0.5 && ABS(xr) < 0.5 && ABS(xc) < 0.5)//如果当前偏移量绝对值中的每个值均小于0.5,退出迭代
break;
c += cvRound(xc);//计算行坐标,cvRound 为四舍五入.
r += cvRound(xr);
intvl += cvRound(xi);
if (intvl < 1 ||//不在计算的图像层中
intvl > intvls ||//高斯差分每组的层数为intvls
c < SIFT_IMG_BORDER ||//靠近图像边缘5个像素的区域不做检测,SIFT_IMG_BORDER=5,
r < SIFT_IMG_BORDER ||
c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER ||//靠近图像边缘5个像素的区域不做检测
r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER)
{
return NULL;
}
i++;//迭代计数
}
#
static void interp_step(IplImage*** dog_pyr, int octv, int intvl, int r, int c,double* xi, double* xr, double* xc)
{
CvMat* dD, *H, *H_inv, X;
double x[3] = { 0 };
dD = deriv_3D(dog_pyr, octv, intvl, r, c);//一阶偏导数
H = hessian_3D(dog_pyr, octv, intvl, r, c);//Hessian 矩阵即二阶导数组成的矩阵
H_inv = cvCreateMat(3, 3, CV_64FC1);
cvInvert(H, H_inv, CV_SVD);//求Hessian矩阵的逆矩阵
cvInitMatHeader(&X, 3, 1, CV_64FC1, x, CV_AUTOSTEP);
cvGEMM(H_inv, dD, -1, NULL, 0, &X, 0); //cvGEMM为矩阵乘法,//第一个矩阵的系数;//H_inv、dD第一二个矩阵//-1矩阵前的常数//X结果矩阵
cvReleaseMat(&dD);
cvReleaseMat(&H);
cvReleaseMat(&H_inv);
*xi = x[2];
*xr = x[1];
*xc = x[0];
}
7. 求取特征点的主方向
经过上面的步骤已经找到了在不同尺度下都存在的特征点,为了实现图像旋转不变性,需要给特征点的方向进行赋值. 利用特征点邻域像素的梯度分布特性来确定其方向参数,再利用图像的梯度直方图求取关键点局部结构的稳定方向.
找到了特征点,也就可以得到该特征点的尺度$\sigma$,也就可以得到特征点所在的尺度图像
$$ L(x, y) = G(x, y, \sigma) * I(x, y) $$
计算以特征点为中心,以 $3 \times 1.5\sigma$ 为半径的区域图像的幅角和幅值,每个点 $L(x, y)$ 的梯度的模 $m(x, y)$ 以及方向 $\theta (x, y)$ 可通过如下公式求得:
$$ m(x, y) = \sqrt{[L(x+1, y) - L(x-1, y)]^2 + [L(x, y+1) - L(x, y-1)]^2} $$
$$ \theta(x, y) = \text{arctan} \frac{L(x, y+1) - L(x, y-1)}{L(x+1, y) - L(x-1, y)} $$
计算得到梯度方向后,就要使用直方图统计特征点邻域内像素对应的梯度方向和幅值. 梯度方向的直方图的横轴是梯度方向的角度(梯度方向的范围是0到360度,直方图每36度一个柱共10个柱,或者没45度一个柱共8个柱),纵轴是梯度方向对应梯度幅值的累加,在直方图的峰值就是特征点的主方向.
在完成关键点的梯度计算后,使用直方图统计领域内像素的梯度和方向. 梯度直方图将0~360度的方向范围分为36个柱(bins),其中每柱10度. 如图所示,直方图的峰值方向代表了关键点的主方向,(为简化,图中只画了八个方向的直方图).
在Lowe的论文还提到了使用高斯函数对直方图进行平滑以增强特征点近的邻域点对关键点方向的作用,并减少突变的影响. 为了得到更精确的方向,通常还可以对离散的梯度直方图进行插值拟合.
对梯度方向直方图的平滑处理:Opencv 采用的平滑公式为:
$$ H(i) = \frac{h(i-2) + h(i+2)}{16} + \frac{4 \times (h(i-1) + h(i+1))}{16} + \frac{6 \times h(i)}{16} $$
其中,$i \in [0, 35]$,$h$ 和 $H$ 分别表示平滑前和平滑后的直方图. 由于角度是循环的,即 $0^o = 360^o$,如果出现 $h(i)$,$i$ 超出了 [0, 35] 的范围,则可以通过圆周循环的方法找到它所对应的、在 $0^o = 360^o$ 之间的值,如 $h(-1) = h(35)$.
具体而言,关键点的方向可以由和主峰值最近的三个柱值通过抛物线插值得到.
梯度直方图抛物线插值:
在梯度直方图中,当存在一个相当于主峰值80%能量的柱值时,则可以将这个方向认为是该特征点辅助方向. 所以,一个特征点可能检测到多个方向(也可以理解为,一个特征点可能产生多个坐标、尺度相同,但是方向不同的特征点). Lowe在论文中指出: 15%的关键点具有多方向,而且这些点对匹配的稳定性很关键.
直方图峰值代表该关键点邻域内图像梯度的主方向,当存在另一个相当于主峰值 80% 能量的峰值时,则认为这个方向是该关键点的辅方向. 所以一个关键点可能检测得到多个方向,这可以增强匹配的鲁棒性.
得到特征点的主方向后,对于每个特征点可以得到三个信息$(x, y, \sigma, \theta)$,即位置、尺度和方向. 由此可以确定一个SIFT特征区域,一个SIFT特征区域由三个值表示,中心表示特征点位置,半径表示关键点的尺度,箭头表示主方向. 具有多个方向的关键点可以被复制成多份,然后将方向值分别赋给复制后的特征点,一个特征点就产生了多个坐标、尺度相等,但是方向不同的特征点.
8. 生成特征描述
通过以上的步骤已经找到了SIFT特征点位置、尺度和方向信息,下面就需要使用一组向量来描述关键点也就是生成特征点描述子,这个描述符不只包含特征点,也含有特征点周围对其有贡献的像素点. 描述子应具有较高的独立性,以保证匹配率.
特征描述符的生成大致有三个步骤:
- 校正旋转主方向,确保旋转不变性.
- 生成描述子,最终形成一个128维的特征向量
- 归一化处理,将特征向量长度进行归一化处理,进一步去除光照的影响.
为了保证特征矢量的旋转不变性,要以特征点为中心,在附近邻域内将坐标轴旋转 $\theta$ (特征点的主方向)角度,即将坐标轴旋转为特征点的主方向.
旋转后邻域内像素的新坐标为:
$$ \left[ \begin{matrix} x' \\ y' \end{matrix} \right] = \left[ \begin{matrix} cos \theta & -sin \theta \\ sin \theta & cos \theta \end{matrix} \right] \left[ \begin{matrix} x \\ y \end{matrix} \right] $$
旋转后以主方向为中心取 8x8 的窗口. 下图所示,左图的中央为当前关键点的位置,每个小格代表为关键点邻域所在尺度空间的一个像素,求取每个像素的梯度幅值与梯度方向,箭头方向代表该像素的梯度方向,长度代表梯度幅值,然后利用高斯窗口对其进行加权运算. 最后在每个 4×4 的小块上绘制8个方向的梯度直方图,计算每个梯度方向的累加值,即可形成一个种子点,如图所示. 每个特征点由4个种子点组成,每个种子点有8个方向的向量信息. 这种邻域方向性信息联合增强了算法的抗噪声能力,同时对于含有定位误差的特征匹配也提供了比较理性的容错性.
与求主方向不同,此时每个种子区域的梯度直方图在0-360之间划分为8个方向区间,每个区间为45度,即每个种子点有8个方向的梯度强度信息.
在实际的计算过程中,为了增强匹配的稳健性,Lowe建议:对每个关键点使用4×4共16个种子点来描述,这样一个关键点就可以产生128维(4x4x8)的SIFT特征向量.
通过对特征点周围的像素进行分块,计算块内梯度直方图,生成具有独特性的向量,这个向量是该区域图像信息的一种抽象,具有唯一性.
特征向量形成后,为了去除光照变化的影响,需要对它们进行归一化处理.
主要处理概括如图:
9. SIFT 特征匹配
在 SIFT 特征提取完成以后,即可生成图像的描述子.
假设两幅图像 A、B 的 SIFT 描述子分别为 k1*128 维、k2*128 ,可以将两图的各个尺度的描述子进行匹配,如果两个特征点匹配上 128 维,则表示其匹配上了.
当两幅图像的SIFT特征向量生成后,下一步我们采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量. 取图像 A 中的某个关键点,并找出其与图像 B 中欧式距离最近的前两个关键点,在这两个关键点中,如果最近的距离除以次近的距离少于某个比例阈值,则接受这一对匹配点. 降低这个比例阈值,SIFT匹配点数目会减少,但更加稳定.
为了排除因为图像遮挡和背景混乱而产生的无匹配关系的关键点, Lowe提出了比较最近邻距离与次近邻距离的方法, 距离比率ratio小于某个阈值的认为是正确匹配. 因为对于错误匹配,由于特征空间的高维性, 相似的距离可能有大量其他的错误匹配,从而它的ratio值比较高. Lowe推荐ratio的阈值为0.8. 但对大量任意存在尺度、旋转和亮度变化的两幅图片进行匹配,结果表明ratio取值在0. 4~0. 6之间最佳,小于0. 4的很少有匹配点,大于0. 6的则存在大量错误匹配点. 故,建议ratio的取值原则如下:
ratio=0. 4 对于准确度要求高的匹配;
ratio=0. 6 对于匹配点数目要求比较多的匹配;
ratio=0. 5 一般情况下.
也可按如下原则:当最近邻距离<200时 ratio=0. 6,反之ratio=0. 4. ratio的取值策略能排分错误匹配点.
10. OpenCV 实现
#!/usr/bin/python3
#!--*-- coding: utf-8 --*--
import cv2
#
imgfile1 = 'test1.jpg'
imgfile2 = 'test2.jpg'
img1 = cv2.imread(imgfile1)
img2 = cv2.imread(imgfile2)
#
sift = cv2.xfeatures2d.SIFT_create()
# 计算关键点和描述子
# kp为关键点keypoints, des为描述子descriptors
kp1, des1 = sift.detectAndCompute(img1, None)
kp2, des2 = sift.detectAndCompute(img2, None)
# 关键点可视化
img3 = cv2.drawKeypoints(img1, kp1, img1, color=(0, 255, 255))
img4 = cv2.drawKeypoints(img2, kp2, img2, color=(0, 255, 255))
# M1:匹配方法一
# 描述子匹配
index_params = dict(algorithm=1, trees=6)
search_params = dict(checks=50)
flann = cv2.FlannBasedMatcher(index_params, search_params)
# knn 计算描述子的匹配
matche = flann.knnMatch(des1, des2, k=2)
matchesMask = [[0, 0] for i in range(len(matche))]
# 匹配结果可视化
result = []
for m, n in matche:
if m.distance < 0.6 * n.distance:
result.append([m])
img5 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, matche, None, flags=2)
cv2.imshow("SIFTMatchResult", img5)
cv2.waitKey(0)
## M2: 匹配方法二
## K近邻算法求取在空间中距离最近的K个数据点,并将这些数据点归为一类
#matcher = cv2.BFMatcher()
#raw_matches = matcher.knnMatch(des1, des2, k = 2)
#good_matches = []
#for m1, m2 in raw_matches:
# # 如果最接近和次接近的比值大于一个既定的值,那么保留这个最接近的值,认为它和其匹配的点为good_match
# if m1.distance < ratio * m2.distance:
# good_matches.append([m1])
#img5 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, matche, None, flags=2)
11. SIFT 的缺点
SIFT在图像的不变特征提取方面拥有无与伦比的优势,但并不完美,仍然存在:
- 实时性不高.
- 有时特征点较少.
- 对边缘光滑的目标无法准确提取特征点,等.
如图所示,对模糊的图像和边缘平滑的图像,检测出的特征点过少,对圆更是无能为力. 近来不断有人改进,其中最著名的有SURF和CSIFT.