Rob Hess sift源码详解（四）

这节为关键点特征描述

static void compute_descriptors( CvSeq* features, IplImage*** gauss_pyr, int d,int n )
{
  struct feature* feat;
  struct detection_data* ddata;
  double*** hist;
  int i, k = features->total;

  for( i = 0; i < k; i++ )
    {
      feat = CV_GET_SEQ_ELEM( struct feature, features, i );
      ddata = feat_detection_data( feat );
      hist = descr_hist( gauss_pyr[ddata->octv][ddata->intvl], ddata->r,
			 ddata->c, feat->ori, ddata->scl_octv, d, n );
      hist_to_descr( hist, d, n, feat );
      release_descr_hist( &hist, d );
    }
}
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;//旋转后的采样点，-0.5意思为向下取整
	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;//高斯加权分布
	      w = exp( -(c_rot * c_rot + r_rot * r_rot) / exp_denom );//计算权值
	      interp_hist_entry( hist, rbin, cbin, obin, grad_mag * w, d, n );
	    }
      }

  return hist;
}
static void interp_hist_entry( double*** hist, double rbin, double cbin,double obin, double mag, int d, int n )
{
  double d_r, d_c, d_o, v_r, v_c, v_o;
  double** row, * h;
  int r0, c0, o0, rb, cb, ob, r, c, o;

  r0 = cvFloor( rbin );
  c0 = cvFloor( cbin );
  o0 = cvFloor( obin );
  d_r = rbin - r0;
  d_c = cbin - c0;
  d_o = obin - o0;

  /*
    The entry is distributed into up to 8 bins.  Each entry into a bin
    is multiplied by a weight of 1 - d for each dimension, where d is the
    distance from the center value of the bin measured in bin units.
  */
  //计算子区域中采样点对每个种子点的贡献，最终累加在每个方向上

  for( r = 0; r <= 1; r++ )
    {
      rb = r0 + r;
      if( rb >= 0  &&  rb < d )
	{
	  v_r = mag * ( ( r == 0 )? 1.0 - d_r : d_r );
	  row = hist[rb];
	  for( c = 0; c <= 1; c++ )
	    {
	      cb = c0 + c;
	      if( cb >= 0  &&  cb < d )
		{
		  v_c = v_r * ( ( c == 0 )? 1.0 - d_c : d_c );
		  h = row[cb];
		  for( o = 0; o <= 1; o++ )
		    {
		      ob = ( o0 + o ) % n;
		      v_o = v_c * ( ( o == 0 )? 1.0 - d_o : d_o );
		      h[ob] += v_o;
		    }
		}
	    }
	}
    }
}
static void hist_to_descr( double*** hist, int d, int n, struct feature* feat )
{
  int int_val, i, r, c, o, k = 0;

  for( r = 0; r < d; r++ )
    for( c = 0; c < d; c++ )
      for( o = 0; o < n; o++ )
	feat->descr[k++] = hist[r][c][o];

  feat->d = k;
  normalize_descr( feat );//归一化特征向量

  for( i = 0; i < k; i++ )//遍历特征描述符，将超过SIFT_DESCR_MAG_THR的元素赋值为SIFT_DESCR_MAG_THR
      if( feat->descr[i] > SIFT_DESCR_MAG_THR )
         feat->descr[i] = SIFT_DESCR_MAG_THR;
  normalize_descr( feat );//归一化特征向量

  /* convert floating-point descriptor to integer valued descriptor */
  for( i = 0; i < k; i++ )//将浮点梯度值转化为整数
    {
      int_val = SIFT_INT_DESCR_FCTR * feat->descr[i];
      feat->descr[i] = MIN( 255, int_val );
    }
}
static void normalize_descr( struct feature* feat )//为去除光照变化的影响，对其进行归一化处理
{
  double cur, len_inv, len_sq = 0.0;
  int i, d = feat->d;

  for( i = 0; i < d; i++ )
    {
      cur = feat->descr[i];
      len_sq += cur*cur;
    }
  len_inv = 1.0 / sqrt( len_sq );
  for( i = 0; i < d; i++ )
    feat->descr[i] *= len_inv;
}
static int feature_cmp( void* feat1, void* feat2, void* param )
{
//比较特征点的尺度，该函数在sift算法最后用到，按照特征点的尺度对特征向量进行排序
//cvSeqSort( features,(CvCmpFunc)feature_cmp, NULL );  
  struct feature* f1 = (struct feature*) feat1;
  struct feature* f2 = (struct feature*) feat2;

  if( f1->scl < f2->scl )
    return 1;
  if( f1->scl > f2->scl )
    return -1;
  return 0;
}
static void release_descr_hist( double**** hist, int d )//释放描述符直方图
{
  int i, j;

  for( i = 0; i < d; i++)
    {
      for( j = 0; j < d; j++ )
	free( (*hist)[i][j] );
      free( (*hist)[i] );
    }
  free( *hist );
  *hist = NULL;
}
static void release_pyr( IplImage**** pyr, int octvs, int n )//释放尺度空间金字塔
{
  int i, j;
  for( i = 0; i < octvs; i++ )
    {
      for( j = 0; j < n; j++ )
	cvReleaseImage( &(*pyr)[i][j] );
      free( (*pyr)[i] );
    }
  free( *pyr );
  *pyr = NULL;
}

终