文章链接：JiauZhang-CSDN博客
机器感知
一个专注于SLAM、机器视觉、Linux 等相关技术文章分享的公众号
从论文摘要中可以知道，SIFT算法的第一步就是进行尺度空间的极值检测，这是后续所有工作的基石。构建尺度空间使用的是差分高斯图像，这也是本论文的一大亮点，利用DoG来近似LoG，简化了计算（虽然计算量上仍然很大^_^)。
为什么可以用DoG来近似LoG，论文中给出了详细的证明过程，原文如下：
所构建的差分高斯图像的生成过程如下图所示， 先生成不同尺度的（即不同sigma）高斯模糊图像，然后相邻高斯模糊图像进行相减，得到高斯差分图像DoG。
生成高斯图像代码如下：其中高斯模糊尺度规则如如下：
void SIFT::buildGaussianPyramid( const Mat& base, vector<Mat>& pyr, int nOctaves ) const
{
//每个octave所需的sigma值
vector<double> sig(nOctaveLayers + 3);
//高斯金字塔总层数
pyr.resize(nOctaves*(nOctaveLayers + 3));
// precompute Gaussian sigmas using the following formula:
// \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2
sig[0] = sigma;
double k = pow( 2., 1. / nOctaveLayers );
//生成对应sigma的值，结合下一个for循环可以看出，不同的octave相同层共用同样的sigma值
for( int i = 1; i < nOctaveLayers + 3; i++ )
{
double sig_prev = pow(k, (double)(i-1))*sigma;
double sig_total = sig_prev*k;
sig[i] = std::sqrt(sig_total*sig_total - sig_prev*sig_prev);
}
for( int o = 0; o < nOctaves; o++ )
{
for( int i = 0; i < nOctaveLayers + 3; i++ )
{
Mat& dst = pyr[o*(nOctaveLayers + 3) + i];
if( o == 0 && i == 0 )
dst = base;
// base of new octave is halved image from end of previous octave
//每个octave的第一层使用上一个octave的倒数第三个进行减半缩放得到的
else if( i == 0 )
{
const Mat& src = pyr[(o-1)*(nOctaveLayers + 3) + nOctaveLayers];
resize(src, dst, Size(src.cols/2, src.rows/2),
0, 0, INTER_NEAREST);
}
//其他层根据相应的sigma进行高斯模糊得到即可
else
{
const Mat& src = pyr[o*(nOctaveLayers + 3) + i-1];
GaussianBlur(src, dst, Size(), sig[i], sig[i]);
}
}
}
}接下来就是求高斯差分图像，代码如下：
//非常简单，对应图像相减即可得到
void SIFT::buildDoGPyramid( const vector<Mat>& gpyr, vector<Mat>& dogpyr ) const
{
int nOctaves = (int)gpyr.size()/(nOctaveLayers + 3);
dogpyr.resize( nOctaves*(nOctaveLayers + 2) );
for( int o = 0; o < nOctaves; o++ )
{
for( int i = 0; i < nOctaveLayers + 2; i++ )
{
const Mat& src1 = gpyr[o*(nOctaveLayers + 3) + i];
const Mat& src2 = gpyr[o*(nOctaveLayers + 3) + i + 1];
Mat& dst = dogpyr[o*(nOctaveLayers + 2) + i];
subtract(src2, src1, dst, noArray(), DataType<sift_wt>::type);
}
}
}得到高斯差分图像之后，接下来就是进行极值检测，采用的方法是：对同一个octave中的DoG图片进行所谓的“非最大值抑制”的方法进行筛选，如下图所示，“x”代表待检测点，其与周围3*9-1=26个点进行比较，当其比任何一个都大或比任何一个都小时就进入下一轮筛选。
void SIFT::findScaleSpaceExtrema( const vector<Mat>& gauss_pyr, const vector<Mat>& dog_pyr,
vector<KeyPoint>& keypoints ) const
{
int nOctaves = (int)gauss_pyr.size()/(nOctaveLayers + 3);
int threshold = cvFloor(0.5 * contrastThreshold / nOctaveLayers * 255 * SIFT_FIXPT_SCALE);
const int n = SIFT_ORI_HIST_BINS;
float hist[n];
KeyPoint kpt;
keypoints.clear();
for( int o = 0; o < nOctaves; o++ )
for( int i = 1; i <= nOctaveLayers; i++ )
{
int idx = o*(nOctaveLayers+2)+i;
//当前图像
const Mat& img = dog_pyr[idx];
//前一幅图像
const Mat& prev = dog_pyr[idx-1];
//后一幅图像
const Mat& next = dog_pyr[idx+1];
int step = (int)img.step1();
int rows = img.rows, cols = img.cols;
for( int r = SIFT_IMG_BORDER; r < rows-SIFT_IMG_BORDER; r++)
{
const sift_wt* currptr = img.ptr<sift_wt>(r);
const sift_wt* prevptr = prev.ptr<sift_wt>(r);
const sift_wt* nextptr = next.ptr<sift_wt>(r);
for( int c = SIFT_IMG_BORDER; c < cols-SIFT_IMG_BORDER; c++)
{
sift_wt val = currptr[c];
// find local extrema with pixel accuracy
//这个if语句就是在进行找最大或最小值
if( std::abs(val) > threshold &&
((val > 0 && val >= currptr[c-1] && val >= currptr[c+1] &&
val >= currptr[c-step-1] && val >= currptr[c-step] && val >= currptr[c-step+1] &&
val >= currptr[c+step-1] && val >= currptr[c+step] && val >= currptr[c+step+1] &&
val >= nextptr[c] && val >= nextptr[c-1] && val >= nextptr[c+1] &&
val >= nextptr[c-step-1] && val >= nextptr[c-step] && val >= nextptr[c-step+1] &&
val >= nextptr[c+step-1] && val >= nextptr[c+step] && val >= nextptr[c+step+1] &&
val >= prevptr[c] && val >= prevptr[c-1] && val >= prevptr[c+1] &&
val >= prevptr[c-step-1] && val >= prevptr[c-step] && val >= prevptr[c-step+1] &&
val >= prevptr[c+step-1] && val >= prevptr[c+step] && val >= prevptr[c+step+1]) ||
(val < 0 && val <= currptr[c-1] && val <= currptr[c+1] &&
val <= currptr[c-step-1] && val <= currptr[c-step] && val <= currptr[c-step+1] &&
val <= currptr[c+step-1] && val <= currptr[c+step] && val <= currptr[c+step+1] &&
val <= nextptr[c] && val <= nextptr[c-1] && val <= nextptr[c+1] &&
val <= nextptr[c-step-1] && val <= nextptr[c-step] && val <= nextptr[c-step+1] &&
val <= nextptr[c+step-1] && val <= nextptr[c+step] && val <= nextptr[c+step+1] &&
val <= prevptr[c] && val <= prevptr[c-1] && val <= prevptr[c+1] &&
val <= prevptr[c-step-1] && val <= prevptr[c-step] && val <= prevptr[c-step+1] &&
val <= prevptr[c+step-1] && val <= prevptr[c+step] && val <= prevptr[c+step+1])))
{
int r1 = r, c1 = c, layer = i;
//找到之后进行下一步筛选，不符合的删除，符合的进入下一步
if( !adjustLocalExtrema(dog_pyr, kpt, o, layer, r1, c1,
nOctaveLayers, (float)contrastThreshold,
(float)edgeThreshold, (float)sigma) )
continue;
//以下为关键点主方向定位，下篇文章详解
.....
}
}
}
}
}
}
}经过以上步骤之后，接下来就是对这些点进行精确定位，即关键点定位。
static bool adjustLocalExtrema( const vector<Mat>& dog_pyr, KeyPoint& kpt, int octv,
int& layer, int& r, int& c, int nOctaveLayers,
float contrastThreshold, float edgeThreshold, float sigma )
{
const float img_scale = 1.f/(255*SIFT_FIXPT_SCALE);
const float deriv_scale = img_scale*0.5f;
const float second_deriv_scale = img_scale;
const float cross_deriv_scale = img_scale*0.25f;
float xi=0, xr=0, xc=0, contr=0;
int i = 0;
//此for循环是在对上一步得到的极值点进行精确定位
//所用方法为：在极值点处进行二阶泰勒展开
for( ; i < SIFT_MAX_INTERP_STEPS; i++ )
{
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);
//利用LU分解方法求解Hx=dD的解，即相对于极值点的各坐标分量的偏移量
Vec3f X = H.solve(dD, DECOMP_LU);
//因为真正的解应该是求Hx=-dD的解，所以此处赋值会取反
xi = -X[2];
xr = -X[1];
xc = -X[0];
//当求解值没有偏移出极值点0.5，即1/2时，直接结束循环
if( std::abs(xi) < 0.5f && std::abs(xr) < 0.5f && std::abs(xc) < 0.5f )
break;
//如果找到的值明显错误，则直接返回false
if( std::abs(xi) > (float)(INT_MAX/3) ||
std::abs(xr) > (float)(INT_MAX/3) ||
std::abs(xc) > (float)(INT_MAX/3) )
return false;
//程序走到这里说明各偏移绝对值都＞=0.5，cvRound进行四舍五入更新极值点坐标
//然后继续下一轮计算
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;
//找到了符合上述要求的极值点坐标，进行下一轮筛选
{
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为此点的二阶泰勒展开在所求得点出的函数值
contr = img.at<sift_wt>(r, c)*img_scale + t * 0.5f;
//检查对比度是否符合要求，不符合的点删除，返回false
if( std::abs( contr ) * nOctaveLayers < contrastThreshold )
return false;
// principal curvatures are computed using the trace and det of Hessian
//利用Hessian矩阵去除边缘响应，详细内容见论文公式
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;
if( det <= 0 || tr*tr*edgeThreshold >= (edgeThreshold + 1)*(edgeThreshold + 1)*det )
return false;
}
kpt.pt.x = (c + xc) * (1 << octv);
kpt.pt.y = (r + xr) * (1 << octv);
kpt.octave = octv + (layer << 8) + (cvRound((xi + 0.5)*255) << 16);
kpt.size = sigma*powf(2.f, (layer + xi) / nOctaveLayers)*(1 << octv)*2;
kpt.response = std::abs(contr);
return true;
}
514
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
