SIFT的原理已經有很多大牛的部落格上做了解析，本文重點將以Rob Hess等人用C實現的程式碼做解析，結合程式碼SIFT原理會更容易理解。一些難理解點的用了☆標註。

      歡迎大家批評指正！

轉載請註明出處：http://blog.csdn.net/luoshixian099/article/details/47377611

     SIFT（Scale-invariant feature transform）即尺度不變特徵轉換，提取的區域性特徵點具有尺度不變性，且對於旋轉，亮度，噪聲等有很高的穩定性。

下圖中，涉及到影象的旋轉，仿射，光照等變化，SIFT演算法依然有很好的匹配效果。

SIFT特徵點提取

 1.影象預處理

 2.構建高斯金字塔（不同尺度下的影象）

 3.生成DOG尺度空間

 4.關鍵點搜尋與定位

 5.計算特徵點所在的尺度

 6.為特徵點分配方向角

 7.構建特徵描述子

/**
   Finds SIFT features in an image using user-specified parameter values.  All
   detected features are stored in the array pointed to by \a feat.
*/
int _sift_features( IplImage* img, struct feature** feat, int intvls,
		    double sigma, double contr_thr, int curv_thr,
		    int img_dbl, int descr_width, int descr_hist_bins )
{
  IplImage* init_img;
  IplImage*** gauss_pyr, *** dog_pyr;
  CvMemStorage* storage;
  CvSeq* features;
  int octvs, i, n = 0;
  
  /* check arguments */
  if( ! img )
    fatal_error( "NULL pointer error, %s, line %d",  __FILE__, __LINE__ );
  if( ! feat )
    fatal_error( "NULL pointer error, %s, line %d",  __FILE__, __LINE__ );

  /* build scale space pyramid; smallest dimension of top level is ~4 pixels */
  init_img = create_init_img( img, img_dbl, sigma );                            //對進行圖片預處理       
  octvs = log( MIN( init_img->width, init_img->height ) ) / log(2) - 2;  //計算高斯金字塔的組數(octave),同時保證頂層至少有4個畫素點
  gauss_pyr = build_gauss_pyr( init_img, octvs, intvls, sigma );  //建立高斯金字塔
  dog_pyr = build_dog_pyr( gauss_pyr, octvs, intvls );   //DOG尺度空間
  
  storage = cvCreateMemStorage( 0 );
  features = scale_space_extrema( dog_pyr, octvs, intvls, contr_thr,   //極值點檢測，並去除不穩定特徵點
				  curv_thr, storage );
  calc_feature_scales( features, sigma, intvls );                      //計算特徵點所在的尺度
  if( img_dbl )
    adjust_for_img_dbl( features );                                       //如果影象初始被擴大了2倍，所有座標與尺度要除以2
  calc_feature_oris( features, gauss_pyr );                               //計算特徵點所在尺度內的方向角
  compute_descriptors( features, gauss_pyr, descr_width, descr_hist_bins );//計算特徵描述子 128維向量

  /* sort features by decreasing scale and move from CvSeq to array */
  cvSeqSort( features, (CvCmpFunc)feature_cmp, NULL );   //對特徵點按尺度排序
  n = features->total;
  *feat = calloc( n, sizeof(struct feature) );
  *feat = cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ );       
  for( i = 0; i < n; i++ )
    {
      free( (*feat)[i].feature_data );
      (*feat)[i].feature_data = NULL;
    }
  
  cvReleaseMemStorage( &storage );
  cvReleaseImage( &init_img );
  release_pyr( &gauss_pyr, octvs, intvls + 3 );
  release_pyr( &dog_pyr, octvs, intvls + 2 );
  return n;
}
 

  —————————————————————————————————————————————————————

 1.影象預處理

/************************ Functions prototyped here **************************/

/*
  Converts an image to 8-bit grayscale and Gaussian-smooths it.  The image is
  optionally doubled in size prior to smoothing.

  @param img input image
  @param img_dbl if true, image is doubled in size prior to smoothing
  @param sigma total std of Gaussian smoothing
*/
static IplImage* create_init_img( IplImage* img, int img_dbl, double sigma )
{
  IplImage* gray, * dbl;
  double sig_diff;

  gray = convert_to_gray32( img );   //轉換為32位灰度圖
  if( img_dbl )                                  // 影象被放大二倍
    {
      sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4 );   //  sigma = 1.6 , SIFT_INIT_SIGMA = 0.5  lowe認為影象在尺度0.5下最清晰
      dbl = cvCreateImage( cvSize( img->width*2, img->height*2 ),
			   IPL_DEPTH_32F, 1 );
      cvResize( gray, dbl, CV_INTER_CUBIC );  //雙三次插值方法 放大影象
      cvSmooth( dbl, dbl, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff );     //高斯平滑
      cvReleaseImage( &gray );
      return dbl;
    }
  else
    {
      sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA );
      cvSmooth( gray, gray, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff ); // 高斯平滑
      return gray;
    }
}
 

2.構建高斯金字塔

其中二維高斯核的計算為             。

                        

                                         

            

Eg. 在第一層上為了得到kσ的高斯模糊影象，可以在原圖上直接採用kσ平滑，也可以在上一層影象上（已被σ平滑）的影象上採用平滑因子為平滑影象，效果是一樣的。

/*
  Builds Gaussian scale space pyramid from an image
  @param base base image of the pyramid
  @param octvs number of octaves of scale space
  @param intvls number of intervals per octave
  @param sigma amount of Gaussian smoothing per octave

  @return Returns a Gaussian scale space pyramid as an octvs x (intvls + 3)
    array 
	
給定組數(octave)和層數(intvls)，以及初始平滑係數sigma,構建高斯金字塔
返回的每組中層數為intvls+3
*/
static IplImage*** build_gauss_pyr( IplImage* base, int octvs,
			     int intvls, double sigma )
{
  IplImage*** gauss_pyr;
  const int _intvls = intvls;             // lowe 採用了每組層數(intvls)為 3
 // double  sig_total, sig_prev;
	 double  k;
  int i, o;
  double *sig = (double *)malloc(sizeof(int)*(_intvls+3));  //儲存每組的高斯平滑因子，每組對應的平滑因子都相同

  gauss_pyr = calloc( octvs, sizeof( IplImage** ) );              
  for( i = 0; i < octvs; i++ )
    gauss_pyr[i] = calloc( intvls + 3, sizeof( IplImage *) );

  /*
    precompute Gaussian sigmas using the following formula:

    \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2

    sig[i] is the incremental sigma value needed to compute 
    the actual sigma of level i. Keeping track of incremental
    sigmas vs. total sigmas keeps the gaussian kernel small.
  */
  k = pow( 2.0, 1.0 / intvls );                 // k = 2^(1/S)
  sig[0] = sigma;
  sig[1] = sigma * sqrt( k*k- 1 );
  for (i = 2; i < intvls + 3; i++)
      sig[i] = sig[i-1] * k;                       //每組對應的平滑因子為 σ ,  sqrt(k^2 -1)* σ, sqrt(k^2 -1)* kσ , ...

  for( o = 0; o < octvs; o++ )
    for( i = 0; i < intvls + 3; i++ )
      {
	if( o == 0  &&  i == 0 )
	  gauss_pyr[o][i] = cvCloneImage(base);                       //第一組，第一層為原圖

	/* base of new octvave is halved image from end of previous octave */
	else if( i == 0 )
	  gauss_pyr[o][i] = downsample( gauss_pyr[o-1][intvls] );  //第一層影象由上一層倒數第三張隔點取樣得到
	  
	/* blur the current octave's last image to create the next one */
	else
	  {
	    gauss_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i-1]),
					     IPL_DEPTH_32F, 1 );
	    cvSmooth( gauss_pyr[o][i-1], gauss_pyr[o][i],
		      CV_GAUSSIAN, 0, 0, sig[i], sig[i] );                       //高斯平滑
	  }
      }
     
   
  return gauss_pyr;
}

3.生成DOG尺度空間

Lindeberg發現高斯差分函式（Difference of Gaussian ，簡稱DOG運算元）與尺度歸一化的高斯拉普拉斯函式非常近似，且

static IplImage*** build_dog_pyr( IplImage*** gauss_pyr, int octvs, int intvls )
{
  IplImage*** dog_pyr;
  int i, o;

  dog_pyr = calloc( octvs, sizeof( IplImage** ) );
  for( i = 0; i < octvs; i++ )
    dog_pyr[i] = calloc( intvls + 2, sizeof(IplImage*) );

  for( o = 0; o < octvs; o++ )
    for( i = 0; i < intvls + 2; i++ )
      {
	dog_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i]),
				       IPL_DEPTH_32F, 1 );
	cvSub( gauss_pyr[o][i+1], gauss_pyr[o][i], dog_pyr[o][i], NULL );   //相鄰兩層影象相減，結果放在dog_pyr陣列內
      }

  return dog_pyr;
}

 4.關鍵點搜尋與定位

/*
  Detects features at extrema in DoG scale space.  Bad features are discarded
  based on contrast and ratio of principal curvatures.

  @return Returns an array of detected features whose scales, orientations,
    and descriptors are yet to be determined.
*/
static CvSeq* scale_space_extrema( IplImage*** dog_pyr, int octvs, int intvls,
				   double contr_thr, int curv_thr,
				   CvMemStorage* storage )
{
  CvSeq* features;
  double prelim_contr_thr = 0.5 * contr_thr / intvls; //極值比較前的閾值處理
  struct feature* feat;
  struct detection_data* ddata;
  int o, i, r, c;

  features = cvCreateSeq( 0, sizeof(CvSeq), sizeof(struct feature), storage );
  for( o = 0; o < octvs; o++ )                     //對DOG尺度空間上，遍歷從第二層影象開始到倒數第二層影象上，每個畫素點
    for( i = 1; i <= intvls; i++ )
      for(r = SIFT_IMG_BORDER; r < dog_pyr[o][0]->height-SIFT_IMG_BORDER; r++)        
	for(c = SIFT_IMG_BORDER; c < dog_pyr[o][0]->width-SIFT_IMG_BORDER; c++)
	  /* perform preliminary check on contrast */
	  if( ABS( pixval32f( dog_pyr[o][i], r, c ) ) > prelim_contr_thr )    // 排除畫素值小於閾值prelim_contr_thr的點，提高穩定性
	    if( is_extremum( dog_pyr, o, i, r, c ) )             //與周圍26個畫素值比較，是否極大值或者極小值點
	      {
		feat = interp_extremum(dog_pyr, o, i, r, c, intvls, contr_thr); //插值處理，找到準確的特徵點座標
		if( feat )
		  {
		    ddata = feat_detection_data( feat );
		    if( ! is_too_edge_like( dog_pyr[ddata->octv][ddata->intvl],   //根據Hessian矩陣 判斷是否為邊緣上的點
					    ddata->r, ddata->c, curv_thr ) )
		      {
			cvSeqPush( features, feat );          //是特徵點進入特徵點序列
		      }
		    else
		      free( ddata );
		    free( feat );
		  }
	      }
  
  return features;
}

尋找極值點

  

/*
  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;
}

準確定位特徵點

                                      ，其中

                                           

                                        

是一個三維向量，向量在任何一個方向上的偏移量大於0.5時，意味著已經偏離了原畫素點，這樣的特徵座標位置需要更新或者繼續插值計算。演算法實現過程中，為了保證插值能夠收斂，設定了最大插值次數（lowe 設定了5次）。同時當時（本文閾值採用了0.04/S） ，特徵點才被保留，因為響應值過小的點，容易受噪聲的干擾而不穩定。

/*
  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;
}

刪除邊緣效應

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;
}

 5.計算特徵點對應的尺度

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 儲存特徵點在組內的尺度，用來下面計算方向角
    }
}

 6.為特徵點分配方向角

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

建立特徵點鄰域內的直方圖

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;  //返回累加完成的直方圖
}

平滑直方圖

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;
    }
}

複製特徵點

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 );
	}
    }
}

 7.構建特徵描述子

  

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;
}

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++ )
    if( feat->descr[i] > SIFT_DESCR_MAG_THR )   //設定門限，門限為0.2
      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 );
    }
}

------------------------------------------------------------------------------------

在此非常感謝CSDN上幾位影象上的大牛，我也是通過他們的文章去學習研究的，本文也是參考了他們的文章才寫成！

推薦看大牛們的文章，原理寫的很好！

http://blog.csdn.net/abcjennifer/article/details/7639681

http://blog.csdn.net/zddblog/article/details/7521424

http://blog.csdn.net/chen825919148/article/details/7685952

http://blog.csdn.net/xiaowei_cqu/article/details/8069548