转载来自：http://blog.csdn.net/laobai1015/article/details/51683076

本文将OpenCV中的RANSAC代码全部挑选出来，进行分析和讲解，以便大家更好的理解RANSAC算法。代码我都试过，可以直接运行。

在计算机视觉和图像处理等很多领域，都需要用到RANSAC算法。openCV中也有封装好的RANSAC算法，以便于人们使用。关于RANSAC算法的一些应用，可以看我的另一篇博客：

利用SIFT和RANSAC算法（openCV框架）实现物体的检测与定位，并求出变换矩阵（findFundamentalMat和findHomography的比较）

但是前几天师弟在使用openCV自带的RANSAC算法时，发现实验的运行时间并不会随着输入数据的增加而增加，感觉和理论上的不太相符。所以我就花了点时间，把openCV中关于RANSAC的源代码全部复制出来研究了一下。以便我们更加清晰的了解RANSAC算法的实际运行过程。

首先看两个类

[cpp] view plain copy    
  1. //模型估计的基类，提供了估计矩阵的各种虚函数
  2. //置信度设为0。99 循环次数设置为了2000
  3. class CvModelEstimator2
  4. {
  5. public:
  6. CvModelEstimator2(int _modelPoints, CvSize _modelSize, int _maxBasicSolutions);
  7. virtual ~CvModelEstimator2();
  8. virtual int runKernel( const CvMat* m1, const CvMat* m2, CvMat* model )=0;
  9. //virtual bool runLMeDS( const CvMat* m1, const CvMat* m2, CvMat* model,
  10. // CvMat* mask, double confidence=0.99, int maxIters=2000 );
  11. virtual bool runRANSAC( const CvMat* m1, const CvMat* m2, CvMat* model,
  12. CvMat* mask, double threshold,
  13. double confidence=0.99, int maxIters=2000 );
  14. virtual bool refine( const CvMat*, const CvMat*, CvMat*, int ) { return true; }
  15. //virtual void setSeed( int64 seed );
  16. protected:
  17. virtual void computeReprojError( const CvMat* m1, const CvMat* m2,
  18. const CvMat* model, CvMat* error ) = 0;
  19. virtual int findInliers( const CvMat* m1, const CvMat* m2,
  20. const CvMat* model, CvMat* error,
  21. CvMat* mask, double threshold );
  22. virtual bool getSubset( const CvMat* m1, const CvMat* m2,
  23. CvMat* ms1, CvMat* ms2, int maxAttempts=1000 );
  24. virtual bool checkSubset( const CvMat* ms1, int count );
  25. CvRNG rng;
  26. int modelPoints;
  27. CvSize modelSize;
  28. int maxBasicSolutions;
  29. bool checkPartialSubsets;
  30. };
  31. //单应矩阵估计的子类
  32. class CvHomographyEstimator : public CvModelEstimator2
  33. {
  34. public:
  35. CvHomographyEstimator( int modelPoints );
  36. virtual int runKernel( const CvMat* m1, const CvMat* m2, CvMat* model );
  37. virtual bool refine( const CvMat* m1, const CvMat* m2,
  38. CvMat* model, int maxIters );
  39. protected:
  40. virtual void computeReprojError( const CvMat* m1, const CvMat* m2,
  41. const CvMat* model, CvMat* error );
  42. };

上面的两个类中，CvModelEstimator2是一个基类，从名字就可以看出，这个类是用来估计模型的。可以看到里面提供了许多虚函数，这些函数有许多，比如runRANSAC是利用RANSAC方法计算单应矩阵，而runLMeDS是利用LMeDS方法计算单应矩阵，我们这里仅仅讲解RANSAC方法，所以其他不需要的内容我就直接注释掉了

CvHomographyEstimator继承自CvModelEstimator2，同样的，从名字也就可以看出，这个类使用来估计单应矩阵的。

接下来是两个类的构造函数和析构函数，这个没啥好说的了，基本都是默认的。

[html] view plain copy    
  1. <pre name="code" class="cpp">//构造函数
  2. CvModelEstimator2::CvModelEstimator2(int _modelPoints, CvSize _modelSize, int _maxBasicSolutions)
  3. {
  4. modelPoints = _modelPoints;
  5. modelSize = _modelSize;
  6. maxBasicSolutions = _maxBasicSolutions;
  7. checkPartialSubsets = true;
  8. rng = cvRNG(-1);
  9. }
  10. //析构函数
  11. CvModelEstimator2::~CvModelEstimator2()
  12. {
  13. }
  14. CvHomographyEstimator::CvHomographyEstimator(int _modelPoints)
  15. : CvModelEstimator2(_modelPoints, cvSize(3,3), 1)
  16. {
  17. assert( _modelPoints == 4 || _modelPoints == 5 );
  18. checkPartialSubsets = false;
  19. }

接下来到重点了。runRANSAC方法就是通过RANSAC来计算矩阵

[html] view plain copy    
  1. <pre name="code" class="cpp">bool CvModelEstimator2::runRANSAC( const CvMat* m1, const CvMat* m2, CvMat* model,
  2. CvMat* mask0, double reprojThreshold,
  3. double confidence, int maxIters )
  4. {
  5. bool result = false;
  6. cv::Ptr<CvMat> mask = cvCloneMat(mask0);   //标记矩阵，标记内点和外点
  7. cv::Ptr<CvMat> models, err, tmask;
  8. cv::Ptr<CvMat> ms1, ms2;
  9. int iter, niters = maxIters;   //这是迭代次数，默认最大的迭代次数为2000次
  10. int count = m1->rows*m1->cols, maxGoodCount = 0;
  11. CV_Assert( CV_ARE_SIZES_EQ(m1, m2) && CV_ARE_SIZES_EQ(m1, mask) );
  12. if( count < modelPoints )  //使用RANSAC算法时，modelPoints为4
  13. return false;
  14. models = cvCreateMat( modelSize.height*maxBasicSolutions, modelSize.width, CV_64FC1 );
  15. err = cvCreateMat( 1, count, CV_32FC1 );
  16. tmask = cvCreateMat( 1, count, CV_8UC1 );
  17. if( count > modelPoints )
  18. {
  19. ms1 = cvCreateMat( 1, modelPoints, m1->type );
  20. ms2 = cvCreateMat( 1, modelPoints, m2->type );
  21. }
  22. else
  23. {
  24. niters = 1;
  25. ms1 = cvCloneMat(m1);
  26. ms2 = cvCloneMat(m2);
  27. }
  28. for( iter = 0; iter < niters; iter++ )
  29. {
  30. int i, goodCount, nmodels;
  31. if( count > modelPoints )
  32. {
  33. bool found = getSubset( m1, m2, ms1, ms2, 300 );//调用函数，300是循环次数，这个函数
  34. if( !found )                                    //就是为了从序列中随机选取4组，以便
  35. {                                               //以便下一步计算单应矩阵
  36. if( iter == 0 )
  37. return false;
  38. break;
  39. }
  40. }
  41. printf("------");
  42. nmodels = runKernel( ms1, ms2, models );//这个函数是通过给定的4组序列计算出矩阵
  43. if( nmodels <= 0 )
  44. continue;
  45. for( i = 0; i < nmodels; i++ )
  46. {
  47. CvMat model_i;
  48. cvGetRows( models, &model_i, i*modelSize.height, (i+1)*modelSize.height );
  49. goodCount = findInliers( m1, m2, &model_i, err, tmask, reprojThreshold );
  50. //输出看看一共循环了多少次
  51. printf("%5d %5d %5d %5d\n",iter,niters,goodCount,maxGoodCount);
  52. if( goodCount > MAX(maxGoodCount, modelPoints-1) )
  53. {
  54. std::swap(tmask, mask);
  55. cvCopy( &model_i, model );
  56. maxGoodCount = goodCount;
  57. //循环的次数会发生变化，原来原因在这里
  58. niters = cvRANSACUpdateNumIters( confidence,
  59. (double)(count - goodCount)/count, modelPoints, niters );
  60. }
  61. }
  62. }
  63. //printf("RANSAC算法实际循环了%d次\n",niters);
  64. if( maxGoodCount > 0 )
  65. {
  66. if( mask != mask0 )
  67. cvCopy( mask, mask0 );
  68. result = true;
  69. }
  70. return result;
  71. }

在这个函数参数中，输入的m1和m2是两个对应的序列，这两组序列的每一对数据一一匹配，其中既有正确的匹配，也有错误的匹配，正确的可以称为内点，错误的称为外点，RANSAC方法就是从这些包含错误匹配的数据中，分离出正确的匹配，并且求得单应矩阵。model就是我们需要求解的单应矩阵，mask我们可以称为标记矩阵，他和m1，m2的长度一样，当一个m1和m2中的点为内点时，mask相应的标记为1，反之为0，说白了，通过mask我们最终可以知道序列中哪些是内点，哪些是外点。reprojThreshold为阈值，当某一个匹配与估计的假设小于阈值时，则被认为是一个内点，这个阈值，openCV默认给的是3，后期使用的时候自己也可以修改。confidence为置信度，其实也就是人为的规定了一个数值，这个数值可以大致表示RANSAC结果的准确性，这个具体有啥用后面咱们再说。这个值初始时被设置为0.995. maxIters为初始迭代次数，RANSAC算法核心就是不断的迭代，这个值就是迭代的次数，默认设为了2000

这个函数的前期，主要是设置了一些变量然后赋初值，然后转换相应的格式等等。最关键的部分，是那个for循环。我们把这个for循环单独拿出来分析一下。代码如下。

[cpp] view plain copy    
  1. for( iter = 0; iter < niters; iter++ )
  2. {
  3. int i, goodCount, nmodels;
  4. if( count > modelPoints )
  5. {
  6. bool found = getSubset( m1, m2, ms1, ms2, 300 );//调用函数，300是循环次数，这个函数
  7. if( !found )                                    //就是为了从序列中随机选取4组，以便
  8. {                                               //以便下一步计算单应矩阵
  9. if( iter == 0 )
  10. return false;
  11. break;
  12. }
  13. }
  14. nmodels = runKernel( ms1, ms2, models );//这个函数是通过给定的4组序列计算出矩阵
  15. if( nmodels <= 0 )
  16. continue;
  17. for( i = 0; i < nmodels; i++ )
  18. {
  19. CvMat model_i;
  20. cvGetRows( models, &model_i, i*modelSize.height, (i+1)*modelSize.height );
  21. goodCount = findInliers( m1, m2, &model_i, err, tmask, reprojThreshold );
  22. //输出看看一共循环了多少次
  23. printf("%5d %5d %5d %5d\n",iter,niters,goodCount,maxGoodCount);
  24. if( goodCount > MAX(maxGoodCount, modelPoints-1) )
  25. {
  26. std::swap(tmask, mask);
  27. cvCopy( &model_i, model );
  28. maxGoodCount = goodCount;
  29. //循环的次数会发生变化，原来原因在这里
  30. niters = cvRANSACUpdateNumIters( confidence,
  31. (double)(count - goodCount)/count, modelPoints, niters );
  32. }
  33. }
  34. }

niters最初的值为2000，这就是初始时的RANSAC算法的循环次数，getSubset（）函数是从一组对应的序列中随机的选出4组（因为要想计算出一个3X3的矩阵，至少需要4组对应的坐标），m1和m2是我们输入序列，ms1和ms2是随机选出的对应的4组匹配。

随机的选出4组匹配后，就应该根据这4个匹配计算相应的矩阵，所以函数runKernel（）就是根据4组匹配计算矩阵，参数里的models就是得到的矩阵。这个矩阵只是一个假设，为了验证这个假设，需要用其他的点去计算，看看其他的点是内点还是外点。

findInliers（）函数就是用来计算内点的。利用前面得到的矩阵，把所有的序列带入，计算得出哪些是内点，哪些是外点，函数的返回值为goodCount，就是此次计算的内点的个数。函数中还有一个值为maxGoodCount，每次循环的内点个数的最大值保存在这个值中，一个估计的矩阵如果有越多的内点，那么这个矩阵就越有可能是正确的。所以计算内点个数以后，紧接着判断一下goodCount和maxGoodCount的大小关系，如果goodCount>maxGoodCount，则把goodCount赋值给maxGoodCount。赋值之后的一行代码非常关键，我们单独拿出来说一下，代码如下：

[cpp] view plain copy    
  1. niters = cvRANSACUpdateNumIters( confidence,
  2. (double)(count - goodCount)/count, modelPoints, niters );

niters本来是迭代的次数，也就是循环的次数。但是通过这行代码我们发现，每次循环后，都会对niters这个值进行更新，也就是每次循环后都会改变循环的总次数。cvRANSACUpdateNumIters（）函数利用confidence（置信度）count（总匹配个数）goodCount（当前内点个数）niters（当前的总迭代次数）这几个参数，来动态的改变总迭代次数的大小。该函数的中心思想就是当内点占的比例较多时，那么很有可能已经找到了正确的估计，所以就适当的减少迭代次数来节省时间。这个迭代次数的减少是以指数形式减少的，所以节省的时间开销也是非常的可观。因此最初设计的2000的迭代次数，可能最终的迭代次数只有几十。同样的，如果你自己一开始把迭代次数设置成10000或者更大，进过几次迭代后，niters又会变得非常小了。所以初始时的niters设置的再大，其实对最终的运行时间也没什么影响。我用我自己的程序简答试了一下，无论初值设为2000,10000,20000，最终的迭代次数都变成了58！！！

所以，们现在应该清楚为什么输入数据增加，而算法运行时间不会增加了。openCV的RANSAC算法首先把迭代的次数设置为2000，然后再迭代的过程中，动态的改变总迭代次数，无论输入数据有多少，总的迭代次数不会增加，并且通过4个匹配计算出估计的矩阵这个时间是不变的，通过估计矩阵来计算内点，这方面的增加的时间开销基本上可以忽略。所以导致的最终结果就是，无论输入点有多少，运算时间基本不会有太大变化。

以上就是RANSAC算法的核心代码，其中用到的一些函数，下面一一给出。

1. 转换为齐次左边，看上去很长，但是完成的功能就是把一般的坐标转换成齐次坐标以方便以后的计算

[cpp] view plain copy    
  1. CV_IMPL void cvConvertPointsHomogeneous( const CvMat* src, CvMat* dst )
  2. {
  3. Ptr<CvMat> temp, denom;
  4. int i, s_count, s_dims, d_count, d_dims;
  5. CvMat _src, _dst, _ones;
  6. CvMat* ones = 0;
  7. if( !CV_IS_MAT(src) )
  8. CV_Error( !src ? CV_StsNullPtr : CV_StsBadArg,
  9. "The input parameter is not a valid matrix" );
  10. if( !CV_IS_MAT(dst) )
  11. CV_Error( !dst ? CV_StsNullPtr : CV_StsBadArg,
  12. "The output parameter is not a valid matrix" );
  13. if( src == dst || src->data.ptr == dst->data.ptr )
  14. {
  15. if( src != dst && (!CV_ARE_TYPES_EQ(src, dst) || !CV_ARE_SIZES_EQ(src,dst)) )
  16. CV_Error( CV_StsBadArg, "Invalid inplace operation" );
  17. return;
  18. }
  19. if( src->rows > src->cols )
  20. {
  21. if( !((src->cols > 1) ^ (CV_MAT_CN(src->type) > 1)) )
  22. CV_Error( CV_StsBadSize, "Either the number of channels or columns or rows must be =1" );
  23. s_dims = CV_MAT_CN(src->type)*src->cols;
  24. s_count = src->rows;
  25. }
  26. else
  27. {
  28. if( !((src->rows > 1) ^ (CV_MAT_CN(src->type) > 1)) )
  29. CV_Error( CV_StsBadSize, "Either the number of channels or columns or rows must be =1" );
  30. s_dims = CV_MAT_CN(src->type)*src->rows;
  31. s_count = src->cols;
  32. }
  33. if( src->rows == 1 || src->cols == 1 )
  34. src = cvReshape( src, &_src, 1, s_count );
  35. if( dst->rows > dst->cols )
  36. {
  37. if( !((dst->cols > 1) ^ (CV_MAT_CN(dst->type) > 1)) )
  38. CV_Error( CV_StsBadSize,
  39. "Either the number of channels or columns or rows in the input matrix must be =1" );
  40. d_dims = CV_MAT_CN(dst->type)*dst->cols;
  41. d_count = dst->rows;
  42. }
  43. else
  44. {
  45. if( !((dst->rows > 1) ^ (CV_MAT_CN(dst->type) > 1)) )
  46. CV_Error( CV_StsBadSize,
  47. "Either the number of channels or columns or rows in the output matrix must be =1" );
  48. d_dims = CV_MAT_CN(dst->type)*dst->rows;
  49. d_count = dst->cols;
  50. }
  51. if( dst->rows == 1 || dst->cols == 1 )
  52. dst = cvReshape( dst, &_dst, 1, d_count );
  53. if( s_count != d_count )
  54. CV_Error( CV_StsUnmatchedSizes, "Both matrices must have the same number of points" );
  55. if( CV_MAT_DEPTH(src->type) < CV_32F || CV_MAT_DEPTH(dst->type) < CV_32F )
  56. CV_Error( CV_StsUnsupportedFormat,
  57. "Both matrices must be floating-point (single or double precision)" );
  58. if( s_dims < 2 || s_dims > 4 || d_dims < 2 || d_dims > 4 )
  59. CV_Error( CV_StsOutOfRange,
  60. "Both input and output point dimensionality must be 2, 3 or 4" );
  61. if( s_dims < d_dims - 1 || s_dims > d_dims + 1 )
  62. CV_Error( CV_StsUnmatchedSizes,
  63. "The dimensionalities of input and output point sets differ too much" );
  64. if( s_dims == d_dims - 1 )
  65. {
  66. if( d_count == dst->rows )
  67. {
  68. ones = cvGetSubRect( dst, &_ones, cvRect( s_dims, 0, 1, d_count ));
  69. dst = cvGetSubRect( dst, &_dst, cvRect( 0, 0, s_dims, d_count ));
  70. }
  71. else
  72. {
  73. ones = cvGetSubRect( dst, &_ones, cvRect( 0, s_dims, d_count, 1 ));
  74. dst = cvGetSubRect( dst, &_dst, cvRect( 0, 0, d_count, s_dims ));
  75. }
  76. }
  77. if( s_dims <= d_dims )
  78. {
  79. if( src->rows == dst->rows && src->cols == dst->cols )
  80. {
  81. if( CV_ARE_TYPES_EQ( src, dst ) )
  82. cvCopy( src, dst );
  83. else
  84. cvConvert( src, dst );
  85. }
  86. else
  87. {
  88. if( !CV_ARE_TYPES_EQ( src, dst ))
  89. {
  90. temp = cvCreateMat( src->rows, src->cols, dst->type );
  91. cvConvert( src, temp );
  92. src = temp;
  93. }
  94. cvTranspose( src, dst );
  95. }
  96. if( ones )
  97. cvSet( ones, cvRealScalar(1.) );
  98. }
  99. else
  100. {
  101. int s_plane_stride, s_stride, d_plane_stride, d_stride, elem_size;
  102. if( !CV_ARE_TYPES_EQ( src, dst ))
  103. {
  104. temp = cvCreateMat( src->rows, src->cols, dst->type );
  105. cvConvert( src, temp );
  106. src = temp;
  107. }
  108. elem_size = CV_ELEM_SIZE(src->type);
  109. if( s_count == src->cols )
  110. s_plane_stride = src->step / elem_size, s_stride = 1;
  111. else
  112. s_stride = src->step / elem_size, s_plane_stride = 1;
  113. if( d_count == dst->cols )
  114. d_plane_stride = dst->step / elem_size, d_stride = 1;
  115. else
  116. d_stride = dst->step / elem_size, d_plane_stride = 1;
  117. denom = cvCreateMat( 1, d_count, dst->type );
  118. if( CV_MAT_DEPTH(dst->type) == CV_32F )
  119. {
  120. const float* xs = src->data.fl;
  121. const float* ys = xs + s_plane_stride;
  122. const float* zs = 0;
  123. const float* ws = xs + (s_dims - 1)*s_plane_stride;
  124. float* iw = denom->data.fl;
  125. float* xd = dst->data.fl;
  126. float* yd = xd + d_plane_stride;
  127. float* zd = 0;
  128. if( d_dims == 3 )
  129. {
  130. zs = ys + s_plane_stride;
  131. zd = yd + d_plane_stride;
  132. }
  133. for( i = 0; i < d_count; i++, ws += s_stride )
  134. {
  135. float t = *ws;
  136. iw[i] = fabs((double)t) > FLT_EPSILON ? t : 1.f;
  137. }
  138. cvDiv( 0, denom, denom );
  139. if( d_dims == 3 )
  140. for( i = 0; i < d_count; i++ )
  141. {
  142. float w = iw[i];
  143. float x = *xs * w, y = *ys * w, z = *zs * w;
  144. xs += s_stride; ys += s_stride; zs += s_stride;
  145. *xd = x; *yd = y; *zd = z;
  146. xd += d_stride; yd += d_stride; zd += d_stride;
  147. }
  148. else
  149. for( i = 0; i < d_count; i++ )
  150. {
  151. float w = iw[i];
  152. float x = *xs * w, y = *ys * w;
  153. xs += s_stride; ys += s_stride;
  154. *xd = x; *yd = y;
  155. xd += d_stride; yd += d_stride;
  156. }
  157. }
  158. else
  159. {
  160. const double* xs = src->data.db;
  161. const double* ys = xs + s_plane_stride;
  162. const double* zs = 0;
  163. const double* ws = xs + (s_dims - 1)*s_plane_stride;
  164. double* iw = denom->data.db;
  165. double* xd = dst->data.db;
  166. double* yd = xd + d_plane_stride;
  167. double* zd = 0;
  168. if( d_dims == 3 )
  169. {
  170. zs = ys + s_plane_stride;
  171. zd = yd + d_plane_stride;
  172. }
  173. for( i = 0; i < d_count; i++, ws += s_stride )
  174. {
  175. double t = *ws;
  176. iw[i] = fabs(t) > DBL_EPSILON ? t : 1.;
  177. }
  178. cvDiv( 0, denom, denom );
  179. if( d_dims == 3 )
  180. for( i = 0; i < d_count; i++ )
  181. {
  182. double w = iw[i];
  183. double x = *xs * w, y = *ys * w, z = *zs * w;
  184. xs += s_stride; ys += s_stride; zs += s_stride;
  185. *xd = x; *yd = y; *zd = z;
  186. xd += d_stride; yd += d_stride; zd += d_stride;
  187. }
  188. else
  189. for( i = 0; i < d_count; i++ )
  190. {
  191. double w = iw[i];
  192. double x = *xs * w, y = *ys * w;
  193. xs += s_stride; ys += s_stride;
  194. *xd = x; *yd = y;
  195. xd += d_stride; yd += d_stride;
  196. }
  197. }
  198. }
  199. }

2. 对迭代值进行更新的函数。这个函数就是对总的迭代次数进行更新，从中可以看到，迭代值以指数形式减少。最初的为2000的迭代次数，有的时候可能经过不断的更新，最终结果成了几十了。

[cpp] view plain copy    
  1. CV_IMPL int
  2. cvRANSACUpdateNumIters( double p, double ep,
  3. int model_points, int max_iters )
  4. {
  5. if( model_points <= 0 )
  6. CV_Error( CV_StsOutOfRange, "the number of model points should be positive" );
  7. p = MAX(p, 0.);
  8. p = MIN(p, 1.);
  9. ep = MAX(ep, 0.);
  10. ep = MIN(ep, 1.);
  11. // avoid inf's & nan's
  12. double num = MAX(1. - p, DBL_MIN);
  13. double denom = 1. - pow(1. - ep,model_points);
  14. if( denom < DBL_MIN )
  15. return 0;
  16. num = log(num);
  17. denom = log(denom);
  18. return denom >= 0 || -num >= max_iters*(-denom) ?
  19. max_iters : cvRound(num/denom);
  20. }

3. 通过4个匹配，计算单应矩阵，就是给你了4个匹配，你把和这四个匹配相符的矩阵计算出来

[cpp] view plain copy    
  1. //通过四个匹配，计算符合要求的单应矩阵
  2. int CvHomographyEstimator::runKernel( const CvMat* m1, const CvMat* m2, CvMat* H )
  3. {
  4. int i, count = m1->rows*m1->cols;
  5. const CvPoint2D64f* M = (const CvPoint2D64f*)m1->data.ptr;
  6. const CvPoint2D64f* m = (const CvPoint2D64f*)m2->data.ptr;
  7. double LtL[9][9], W[9][1], V[9][9];
  8. CvMat _LtL = cvMat( 9, 9, CV_64F, LtL );
  9. CvMat matW = cvMat( 9, 1, CV_64F, W );
  10. CvMat matV = cvMat( 9, 9, CV_64F, V );
  11. CvMat _H0 = cvMat( 3, 3, CV_64F, V[8] );
  12. CvMat _Htemp = cvMat( 3, 3, CV_64F, V[7] );
  13. CvPoint2D64f cM={0,0}, cm={0,0}, sM={0,0}, sm={0,0};
  14. for( i = 0; i < count; i++ )
  15. {
  16. cm.x += m[i].x; cm.y += m[i].y;
  17. cM.x += M[i].x; cM.y += M[i].y;
  18. }
  19. cm.x /= count; cm.y /= count;
  20. cM.x /= count; cM.y /= count;
  21. for( i = 0; i < count; i++ )
  22. {
  23. sm.x += fabs(m[i].x - cm.x);
  24. sm.y += fabs(m[i].y - cm.y);
  25. sM.x += fabs(M[i].x - cM.x);
  26. sM.y += fabs(M[i].y - cM.y);
  27. }
  28. if( fabs(sm.x) < DBL_EPSILON || fabs(sm.y) < DBL_EPSILON ||
  29. fabs(sM.x) < DBL_EPSILON || fabs(sM.y) < DBL_EPSILON )
  30. return 0;
  31. sm.x = count/sm.x; sm.y = count/sm.y;
  32. sM.x = count/sM.x; sM.y = count/sM.y;
  33. double invHnorm[9] = { 1./sm.x, 0, cm.x, 0, 1./sm.y, cm.y, 0, 0, 1 };
  34. double Hnorm2[9] = { sM.x, 0, -cM.x*sM.x, 0, sM.y, -cM.y*sM.y, 0, 0, 1 };
  35. CvMat _invHnorm = cvMat( 3, 3, CV_64FC1, invHnorm );
  36. CvMat _Hnorm2 = cvMat( 3, 3, CV_64FC1, Hnorm2 );
  37. cvZero( &_LtL );
  38. for( i = 0; i < count; i++ )
  39. {
  40. double x = (m[i].x - cm.x)*sm.x, y = (m[i].y - cm.y)*sm.y;
  41. double X = (M[i].x - cM.x)*sM.x, Y = (M[i].y - cM.y)*sM.y;
  42. double Lx[] = { X, Y, 1, 0, 0, 0, -x*X, -x*Y, -x };
  43. double Ly[] = { 0, 0, 0, X, Y, 1, -y*X, -y*Y, -y };
  44. int j, k;
  45. for( j = 0; j < 9; j++ )
  46. for( k = j; k < 9; k++ )
  47. LtL[j][k] += Lx[j]*Lx[k] + Ly[j]*Ly[k];
  48. }
  49. cvCompleteSymm( &_LtL );
  50. //cvSVD( &_LtL, &matW, 0, &matV, CV_SVD_MODIFY_A + CV_SVD_V_T );
  51. cvEigenVV( &_LtL, &matV, &matW );
  52. cvMatMul( &_invHnorm, &_H0, &_Htemp );
  53. cvMatMul( &_Htemp, &_Hnorm2, &_H0 );
  54. cvConvertScale( &_H0, H, 1./_H0.data.db[8] );
  55. return 1;
  56. }

4. 给定输入序列后，从中随机的选出4对匹配

[cpp] view plain copy    
  1. bool CvModelEstimator2::getSubset( const CvMat* m1, const CvMat* m2,
  2. CvMat* ms1, CvMat* ms2, int maxAttempts )  //maxAttempts被设为300
  3. {
  4. cv::AutoBuffer<int> _idx(modelPoints);
  5. int* idx = _idx;
  6. int i = 0, j, k, idx_i, iters = 0;
  7. int type = CV_MAT_TYPE(m1->type), elemSize = CV_ELEM_SIZE(type);
  8. const int *m1ptr = m1->data.i, *m2ptr = m2->data.i;
  9. int *ms1ptr = ms1->data.i, *ms2ptr = ms2->data.i;
  10. int count = m1->cols*m1->rows;
  11. assert( CV_IS_MAT_CONT(m1->type & m2->type) && (elemSize % sizeof(int) == 0) );
  12. elemSize /= sizeof(int);
  13. for(; iters < maxAttempts; iters++)
  14. {
  15. for( i = 0; i < modelPoints && iters < maxAttempts; )
  16. {
  17. idx[i] = idx_i = cvRandInt(&rng) % count;    //产生count以内的随机数，count是序列长度
  18. for( j = 0; j < i; j++ )                    //保证产生的随机数没有重复的
  19. if( idx_i == idx[j] )
  20. break;
  21. if( j < i )
  22. continue;
  23. for( k = 0; k < elemSize; k++ )
  24. {
  25. ms1ptr[i*elemSize + k] = m1ptr[idx_i*elemSize + k];   //把随机产生的数给了ms1和ms2
  26. ms2ptr[i*elemSize + k] = m2ptr[idx_i*elemSize + k];
  27. }
  28. if( checkPartialSubsets && (!checkSubset( ms1, i+1 ) || !checkSubset( ms2, i+1 )))   //调用函数checkSubset
  29. {
  30. iters++;
  31. continue;
  32. }
  33. i++;
  34. }
  35. if( !checkPartialSubsets && i == modelPoints &&
  36. (!checkSubset( ms1, i ) || !checkSubset( ms2, i )))
  37. continue;
  38. break;
  39. }
  40. return i == modelPoints && iters < maxAttempts;
  41. }

5. 对生成的4组匹配进行检验，观察其是否合乎要求。

[cpp] view plain copy    
  1. bool CvModelEstimator2::checkSubset( const CvMat* m, int count )
  2. {
  3. int j, k, i, i0, i1;
  4. CvPoint2D64f* ptr = (CvPoint2D64f*)m->data.ptr;
  5. assert( CV_MAT_TYPE(m->type) == CV_64FC2 );
  6. if( checkPartialSubsets )
  7. i0 = i1 = count - 1;
  8. else
  9. i0 = 0, i1 = count - 1;
  10. for( i = i0; i <= i1; i++ )
  11. {
  12. // check that the i-th selected point does not belong
  13. // to a line connecting some previously selected points
  14. for( j = 0; j < i; j++ )
  15. {
  16. double dx1 = ptr[j].x - ptr[i].x;
  17. double dy1 = ptr[j].y - ptr[i].y;
  18. for( k = 0; k < j; k++ )
  19. {
  20. double dx2 = ptr[k].x - ptr[i].x;
  21. double dy2 = ptr[k].y - ptr[i].y;
  22. if( fabs(dx2*dy1 - dy2*dx1) <= FLT_EPSILON*(fabs(dx1) + fabs(dy1) + fabs(dx2) + fabs(dy2)))
  23. break;
  24. }
  25. if( k < j )
  26. break;
  27. }
  28. if( j < i )
  29. break;
  30. }
  31. return i >= i1;
  32. }

6. 计算内点的个数并且标记序列中哪些点是内点。

[cpp] view plain copy    
  1. int CvModelEstimator2::findInliers( const CvMat* m1, const CvMat* m2,
  2. const CvMat* model, CvMat* _err,
  3. CvMat* _mask, double threshold )
  4. {
  5. int i, count = _err->rows*_err->cols, goodCount = 0;
  6. const float* err = _err->data.fl;
  7. uchar* mask = _mask->data.ptr;
  8. computeReprojError( m1, m2, model, _err );  //_err里面是计算后的矩阵的大小，用于与阈值比较
  9. threshold *= threshold;
  10. for( i = 0; i < count; i++ )
  11. goodCount += mask[i] = err[i] <= threshold;//goodCount为计算出的内点的个数
  12. return goodCount;
  13. }

7.上面的函数调用的一些函数，这些函数不难，所以下面相应的列举一下

[cpp] view plain copy    
  1. bool CvHomographyEstimator::refine( const CvMat* m1, const CvMat* m2, CvMat* model, int maxIters )
  2. {
  3. CvLevMarq solver(8, 0, cvTermCriteria(CV_TERMCRIT_ITER+CV_TERMCRIT_EPS, maxIters, DBL_EPSILON));
  4. int i, j, k, count = m1->rows*m1->cols;
  5. const CvPoint2D64f* M = (const CvPoint2D64f*)m1->data.ptr;
  6. const CvPoint2D64f* m = (const CvPoint2D64f*)m2->data.ptr;
  7. CvMat modelPart = cvMat( solver.param->rows, solver.param->cols, model->type, model->data.ptr );
  8. cvCopy( &modelPart, solver.param );
  9. for(;;)
  10. {
  11. const CvMat* _param = 0;
  12. CvMat *_JtJ = 0, *_JtErr = 0;
  13. double* _errNorm = 0;
  14. if( !solver.updateAlt( _param, _JtJ, _JtErr, _errNorm ))
  15. break;
  16. for( i = 0; i < count; i++ )
  17. {
  18. const double* h = _param->data.db;
  19. double Mx = M[i].x, My = M[i].y;
  20. double ww = h[6]*Mx + h[7]*My + 1.;
  21. ww = fabs(ww) > DBL_EPSILON ? 1./ww : 0;
  22. double _xi = (h[0]*Mx + h[1]*My + h[2])*ww;
  23. double _yi = (h[3]*Mx + h[4]*My + h[5])*ww;
  24. double err[] = { _xi - m[i].x, _yi - m[i].y };
  25. if( _JtJ || _JtErr )
  26. {
  27. double J[][8] =
  28. {
  29. { Mx*ww, My*ww, ww, 0, 0, 0, -Mx*ww*_xi, -My*ww*_xi },
  30. { 0, 0, 0, Mx*ww, My*ww, ww, -Mx*ww*_yi, -My*ww*_yi }
  31. };
  32. for( j = 0; j < 8; j++ )
  33. {
  34. for( k = j; k < 8; k++ )
  35. _JtJ->data.db[j*8+k] += J[0][j]*J[0][k] + J[1][j]*J[1][k];
  36. _JtErr->data.db[j] += J[0][j]*err[0] + J[1][j]*err[1];
  37. }
  38. }
  39. if( _errNorm )
  40. *_errNorm += err[0]*err[0] + err[1]*err[1];
  41. }
  42. }
  43. cvCopy( solver.param, &modelPart );
  44. return true;
  45. }
  46. void CvHomographyEstimator::computeReprojError( const CvMat* m1, const CvMat* m2,
  47. const CvMat* model, CvMat* _err )
  48. {
  49. int i, count = m1->rows*m1->cols;
  50. const CvPoint2D64f* M = (const CvPoint2D64f*)m1->data.ptr;
  51. const CvPoint2D64f* m = (const CvPoint2D64f*)m2->data.ptr;
  52. const double* H = model->data.db;
  53. float* err = _err->data.fl;
  54. for( i = 0; i < count; i++ )
  55. {
  56. double ww = 1./(H[6]*M[i].x + H[7]*M[i].y + 1.);
  57. double dx = (H[0]*M[i].x + H[1]*M[i].y + H[2])*ww - m[i].x;
  58. double dy = (H[3]*M[i].x + H[4]*M[i].y + H[5])*ww - m[i].y;
  59. err[i] = (float)(dx*dx + dy*dy);
  60. }
  61. }

8，最后一部分是比较关键的。就是FindHomography函数本身。这个函数又去调用了cvFindHomography函数，估计就是openCV不同版本的函数吧，其实现的功能和思想都是一样的。这个函数内部基本上也就是做一些判断防止溢出，排查错误，检验变量以及变换格式等辅助性的内容，真正的方法性质的代码还是在上面的提到的CvHomographyEstimator类中。

[cpp] view plain copy    
  1. cv::Mat findHomography( InputArray _points1, InputArray _points2,
  2. int method, double ransacReprojThreshold, OutputArray _mask)
  3. {
  4. Mat points1 = _points1.getMat(), points2 = _points2.getMat();
  5. int npoints = points1.checkVector(2);//返回矩阵的序列个数
  6. CV_Assert( npoints >= 0 && points2.checkVector(2) == npoints &&
  7. points1.type() == points2.type());  //检验初始条件是否正确
  8. Mat H(3, 3, CV_64F);
  9. CvMat _pt1 = points1, _pt2 = points2;
  10. CvMat matH = H, c_mask, *p_mask = 0;
  11. if( _mask.needed() )
  12. {
  13. _mask.create(npoints, 1, CV_8U, -1, true);
  14. p_mask = &(c_mask = _mask.getMat());
  15. }
  16. bool ok = cvFindHomography( &_pt1, &_pt2, &matH, method, ransacReprojThreshold, p_mask ) > 0;  //函数调用
  17. if( !ok )
  18. H = Scalar(0);
  19. return H;
  20. }
  21. CV_IMPL int
  22. cvFindHomography( const CvMat* objectPoints, const CvMat* imagePoints,
  23. CvMat* __H, int method, double ransacReprojThreshold,
  24. CvMat* mask )
  25. {
  26. const double confidence = 0.995;
  27. const int maxIters = 2000;                           //修改这里来修改迭代次数
  28. const double defaultRANSACReprojThreshold = 3;
  29. bool result = false;
  30. Ptr<CvMat> m, M, tempMask;
  31. double H[9];
  32. CvMat matH = cvMat( 3, 3, CV_64FC1, H );    //这就是单应矩阵，矩阵初始化
  33. int count;
  34. CV_Assert( CV_IS_MAT(imagePoints) && CV_IS_MAT(objectPoints) );
  35. count = MAX(imagePoints->cols, imagePoints->rows);    //序列个数
  36. CV_Assert( count >= 4 );
  37. if( ransacReprojThreshold <= 0 )
  38. ransacReprojThreshold = defaultRANSACReprojThreshold;
  39. m = cvCreateMat( 1, count, CV_64FC2 );
  40. cvConvertPointsHomogeneous( imagePoints, m );  //转换齐次坐标
  41. M = cvCreateMat( 1, count, CV_64FC2 );
  42. cvConvertPointsHomogeneous( objectPoints, M );
  43. if( mask )
  44. {
  45. CV_Assert( CV_IS_MASK_ARR(mask) && CV_IS_MAT_CONT(mask->type) &&
  46. (mask->rows == 1 || mask->cols == 1) &&
  47. mask->rows*mask->cols == count );
  48. }
  49. if( mask || count > 4 )
  50. tempMask = cvCreateMat( 1, count, CV_8U );
  51. if( !tempMask.empty() )
  52. cvSet( tempMask, cvScalarAll(1.) );
  53. CvHomographyEstimator estimator( MIN(count, 4) );   //参数是一个小于等于4的值，只有大于4，才能用RANSAC计算
  54. if( count == 4 )
  55. method = 0;
  56. if( method == CV_LMEDS )
  57. //result = estimator.runLMeDS( M, m, &matH, tempMask, confidence, maxIters );
  58. printf("");
  59. else if( method == CV_RANSAC )
  60. result = estimator.runRANSAC( M, m, &matH, tempMask, ransacReprojThreshold, confidence, maxIters);
  61. else
  62. result = estimator.runKernel( M, m, &matH ) > 0;
  63. if( result && count > 4 )
  64. {
  65. icvCompressPoints( (CvPoint2D64f*)M->data.ptr, tempMask->data.ptr, 1, count );  //压缩，使序列紧凑
  66. count = icvCompressPoints( (CvPoint2D64f*)m->data.ptr, tempMask->data.ptr, 1, count );
  67. M->cols = m->cols = count;    //筛选过后，这个count是内点的个数
  68. if( method == CV_RANSAC )
  69. estimator.runKernel( M, m, &matH );  //重新计算最终的单应矩阵，matH
  70. estimator.refine( M, m, &matH, 10 );
  71. }
  72. if( result )
  73. cvConvert( &matH, __H );
  74. if( mask && tempMask )
  75. {
  76. if( CV_ARE_SIZES_EQ(mask, tempMask) )    //复制这个矩阵
  77. cvCopy( tempMask, mask );
  78. else
  79. cvTranspose( tempMask, mask );        //行列调换的 复制这个矩阵
  80. }
  81. return (int)result;
  82. }

OpenCV中的RANSAC详解相关推荐

  1. OpenCV 中Iplimage结构详解

    Iplimage数据结构 主要困扰我许久的就是其中的widthStep与width*nChannels是否相等,事实上我们可以在源码opencv\modules\core\src\array.cpp中 ...

  2. OpenCV 中cv2.threshold详解，（大白话版）

    cv2.threshold (src, thresh, maxval, type) cv2.threshold (源图片, 阈值, 填充色, 阈值类型) 函数含义请看这篇博客:https://blog ...

  3. OpenCV中的waitkey()详解

    函数原型: int waitkey(int delay = 0) 其中delay表示延迟时间. 1) 当delay <= 0时,表示如果没有键盘时间,一直处于等待状态.这一般用于程序暂停. 2) ...

  4. Python-Matplotlib可视化（番外篇）——Matplotlib中的事件处理详解与实战

    Python-Matplotlib可视化(番外篇)--Matplotlib中的事件处理详解与实战 前言 事件连接 事件属性 实战1:直方图中矩形的拖拽 实战2:鼠标进入和离开 相关链接与参考 前言 在 ...

  5. ALSA声卡驱动中的DAPM详解之四：在驱动程序中初始化并注册widget和route

    前几篇文章我们从dapm的数据结构入手,了解了代表音频控件的widget,代表连接路径的route以及用于连接两个widget的path.之前都是一些概念的讲解以及对数据结构中各个字段的说明,从本章开 ...

  6. Asp.net中GridView使用详解(引)【转】

    Asp.net中GridView使用详解(引) GridView无代码分页排序 GridView选中,编辑,取消,删除 GridView正反双向排序 GridView和下拉菜单DropDownList ...

  7. Linux中iptraf命令详解(IP局域网监控工具)

    2019独角兽企业重金招聘Python工程师标准>>> Linux中iptraf命令详解(IP局域网监控工具) 发布时间:2017-12-27 20:46:03   作者:佚名    ...

  8. ArcGIS Engine中的Symbols详解

    转自原文 ArcGIS Engine中的Symbols详解 本文由本人翻译ESRI官方帮助文档.尊重劳动成果,转载请注明来源. Symbols ArcObjects用了三种类型的Symbol(符号样式 ...

  9. js路由在php上面使用,React中路由使用详解

    这次给大家带来React中路由使用详解,React中路由使用的注意事项有哪些,下面就是实战案例,一起来看一下. 路由 通过 URL 映射到对应的功能实现,React 的路由使用要先引入 react-r ...

最新文章

  1. 自定义一个安全的rm指令
  2. 定制键盘输入处理(1503)
  3. xcode symbol(s) not found for architecture i386错误解决方法
  4. 数据库连接池和Tomcat连接池的配置问题
  5. wps的高亮怎么用_怎样在WPS上实现代码语法高亮
  6. 基于python的图形化邮件发送程序（支持添加附件）
  7. 实现定时中断_EPIT 定时器，仅需做到如下几步，即可轻松配置使用
  8. savexml php,PHP DOMDocument saveXML()用法及代码示例
  9. STEAM 97%好评，体验堪比《杀戮尖塔》，为什么玩家说这是2020年上半年最超值的游戏？
  10. 黄聪： 50 个 Bootstrap 插件
  11. 不安全的反序列化_CVE202027131 思科安全管理器反序列化漏洞 POC
  12. 设计模式1【续】：动态设定行为
  13. php处理excel类,30 个 PHP 的 Excel 处理类
  14. Docker run centos 内部使用systemctl 启动服务的方法
  15. 补单平台哪个靠谱 天猫补单哪个安全
  16. netsetman使用教程_人性化的IP切换工具—NetSetMan
  17. NanoPi NEO3上手日记第一天——把玩＆刷固件
  18. python 图像相似度；用0-1矩阵表示两幅图像的相似度
  19. 真假内推？直拿offer？别被无良中介给骗了
  20. DAS、NAS、SAN简介以及区别

热门文章

  1. 【一周头条盘点】中国软件网（2018.6.25~2018.6.29）
  2. 如何成为一个自律的人？这5本书里藏着你想要的答案
  3. nosql之mongodb的数据库操作+集合的插入和更新操作
  4. 信号波峰波谷二阶差分识别算法
  5. 不识字也能翻译：谷歌AI直接用音频翻音频，不用先转文本
  6. github Pull请求(Pull request)
  7. Delphi字符串转换成代码
  8. 程序猿开发大牛成长记 | 龙果社区有奖征文活动
  9. PyCharm的使用
  10. 基于单片机的室内智能照明系统设计(#0419)