SIFT四部曲之——极值检测和定位

2021年9月20日 8点热度 0条评论 来源: 晨凫追风

版权声明:本文为博主原创文章,未经博主允许不得转载。博客不用于商业活动,博主对博客的使用,拥有最终解释权

本文为原创作品,未经本人同意,禁止转载,禁止用于商业用途!本人对博客使用拥有最终解释权

欢迎关注我的网上图书室:晨凫追风 和 微信公众号:青春当追风

 

 

这篇博客要做的是对sift的第二个步骤进行介绍和总结。然后对opencv中的源码进行解释,因为对这一部分的理解主要是根据Opencv中的源码反推的,感谢开源的作者们!

Ok言归正传!

Sift的第二个步骤主要是在图像中找到那些尺度不变性的点,这里叫做特征点。这些特征点是由DOG金字塔中的图像通过一系列的筛选规则得到的,所以这一步骤就是为了介绍这些筛选规则,以及用这些筛选规则怎么选出代表。

首先在前面介绍一些数学的知识吧,这些知识是经过推导得到的一些结论,所以记住它们(微笑脸)。

1.有限差分法求导(这里没有推导,也不太严谨,不用怕

这个知识主要是介绍对于图像中函数的求导法则,介绍的是二元函数的求导。当然在尺度空间中,像素值f是坐标(x,y)和尺度σ的三元函数,它的求导法则可以类比。看图1和公式(1-5)这里的求导其实可以理解成,每一个像素值对哪一个方向求导,则与其他方向上的无关,就是把其他方向看成一个常数来理解。在同一层的图像中,与尺度的方向无关,所以把尺度看成一个常数。

图1

推广到三元函数,就是增加σ轴向的求导如图2,这里把y方向上的看成无关,来求导尺度轴和x方向的关系。

图2

以此类推y和σ方向上的关系:

图3

2.三阶矩阵求逆的公式

该矩阵存在逆矩阵,则它的行列式不等0,即:

所以

    

至此所要用到的一部分数学知识就讲到这里。接下来开始讲选举条例!

1.第0轮推举(阈值检测)

由于我们不能保证图片就是完全没有噪声的,所以在推选之前我们要排除一些对比度比较低的点,这些点很不稳定,容易受到干扰,本着宁缺毋滥的原则,先把不稳定因子排除。这个阈值一般是自己设定的,在opencv中设定它为:0.5*T/S,这里的T lowe 定义它为0.04,s是该点所在的层。拿到这个入场券的同学,进入下一轮选举

2.第一轮选举(极值检测)

在上一篇博客里讲到了,某某某证明了,高斯差分金字塔里的极值点比一般的(角点,Hessian,梯度函数)这些方法得到的点性质稳定。所以很自然的就要用到极值点来做候选人。在sift算法里面极值点的检测方法是(图5):

图5

对图像中的每一个点进行遍历,判断每一个点是否是极值,判断的标准是,把这个点与相邻的层,上下层9+9=18个点和同一层中8个点,共18+8=26个点进行比较,看这个点是否是最大值,还是最小值,若是,则保留下来,当做候选人。

3.第二轮选举(极值点精确定位)

在sift算法中,把尺度看成连续的,而在前面的极值检测中,只是把它当成离散点来计算,于是极值点的检测就会出现图6的情况,这怎么能符合Lowe同学严谨的风格呢?

图6

所以在这一轮的选举中,我们的目标就是找到精确的极值点(成仙之前当然要经过考验)其实这里的精确点的位置,是一个亚像素级别的概念,比如一个数的个位十位百位都确定了,要更加精细的去找这个数的小数点后的那个数据,就用插值来计算。

这里又是一系列的数学公式(打公式好麻烦呀!哭!)数学是根本呀!大一高数学不好,本科兵败如山倒!

言归正传!这里我们需要高数和线性代数的一些知识。高数中提过,函数f(x)在x0处可以进行泰勒展开。于是我们也在这个检测到的极值点处处,对它进行泰勒级数展开(这里只展开到二阶项,高次项砍掉),以此来拟合这个三维二次函数。于是得到

把上面的式子写成矢量的形式:

这真是个伟大的公式!,好了我们拟合出来了该点附近的函数了,表示拟合之后连续空间下面的插值点坐标,设表示:相对于插值中心的偏移量,这下公式(13)可以用偏移量表示:

(14)

该咋求极值呢?当然是求导啦!导数为0的那个点(排除驻点的情况,非要考虑驻点的话,不好意思,我也不会,微笑脸)于是得到下面的求导式子:

(15)

让导数等于0,就可以得到极值点下的相对于插值中心的偏移量:

(16)

把式(16)代入式(14)可以得到该极值点下的极值为:

(17)

Ok,其实到这里,已经把精确的极值点找到了,但是真的是这样的吗?答案当然不是。上述的知识只是告诉我们找到偏移量的一个大的方向,真正的细节并没有说明白!

这里该注意几个问题:

  1. 当求出来的偏移量很大的时候,这时就表明精确的极值点已经完全偏离了离散比较得到的极值点,这时候就得删除它
  2. 当求出来的偏移量大于0.5,(只要x,y,和σ任意一个量大于0.5)就表示这个插值点偏离了插值中心,这时候就应该改变插值中心,继续使用上述的泰勒展开进行拟合,一直到插值偏移量小于0.5,(指的是三个量的偏移都小于0.5)
  3. 当然在第二个问题中,它的解决办法是一个迭代的过程,当然它也有自己的迭代停止条件,但是如果要迭代很多次才能找到那个精确点,那说明这个点本身就有自己不稳定的成分,(与问题1类似)所以就要设置一个迭代次数的限制阈值,超过这个阈值,迭代停止!(不是没给你机会,给你机会,你不珍惜!微笑脸)这个点也得出局!

4.第三轮选举(低对比度筛选)

在上一步的检测中我们得到了精确点的位置了,也得到精确点的值了,这时候又有一个问题来了,精确点的值如果很小,那很大程度上是不稳定的点,于是很遗憾,这些精确点应该出局!所以这一轮的标准是:

(18)

论文中T=0.04,s表示处于该组的第几层。

5.第四轮筛选(消除边缘效应)

首先恭喜存活下来的选手,来到了最后一轮的筛选了,过了这一关就可以晋级特征点行列了!

在DOG函数中,它存在比较强的边缘效应,而当特征点在边缘的时候,这些点就会很不稳定,(1是因为这些点比较难定位,他们具有定位的歧义性。2是因为这些点容易受到噪声的干扰而变得不稳定)所以呢,我们应该把这些披着羊皮的狼(边缘效应很强的点)找出来。啊!又来数学公式!推导!这里的方法和Harris角点检测算法的原理相似,即DOG函数欠佳的峰值点附近,在横跨边缘的方向上有较大的主曲率,而在垂直于边缘的方向上具有比较小的主曲率,这个主曲率可以通过Hessian矩阵得到:

(19)

其中 分别表示对DOG图像的像素在x轴和y轴方向上的二阶导数和二阶混合偏导数。由于该像素点的主曲率和H的特征值成正比,该矩阵的特征值代表着x轴和y轴方向上的梯度。但是为了避免求解矩阵的特征值我们只要知道这两个值的比例就可以得到该点的主曲率。

引入两个量

(20)

这里我们先删除掉那些行列式为负数的点,即,因为如果像素的曲率有不同的符号,则该点肯定不会是特征点。

接着设并且,其中,于是得到:

(21)

上面式子的结果只和两个特征值的比例有关,只有在两个特征值相等的时候,式(21)才最小,随着的增大,该式越大,说明两个特征值的比值就越大,即在某一方向上的梯度值就越大,而另一方向上的梯度越小,这便是边缘所符合的情况。所以为了剔除这些边缘点,我们让值小于一定的阈值,因此为了检测主曲率是否在某个阈值之下,只需检测

(22)

对于不满足式22的点(用拟合的精确点来计算),只有一个原则,删除!当然这里需要给出Lowe所建议的值,为10.

 

 

Ok!恭喜!通过上述选举的同学,满级!

最后附上OPENCV关于这一部分的源码解释

//在DOG金字塔内找到极值点的函数


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( 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;
                        float scl_octv = kpt.size*0.5f/(1 << o);
                        float omax = calcOrientationHist(gauss_pyr[o*(nOctaveLayers+3) + layer],
                                                         Point(c1, r1),
                                                         cvRound(SIFT_ORI_RADIUS * scl_octv),
                                                         SIFT_ORI_SIG_FCTR * scl_octv,
                                                         hist, n);
                        float mag_thr = (float)(omax * SIFT_ORI_PEAK_RATIO);
                        for( int j = 0; j < n; j++ )
                        {
                            int l = j > 0 ? j - 1 : n - 1;
                            int r2 = j < n-1 ? j + 1 : 0;

                            if( hist[j] > hist[l]  &&  hist[j] > hist[r2]  &&  hist[j] >= mag_thr )
                            {
                                float bin = j + 0.5f * (hist[l]-hist[r2]) / (hist[l] - 2*hist[j] + hist[r2]);
                                bin = bin < 0 ? n + bin : bin >= n ? bin - n : bin;
                                kpt.angle = 360.f - (float)((360.f/n) * bin);
                                if(std::abs(kpt.angle - 360.f) < FLT_EPSILON)
                                    kpt.angle = 0.f;
                                keypoints.push_back(kpt);
                            }
                        }
                    }
                }
            }
        }
}

//精确点定位函数

//dog_pyr为DOG金字塔,kpt为特征点,octv和layer为极值点所在的组和层,r和c为极值点坐标

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);//先进行归一化处理,这里的插值函数都为[0,1]
    const float deriv_scale = img_scale*0.5f;              //公式(6)中分母的倒数,这里h=1
    const float second_deriv_scale = img_scale;          //公式(7)中分母的倒数,这里h=1
    const float cross_deriv_scale = img_scale*0.25f;   //公式(8)中分母的倒数,这里h=1

    float xi=0, xr=0, xc=0, contr=0;
    int i = 0;

    for( ; i < SIFT_MAX_INTERP_STEPS; i++ )            //SIFT_MAX_INTERP_STEPS表示迭代的次数,为5次
    {
        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];                  //金字塔上层图像(上方)

		
		//变量dD是对X的一阶偏导数公式1,2,6
		
        Vec3f dD((img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1))*deriv_scale,           //对x的一阶偏导
                 (img.at<sift_wt>(r+1, c) - img.at<sift_wt>(r-1, c))*deriv_scale,              //对y的一阶偏导
                 (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;

					 
					 //公式13中的f对向量X求二阶导,展开式见公式12
        Matx33f H(dxx, dxy, dxs,
                  dxy, dyy, dys,
                  dxs, dys, dss);

        Vec3f X = H.solve(dD, DECOMP_LU);   //公式16
		
		
//上面式子和公式16差一个负号,这里补上

        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;

    {
        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;    //对比度检测公式18
        if( std::abs( contr ) * nOctaveLayers < contrastThreshold )
            return false;

			
			
			//第4轮筛选
        // principal curvatures are computed using the trace and det of 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;
}

下面是matlab代码,这部分代码是承接上一篇博客的下一部分

% 下一步是查找差分高斯金字塔中的局部极值,并通过曲率和照度进行检验
curvature_threshold = ((curvature_threshold + 1)^2)/curvature_threshold;


% 二阶微分核
xx = [ 1 -2  1 ];
yy = xx';
xy = [ 1 0 -1; 0 0 0; -1 0 1 ]/4;


raw_keypoints = [];%%原始特征点
contrast_keypoints = [];%%对比度
curve_keypoints = [];%%曲率


% 在高斯金字塔中查找局部极值
if interactive >= 1
   fprintf( 2, 'Locating keypoints...\n' );
end
tic;
loc = cell(size(DOG_pyr));
for octave = 1:octaves
   if interactive >= 1
      fprintf( 2, '\tProcessing octave %d\n', octave );
   end
   for interval = 2:(intervals+1)
      keypoint_count = 0;
      contrast_mask = abs(DOG_pyr{octave}(:,:,interval)) >= contrast_threshold;
      loc{octave,interval} = zeros(size(DOG_pyr{octave}(:,:,interval)));
      if exist('corrsep') == 3
         edge = 1;
      else         
         edge = ceil(filter_size(octave,interval)/2);
      end      
      for y=(1+edge):(size(DOG_pyr{octave}(:,:,interval),1)-edge)
          y_size = size(DOG_pyr{octave}(:,:,interval),1)-edge;
         for x=(1+edge):(size(DOG_pyr{octave}(:,:,interval),2)-edge)
            x_size = size(DOG_pyr{octave}(:,:,interval),2)-edge;
            if object_mask(round(y*subsample(octave)),round(x*subsample(octave))) == 1 
               
               
               if( (interactive >= 2) | (contrast_mask(y,x) == 1) ) 
                  
               % 通过空间核尺度检测最大值和最小值
                  tmp = DOG_pyr{octave}((y-1):(y+1),(x-1):(x+1),(interval-1):(interval+1)); %%取出该位置的上下3层,共27个点 
                  pt_val = tmp(2,2,2);            %%该点是中间的待检测点第二层第二排第二列
                  if( (pt_val == min(tmp(:))) | (pt_val == max(tmp(:))) )
%                       % 存储极值点                                           
%                      raw_keypoints = [raw_keypoints; x*subsample(octave) y*subsample(octave)]; %%乘以采样的数,还原到原始图片的坐标
 
%%迭代插值
                         for i=1:interactive
                             Ddx = (DOG_pyr{octave}(y,(x+1),interval)-DOG_pyr{octave}(y,(x-1),interval))*0.5;
                             Ddy = (DOG_pyr{octave}(y+1,x,interval)-DOG_pyr{octave}(y-1,x,interval))*0.5;
                             Ddsigama = (DOG_pyr{octave}(y,x,interval+1)-DOG_pyr{octave}(y,x,interval-1))*0.5;
                             
                             
                             v2 = DOG_pyr{octave}(y,x,interval)*2;
                             dxx = (DOG_pyr{octave}(y,(x+1),interval)-DOG_pyr{octave}(y,(x-1),interval) - v2);
                             dyy = (DOG_pyr{octave}(y+1,x,interval)-DOG_pyr{octave}(y-1,x,interval) - v2);
                             dss = (DOG_pyr{octave}(y,x,interval+1)-DOG_pyr{octave}(y,x,interval-1)-v2);
                             
                             
                             dxy = (DOG_pyr{octave}(y+1,x+1,interval)-DOG_pyr{octave}(y+1,x-1,interval)-DOG_pyr{octave}(y-1,x+1,interval)+DOG_pyr{octave}(y-1,x-1,interval))*0.25;
                             dxs = (DOG_pyr{octave}(y,x+1,interval+1)-DOG_pyr{octave}(y,x-1,interval+1)-DOG_pyr{octave}(y,x+1,interval-1)+DOG_pyr{octave}(y,x-1,interval-1))*0.25;
                             dys = (DOG_pyr{octave}(y+1,x,interval+1)-DOG_pyr{octave}(y-1,x,interval+1)-DOG_pyr{octave}(y+1,x,interval-1)+DOG_pyr{octave}(y-1,x,interval-1))*0.25;
                             HH = [dxx,dxy,dxs;dxy,dyy,dys;dxs,dys,dss];
                         
                             D = [Ddx,Ddy,Ddsigama];
                             X = inv(HH)*D';
                             xi = -X(3);
                             xr = -X(2);
                             xc = -X(1);
                             if (abs(xi)<0.5&abs(xr)<0.5&abs(xc)<0.5)
                                 break;
                             end
                              if (abs(xi)>2147483647|abs(xr)>2147483647|abs(xc)>2147483647)
                                  i = interactive + 1;
                                 break;
                              end
                             x = x+ round(xc);
                             y= y + round(xr);
                             interval = interval + round(xi);
                             if(interval<1|interval>(intervals+1)|x<edge|x>x_size|y<edge|y>y_size)
                                 i = interactive + 1;
                             break;
                             end                             
                             end
                         
                      if(i<interactive)
                          Ddx = (DOG_pyr{octave}(y,(x+1),interval)-DOG_pyr{octave}(y,(x-1),interval))*0.5;
                          Ddy = (DOG_pyr{octave}(y+1,x,interval)-DOG_pyr{octave}(y-1,x,interval))*0.5;
                          Ddsigama = (DOG_pyr{octave}(y,x,interval+1)-DOG_pyr{octave}(y,x,interval-1))*0.5;
                          D = [Ddx,Ddy,Ddsigama];
                          t = (D) * (X) ;
                          contr = DOG_pyr{octave}(y,x,interval) + 0.5*t;
                                        

                      % 存储对灰度大于对比度阈值的点的坐标
                     if (abs(contr) >= contrast_threshold)
                        raw_keypoints = [raw_keypoints; x*subsample(octave) y*subsample(octave)]; %%乘以采样的数,还原到原始图片的坐标
                        contrast_keypoints = [contrast_keypoints; raw_keypoints(end,:)];%%存储的是最后一个的值,即最新检测到的值
                        
                        % 计算局部极值的Hessian矩阵
                        Dxx = sum(DOG_pyr{octave}(y,x-1:x+1,interval) .* xx);
                        Dyy = sum(DOG_pyr{octave}(y-1:y+1,x,interval) .* yy);
                        Dxy = sum(sum(DOG_pyr{octave}(y-1:y+1,x-1:x+1,interval) .* xy));
                        
                        % 计算Hessian矩阵的直迹和行列式.
                        Tr_H = Dxx + Dyy;
                        Det_H = Dxx*Dyy - Dxy^2;
                        

                        % 计算主曲率.
                        curvature_ratio = (Tr_H^2)/Det_H;
                        
                        if ((Det_H >= 0) & (curvature_ratio < curvature_threshold))

                           % 存储主曲率小于阈值的的极值点的坐标(非边缘点)
                           curve_keypoints = [curve_keypoints; raw_keypoints(end,:)];
                           
                           % 将该点的位置的坐标设为1,并计算点的数量.
                           loc{octave,interval}(y,x) = 1;
                           keypoint_count = keypoint_count + 1;
                        end
                     end
                     end                  
                  end
               end               
            end
         end         
      end
      if interactive >= 1
         fprintf( 2, '\t\t%d keypoints found on interval %d\n', keypoint_count, interval );
      end
   end
end
keypoint_time = toc;
if interactive >= 1
   fprintf( 2, 'Keypoint location time %.2f seconds.\n', keypoint_time );
end   

% 在交互模式下显示特征点检测的结果.
if interactive >= 2
   fig = figure;
   clf;
   imshow(im);
   hold on;
   plot(raw_keypoints(:,1),raw_keypoints(:,2),'y+');
%    resizeImageFig( fig, size(im), 2 );
   fprintf( 2, 'DOG extrema (2x scale).\nPress any key to continue.\n' );
%    pause;
%    close(fig);
   fig = figure;
%    clf;
   imshow(im);
   hold on;
   plot(contrast_keypoints(:,1),contrast_keypoints(:,2),'y+');
%    resizeImageFig( fig, size(im), 2 );
   fprintf( 2, 'Keypoints after removing low contrast extrema (2x scale).\nPress any key to continue.\n' );
%    pause;
%    close(fig);
   fig = figure;
%    clf;
   imshow(im);
   hold on;
   plot(curve_keypoints(:,1),curve_keypoints(:,2),'y+');
%    resizeImageFig( fig, size(im), 2 );
   fprintf( 2, 'Keypoints after removing edge points using principal curvature filtering (2x scale).\nPress any key to continue.\n' );
%    pause;
%    close(fig);  
end
clear raw_keypoints contrast_keypoints curve_keypoints

经过前面三轮选举得到:

经过第四轮选举得到:

本文为原创作品,未经本人同意,禁止转载,禁止用于商业用途!本人对博客使用拥有最终解释权


更多内容请看下一篇



    原文作者:晨凫追风
    原文地址: https://blog.csdn.net/hit2015spring/article/details/52972890
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系管理员进行删除。