Sift：（第四步）为关键点构造描述符

zzz_zzzz_

已于 2023-09-08 23:17:22 修改

阅读量692

点赞数 2

分类专栏： sift 文章标签：人工智能

于 2023-08-17 07:33:54 首次发布

本文链接：https://blog.youkuaiyun.com/zzz_zzzz_/article/details/132270450

版权

sift 专栏收录该内容

7 篇文章

订阅专栏

上一步：sift：（第三步）为关键点赋予方向_zzz_zzzz_的博客-优快云博客

我觉得最难理解的是这个第四步，就是看似简单，细究的话，还需要好好想想。

先说目的：通过前面的三大步骤，我们找到了关键点并且知道了它们的位置、尺度以及方向信息。下面就是要利用这些信息来构造一种特征描述方法，进而对关键点进行描述。描述是为了后面特征点匹配的任务。剧透一下，这个描述符就是一个【128维的向量】，叫做关键点描述子（Sift描述子），其中包含了特征点周围像素的信息，实际上是对关键点周围区域进行的信息描述。

看看具体怎么做的？

旋转到主方向

为了满足旋转不变性，要先将该关键点的局部区域旋转到其主方向上。旋转公式如下：

推导如下：

如果旋转的中心是原心，而图像又是关于旋转中心（原心）全对称的，所以，让图像旋转 $\theta$ 度后图像的新坐标等同于将两个坐标轴旋转 $\theta$ 度。示例如下：

我们要将关键点附近一定区域内的像素进行旋转，那么再来看这个区域的大小应该选取多少？

这个是答案：

解释如下：

描述子中要涵盖关键点周围一定区域内像素的信息，先把这个区域分为d*d个子区域，建议d取4，也就是4*4的子区域。每个子区域的大小选择经验值3* $\delta _{L}$ （这里的 $\delta _{L}$ 和第三步中的高斯模糊的系数1.5* $\delta _{L}$ 中的 $\delta _{L}$ 是一样的，都是 $\delta_{0} *2 ^{\tfrac{L}{S}}$ ，s是要从多少张图片中提取特征的值，L是关键点所在组的层数)。那么也就是这个区域的像素宽度为 $3\delta _{L}*d$ 。

像素就是一个小正方形块，这一个小块的值全都是同一个数。由于正方形中心计算更复杂，我们将正方形左下点看作像素位置。关键点并不一定是在像素位置，因为对关键点进行了精确定位。因此，其他像素位置的采样也不一定是在像素上，所以需要插值。进行双线性插值，所以，旋转的区域扩大，像素宽度为 $3\delta _{L}*(d+1)$ 。

又由于图像需要旋转，以旋转中可能出现的最大边长为需要进行旋转的区域范围，像素宽度为 $3\delta _{L}*(d+1)*\sqrt{2}$ 。所以，半径就变成了 $\frac{3\delta _{L}*(d+1)*\sqrt{2}}{2}$ 。

确定了区域之后，这个区域内的采样点都要进行旋转，也就是：

根据旋转的要求，我们需要将更大范围的像素进行旋转，但是我们生成描述子向量的采样范围其实就是d*d范围，旋转取更大范围只是为了让采样范围都包含其中。

确定采样点

旋转的范围比采样的范围大，所以要进一步确定采样点所处的范围。

这里的公式是：

先不看-0.5，运算公式将坐标都计算在[0,0]-[d,d]范围内。前面部分差不多是先缩放，再平移的。我觉得从实例化的角度能比较好的理解：

至于再减去-0.5，是因为我们将像素左下角位置等同于该点的像素位置，所以经常会将当前点偏移0.5个像素。如下图：

如果计算后处于[0,0]-[4,4]之间的，认为是采样点。

128维向量？

将区域内每个采样点像素的梯度大小和方向进行采样。采样后，根据主方向计算旋转后新的梯度方向：grad_ori -= ori（ori为主方向），对于这个式子的解释如下：

由旋转后的梯度方向可以将对应方向划分到8个bin中，公式为：obin = grad_ori * bins_per_rad （ grad_ori为旋转后的梯度方向，bins_per_rad 为8/(2*CV_PI)

此外，还要进行高斯滤波，离关键点近的权重大，远的权重小。计算权重公式为：

$2*(0.5d)^{2}$ 等于 $0.5*d^{2}$

然后，对每个子区域内的梯度方向进行统计（也就像投票那样，把高斯模糊后的梯度幅值加和到对应的方向中），分为8个方向。这样以来，Sift描述子就是4*4*8=128维的向量了。

三线性插值

因为采样点不一定是在像素位置，而是更精确的位置。在构建Sift描述子时，要将该采样点的值按照一定规则的比重“投射”到像素位置上。

双线性插值：

再举一个双线性插值的例子：现有一个点（x整+x小，y整+y小）x整为1，x小为0.25，y整为1,y小为0.75，即（1+0.25，1+0.75）。被（1，1）（1，2）（2，1）（2，2）包围。（2，1）这个点对应贡献值为f(2,1)*0.25*(1-0.25)。

三线性插值也是同理。

上例都是用周围点求其中的一个精确点，而Sift相当于“反其道而行之”，是将一个点的值按比重分配到周围点上。

上图左表示了向量维数是4*4*8。【形象化表示】上图右，实际采样点的位置是绿色的（比较的精确），要将该点的值插值到它周围8个点中。

看懂这个代码就明白了

for( r = 0; r <= 1; r++ ) //r=0,r=1
{
rb = r0 + r;//rb取整，值范围：x整+0（rb:-1 ~ 3）,x整+1（rb:0 ~ 4）
if( rb >= 0 && rb < d )//范围之外无需插值，指的是-1和4
   {
   v_r = mag * ( ( r == 0 )? 1.0 - d_r : d_r );
   row = hist[rb];//指针
   for( c = 0; c <= 1; c++ ) //c=0,c=1
   {
   cb = c0 + c;//cb取整，值范围：y整+0（cb:-1 ~ 3）,y整+1（cb:0 ~ 4）
   if( cb >= 0 && cb < d )//范围之外无需插值，指的是-1和4
       {
       v_c = v_r * ( ( c == 0 )? 1.0 - d_c : d_c );
       h = row[cb];//指针 hist[rb][cb]
       for( o = 0; o <= 1; o++ )//o=0,o=1
       {
       ob = ( o0 + o ) % n;//ob取整，值范围：0-7
       v_o = v_c * ( ( o == 0 )? 1.0 - d_o : d_o );
       h[ob] += v_o;//hist[rb][cb][ob]+=v_o
       }
       }
   }
   }
}

光照不变性

Lowe,D 的SIFT代码

调用大框架

//第四步：特征描述符的构建的大框架 
/*
  Computes feature descriptors for features in an array.  Based on Section 6
  of Lowe's paper.

  @param features array of features（特征点序列） 
  @param gauss_pyr Gaussian scale space pyramid（高斯金字塔） 
  @param d width of 2D array of orientation histograms（d为4） 
  @param n number of bins per orientation histogram（n为8） 
*/
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;//k为特征点的总个数 

  for( i = 0; i < k; i++ )
    {
      feat = CV_GET_SEQ_ELEM( struct feature, features, i );//将特征点序列中的第i个特征赋值给feat 
      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 );//调用函数descr_hist,对一个指定的特征点构造描述符 
      hist_to_descr( hist, d, n, feat );//调用函数hist_to_descr，光照不变性 
      release_descr_hist( &hist, d );//调用函数release_descr_hist，释放 
    }
}

旋转到主方向，调用三线性插值计算4*4*8维描述符

//对指定的特征点构造描述符 
/*
  Computes the 2D array of orientation histograms that form the feature
  descriptor.  Based on Section 6.1 of Lowe's paper.

  @param img image used in descriptor computation（关键点所在高斯金字塔对应的图像） 
  @param r row coord of center of orientation histogram array（关键点x坐标） 
  @param c column coord of center of orientation histogram array（关键点y坐标）
  @param ori canonical orientation of feature whose descr is being computed（关键点的方向） 
  @param scl scale relative to img of feature whose descr is being computed（ddata->scl_octv，高斯金字塔层数对应的高斯系数） 
  @param d width of 2d array of orientation histograms（d为4） 
  @param n bins per orientation histogram（n为8）

  @return Returns a d x d array of n-bin orientation histograms.
*/
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分配空间（4*4*8=128位向量） 
  //4*4*8
  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;// 8/2*CV_PI 
  exp_denom = d * d * 0.5;//计算权重需要用到 
  hist_width = SIFT_DESCR_SCL_FCTR * scl;//子区域的像素宽度 
  //根据旋转的要求，我们需要将更大范围的像素进行旋转，但是我们生成描述子向量的采样范围其实就是d*d范围，旋转取更大范围知识为了让采样范围都包含其中。 
  radius = hist_width * sqrt(2) * ( d + 1.0 ) * 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域（-0.5~3.5） 
	cbin = c_rot + d / 2 - 0.5;//当前采样点合理的y域（-0.5~3.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 )//控制调整后的梯度方向要在0-2*CV_PI这个范围内 
		      grad_ori += PI2;
	      while( grad_ori >= PI2 )
		      grad_ori -= PI2;
	      
	      obin = grad_ori * bins_per_rad;//grad_ori*(8/2*CV_PI)，obin在0~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 );//调用函数interp_hist_entry，三线性插值 
	    }
      }

  return hist;
}

三线性插值

//三线性插值，将一个采样点的值（赋值*权重）进行贡献分配 
/*
  Interpolates an entry into the array of orientation histograms that form
  the feature descriptor.

  @param hist 2D array of orientation histograms（hist:4*4*8的Sift描述子） 
  @param rbin sub-bin row coordinate of entry（rbin:0-4之间） 
  @param cbin sub-bin column coordinate of entry（cbin:0-4之间） 
  @param obin sub-bin orientation coordinate of entry（obin:0-8之间） 
  @param mag size of entry（赋值*权重） 
  @param d width of 2D array of orientation histograms（d为4） 
  @param n number of bins per orientation histogram（n为8） 
*/
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.
  */
  //三线性插值 hist[4][4][4]
  //rbin:-0.5 - 3.5  r0:-1~3
  //cbin:-0.5 - 3.5  c0:-1~3
  //obin:0-8  o0:0-7
  for( r = 0; r <= 1; r++ ) //r=0,r=1
    {
      rb = r0 + r;//rb取整，值范围：x整+0（rb:-1 ~ 3）,x整+1（rb:0 ~ 4） 
      if( rb >= 0  &&  rb < d )//范围之外无需插值，指的是-1和4 
	{
	  v_r = mag * ( ( r == 0 )? 1.0 - d_r : d_r );
	  row = hist[rb];//指针 
	  for( c = 0; c <= 1; c++ ) //c=0,c=1
	    {
	      cb = c0 + c;//cb取整，值范围：y整+0（cb:-1 ~ 3）,y整+1（cb:0 ~ 4）
	      if( cb >= 0  &&  cb < d )//范围之外无需插值，指的是-1和4 
		{
		  v_c = v_r * ( ( c == 0 )? 1.0 - d_c : d_c );
		  h = row[cb];//指针 hist[rb][cb]
		  for( o = 0; o <= 1; o++ )//o=0,o=1
		    {
		      ob = ( o0 + o ) % n;//ob取整，值范围：0-7 
		      v_o = v_c * ( ( o == 0 )? 1.0 - d_o : d_o );
		      h[ob] += v_o;//hist[rb][cb][ob]+=v_o
		    }
		}
	    }
	}
    }
}

两次归一化处理，满足光照不变性（不懂的是，为什么转为uchar时，乘以512.0）

//光照不变性 
/*
  Converts the 2D array of orientation histograms into a feature's descriptor
  vector.
  
  @param hist 2D array of orientation histograms（Sift描述子） 
  @param d width of hist（d为4） 
  @param n bins per histogram（n为8） 
  @param feat feature into which to store descriptor（特征点） 
*/
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]; //4*4*8 

  feat->d = k;//128 
  normalize_descr( feat );//调用函数normalize_descr，进行单位化操作 
  for( i = 0; i < k; i++ )
    if( feat->descr[i] > SIFT_DESCR_MAG_THR )//SIFT_DESCR_MAG_THR为0.2（阈值） 
      feat->descr[i] = SIFT_DESCR_MAG_THR;//大于0.2的用0.2替代 
  normalize_descr( feat );//调用函数normalize_descr，再进行单位化操作 

  /* convert floating-point descriptor to integer valued descriptor */
  for( i = 0; i < k; i++ )
    {
      int_val = SIFT_INT_DESCR_FCTR * feat->descr[i];//SIFT_INT_DESCR_FCTR为512.0，作为浮点型转换为整型时所用到的系数 
      feat->descr[i] = MIN( 255, int_val );//最后要将数据保存再uchar类型的存储结构中，因此乘以这个系数 
    }
}

归一化操作

//Sift特征子的单位化 
/*
  Normalizes a feature's descriptor vector to unitl length

  @param feat feature（特征值序列） 
*/
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;//代公式 
}

释放描述子空间

//释放描述符空间 
/*
  De-allocates memory held by a descriptor histogram

  @param hist pointer to a 2D array of orientation histograms
  @param d width of hist
*/
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;//最后，将指针置为空 
}