sift：（第三步）为关键点赋予方向

上一步：Sift：(第二步）关键点位置的确定_zzz_zzzz_的博客-CSDN博客

目的：为了使关键点描述符具有旋转不变性，需要利用图像的局部特征为每一个关键点分配一个基准方向。

第二步已经找到了关键点，并在此之前已经确定了各个关键点的属性值。下面要做的是，为关键点赋予方向。首先，找到关键点所在尺度相近的高斯金字塔图像。然后在关键点坐标位置附近的一定邻域内（如上图左），计算各像素点的梯度幅值和梯度方向。并以关键点位置为核心，做1.5*的高斯滤波（使得关键点位置的赋值占比重大，距离关键点越远的位置的梯度赋值占比越小），这样就得到了经高斯滤波之后的赋值。至于，这个领域有多大，根据高斯滤波的3原则，领域的大小为3*1.5*。（注意，这里的sigma为*）

将我们计算得到的梯度方向（0~360度），分为36个柱（bins)，每柱10度的区间范围。关键点邻域内的每一个点都有一个值，也都有一个bin。然后就像投票一样，每个点的值加到其对应的bin中，得到上图右边的直方图（图中仅画了8个方向的直方图，实际上是36个）。

全都“投票”之后，选择对应最大值的方向为该关键点的主方向。如果另有方向的值大于主方向值的80%，那么选择该方向为该关键点的辅方向。

以上为该步骤的理论知识。

Lowe,D 的SIFT代码 

对于每一个特征点，先把它从特征值序列features中弹出，计算得它的36个方向的梯度直方图（hist[36]) ，对该梯度直方图进行平滑。再使用“打擂台”算法得到hist[36]中的最大值，也就是主方向的值。然后，将hist[36]中大于等于【80%*主方向的值】视做主方向/辅方向，分别再压入features序列中。

对于features前面弹出，待增加方向后，再压入。

//第三步：为关键点赋予方向的大框架 
/*
  Computes a canonical orientation for each image feature in an array.  Based
  on Section 5 of Lowe's paper.  This function adds features to the array when
  there is more than one dominant orientation at a given feature location.

  @param features an array of image features（关键点/特征点序列） 
  @param gauss_pyr Gaussian scale space pyramid（高斯金字塔） 
*/
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;//n是关键点序列中关键点的总数 

  for( i = 0; i < n; i++ )//对于每一个关键点的循环操作 
    {
      feat = malloc( sizeof( struct feature ) );//分配内存空间 
      cvSeqPopFront( features, feat );//将featutes序列头部元素弹出，放入feat 
      ddata = feat_detection_data( feat );
      //得到梯度直方图 
      hist = ori_hist( gauss_pyr[ddata->octv][ddata->intvl],
		       ddata->r, ddata->c, SIFT_ORI_HIST_BINS,
		       cvRound( SIFT_ORI_RADIUS * ddata->scl_octv ),
		       SIFT_ORI_SIG_FCTR * ddata->scl_octv );//调用函数ori_hist
	  //平滑梯度直方图 
      for( j = 0; j < SIFT_ORI_SMOOTH_PASSES; j++ )  //SIFT_ORI_SMOOTH_PASSES为2 
	        smooth_ori_hist( hist, SIFT_ORI_HIST_BINS );//调用函数smooth_ori_hist 
	  //得到主方向 
      omax = dominant_ori( hist, SIFT_ORI_HIST_BINS );//调用函数dominant_ori，返回36个bin中的最大幅值 
      //添加主方向，辅方向 
      add_good_ori_features( features, hist, SIFT_ORI_HIST_BINS,
			     omax * SIFT_ORI_PEAK_RATIO, feat );//调用函数add_good_ori_features
      free( ddata );
      free( feat );
      free( hist );
    }
}

生成梯度直方图，实际上是hist[36]，里面存的是36个大致方向的赋值（已平滑）之和

//生成梯度直方图 
/*
  Computes a gradient orientation histogram at a specified pixel.

  @param img image（ 关键点所在的对应高斯金字塔的图片）gauss_pyr[ddata->octv][ddata->intvl]
  @param r pixel row（关键点的x坐标） ddata->r
  @param c pixel col（关键点的y坐标） ddata->c
  @param n number of histogram bins（分为36个方向，梯度直方图的36个柱子）SIFT_ORI_HIST_BINS
  注意：ddata->scl_octv = sigma * pow( 2.0, intvl / intvls );
  @param rad radius of region over which histogram is computed（邻域：3*1.5倍的对应尺度） cvRound( SIFT_ORI_RADIUS * ddata->scl_octv )
  @param sigma std for Gaussian weighting of histogram entries（对梯度幅值进行高斯滤波的系数：1.5倍的对应尺度）SIFT_ORI_SIG_FCTR * ddata->scl_octv

  @return Returns an n-element array containing an orientation histogram
    representing orientations between 0 and 2 PI.
*/
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;//2*sigma^2 
  for( i = -rad; i <= rad; i++ )
    for( j = -rad; j <= rad; j++ )//在关键点的领域半径内 
      if( calc_grad_mag_ori( img, r + i, c + j, &mag, &ori ) )//调用函数calc_grad_mag_ori，求梯度赋值和方向，求出返回1，用指针带回 
	{
	  w = exp( -( i*i + j*j ) / exp_denom );//生成领域内像素的高斯权重 
	  bin = cvRound( n * ( ori + CV_PI ) / PI2 );//求当前梯度方向对应的bin 
	  bin = ( bin < n )? bin : 0;
	  hist[bin] += w * mag;//相当于“投票” 
	}

  return hist;
}

计算当前像素点的梯度幅值和方向

//求梯度方向和梯度赋值 
/*
  Calculates the gradient magnitude and orientation at a given pixel.

  @param img image(关键点所在图片） 
  @param r pixel row（关键点的x坐标）
  @param c pixel col（关键点的y坐标）
  @param mag output as gradient magnitude at pixel (r,c) （梯度赋值） 
  @param ori output as gradient orientation at pixel (r,c)（梯度方向） 

  @return Returns 1 if the specified pixel is a valid one and sets mag and
    ori accordingly; otherwise returns 0
*/
static int calc_grad_mag_ori( IplImage* img, int r, int c, double* mag,
			      double* ori )
{
  double dx, dy;

  if( r > 0  &&  r < img->height - 1  &&  c > 0  &&  c < img->width - 1 )
    {
      dx = pixval32f( img, r, c+1 ) - pixval32f( img, r, c-1 );
      dy = pixval32f( img, r-1, c ) - pixval32f( img, r+1, c );
      *mag = sqrt( dx*dx + dy*dy );//梯度赋值 
      *ori = atan2( dy, dx );//梯度方向，第一象限(0,PI/2) 第二象限(PI/2,PI) 第三象限（-PI,-PI/2) 第四象限(-PI/2,0) 
      return 1;
    }

  else
    return 0;
}

平滑梯度直方图

//平滑梯度直方图
//平滑目的：防止某个方向因为噪声干扰等原因突变 
/*
  Gaussian smooths an orientation histogram.

  @param hist an orientation histogram（梯度直方图） 
  @param n number of bins（36） 
*/
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] );//25%的前一个，50%自身，25%后一个
	  // hist[0]=0.25*hist[35]+0.5*hist[0]+0.25*hist[1]
      prev = tmp;
    }
}

获取主方向的对应hist中的值

//主方向对应值的获取 
//应用“打擂台”算法得最大值 
/*
  Finds the magnitude of the dominant orientation in a histogram

  @param hist an orientation histogram（平滑后的梯度直方图） 
  @param n number of bins（36） 

  @return Returns the value of the largest bin in hist
*/
static double dominant_ori( double* hist, int n )
{
  double omax;
  int maxbin, i;

  omax = hist[0];
  maxbin = 0;
  for( i = 1; i < n; i++ )
    if( hist[i] > omax )
      {
	      omax = hist[i];
	      maxbin = i;
      }
  return omax;
}

添加特征点的主方向/辅方向，并将其添加到features中 

如果一个特征点有一个主方向，一个辅方向，那么会有两个new_feature被压入features序列

//添加主方向、辅方向 
/*
  Adds features to an array for every orientation in a histogram greater than
  a specified threshold.

  @param features new features are added to the end of this array（特征点序列） 
  @param hist orientation histogram（梯度直方图） 
  @param n number of bins in hist（36） 
  @param mag_thr new features are added for entries in hist greater than this（主方向的值*0.8） 
  @param feat new features are clones of this with different orientations（当前特征点：已经 
*/
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;//前一个bin，当i=0 l=36；i=1 l=0；i=2 l=1;... 
      r = ( i + 1 ) % n;//后一个bin，当i=0 r=1；i=1 r=2;...i=35 r=0 
      
      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 = ( bin < 0 )? n + bin : ( bin >= n )? bin - n : bin;//bin控制在0-36之间 
	  new_feat = clone_feature( feat );//调用函数clone_feature 
	  new_feat->ori = ( ( PI2 * bin ) / n ) - CV_PI;//根据bin求梯度方向ori，bin = cvRound( n * ( ori + CV_PI ) / PI2 )反过来 
	  cvSeqPush( features, new_feat );//把new_feature压入关键点序列（前面有弹出） 
	  free( new_feat );//释放空间 
	}
    }
}

抛物线插值 

公式:

//区间变为值 
//原先的方向是一个角度区间，抛物线插值后，用顶点的横坐标作为主方向的值 
/*
  Interpolates a histogram peak from left, center, and right values
*/
#define interp_hist_peak( l, c, r ) ( 0.5 * ((l)-(r)) / ((l) - 2.0*(c) + (r)) )//代入抛物线插值的公式 

特征点复制 

//特征点的复制 
/*
  Makes a deep copy of a feature

  @param feat feature to be cloned（当前特征点） 

  @return Returns a deep copy of feat
*/
static struct feature* clone_feature( struct feature* feat )
{
  struct feature* new_feat;
  struct detection_data* ddata;

  new_feat = new_feature();
  ddata = feat_detection_data( new_feat );
  memcpy( new_feat, feat, sizeof( struct feature ) );//分配空间，完成复制 
  memcpy( ddata, feat_detection_data(feat), sizeof( struct detection_data ) );//分配空间，完成复制 
  new_feat->feature_data = ddata;

  return new_feat;
}

下一步：Sift：（第四步）为关键点构造描述符_zzz_zzzz_的博客-CSDN博客