【特征匹配】SIFT原理与C源代码剖析(3)

4.1

在DOG尺度空间上。每组有S+2层图像。每一组都从第二层開始每个像素点都要与它相邻的像素点比較，看是否比它在图像域或尺度域的全部点的值大或者小。

与它同尺度的相邻像素点有8个，上下相邻尺度的点共同拥有2×9=18。共同拥有26个像素点。也就在一个3×3的立方体内进行。搜索的过程是第二层開始到倒数第二层结束，共检測了octave组。每组S层。

/*
  Determines whether a pixel is a scale-space extremum by comparing it to it's
  3x3x3 pixel neighborhood.
*/
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;
  /* 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;
}

4.2

以上的极值点搜索是在离散空间进行的，极值点不真正意义上的极值点。通过对空间尺度函数拟合。能够得到亚像素级像素点坐标。

尺度空间的Taylor展开式：

118.98.232.202:8080@http$6&141,1016,1016#欧洲 unknown。bikon 1016-17x28/18x28/19x28/20x28、bikon 1016-22x36/24x36/25x36/28x36、bikon 1016-30x48/32x48/35x48/38x48。218.246.181.241:843 http 79,1016,1016 黑龙江省大庆市 中基电信石油建设公司。

求导并令其为0，得到亚像素级：

相应的函数值为：

是一个三维矢量，矢量在不论什么一个方向上的偏移量大于0.5时，意味着已经偏离了原像素点，这种特征坐标位置须要更新或者继续插值计算。算法实现过程中，为了保证插值可以收敛，设置了最大插值次数（lowe 设置了5次）。

同一时候当时（本文阈值採用了0.04/S） ，特征点才被保留，由于响应值过小的点。easy受噪声的干扰而不稳定。

对离散空间进行函数拟合（插值）：

/*
  Performs one step of extremum interpolation.  Based on Eqn. (3) in Lowe's
  paper.
  r,c 为特征点位置，xi,xr,xc,保存三个方向的偏移量
*/
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 );    // 计算3维空间的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];
}

/*
  Interpolates a scale-space extremum's location and scale to subpixel
  accuracy to form an image feature. 
*/
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 )   //在最大迭代次数范围内进行
    {
      interp_step( dog_pyr, octv, intvl, r, c, &xi, &xr, &xc );          //插值后得到的三个方向的偏移量(xi,xr,xc)
      if( ABS( xi ) < 0.5  &&  ABS( xr ) < 0.5  &&  ABS( xc ) < 0.5 )
	break;
      
      c += cvRound( xc );    //更新位置
      r += cvRound( xr );
      intvl += cvRound( xi );
      
      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 )   //小于阈值（0.04/S）的点。则丢弃
    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;
}

为了得到稳定的特征点。要删除掉落在图像边缘上的点。

一个落在边缘上的点。能够依据主曲率计算推断。主曲率能够通过2维的 Hessian矩阵求出；

在边缘上的点，必然使得Hessian矩阵的两个特征值相差比較大。而特征值与矩阵元素有下面关系；

令α=rβ ，所以有：

我们能够推断上述公式的比值大小，大于阈值（lowe採用 r =10）的点排除。

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);                                                                             //计算Hessian 矩阵内的4个元素值
  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;                          //矩阵的迹
  det = dxx * dyy - dxy * dxy;     //矩阵的值
  /* 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 = 10
    return 0;
  return 1;
}