在上一章节中,我们介绍了如何对图像模拟不同尺度空间,以在尺度空间下检测稳定特征点。构造尺度空间的方法是通过构造高斯金字塔和DoG实现的。那么在本章节中,我们将讨论如何在构造好的DoG空间内检测关键点,以及筛选出不稳定的关键点。
4.检测并筛选关键点(keypoints)
4.1 检测尺度空间极值点
计算DoG空间的极值点(局部最大值或最小),是通过这样的方法实现的:每一个像素点需要与同一层内相邻的8个像素点以及相邻上下两层的
9∗2=18
个,共26个像素点比较大小。如果当前像素点比这26个像素点都要大或者都要小,则当前像素点为极值点。如下图所示,其中x号像素表示当前像素,绿色圈表示像素为相邻比较像素,
4.2确定检测点的精确位置
实际上,经过这样步骤检测的极值点是离散的,可能并不是真实的极值点,这可以从下面这个示意图看出,真正的极值点是位于检测的离散极值点附近的。所以为了得到关键点的准确位置,Lowe是通过对关键点附近的离散点进行拟合,然后求极值,从而得到关键点的准确位置的。
为了说明对离散极值点拟合逼近真实极值点的过程,我们首先举出一个一维函数插值的例子。
假设有,
f(−1)=1,f(0)=6,f(1)=5
, 求
f(x)
在区间 [-1,1]上的最大值。
f(x)
在
f(0)
附近的泰勒展开式为,
另外有, f(x) 在 x 处的导数可以写成离散的形式为,
所以有,
对函数 f(x) 求导得到极值点的坐标为 x^=13 , 则极值点的坐标为 f(x^)=613
同样的道理,我们对检测出的离散关键点进行拟合求极值,假设在尺度
σ
的图像空间
D(x,y)
处检测到一个极值点(关键点),那么这个极值点的空间坐标为
(x,y,σ)
,通过上面分析我们知道,真实的极值点有可能落在这个离散极值点的附近。设
(x,y,σ)
附近的点相对于该极值点的偏移量为
(Δx,Δy,Δσ)
,那么
D(Δx,Δy,Δσ)
可以表示为点
(x,y,σ)
处的泰勒展开式,
X 表示偏移量
多次迭代,每次更新极值点的位置, Lowe算法中最多迭代5次,最终得到候选极值点的精确位置 x^ 。
将极值点的精确位置 x^ 带入 D ,得到极值点的像素值。
如果
4.3 消除边缘响应
DoG对图像中的边缘具有较强的响应,而一旦特征点落在边缘上,这些点就是不稳定的检测点,所以需要将这些边缘点去掉。
而在横跨边缘处,DoG往往具有较大的主曲率,而在垂直边缘的方向具有较小的主曲率,所以我们可以考虑,如果检查最大主曲率与最小主曲率的比值,如果比值大于某一阈值,就认为是边缘点,并剔除该关键点,就可以将边缘不稳定的关键点排除掉。
海塞矩阵的特征值是和D的主曲率是成正比的,即计算主曲率可以转换为计算海塞矩阵。
D 可以通过邻近像素点求差分得到。设H的最大特征值和最小特征值分别为
如果 γ 表示最大特征值和最小特征值之间的比例,即 α=γβ ,那么有,
当两个特征值相等时, (γ+1)2γ 有最小值,而随着两者比例 γ 的增加, (γ+1)2γ 也增加。所以想要确定D的主曲率是否大于某一阈值 γ^ , 只需要检查下面的式子是否成立,
在Lowe算法中,给定的阈值 γ^ 是10,也就是说当最大主曲率与最小主曲率之间的比值大于10时,该关键点被删除。
经过上面三个步骤,即首先通过比较像素点与同层8邻域像素点和相邻上下两层像素点确定候选关键点;然后通过曲线拟合确定关键点的准确位置,并删除对比度较低的极值点;最后通过判断最大主曲率与最小主曲率之间的比值,删除边缘DoG相应较强的关键点。最终得到的关键点用于SIFT特征计算。
4.4 C++源码解读
1. 曲线拟合精确确定关键点的坐标
//迭代插值过程确定关键点精确坐标
for( ; i < SIFT_MAX_INTERP_STEPS; i++ )
{
//octv当前关键点的组数;layer当前关键点的层数;nOctaveLayers每组的总层数
//dog_pyr:DoG金字塔
int idx = octv*(nOctaveLayers+2) + layer;
const Mat& img = dog_pyr[idx];
const Mat& prev = dog_pyr[idx-1];
const Mat& next = dog_pyr[idx+1];
Vec3f dD((img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1))*deriv_scale,
(img.at<sift_wt>(r+1, c) - img.at<sift_wt>(r-1, c))*deriv_scale,
(next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c))*deriv_scale);
float v2 = (float)img.at<sift_wt>(r, c)*2;
float dxx = (img.at<sift_wt>(r, c+1) + img.at<sift_wt>(r, c-1) - v2)*second_deriv_scale;
float dyy = (img.at<sift_wt>(r+1, c) + img.at<sift_wt>(r-1, c) - v2)*second_deriv_scale;
float dss = (next.at<sift_wt>(r, c) + prev.at<sift_wt>(r, c) - v2)*second_deriv_scale;
float dxy = (img.at<sift_wt>(r+1, c+1) - img.at<sift_wt>(r+1, c-1) -
img.at<sift_wt>(r-1, c+1) + img.at<sift_wt>(r-1, c-1))*cross_deriv_scale;
float dxs = (next.at<sift_wt>(r, c+1) - next.at<sift_wt>(r, c-1) -
prev.at<sift_wt>(r, c+1) + prev.at<sift_wt>(r, c-1))*cross_deriv_scale;
float dys = (next.at<sift_wt>(r+1, c) - next.at<sift_wt>(r-1, c) -
prev.at<sift_wt>(r+1, c) + prev.at<sift_wt>(r-1, c))*cross_deriv_scale;
Matx33f H(dxx, dxy, dxs,
dxy, dyy, dys,
dxs, dys, dss);
Vec3f X = H.solve(dD, DECOMP_LU);
//关键步骤-》
//插值求得的偏移极值点坐标
xi = -X[2];
xr = -X[1];
xc = -X[0];
//如果求得的偏移极值点大于0.5,则认为真实的极值点是位于其他的采样点的,此时需要更新极值点的坐标
if( std::abs(xi) < 0.5f && std::abs(xr) < 0.5f && std::abs(xc) < 0.5f )
break;
if( std::abs(xi) > (float)(INT_MAX/3) ||
std::abs(xr) > (float)(INT_MAX/3) ||
std::abs(xc) > (float)(INT_MAX/3) )
return false;
//根据求得的偏移极值点更新当前极值点的坐标
c += cvRound(xc);
r += cvRound(xr);
layer += cvRound(xi);
if( layer < 1 || layer > nOctaveLayers ||
c < SIFT_IMG_BORDER || c >= img.cols - SIFT_IMG_BORDER ||
r < SIFT_IMG_BORDER || r >= img.rows - SIFT_IMG_BORDER )
return false;
}
// ensure convergence of interpolation
if( i >= SIFT_MAX_INTERP_STEPS )
return false;
2.判断最大主曲率与最小主曲率比值来判断是否是边缘
{
//octv当前关键点的组数;layer当前关键点的层数;nOctaveLayers每组的总层数
//dog_pyr:DoG金字塔
int idx = octv*(nOctaveLayers+2) + layer;
const Mat& img = dog_pyr[idx];
const Mat& prev = dog_pyr[idx-1];
const Mat& next = dog_pyr[idx+1];
Matx31f dD((img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1))*deriv_scale,
(img.at<sift_wt>(r+1, c) - img.at<sift_wt>(r-1, c))*deriv_scale,
(next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c))*deriv_scale);
float t = dD.dot(Matx31f(xc, xr, xi));
contr = img.at<sift_wt>(r, c)*img_scale + t * 0.5f;
if( std::abs( contr ) * nOctaveLayers < contrastThreshold )
return false;
// 利用海塞矩阵的迹和行列式计算D主曲率的比值
float v2 = img.at<sift_wt>(r, c)*2.f;
float dxx = (img.at<sift_wt>(r, c+1) + img.at<sift_wt>(r, c-1) - v2)*second_deriv_scale;
float dyy = (img.at<sift_wt>(r+1, c) + img.at<sift_wt>(r-1, c) - v2)*second_deriv_scale;
float dxy = (img.at<sift_wt>(r+1, c+1) - img.at<sift_wt>(r+1, c-1) -
img.at<sift_wt>(r-1, c+1) + img.at<sift_wt>(r-1, c-1)) * cross_deriv_scale;
float tr = dxx + dyy;//迹
float det = dxx * dyy - dxy * dxy;//行列式
//判断主曲率的比例是否小于阈值edgeThreshold,如果不是,则认为是边缘关键点,删除该关键点
if( det <= 0 || tr*tr*edgeThreshold >= (edgeThreshold + 1)*(edgeThreshold + 1)*det )
return false;
}