计算机视觉
图像处理

【OpenCV】SIFT原理与源码分析:关键点搜索与定位

文章目录

《SIFT原理与源码分析》系列文章索引:http://www.cvvision.cn/2178.html

由前一步《DoG尺度空间构造》,我们得到了DoG高斯差分金字塔:

如上图的金字塔,高斯尺度空间金字塔中每组有五层不同尺度图像,相邻两层相减得到四层DoG结果。关键点搜索就在这四层DoG图像上寻找局部极值点。

DoG局部极值点

寻 找DoG极值点时,每一个像素点和它所有的相邻点比较,当其大于(或小于)它的图像域和尺度域的所有相邻点时,即为极值点。如下图所示,比较的范围是个 3×3的立方体:中间的检测点和它同尺度的8个相邻点,以及和上下相邻尺度对应的9×2个点——共26个点比较,以确保在尺度空间和二维图像空间都检测到 极值点。

在一组中,搜索从每组的第二层开始,以第二层为当前层,第一层和第三层分别作为立方体的的上下层;搜索完成后再以第三层为当前层做同样的搜索。所以每层的点搜索两次。通常我们将组Octaves索引以-1开始,则在比较时牺牲了-1组的第0层和第N组的最高层

高斯金字塔,DoG图像及极值计算的相互关系如上图所示。

关键点精确定位

以上极值点的搜索是在离散空间进行搜索的,由下图可以看到,在离散空间找到的极值点不一定是真正意义上的极值点。可以通过对尺度空间DoG函数进行曲线拟合寻找极值点来减小这种误差。
利用DoG函数在尺度空间的Taylor展开式:
则极值点为:
程序中还除去了极值小于0.04的点。如下所示:
  1. // Detects features at extrema in DoG scale space.  Bad features are discarded
  2. // based on contrast and ratio of principal curvatures.
  3. // 在DoG尺度空间寻特征点(极值点)
  4. void SIFT::findScaleSpaceExtrema( const vector<Mat>& gauss_pyr, const vector<Mat>& dog_pyr,
  5.                                   vector<KeyPoint>& keypoints ) const
  6. {
  7.     int nOctaves = (int)gauss_pyr.size()/(nOctaveLayers + 3);
  8.     // The contrast threshold used to filter out weak features in semi-uniform
  9.     // (low-contrast) regions. The larger the threshold, the less features are produced by the detector.
  10.     // 过滤掉弱特征的阈值 contrastThreshold默认为0.04
  11.     int threshold = cvFloor(0.5 * contrastThreshold / nOctaveLayers * 255 * SIFT_FIXPT_SCALE);
  12.     const int n = SIFT_ORI_HIST_BINS; //36
  13.     float hist[n];
  14.     KeyPoint kpt;
  15.     keypoints.clear();
  16.     forint o = 0; o < nOctaves; o++ )
  17.         forint i = 1; i <= nOctaveLayers; i++ )
  18.         {
  19.             int idx = o*(nOctaveLayers+2)+i;
  20.             const Mat& img = dog_pyr[idx];
  21.             const Mat& prev = dog_pyr[idx-1];
  22.             const Mat& next = dog_pyr[idx+1];
  23.             int step = (int)img.step1();
  24.             int rows = img.rows, cols = img.cols;
  25.             forint r = SIFT_IMG_BORDER; r < rows-SIFT_IMG_BORDER; r++)
  26.             {
  27.                 const short* currptr = img.ptr<short>(r);
  28.                 const short* prevptr = prev.ptr<short>(r);
  29.                 const short* nextptr = next.ptr<short>(r);
  30.                 forint c = SIFT_IMG_BORDER; c < cols-SIFT_IMG_BORDER; c++)
  31.                 {
  32.                     int val = currptr[c];
  33.                     // find local extrema with pixel accuracy
  34.                     // 寻找局部极值点,DoG中每个点与其所在的立方体周围的26个点比较
  35.                     // if (val比所有都大 或者 val比所有都小)
  36.                     if( std::abs(val) > threshold &&
  37.                        ((val > 0 && val >= currptr[c-1] && val >= currptr[c+1] &&
  38.                          val >= currptr[c-step-1] && val >= currptr[c-step] &&
  39.                          val >= currptr[c-step+1] && val >= currptr[c+step-1] &&
  40.                          val >= currptr[c+step] && val >= currptr[c+step+1] &&
  41.                          val >= nextptr[c] && val >= nextptr[c-1] &&
  42.                          val >= nextptr[c+1] && val >= nextptr[c-step-1] &&
  43.                          val >= nextptr[c-step] && val >= nextptr[c-step+1] &&
  44.                          val >= nextptr[c+step-1] && val >= nextptr[c+step] &&
  45.                          val >= nextptr[c+step+1] && val >= prevptr[c] &&
  46.                          val >= prevptr[c-1] && val >= prevptr[c+1] &&
  47.                          val >= prevptr[c-step-1] && val >= prevptr[c-step] &&
  48.                          val >= prevptr[c-step+1] && val >= prevptr[c+step-1] &&
  49.                          val >= prevptr[c+step] && val >= prevptr[c+step+1]) ||
  50.                         (val < 0 && val <= currptr[c-1] && val <= currptr[c+1] &&
  51.                          val <= currptr[c-step-1] && val <= currptr[c-step] &&
  52.                          val <= currptr[c-step+1] && val <= currptr[c+step-1] &&
  53.                          val <= currptr[c+step] && val <= currptr[c+step+1] &&
  54.                          val <= nextptr[c] && val <= nextptr[c-1] &&
  55.                          val <= nextptr[c+1] && val <= nextptr[c-step-1] &&
  56.                          val <= nextptr[c-step] && val <= nextptr[c-step+1] &&
  57.                          val <= nextptr[c+step-1] && val <= nextptr[c+step] &&
  58.                          val <= nextptr[c+step+1] && val <= prevptr[c] &&
  59.                          val <= prevptr[c-1] && val <= prevptr[c+1] &&
  60.                          val <= prevptr[c-step-1] && val <= prevptr[c-step] &&
  61.                          val <= prevptr[c-step+1] && val <= prevptr[c+step-1] &&
  62.                          val <= prevptr[c+step] && val <= prevptr[c+step+1])))
  63.                     {
  64.                         int r1 = r, c1 = c, layer = i;
  65.                         // 关键点精确定位
  66.                         if( !adjustLocalExtrema(dog_pyr, kpt, o, layer, r1, c1,
  67.                                                 nOctaveLayers, (float)contrastThreshold,
  68.                                                 (float)edgeThreshold, (float)sigma) )
  69.                             continue;
  70.                         float scl_octv = kpt.size*0.5f/(1 << o);
  71.                         // 计算梯度直方图
  72.                         float omax = calcOrientationHist(
  73.                             gauss_pyr[o*(nOctaveLayers+3) + layer],
  74.                             Point(c1, r1),
  75.                             cvRound(SIFT_ORI_RADIUS * scl_octv),
  76.                             SIFT_ORI_SIG_FCTR * scl_octv,
  77.                             hist, n);
  78.                         float mag_thr = (float)(omax * SIFT_ORI_PEAK_RATIO);
  79.                         forint j = 0; j < n; j++ )
  80.                         {
  81.                             int l = j > 0 ? j – 1 : n – 1;
  82.                             int r2 = j < n-1 ? j + 1 : 0;
  83.                             if( hist[j] > hist[l]  &&  hist[j] > hist[r2]  &&  hist[j] >= mag_thr )
  84.                             {
  85.                                 float bin = j + 0.5f * (hist[l]-hist[r2]) /
  86.                                 (hist[l] – 2*hist[j] + hist[r2]);
  87.                                 bin = bin < 0 ? n + bin : bin >= n ? bin – n : bin;
  88.                                 kpt.angle = (float)((360.f/n) * bin);
  89.                                 keypoints.push_back(kpt);
  90.                             }
  91.                         }
  92.                     }
  93.                 }
  94.             }
  95.         }
  96. }

删除边缘效应

除了DoG响应较低的点,还有一些响应较强的点也不是稳定的特征点。DoG对图像中的边缘有较强的响应值,所以落在图像边缘的点也不是稳定的特征点。
一个平坦的DoG响应峰值在横跨边缘的地方有较大的主曲率,而在垂直边缘的地方有较小的主曲率。主曲率可以通过2×2的Hessian矩阵H求出:
D值可以通过求临近点差分得到。H的特征值与D的主曲率成正比,具体可参见Harris角点检测算法。
为了避免求具体的值,我们可以通过H将特征值的比例表示出来。令为最大特征值,为最小特征值,那么:
Tr(H)表示矩阵H的迹,Det(H)表示H的行列式。
表示最大特征值与最小特征值的比值,则有:
上式与两个特征值的比例有关。随着主曲率比值的增加,也会增加。我们只需要去掉比率大于一定值的特征点。Lowe论文中去掉r=10的点。
  1. // Interpolates a scale-space extremum’s location and scale to subpixel
  2. // accuracy to form an image feature.  Rejects features with low contrast.
  3. // Based on Section 4 of Lowe’s paper.
  4. // 特征点精确定位
  5. static bool adjustLocalExtrema( const vector<Mat>& dog_pyr, KeyPoint& kpt, int octv,
  6.                                 int& layer, int& r, int& c, int nOctaveLayers,
  7.                                 float contrastThreshold, float edgeThreshold, float sigma )
  8. {
  9.     const float img_scale = 1.f/(255*SIFT_FIXPT_SCALE);
  10.     const float deriv_scale = img_scale*0.5f;
  11.     const float second_deriv_scale = img_scale;
  12.     const float cross_deriv_scale = img_scale*0.25f;
  13.     float xi=0, xr=0, xc=0, contr;
  14.     int i = 0;
  15.     //三维子像元插值
  16.     for( ; i < SIFT_MAX_INTERP_STEPS; i++ )
  17.     {
  18.         int idx = octv*(nOctaveLayers+2) + layer;
  19.         const Mat& img = dog_pyr[idx];
  20.         const Mat& prev = dog_pyr[idx-1];
  21.         const Mat& next = dog_pyr[idx+1];
  22.         Vec3f dD((img.at<short>(r, c+1) – img.at<short>(r, c-1))*deriv_scale,
  23.                  (img.at<short>(r+1, c) – img.at<short>(r-1, c))*deriv_scale,
  24.                  (next.at<short>(r, c) – prev.at<short>(r, c))*deriv_scale);
  25.         float v2 = (float)img.at<short>(r, c)*2;
  26.         float dxx = (img.at<short>(r, c+1) +
  27.                 img.at<short>(r, c-1) – v2)*second_deriv_scale;
  28.         float dyy = (img.at<short>(r+1, c) +
  29.                 img.at<short>(r-1, c) – v2)*second_deriv_scale;
  30.         float dss = (next.at<short>(r, c) +
  31.                 prev.at<short>(r, c) – v2)*second_deriv_scale;
  32.         float dxy = (img.at<short>(r+1, c+1) –
  33.                 img.at<short>(r+1, c-1) – img.at<short>(r-1, c+1) +
  34.                 img.at<short>(r-1, c-1))*cross_deriv_scale;
  35.         float dxs = (next.at<short>(r, c+1) –
  36.                 next.at<short>(r, c-1) – prev.at<short>(r, c+1) +
  37.                 prev.at<short>(r, c-1))*cross_deriv_scale;
  38.         float dys = (next.at<short>(r+1, c) –
  39.                 next.at<short>(r-1, c) – prev.at<short>(r+1, c) +
  40.                 prev.at<short>(r-1, c))*cross_deriv_scale;
  41.         Matx33f H(dxx, dxy, dxs,
  42.                   dxy, dyy, dys,
  43.                   dxs, dys, dss);
  44.         Vec3f X = H.solve(dD, DECOMP_LU);
  45.         xi = -X[2];
  46.         xr = -X[1];
  47.         xc = -X[0];
  48.         if( std::abs( xi ) < 0.5f  &&  std::abs( xr ) < 0.5f  &&  std::abs( xc ) < 0.5f )
  49.             break;
  50.         //将找到的极值点对应成像素(整数)
  51.         c += cvRound( xc );
  52.         r += cvRound( xr );
  53.         layer += cvRound( xi );
  54.         if( layer < 1 || layer > nOctaveLayers ||
  55.            c < SIFT_IMG_BORDER || c >= img.cols – SIFT_IMG_BORDER  ||
  56.            r < SIFT_IMG_BORDER || r >= img.rows – SIFT_IMG_BORDER )
  57.             return false;
  58.     }
  59.     /* ensure convergence of interpolation */
  60.     // SIFT_MAX_INTERP_STEPS:插值最大步数,避免插值不收敛,程序中默认为5
  61.     if( i >= SIFT_MAX_INTERP_STEPS )
  62.         return false;
  63.     {
  64.         int idx = octv*(nOctaveLayers+2) + layer;
  65.         const Mat& img = dog_pyr[idx];
  66.         const Mat& prev = dog_pyr[idx-1];
  67.         const Mat& next = dog_pyr[idx+1];
  68.         Matx31f dD((img.at<short>(r, c+1) – img.at<short>(r, c-1))*deriv_scale,
  69.                    (img.at<short>(r+1, c) – img.at<short>(r-1, c))*deriv_scale,
  70.                    (next.at<short>(r, c) – prev.at<short>(r, c))*deriv_scale);
  71.         float t = dD.dot(Matx31f(xc, xr, xi));
  72.         contr = img.at<short>(r, c)*img_scale + t * 0.5f;
  73.         if( std::abs( contr ) * nOctaveLayers < contrastThreshold )
  74.             return false;
  75.         /* principal curvatures are computed using the trace and det of Hessian */
  76.        //利用Hessian矩阵的迹和行列式计算主曲率的比值
  77.        float v2 = img.at<short>(r, c)*2.f;
  78.         float dxx = (img.at<short>(r, c+1) +
  79.                 img.at<short>(r, c-1) – v2)*second_deriv_scale;
  80.         float dyy = (img.at<short>(r+1, c) +
  81.                 img.at<short>(r-1, c) – v2)*second_deriv_scale;
  82.         float dxy = (img.at<short>(r+1, c+1) –
  83.                 img.at<short>(r+1, c-1) – img.at<short>(r-1, c+1) +
  84.                 img.at<short>(r-1, c-1)) * cross_deriv_scale;
  85.         float tr = dxx + dyy;
  86.         float det = dxx * dyy – dxy * dxy;
  87.         //这里edgeThreshold可以在调用SIFT()时输入;
  88.         //其实代码中定义了 static const float SIFT_CURV_THR = 10.f 可以直接使用
  89.         if( det <= 0 || tr*tr*edgeThreshold >= (edgeThreshold + 1)*(edgeThreshold + 1)*det )
  90.             return false;
  91.     }
  92.     kpt.pt.x = (c + xc) * (1 << octv);
  93.     kpt.pt.y = (r + xr) * (1 << octv);
  94.     kpt.octave = octv + (layer << 8) + (cvRound((xi + 0.5)*255) << 16);
  95.     kpt.size = sigma*powf(2.f, (layer + xi) / nOctaveLayers)*(1 << octv)*2;
  96.     return true;
  97. }

至此,SIFT第二步就完成了。参见《SIFT原理与源码分析

转载注明来源:CV视觉网 » 【OpenCV】SIFT原理与源码分析:关键点搜索与定位

分享到:更多 ()
扫描二维码,给作者 打赏
pay_weixinpay_weixin

请选择你看完该文章的感受:

0不错 0超赞 0无聊 0扯淡 0不解 0路过

评论 3

评论前必须登录!