blue_lg的surf c++源码 解析（四）

［来源］ 达内    ［编辑］ 达内   ［时间］2012-08-30

由于各层大小不一致，顾在比较大小的时候，要统一到统一尺度下，在getResponse()里面有所体现，即先计算两者尺度之比，比如尺度之比为2，最上层当前点为（20，25），那么需要比较的紧邻下层点为（40,50）而不是下层的（20,25）。

（接上文） 

③根据上一步计算每层的每个点的det(Happrox)，判断极值点。

在fasthessian.cpp里面找到getIpoints()，下面开始抽取每组（共octaves组）相邻的3层（共有4层，每组有2个这样层的满足）：                            

　　在判断的时候用到了isExtremum()和interpolateExtremum()子函数，在当前cpp内找到它们，并分析。

　　大家务必注意，由于各层大小不一致，顾在比较大小的时候，要统一到统一尺度下，在getResponse()里面有所体现，即先计算两者尺度之比，比如尺度之比为2，最上层当前点为（20，25），那么需要比较的紧邻下层点为（40,50）而不是下层的（20,25）。再看若当前点为极值点，那么调用 interpolateExtremum():

　　本函数实现功能，利用插值精确计算极值点在原图像中的位置，并保存在ipt中。

      打开interpolateStep，其中deriv3D是求当前点的3的方向上的一阶导，hessian3D是求当前点的3维hessian二阶导，最后计算出3个方向的偏差值xi,xr,xc .

　　下面开始在特征点周围提取特征描述符

④ 提取特征描述符

转到surf.cpp里寻找getDescriptors()，由于upright初始化设置为false（为特征点分配主方向，并旋转找到邻域提取特征描述符），若为true,则直接提取邻域特征描述符。

　　其实两者区别在于提取特征描述符之前判断upright的时候是否需要多调用一句getOrientation()就是了。前者不考虑旋转，两幅图只是尺度有异，而后者还考虑了旋转不变性，更加全面。

 几个地方：   

1.如论文提出的采取r=6的圆域，共计包含109个像素点，大家可以算算，在此不多解释了。

2.像素点的方向由harrx和harry计算得到,相当于梯度公式里面的dx和dy，然后利用getAngle得到θ（分布在0~2∏，而不是简单的tan-1 (harrx/harry)取值在0~∏）。

?
//! Calculate Haar wavelet responses in x direction
inline float SurhaarX(int row, int column, int s)
{
  return BoxIntegral(img, row-s/2, column, s, s/2)
    -1 * BoxIntegral(img, row-s/2, column-s/2, s, s/2);
}
 
//-------------------------------------------------------
 
//! Calculate Haar wavelet responses in y direction
inline float SurhaarY(int row, int column, int s)
{
  return BoxIntegral(img, row, column-s/2, s/2, s)
    -1 * BoxIntegral(img, row-s/2, column-s/2, s/2, s);
}
 
float SurgetAngle(float X, float Y)//计算每个点的方向角
{
  if(X > 0 && Y >= 0)
    return atan(Y/X);
 
  if(X < 0 && Y >= 0)
    return pi - atan(-Y/X);
 
  if(X < 0 && Y < 0)
    return pi + atan(Y/X);
 
  if(X > 0 && Y < 0)
    return 2*pi - atan(-Y/X);
 
  return 0;
}

　　图示如右：                          harr  （黑色代表-1，白色为1，即白色区域积分值减去黑色区域积分值）

 

在getOrienation()里面resx和resy分别是上面计算出来的harrx和harry分别乘上高斯模板系数gauss25。

?
void SurgetOrientation()
{
  ......
   //将像素点按方向角分配到6个宽为60度的区间去，比如说可以将80度分配到ang1=90度，因为90-30<80<90+30
   // loop slides pi/3 window around feature point
  for(ang1 = 0; ang1 < 2*pi;  ang1+=0.15f) {
    ang2 = ( ang1+pi/3.0f > 2*pi ? ang1-5.0f*pi/3.0f : ang1+pi/3.0f);//保证ang1+60 不会大于360度
    sumX = sumY = 0.f;
    for(unsigned int k = 0; k < Ang.size(); ++k) //对于半径为6的圆邻域，这里面的Ang.size()=109
    {
      // get angle from the x-axis of the sample point
      const float & ang = Ang[k];
      // determine whether the point is within the window
      if (ang1 < ang2 && ang1 < ang && ang < ang2)
      {
        sumX+=resX[k]; 
        sumY+=resY[k];
      }
      else if (ang2 < ang1 &&
        ((ang > 0 && ang < ang2) || (ang > ang1 && ang < 2*pi) ))
      {
        sumX+=resX[k]; 
        sumY+=resY[k];
      }
    }
    // if the vector produced from this window is longer than all
    // previous vectors then this forms the new dominant direction
    if (sumX*sumX + sumY*sumY > max) //寻找一个ang1使得角度在此区间的长度最大 ，然后计算出累加的resx和resy，得出的方向角
    {
      // store largest orientation
      max = sumX*sumX + sumY*sumY;
      orientation = getAngle(sumX, sumY);
    }
  }
  // assign orientation of the dominant response vector
  ipt->orientation = orientation;
}

　　好了，终于要开始提取特征描述符了哈~.~

　　其中i和j分别取的值为-8，-3，2，7，很明显i,j确定的邻域为7-（-8）+1=16，16x16的邻域，旋转对应的在原图像的点位置为   

                               xs = fRound(x + ( -jx*scale*si + ix*scale*co));                   ys = fRound(y + ( jx*scale*co + ix*scale*si));

                               co = cos(ipt->orientation);                                                  si = sin(ipt->orientation);                        【ipt->orientation为特征点的方向角】

共有 4x4个子块（9x9=81个像素点），每个子块分别计算了其中16个dx,dy,|dx|,|dy|之和（当然还要考虑高斯滤波权重系数），则最终的特征描述符为4x4x4=64维向量。

 

main.cpp内mainImage函数内部drawIpoints(img, ipts)就不用再做解释了吧。

 

搞了一天，终于搞完了~~

 来个截图~~

谢谢大家~（全文结束，转自cnblogs,达内c++培训整理）

顶

踩

资源下载