Rob Hess sift源码详解（二）

接上回继续，这部分为关键点定位。

static CvSeq* scale_space_extrema( IplImage*** dog_pyr, int octvs, int intvls,   double contr_thr, int curv_thr,   CvMemStorage* storage )//关键点定位{  CvSeq* features;  double prelim_contr_thr = 0.5 * contr_thr / intvls;  struct feature* feat;  struct detection_data* ddata;  int o, i, r, c;  features = cvCreateSeq( 0, sizeof(CvSeq), sizeof(struct feature), storage );  for( o = 0; o < octvs; o++ )    for( i = 1; i <= intvls; i++ )      for(r = SIFT_IMG_BORDER; r < dog_pyr[o][0]->height-SIFT_IMG_BORDER; r++)for(c = SIFT_IMG_BORDER; c < dog_pyr[o][0]->width-SIFT_IMG_BORDER; c++)//加入一个长度为SIFT_IMG_BORDER的边界，目的是排除边界的影响  /* perform preliminary check on contrast */  if( ABS( pixval32f( dog_pyr[o][i], r, c ) ) > prelim_contr_thr )//不理解为什么要大于某个阈值，希望知道者赐教    if( is_extremum( dog_pyr, o, i, r, c ) )//关键点初始定位      {feat = interp_extremum(dog_pyr, o, i, r, c, intvls, contr_thr);//关键点精确定位if( feat )  {    ddata = feat_detection_data( feat );                    //去除边缘响应    if( ! is_too_edge_like( dog_pyr[ddata->octv][ddata->intvl],    ddata->r, ddata->c, curv_thr ) )//阈值判断      {cvSeqPush( features, feat );      }    else      free( ddata );    free( feat );  }      }    return features;}static void interp_step( IplImage*** dog_pyr, int octv, int intvl, int r, int c,double* xi, double* xr, double* xc ){  CvMat* dD, * H, * H_inv, X;  double x[3] = { 0 };    dD = deriv_3D( dog_pyr, octv, intvl, r, c );//求偏导  H = hessian_3D( dog_pyr, octv, intvl, r, c );//计算Hessian矩阵  H_inv = cvCreateMat( 3, 3, CV_64FC1 );  cvInvert( H, H_inv, CV_SVD );//矩阵求逆  cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );  cvGEMM( H_inv, dD, -1, NULL, 0, &X, 0 );    cvReleaseMat( &dD );  cvReleaseMat( &H );  cvReleaseMat( &H_inv );  *xi = x[2];  *xr = x[1];  *xc = x[0];}static int is_extremum( IplImage*** dog_pyr, int octv, int intvl, int r, int c )//关键点初始定位{  double val = pixval32f( dog_pyr[octv][intvl], r, c );  int i, j, k;//为寻找DoG函数的极值点，每一个像素点都要和它所有的相邻点相比较，看其是否比它的图像域和尺度域的相邻点大或者小。包括和它同尺度的8个相邻点和上下相邻尺度对应的9*2个点共计26个点比较  /* check for maximum */  if( val > 0 )    {         for( i = -1; i <= 1; i++ )      for( j = -1; j <= 1; j++ )           for( k = -1; k <= 1; k++ )                 if( val < pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )                      return 0;    }  /* check for minimum */  else    {         for( i = -1; i <= 1; i++ )      for( j = -1; j <= 1; j++ )            for( k = -1; k <= 1; k++ )                  if( val > pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )                       return 0;   }  return 1;}//关键点精确定位static struct feature* interp_extremum( IplImage*** dog_pyr, int octv,int intvl, int r, int c, int intvls,double contr_thr ){  struct feature* feat;  struct detection_data* ddata;  double xi, xr, xc, contr;  int i = 0;    while( i < SIFT_MAX_INTERP_STEPS )//lowe选择迭代SIFT_MAX_INTERP_STEPS次    {      interp_step( dog_pyr, octv, intvl, r, c, &xi, &xr, &xc );//通过三维二次函数来精确定位      if( ABS( xi ) < 0.5  &&  ABS( xr ) < 0.5  &&  ABS( xc ) < 0.5 )      //当相对插值中心的偏移量在任一维度大于0.5时代表插值中心已经偏移到其邻近点上了，故在新位置上反复插值直到收敛break;            //四舍五入叠加到原位置上      c += cvRound( xc );//cols            r += cvRound( xr );//rows            intvl += cvRound( xi );//intvls                  if( intvl < 1  ||            intvl > intvls  ||            c < SIFT_IMG_BORDER  ||            r < SIFT_IMG_BORDER  ||            c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER  ||            r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER )          {               return NULL;          }                  i++;        }        /* ensure convergence of interpolation */      if( i >= SIFT_MAX_INTERP_STEPS )           return NULL;        contr = interp_contr( dog_pyr, octv, intvl, r, c, xi, xr, xc );      if( ABS( contr ) < contr_thr / intvls )//过小的点易受噪声的干扰而变得不稳定，所以将其小于contr_thr / intvls的极值点删除。        return NULL;      feat = new_feature();      ddata = feat_detection_data( feat );  

    //获取特征点的精确位置（原位置上加上拟合的偏移量）以及尺度    feat->img_pt.x = feat->x = ( c + xc ) * pow( 2.0, octv );      feat->img_pt.y = feat->y = ( r + xr ) * pow( 2.0, octv );      ddata->r = r;      ddata->c = c;      ddata->octv = octv;      ddata->intvl = intvl;      ddata->subintvl = xi;      return feat;}static CvMat* deriv_3D( IplImage*** dog_pyr, int octv, int intvl, int r, int c ){  CvMat* dI;  double dx, dy, ds;  //通过有限差分求导代替微分方程中的微分，将连续变量离散化  dx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) - pixval32f( dog_pyr[octv][intvl], r, c-1 ) ) / 2.0;  dy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) - pixval32f( dog_pyr[octv][intvl], r-1, c ) ) / 2.0;  ds = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) - pixval32f( dog_pyr[octv][intvl-1], r, c ) ) / 2.0;    dI = cvCreateMat( 3, 1, CV_64FC1 );  cvmSet( dI, 0, 0, dx );  cvmSet( dI, 1, 0, dy );  cvmSet( dI, 2, 0, ds );  return dI;}static CvMat* hessian_3D( IplImage*** dog_pyr, int octv, int intvl, int r,int c ){  CvMat* H;  double v, dxx, dyy, dss, dxy, dxs, dys;    v = pixval32f( dog_pyr[octv][intvl], r, c );  dxx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) +   pixval32f( dog_pyr[octv][intvl], r, c-1 ) - 2 * v );  dyy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) +  pixval32f( dog_pyr[octv][intvl], r-1, c ) - 2 * v );  dss = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) +  pixval32f( dog_pyr[octv][intvl-1], r, c ) - 2 * v );  dxy = ( pixval32f( dog_pyr[octv][intvl], r+1, c+1 ) -  pixval32f( dog_pyr[octv][intvl], r+1, c-1 ) -  pixval32f( dog_pyr[octv][intvl], r-1, c+1 ) +  pixval32f( dog_pyr[octv][intvl], r-1, c-1 ) ) / 4.0;  dxs = ( pixval32f( dog_pyr[octv][intvl+1], r, c+1 ) -  pixval32f( dog_pyr[octv][intvl+1], r, c-1 ) -  pixval32f( dog_pyr[octv][intvl-1], r, c+1 ) +  pixval32f( dog_pyr[octv][intvl-1], r, c-1 ) ) / 4.0;  dys = ( pixval32f( dog_pyr[octv][intvl+1], r+1, c ) -  pixval32f( dog_pyr[octv][intvl+1], r-1, c ) -  pixval32f( dog_pyr[octv][intvl-1], r+1, c ) +  pixval32f( dog_pyr[octv][intvl-1], r-1, c ) ) / 4.0;    H = cvCreateMat( 3, 3, CV_64FC1 );  cvmSet( H, 0, 0, dxx );  cvmSet( H, 0, 1, dxy );  cvmSet( H, 0, 2, dxs );  cvmSet( H, 1, 0, dxy );  cvmSet( H, 1, 1, dyy );  cvmSet( H, 1, 2, dys );  cvmSet( H, 2, 0, dxs );  cvmSet( H, 2, 1, dys );  cvmSet( H, 2, 2, dss );  /// Ixx  Ixy  Ixs \   //| Ixy  Iyy  Iys |   //\ Ixs  Iys  Iss /  return H;}static double interp_contr( IplImage*** dog_pyr, int octv, int intvl, int r,    int c, double xi, double xr, double xc ){  CvMat* dD, X, T;  double t[1], x[3] = { xc, xr, xi };  cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );  cvInitMatHeader( &T, 1, 1, CV_64FC1, t, CV_AUTOSTEP );  dD = deriv_3D( dog_pyr, octv, intvl, r, c );  cvGEMM( dD, &X, 1, NULL, 0, &T,  CV_GEMM_A_T );  cvReleaseMat( &dD );  return pixval32f( dog_pyr[octv][intvl], r, c ) + t[0] * 0.5;}static struct feature* new_feature( void ){  struct feature* feat;  struct detection_data* ddata;  feat = malloc( sizeof( struct feature ) );  memset( feat, 0, sizeof( struct feature ) );  ddata = malloc( sizeof( struct detection_data ) );  memset( ddata, 0, sizeof( struct detection_data ) );  feat->feature_data = ddata;  feat->type = FEATURE_LOWE;  return feat;}static int is_too_edge_like( IplImage* dog_img, int r, int c, int curv_thr )//消除边缘响应{  double d, dxx, dyy, dxy, tr, det;  /* principal curvatures are computed using the trace and det of Hessian */  d = pixval32f(dog_img, r, c);  dxx = pixval32f( dog_img, r, c+1 ) + pixval32f( dog_img, r, c-1 ) - 2 * d;  dyy = pixval32f( dog_img, r+1, c ) + pixval32f( dog_img, r-1, c ) - 2 * d;  dxy = ( pixval32f(dog_img, r+1, c+1) - pixval32f(dog_img, r+1, c-1) -  pixval32f(dog_img, r-1, c+1) + pixval32f(dog_img, r-1, c-1) ) / 4.0;  tr = dxx + dyy;//矩阵H的对角线元素之和  det = dxx * dyy - dxy * dxy;//矩阵H的行列式  /* negative determinant -> curvatures have different signs; reject feature */  if( det <= 0 )    return 1;  if( tr * tr / det < ( curv_thr + 1.0 )*( curv_thr + 1.0 ) / curv_thr )//检测主曲率是否小于阈值curv_thr    return 0;  return 1;}

未完待续