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

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