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

static void calc_feature_scales( CvSeq* features, double sigma, int intvls )
{
  struct feature* feat;
  struct detection_data* ddata;
  double intvl;
  int i, n;
  n = features->total;
  for( i = 0; i < n; i++ )
    {
      feat = CV_GET_SEQ_ELEM( struct feature, features, i );
      ddata = feat_detection_data( feat );
      intvl = ddata->intvl + ddata->subintvl;                        
      feat->scl = sigma * pow( 2.0, ddata->octv + intvl / intvls );      // feat->scl 保存特征点在整体上尺度
      ddata->scl_octv = sigma * pow( 2.0, intvl / intvls );     // feat->feature_data->scl__octv 保存特征点在组内的尺度，用来以下计算方向角
    }
}

这部分包含：计算邻域内梯度直方图，平滑直方图，复制特征点（有辅方向的特征点）

static void calc_feature_oris( CvSeq* features, IplImage*** gauss_pyr )  
{
  struct feature* feat;
  struct detection_data* ddata;
  double* hist;
  double omax;
  int i, j, n = features->total;
  for( i = 0; i < n; i++ )
    {
      feat = malloc( sizeof( struct feature ) );
      cvSeqPopFront( features, feat );
      ddata = feat_detection_data( feat );
      hist = ori_hist( gauss_pyr[ddata->octv][ddata->intvl],     // 计算邻域内的梯度直方图,邻域半径radius = 3*1.5*sigma;  高斯加权系数= 1.5 *sigma 
		       ddata->r, ddata->c, SIFT_ORI_HIST_BINS,
		       cvRound( SIFT_ORI_RADIUS * ddata->scl_octv ),
		       SIFT_ORI_SIG_FCTR * ddata->scl_octv );
      for( j = 0; j < SIFT_ORI_SMOOTH_PASSES; j++ )
	smooth_ori_hist( hist, SIFT_ORI_HIST_BINS );    // 对直方图平滑
      omax = dominant_ori( hist, SIFT_ORI_HIST_BINS ); // 返回直方图的主方向
      add_good_ori_features( features, hist, SIFT_ORI_HIST_BINS,//大于主方向的80%为辅方向
			     omax * SIFT_ORI_PEAK_RATIO, feat );
      free( ddata );
      free( feat );
      free( hist );
    }
}

6.1

上一步scl_octv保存了每一个特征点所在的组内尺度。以下计算特征点所在尺度内的高斯图像，以3×1.5×scl_octv为半径的区域内的全部像素点的梯度幅值与幅角；

nbins表示在一个胞元（cell）中统计梯度的方向数目，例如nbins=9时，在一个胞元内统计9个方向的梯度直方图，每个方向为180/9=20度。对梯度影像进行直方图统计，并进行概率累积，概率累积值作为 阈值将梯度影像划分为标记像素和非标记像素。hog算法的具体操作就是将图像分割成很多小的连接区域，在每个小的区域之中生产多个方向梯度直方图或像素点的边缘方向，图像所有的这些区域的方向梯度直方图或像素点边缘方向就构成了图像的描述子hogdescriptor，这个描述子即可作为识别该图像物体的特征向量。

下图描写叙述的划分为8个bin,每一个bin的宽度为45的效果图：

其次。每一个被增加直方图的幅值，要进行权重处理，权重也是採用高斯加权函数。当中高斯系数为1.5×scl_octv。通过高斯加权使特征点附近的点有较大的权重，能够弥补部分因没有仿射不变性而产生的不稳定问题；

即每一个bin值按以下的公式累加，mag是幅值，后面为权重；i，j,为偏离特征点距离：

☆程序上能够帮助你理解上面的概念：

static double* ori_hist( IplImage* img, int r, int c, int n, int rad,
			 double sigma )
{
  double* hist;
  double mag, ori, w, exp_denom, PI2 = CV_PI * 2.0;
  int bin, i, j;
  hist = calloc( n, sizeof( double ) );
  exp_denom = 2.0 * sigma * sigma;
  for( i = -rad; i <= rad; i++ )
    for( j = -rad; j <= rad; j++ )
      if( calc_grad_mag_ori( img, r + i, c + j, &mag, &ori ) )    //计算领域像素点的梯度幅值mag与方向ori
	{
	  w = exp( -( i*i + j*j ) / exp_denom );                    //高斯权重
	  bin = cvRound( n * ( ori + CV_PI ) / PI2 );
	  bin = ( bin < n )?
 bin : 0;
	  hist[bin] += w * mag;                             //直方图上累加
	}
  return hist;  //返回累加完毕的直方图
}

lowe建议对直方图进行平滑，降低突变的影响。

static void smooth_ori_hist( double* hist, int n )
{
  double prev, tmp, h0 = hist[0];
  int i;
  prev = hist[n-1];
  for( i = 0; i < n; i++ )
    {
      tmp = hist[i];
      hist[i] = 0.25 * prev + 0.5 * hist[i] + 
	0.25 * ( ( i+1 == n )? h0 : hist[i+1] );
      prev = tmp;
    }
}

使用sift算法进行识别（特征点的提取并用特征向量对特征点描述，接着当前视图的特征向量与目标对象的特征向量进行匹配）。sift/surf采用henssian矩阵获取图像局部最值还是十分稳定的，但是在求主方向阶段太过于依赖局部区域像素的梯度方向，有可能使得找到的主 方向不准确，后面的特征向量提取以及匹配都严重依赖于主方向，即使不大偏差角度也可以造成后面特征匹配的放大误差，从而匹配不成功。具体思路是：使用sift算法进行识别（特征点的提取并用特征向量对特征点描述，接着当前视图的特征向量与目标对象的特征向量进行匹配），根据识别出来的原目标和帧图像匹配关系得到变化矩阵，来显示三维物体（使用opengl来绘制），实现跟踪。

static void add_good_ori_features( CvSeq* features, double* hist, int n,
				   double mag_thr, struct feature* feat )
{
  struct feature* new_feat;
  double bin, PI2 = CV_PI * 2.0;
  int l, r, i;
  for( i = 0; i < n; i++ )                   //检查全部的方向
    {
      l = ( i == 0 )?
 n - 1 : i-1;
      r = ( i + 1 ) % n;
      
      if( hist[i] > hist[l]  &&  hist[i] > hist[r]  &&  hist[i] >= mag_thr ) //推断是不是幅方向
	{
	  bin = i + interp_hist_peak( hist[l], hist[i], hist[r] );     //插值离散处理。取得更精确的方向
	  bin = ( bin < 0 )? n + bin : ( bin >= n )?
 bin - n : bin;
	  new_feat = clone_feature( feat );       //复制特征点
	  new_feat->ori = ( ( PI2 * bin ) / n ) - CV_PI;//  为特征点方向赋值[-180,180]
	  cvSeqPush( features, new_feat );  //
	  free( new_feat );
	}
    }
}

眼下每一个特征点具有属性有位置、方向、尺度三个信息，如今要用一个向量去描写叙述这个特征点，使其具有高度的唯一特征性。

1.lowe採用了把特征点邻域划分成 d×d (lowe建议d=4) 个子区域，然后再统计每一个子区域的方向直方图（8个方向），直方图横轴有8个bin,纵轴为梯度幅值（×权重）的累加。这样描写叙述这个特征点的向量为4×4×8=128维。每一个子区域的宽度建议为3×octv，octv为组内的尺度。考虑到插值问题。特征点的邻域范围边长为3×octv×(d+1)。考虑到旋转问题,邻域的范围边长为3×octv×(d+1)×sqrt(2)。最后半径为：

2.把坐标系旋转到主方向位置。再次统计邻域内全部像素点的梯度幅值与方向，计算所在子区域。并把幅值×权重累加到这个子区域的直方图上。

算法上即统计每一个邻域的方向直方图时。所有是相对于这个特征点的主方向的方向。

假设主方向为30度,某个像素点的梯度方向为50度。这时统计到该子区域直方图上就成了20度。同一时候因为旋转，这时权重也必须是按旋转后的坐标。

计算所在的子区域的位置：

权重按高斯加权函数。系数为描写叙述子宽度的一半，即0.5d：

static double*** descr_hist( IplImage* img, int r, int c, double ori,
			     double scl, int d, int n )
{
  double*** hist;
  double cos_t, sin_t, hist_width, exp_denom, r_rot, c_rot, grad_mag,
    grad_ori, w, rbin, cbin, obin, bins_per_rad, PI2 = 2.0 * CV_PI;
  int radius, i, j;
  hist = calloc( d, sizeof( double** ) );
  for( i = 0; i < d; i++ )
    {
      hist[i] = calloc( d, sizeof( double* ) );
      for( j = 0; j < d; j++ )
	hist[i][j] = calloc( n, sizeof( double ) );
    }
  
  cos_t = cos( ori );
  sin_t = sin( ori );
  bins_per_rad = n / PI2;
  exp_denom = d * d * 0.5;
  hist_width = SIFT_DESCR_SCL_FCTR * scl;
  radius = hist_width * sqrt(2) * ( d + 1.0 ) * 0.5 + 0.5;   //计算邻域范围半径，+0.5为了取得很多其它元素
  for( i = -radius; i <= radius; i++ )
    for( j = -radius; j <= radius; j++ )
      {
	/*
	  Calculate sample's histogram array coords rotated relative to ori.
	  Subtract 0.5 so samples that fall e.g. in the center of row 1 (i.e.
	  r_rot = 1.5) have full weight placed in row 1 after interpolation.
	*/
	c_rot = ( j * cos_t - i * sin_t ) / hist_width;              //
	r_rot = ( j * sin_t + i * cos_t ) / hist_width;
	rbin = r_rot + d / 2 - 0.5;                                        //旋转后相应的x``及y''
	cbin = c_rot + d / 2 - 0.5;
	
	if( rbin > -1.0  &&  rbin < d  &&  cbin > -1.0  &&  cbin < d )
	  if( calc_grad_mag_ori( img, r + i, c + j, &grad_mag, &grad_ori ))       //计算每一个像素点的梯度方向、幅值、
	    {
	      grad_ori -= ori;              //每一个像素点相对于特征点的梯度方向
	      while( grad_ori < 0.0 )
		grad_ori += PI2;
	      while( grad_ori >= PI2 )
		grad_ori -= PI2;
	      
	      obin = grad_ori * bins_per_rad;     //像素梯度方向转换成8个方向
	      w = exp( -(c_rot * c_rot + r_rot * r_rot) / exp_denom );     //每一个子区域内像素点，相应的权重
	      interp_hist_entry( hist, rbin, cbin, obin, grad_mag * w, d, n );   //为了减小突变影响对附近三个bin值三线性插值处理
	    }
      }
  return hist;
}