Ransac算法代码分析和改进

SIFT算法中所用的RANSAC算法源码分析如下，算法的改进（改进是为适合我的应用而修改）包括5个方面：

（1）去除图像中某些噪声的影响，如摄像机的编号等等，在特征点比较少的情况下对算法影响较大；

（2）去除代码中对最小内点数的限制；

（3）增加最小循环次数的限制（保证M矩阵的精度）；

（4）增加最大循环次数的限制（保证循环可终结）；

（5）增加随机循环过程中的M矩阵的正确性的判别（在M矩阵满足某先验条件时，才执行后续操作）；

CvMat* ransac_xform( struct feature* features, int n, int mtype,
     ransac_xform_fn xform_fn, int m, double p_badxform,
     ransac_err_fn err_fn, double err_tol,
struct feature*** inliers, int* n_in )
{
 struct feature** matched, ** sample, ** consensus, ** consensus_max = NULL;
 struct ransac_data* rdata;
 CvPoint2D64f* pts, * mpts;
 CvMat* M = NULL;
 gsl_rng* rng;
 double p, in_frac = RANSAC_INLIER_FRAC_EST;
 int i, nm, in, in_min, in_max = 0, k = 0;
 const int MaxRunCount = 20000;
 const int MinRunCount = 2000;
 int runCount = 0;
 int v;

 //获得所有的匹配点对的个数
 nm = get_matched_features( features, n, mtype, &matched );
 if( nm < m )
 {
  fprintf( stderr, "Warning: not enough matches to compute xform, %s" \
   " line %d\n", __FILE__, __LINE__ );
  goto end;
 }

 
 rng = gsl_rng_alloc( gsl_rng_mt19937 );
 gsl_rng_set( rng, time(NULL) );

 //最小内点的个数
 in_min = calc_min_inliers( nm, m, RANSAC_PROB_BAD_SUPP, p_badxform );
 p = pow( 1.0 - pow( in_frac, m ), k );
 i = 0;
 
 while( p > p_badxform )
 {
  //随机采样m个样本
  sample = draw_ransac_sample( matched, nm, m, rng );
  //提取采样中的对应点对
  extract_corresp_pts( sample, m, mtype, &pts, &mpts );
  //求取采样点对的最小二乘解M矩阵
  M = xform_fn( pts, mpts, m );
  if( ! M )
   goto iteration_end;
  //寻找与随机模型一致的采样点对，err_fn函数定义了一致性的判别条件
  in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);

  //新增加的M矩阵的判别，由于我们应用中不存在缩放变换，故我们可以判断随机估计矩阵是否正确
  v = ValidateM(M);

  //记录满足最大一致性相对应的匹配点的数组和一致性个数
  if( in > in_max && v)
  {
   if( consensus_max )
    free( consensus_max );
   consensus_max = consensus;
   in_max = in;
   in_frac = (double)in_max / nm;

   //p = pow( 1.0 - pow( in_frac, m ), ++k );
  }
  else
   free( consensus );
  cvReleaseMat( &M );

iteration_end:
  release_mem( pts, mpts, sample );
  

  //增加最小迭代次数的限制
  if(runCount > MinRunCount)
  {
   //计算错误率
   p = pow( 1.0 - pow( in_frac, m ), ++k );
  }
  else
  {
   ++k;
  }

  //增加最大迭代次数的限制
  if(runCount > MaxRunCount)
  {
   break;
  }
  

  runCount++;
//    
 }

 printf("%d\n",runCount);
 
 
 if( 1)
 {
  extract_corresp_pts( consensus_max, in_max, mtype, &pts, &mpts );
  M = xform_fn( pts, mpts, in_max );
  in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);
  cvReleaseMat( &M );
  release_mem( pts, mpts, consensus_max );
  extract_corresp_pts( consensus, in, mtype, &pts, &mpts );
  M = xform_fn( pts, mpts, in );
  if( inliers )
  {
   *inliers = consensus;
   consensus = NULL;
  }
  if( n_in )
   *n_in = in;
  release_mem( pts, mpts, consensus );
 }
 else if( consensus_max )
 {
  if( inliers )
   *inliers = NULL;
  if( n_in )
   *n_in = 0;
  free( consensus_max );
 }

 gsl_rng_free( rng );
end:
 for( i = 0; i < nm; i++ )
 {
  rdata = feat_ransac_data( matched[i] );
  matched[i]->feature_data = rdata->orig_feat_data;
  free( rdata );
 }
 free( matched );
 return M;
}