文章链接：https://blog.csdn.net/q_z_r_s

机器感知

一个专注于SLAM、机器视觉、Linux 等相关技术文章分享的公众号

从论文摘要中可以知道，SIFT算法的第一步就是进行尺度空间的极值检测，这是后续所有工作的基石。构建尺度空间使用的是差分高斯图像，这也是本论文的一大亮点，利用DoG来近似LoG，简化了计算（虽然计算量上仍然很大^_^)。

为什么可以用DoG来近似LoG，论文中给出了详细的证明过程，原文如下：

所构建的差分高斯图像的生成过程如下图所示，先生成不同尺度的（即不同sigma）高斯模糊图像，然后相邻高斯模糊图像进行相减，得到高斯差分图像DoG。

生成高斯图像代码如下：其中高斯模糊尺度规则如如下：
$\sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2$

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}^2sig[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;//如果找到的值明显错误，则直接返回falseif( 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 interpolationif( 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;//检查对比度是否符合要求，不符合的点删除，返回falseif( 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;
}

【SIFT算法】极值检测关键点精确定位相关推荐

目标检测_精确定位_2020
Side-Aware Boundary Localization for More Precise Object Detection 论文:https://arxiv.org/pdf/1912.042 ...
一维测量助手（尺寸检测、精确定位）
SIFT、SURF等关键点特征提取算法代码
文章目录 1.关键点特征提取算法 2.SIFT代码(python+opencv) 2.SURF代码(python+opencv) 3.SIFT和SURF的比较 1.关键点特征提取算法特征提取是提取出 ...
【目标识别】SIFT算法理论部分
SIFT算法理论部分目录一.介绍 1.1 SIFT算法 1.2 SIFT特征的获取方法 1.3 图像匹配和识别的方法 1.4 如何提高匹配准确率二.尺度空间极值检测 2.1 高斯模糊 2.2 尺度 ...
无人机航拍图像匹配——SIFT算法实践（含代码）
无人机航拍图像匹配--SIFT算法实践(含代码) 一.摘要二.SIFT算法的原理 1.尺度空间极值检测 &关键点定位尺度不变性&尺度空间高斯金字塔 2.方向分配 3.特征描述 4 ...
图像拼接（一）——SIFT算法新手入门级介绍！！！
前言因为博主毕设和图像拼接有关,所以博主简单地学习了SIFT算法,不知道以后研究生会不会接触图像这块,顺便也为了检验一下自己学的是否透彻,所以在CSDN上写一篇Blog,纪录一下,仅供各位小伙伴参考 ...
SIFT算法原理详解
通过<图像局部不变性特征与描述>学习SIFT,遇到各种Issue,总结了这篇博客和另外九篇博客.感谢关注,希望可以互相学习,不断提升.转载请注明链接:https://www.cnblogs ...
SIFT算法原理（2）-极值点的精确定位
在SIFT解析(一)建立高斯金字塔中,我们得到了高斯差分金字塔: 检测DOG尺度空间极值点 SIFT关键点是由DOG空间的局部极值点组成的.以中心点进行3X3X3的相邻点比较,检测其是否是图像域和尺度 ...
SIFT四部曲之——极值检测和定位
版权声明:本文为博主原创文章,未经博主允许不得转载.博客不用于商业活动,博主对博客的使用,拥有最终解释权本文为原创作品,未经本人同意,禁止转载,禁止用于商业用途!本人对博客使用拥有最终解释权欢迎关 ...

【SIFT算法】极值检测关键点精确定位

机器感知

【SIFT算法】极值检测关键点精确定位相关推荐

最新文章

热门文章