计算机视觉
图像处理

【OpenCV】SIFT原理与源码分析:关键点描述

文章目录

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

由前一篇《方向赋值》,为找到的关键点即SIFT特征点赋了值,包含位置、尺度和方向的信息。接下来的步骤是关键点描述,即用用一组向量将这个关键点描述出来,这个描述子不但包括关键点,也包括关键点周围对其有贡献的像素点。用来作为目标匹配的依据(所以描述子应该有较高的独特性,以保证匹配率),也可使关键点具有更多的不变特性,如光照变化、3D视点变化等。

SIFT描述子h(x,y,θ)是对关键点附近邻域内高斯图像梯度统计的结果,是一个三维矩阵,但通常用一个矢量来表示。矢量通过对三维矩阵按一定规律排列得到。

描述子采样区域

特 征描述子与关键点所在尺度有关,因此对梯度的求取应在特征点对应的高斯图像上进行。将关键点附近划分成d×d个子区域,每个子区域尺寸为mσ个像元 (d=4,m=3,σ为尺特征点的尺度值)。考虑到实际计算时需要双线性插值,故计算的图像区域为mσ(d+1),再考虑旋转,则实际计算的图像区域为,如下图所示:

源码

  1.    Point pt(cvRound(ptf.x), cvRound(ptf.y));
  2. //计算余弦,正弦,CV_PI/180:将角度值转化为幅度值
  3.    float cos_t = cosf(ori*(float)(CV_PI/180));
  4.    float sin_t = sinf(ori*(float)(CV_PI/180));
  5.    float bins_per_rad = n / 360.f;
  6.    float exp_scale = -1.f/(d * d * 0.5f); //d:SIFT_DESCR_WIDTH 4    
  7.    float hist_width = SIFT_DESCR_SCL_FCTR * scl;  // SIFT_DESCR_SCL_FCTR: 3 
  8.                                                // scl: size*0.5f
  9. // 计算图像区域半径mσ(d+1)/2*sqrt(2)
  10. // 1.4142135623730951f 为根号2
  11.    int radius = cvRound(hist_width * 1.4142135623730951f * (d + 1) * 0.5f);
  12.    cos_t /= hist_width;
  13.    sin_t /= hist_width;

区域坐标轴旋转

为了保证特征矢量具有旋转不变性,要以特征点为中心,在附近邻域内旋转θ角,即旋转为特征点的方向。
旋转后区域内采样点新的坐标为:

源码

  1. //计算采样区域点坐标旋转
  2.     for( i = -radius, k = 0; i <= radius; i++ )
  3.         for( j = -radius; j <= radius; j++ )
  4.         {
  5.             /*
  6.              Calculate sample’s histogram array coords rotated relative to ori.
  7.              Subtract 0.5 so samples that fall e.g. in the center of row 1 (i.e.
  8.              r_rot = 1.5) have full weight placed in row 1 after interpolation.
  9.              */
  10.             float c_rot = j * cos_t – i * sin_t;
  11.             float r_rot = j * sin_t + i * cos_t;
  12.             float rbin = r_rot + d/2 – 0.5f;
  13.             float cbin = c_rot + d/2 – 0.5f;
  14.             int r = pt.y + i, c = pt.x + j;
  15.             if( rbin > -1 && rbin < d && cbin > -1 && cbin < d &&
  16.                r > 0 && r < rows – 1 && c > 0 && c < cols – 1 )
  17.             {
  18.                 float dx = (float)(img.at<short>(r, c+1) – img.at<short>(r, c-1));
  19.                 float dy = (float)(img.at<short>(r-1, c) – img.at<short>(r+1, c));
  20.                 X[k] = dx; Y[k] = dy; RBin[k] = rbin; CBin[k] = cbin;
  21.                 W[k] = (c_rot * c_rot + r_rot * r_rot)*exp_scale;
  22.                 k++;
  23.             }
  24.         }

计算采样区域梯度直方图

将旋转后区域划分为d×d个子区域(每个区域间隔为mσ像元),在子区域内计算8个方向的梯度直方图,绘制每个方向梯度方向的累加值,形成一个种子点。
与求主方向不同的是,此时,每个子区域梯度方向直方图将0°~360°划分为8个方向区间,每个区间为45°。即每个种子点有8个方向区间的梯度强度信息。由于存在d×d,即4×4个子区域,所以最终共有4×4×8=128个数据,形成128维SIFT特征矢量。
对特征矢量需要加权处理,加权采用mσd/2的标准高斯函数。为了除去光照变化影响,还有一步归一化处理。

源码

  1. //计算梯度直方图
  2.     for( k = 0; k < len; k++ )
  3.     {
  4.         float rbin = RBin[k], cbin = CBin[k];
  5.         float obin = (Ori[k] – ori)*bins_per_rad;
  6.         float mag = Mag[k]*W[k];
  7.         int r0 = cvFloor( rbin );
  8.         int c0 = cvFloor( cbin );
  9.         int o0 = cvFloor( obin );
  10.         rbin -= r0;
  11.         cbin -= c0;
  12.         obin -= o0;
  13.         //n为SIFT_DESCR_HIST_BINS:8,即将360°分为8个区间
  14.         if( o0 < 0 )
  15.             o0 += n;
  16.         if( o0 >= n )
  17.             o0 -= n;
  18.         // histogram update using tri-linear interpolation
  19.         // 双线性插值
  20.         float v_r1 = mag*rbin, v_r0 = mag – v_r1;
  21.         float v_rc11 = v_r1*cbin, v_rc10 = v_r1 – v_rc11;
  22.         float v_rc01 = v_r0*cbin, v_rc00 = v_r0 – v_rc01;
  23.         float v_rco111 = v_rc11*obin, v_rco110 = v_rc11 – v_rco111;
  24.         float v_rco101 = v_rc10*obin, v_rco100 = v_rc10 – v_rco101;
  25.         float v_rco011 = v_rc01*obin, v_rco010 = v_rc01 – v_rco011;
  26.         float v_rco001 = v_rc00*obin, v_rco000 = v_rc00 – v_rco001;
  27.         int idx = ((r0+1)*(d+2) + c0+1)*(n+2) + o0;
  28.         hist[idx] += v_rco000;
  29.         hist[idx+1] += v_rco001;
  30.         hist[idx+(n+2)] += v_rco010;
  31.         hist[idx+(n+3)] += v_rco011;
  32.         hist[idx+(d+2)*(n+2)] += v_rco100;
  33.         hist[idx+(d+2)*(n+2)+1] += v_rco101;
  34.         hist[idx+(d+3)*(n+2)] += v_rco110;
  35.         hist[idx+(d+3)*(n+2)+1] += v_rco111;
  36.     }

关键点描述源码

  1. // SIFT关键点特征描述
  2. // SIFT描述子是关键点领域高斯图像提取统计结果的一种表示
  3. static void calcSIFTDescriptor( const Mat& img, Point2f ptf, float ori, float scl,
  4.                                int d, int n, float* dst )
  5. {
  6.     Point pt(cvRound(ptf.x), cvRound(ptf.y));
  7.     //计算余弦,正弦,CV_PI/180:将角度值转化为幅度值
  8.     float cos_t = cosf(ori*(float)(CV_PI/180));
  9.     float sin_t = sinf(ori*(float)(CV_PI/180));
  10.     float bins_per_rad = n / 360.f;
  11.     float exp_scale = -1.f/(d * d * 0.5f); //d:SIFT_DESCR_WIDTH 4   
  12.     float hist_width = SIFT_DESCR_SCL_FCTR * scl;  // SIFT_DESCR_SCL_FCTR: 3 
  13.                                                    // scl: size*0.5f
  14.     // 计算图像区域半径mσ(d+1)/2*sqrt(2)
  15.     // 1.4142135623730951f 为根号2
  16.     int radius = cvRound(hist_width * 1.4142135623730951f * (d + 1) * 0.5f);
  17.     cos_t /= hist_width;
  18.     sin_t /= hist_width;
  19.     int i, j, k, len = (radius*2+1)*(radius*2+1), histlen = (d+2)*(d+2)*(n+2);
  20.     int rows = img.rows, cols = img.cols;
  21.     AutoBuffer<float> buf(len*6 + histlen);
  22.     float *X = buf, *Y = X + len, *Mag = Y, *Ori = Mag + len, *W = Ori + len;
  23.     float *RBin = W + len, *CBin = RBin + len, *hist = CBin + len;
  24.     //初始化直方图
  25.     for( i = 0; i < d+2; i++ )
  26.     {
  27.         for( j = 0; j < d+2; j++ )
  28.             for( k = 0; k < n+2; k++ )
  29.                 hist[(i*(d+2) + j)*(n+2) + k] = 0.;
  30.     }
  31.     //计算采样区域点坐标旋转
  32.     for( i = -radius, k = 0; i <= radius; i++ )
  33.         for( j = -radius; j <= radius; j++ )
  34.         {
  35.             /*
  36.              Calculate sample’s histogram array coords rotated relative to ori.
  37.              Subtract 0.5 so samples that fall e.g. in the center of row 1 (i.e.
  38.              r_rot = 1.5) have full weight placed in row 1 after interpolation.
  39.              */
  40.             float c_rot = j * cos_t – i * sin_t;
  41.             float r_rot = j * sin_t + i * cos_t;
  42.             float rbin = r_rot + d/2 – 0.5f;
  43.             float cbin = c_rot + d/2 – 0.5f;
  44.             int r = pt.y + i, c = pt.x + j;
  45.             if( rbin > -1 && rbin < d && cbin > -1 && cbin < d &&
  46.                r > 0 && r < rows – 1 && c > 0 && c < cols – 1 )
  47.             {
  48.                 float dx = (float)(img.at<short>(r, c+1) – img.at<short>(r, c-1));
  49.                 float dy = (float)(img.at<short>(r-1, c) – img.at<short>(r+1, c));
  50.                 X[k] = dx; Y[k] = dy; RBin[k] = rbin; CBin[k] = cbin;
  51.                 W[k] = (c_rot * c_rot + r_rot * r_rot)*exp_scale;
  52.                 k++;
  53.             }
  54.         }
  55.     len = k;
  56.     fastAtan2(Y, X, Ori, len, true);
  57.     magnitude(X, Y, Mag, len);
  58.     exp(W, W, len);
  59.     //计算梯度直方图
  60.     for( k = 0; k < len; k++ )
  61.     {
  62.         float rbin = RBin[k], cbin = CBin[k];
  63.         float obin = (Ori[k] – ori)*bins_per_rad;
  64.         float mag = Mag[k]*W[k];
  65.         int r0 = cvFloor( rbin );
  66.         int c0 = cvFloor( cbin );
  67.         int o0 = cvFloor( obin );
  68.         rbin -= r0;
  69.         cbin -= c0;
  70.         obin -= o0;
  71.         //n为SIFT_DESCR_HIST_BINS:8,即将360°分为8个区间
  72.         if( o0 < 0 )
  73.             o0 += n;
  74.         if( o0 >= n )
  75.             o0 -= n;
  76.         // histogram update using tri-linear interpolation
  77.         // 双线性插值
  78.         float v_r1 = mag*rbin, v_r0 = mag – v_r1;
  79.         float v_rc11 = v_r1*cbin, v_rc10 = v_r1 – v_rc11;
  80.         float v_rc01 = v_r0*cbin, v_rc00 = v_r0 – v_rc01;
  81.         float v_rco111 = v_rc11*obin, v_rco110 = v_rc11 – v_rco111;
  82.         float v_rco101 = v_rc10*obin, v_rco100 = v_rc10 – v_rco101;
  83.         float v_rco011 = v_rc01*obin, v_rco010 = v_rc01 – v_rco011;
  84.         float v_rco001 = v_rc00*obin, v_rco000 = v_rc00 – v_rco001;
  85.         int idx = ((r0+1)*(d+2) + c0+1)*(n+2) + o0;
  86.         hist[idx] += v_rco000;
  87.         hist[idx+1] += v_rco001;
  88.         hist[idx+(n+2)] += v_rco010;
  89.         hist[idx+(n+3)] += v_rco011;
  90.         hist[idx+(d+2)*(n+2)] += v_rco100;
  91.         hist[idx+(d+2)*(n+2)+1] += v_rco101;
  92.         hist[idx+(d+3)*(n+2)] += v_rco110;
  93.         hist[idx+(d+3)*(n+2)+1] += v_rco111;
  94.     }
  95.     // finalize histogram, since the orientation histograms are circular
  96.     // 最后确定直方图,目标方向直方图是圆的
  97.     for( i = 0; i < d; i++ )
  98.         for( j = 0; j < d; j++ )
  99.         {
  100.             int idx = ((i+1)*(d+2) + (j+1))*(n+2);
  101.             hist[idx] += hist[idx+n];
  102.             hist[idx+1] += hist[idx+n+1];
  103.             for( k = 0; k < n; k++ )
  104.                 dst[(i*d + j)*n + k] = hist[idx+k];
  105.         }
  106.     // copy histogram to the descriptor,
  107.     // apply hysteresis thresholding
  108.     // and scale the result, so that it can be easily converted
  109.     // to byte array
  110.     float nrm2 = 0;
  111.     len = d*d*n;
  112.     for( k = 0; k < len; k++ )
  113.         nrm2 += dst[k]*dst[k];
  114.     float thr = std::sqrt(nrm2)*SIFT_DESCR_MAG_THR;
  115.     for( i = 0, nrm2 = 0; i < k; i++ )
  116.     {
  117.         float val = std::min(dst[i], thr);
  118.         dst[i] = val;
  119.         nrm2 += val*val;
  120.     }
  121.     nrm2 = SIFT_INT_DESCR_FCTR/std::max(std::sqrt(nrm2), FLT_EPSILON);
  122.     for( k = 0; k < len; k++ )
  123.     {
  124.         dst[k] = saturate_cast<uchar>(dst[k]*nrm2);
  125.     }
  126. }

至此SIFT描述子生成,SIFT算法也基本完成了~参见《SIFT原理与源码分析

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

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

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

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

评论 3

评论前必须登录!