RANSAC 算法介绍 - Zhening's Blog

RANSAC 算法

背景

在美留学学习传统计算机视觉课程的时候，曾经做过一个非常有趣的项目，通过拍摄几张连续的图片生成一张全景照片（如果用过早期的Google手机应该对这种全景照片的生成过程比较熟悉，并非苹果手机那种连续滑动的拍摄方式）。这个项目涉及到了很多有趣可聊的部分，包括摄像机参数的标定，图像畸变的矫正，SIFT算法标记特征点，特征点的匹配，RANSAC算法寻找合适的运动模型以及图像融合和裁剪，其中SIFT算法和RANSAC算法在整个过程中起到了至关重要的作用，而RANSAC算法又完全由自己的代码实现，所以这里结合计算机视觉背景对RANSAC算法做一个简单的回顾。

介绍

RANSAC是Ransom Sample Consensus的缩写，中文我们一般叫做随机抽样一致算法，是Martin A. Fischler 和Robert C. Bolles于1981年发表在ACM期刊上的论文《Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography》^[1] 中阐述的一种用于将模型拟合到实验数据的新范式。论文中提到说在很大程度上，场景分析（甚至于是大部分的科学研究）实际上是根据一组预定义的模型来解释感测到的数据。而数据集中总是包含正确的数据和异常的数据，异常数据由于错误的测量，错误的假设或者错误的计算而导致的误差会影响我们寻找正确的模型。RANSAC算法可以很好地过滤样本中离群点对于模型估计的影响。

论文中有一段很小的篇幅描述了RANSAC算法的实现过程。

Given a model that requires a minimum of n data points to instantiate its free parameters, and a set of data points P such that the number of points in P is greater than n [#(P) >= n], randomly select a subset S1 of n data points from P and instantiate the model. Use the instantiated model M1 to determine the subset S1* of points in P that are within some error tolerance of M1. The set S1* is called the consensus set of S1. If #(S1*) is greater than some threshold t, which is a function of the estimate of the number of gross errors in P, use S1* to compute (possibly using least squares) a new model M1*. If # (S1*) is less than t, randomly select a new subset S2 and repeat the above process. If, after some predetermined number of trials, no consensus set with t or more members has been found, either solve the model with the largest consensus set found, or terminate in failure.

翻译一下，大致可以分成以下几个步骤。

随机选择一个属于数据集P的子集S1，S1有n个元素（n小于P的总元素数，但n的个数能唯一确定一个模型），用S1中的数据初始化一个模型M1
用这个模型M1，在一定的容错范围内，计算得到另一个子集S1* (内点集)
如果S1*中的元素大于某个阈值，我们就直接使用这个模型（注：通常不这么做，不如设置阈值，直接进入第四布）
如果S1*中的元素小于某个阈值，我们就重新随机选取一个子集S2，重复上面的步骤，直到完成我们设定的循环次数。取内点集元素最多的对应模型当做我们的结果。

RANSAC算法其实不难理解，通俗来讲就是在一定的条件下随机生成一个模型，然后通过类似选举机制，选出票数最多的模型。核心思想在于它的随机性和假设性，平等对待每一个数据点，属于一种不确定的算法，有一定的概率得出一个合理的结果。提高模型的准确性程度依赖于循环次数和允许的容错范围，这两个参数在整个过程中需要人为选取。

结合项目

在这个项目里面，我们已经通过SIFT算法获得了两张相邻图像的特征点，并且通过特征点匹配算法获得了对应的特征点对。RANSAC算法在这里起到的作用就是使用这些特征点对，过滤掉那些明显错误的特征点对，然后生成一个图像的运动模型（这里我们只考虑两个自由度的运动，即只在X和Y轴方向上移动图像，不考虑图像旋转和尺寸变化）。后续的步骤就可以使用这个运动模型使两张相邻图像重叠在一起，在通过图像的融合算法（Alpha-Blending等）生成一张完整的图片。

图片占位

对应到RANSAC算法的实现，可以分为以下几部分。

随机选取一个特征点对（两个特征点就可以计算出一个运动模型，只有两个自由度），由这个个特征点对计算出一个临时的运动模型。
计算当前运动模型和其他特征点对之间欧拉距离的差值，如果小于我们设定的阈值就判断这个特征点同意当前模型，记录下有多少特征点同意当前模型
循环上面两步，不断更新记录获得对多特征点同意的模型，得到我们最终的运动模型。

下面看代码再来理解一下。

/****************************************
* alignPair:
*	输入:
*		f1, f2: 		两个源特征点集
*		matches: 		特征点对 （保存对应的f1和f2的id）
*		m: 				运动模型
*		f: 				相机参数
*		nRANSAC: 		循环次数
*		RANSACthresh: 	阈值
*		M: 				图形学变换矩阵 (output)
*	逻辑:
*		repeat for nRANSAC iterations:
*			choose a minimal set of feature matches
*			estimate the transformation implied by these matches
*			count the number of inliers
*		for the transformation with the maximum number of inliers,
*		compute the least squares motion estimate using the inliers,
*		and store it in M
*/
int alignPair(const FeatureSet &f1, const FeatureSet &f2,
              const vector<FeatureMatch> &matches, MotionModel m, float f,
              int nRANSAC, double RANSACthresh, CTransform3x3& M)
{
    // 初始化随机种子
	srand((unsigned int)(time(NULL)));
	
	int maxInliersCount = 0;
	int match_id = 0;
	CTransform3x3 tempM;
	vector<int> inliers;
	vector<int> tempInliers;
	
	for(int i = 0;i < nRANSAC;i++){ // 循环次数需要我们设定 tunning
        // 随机选取一个特征点对
		match_id = rand()%(matches.size());
		int id1 = matches[match_id].id1 - 1;  // f1的id
		int id2 = matches[match_id].id2 - 2;  // f2的id
        // 计算出一个临时的运动模型
		int u = f2[id2].x - f1[id1].x;  // X 轴方向的偏移量
		int v = f2[id2].y - f1[id1].y;  // Y 轴方向的偏移量
		tempM[0][2] = u;
		tempM[1][2] = v;
		
        // 计算同意该模型的特征点对数量
		int inliersCount = countInliers(f1, f2, matches, m, f, tempM, RANSACthresh, tempInliers);
		// 更新同意该模型的特征点对数量的最大值
		if (inliersCount > maxInliersCount) {
			inliers = tempInliers;
			maxInliersCount = inliersCount;
		}
	}
    // 计算图形学变换矩阵
	leastSquaresFit(f1, f2, matches, m, f, inliers,M);
    return 0;
}
/****************************************
* countInliers:
*	输入:
*		f1, f2: 		两个源特征点集
*		matches: 		特征点对 （保存对应的f1和f2的id）
*		m: 				运动模型
*		f: 				相机参数
*		RANSACthresh: 	阈值
*		M: 				图形学变换矩阵 (output)
*		inliers:  		inlier feature IDs
*	逻辑:
*		transform the matched features in f1 by M
*
*		count the number of matching features for which the transformed
*		feature f1[id1-1] is within SSD distance RANSACthresh of its match
*		f2[id2-1]
*
*		store the indices of these matches in inliers
*
*		
*/
int countInliers(const FeatureSet &f1, const FeatureSet &f2,
                 const vector<FeatureMatch> &matches, MotionModel m, float f,
                 CTransform3x3 M, double RANSACthresh, vector<int> &inliers)
{
    inliers.clear();
    int count = 0;
	
    // 循环所有特征点对
    for (unsigned int i=0; i<(int) matches.size(); i++) {
        // determine if the ith matched feature f1[id1-1], when transformed by M,
        // is within RANSACthresh of its match in f2
        //
        // if so, increment count and append i to inliers
        //

		//feature index 
		int id1 = matches[i].id1 - 1;
		int id2 = matches[i].id2 - 1;
		//compute original image add translation
		int xTranslation = f1[id1].x + M[0][2];
		int yTranslation = f1[id1].y + M[1][2];
		//compute distance error with align image
		int xDistance = f2[id2].x - xTranslation;
		int yDistance = f2[id2].y - yTranslation;
		double errorDistance = pow(xDistance,2) + pow(yDistance,2);
		
		//if error is small take it as inlier
		if(errorDistance < pow(RANSACthresh,2)){ // 与阈值进行比较 tunning
			count++;
			inliers.push_back(i);
		}
    }

    return count;
}

/******************这部分与RANSAC无关**********************
* leastSquaresFit:
*	INPUT:
*		f1, f2: 		两个源特征点集
*		matches: 		特征点对 （保存对应的f1和f2的id）
*		m: 				运动模型
*		f: 				相机参数
*		M: 				图形学变换矩阵 (output)
*		inliers:  		inlier feature IDs
*	OUTPUT:
*		compute the transformation from f1 to f2 using only the inliers
*		and return it in M
*/
int leastSquaresFit(const FeatureSet &f1, const FeatureSet &f2,
                    const vector<FeatureMatch> &matches, MotionModel m, float f,
                    const vector<int> &inliers, CTransform3x3& M)
{
    // the transformation is a translation and
    // only has two degrees of freedom
    //
    // therefore, we simply compute the average translation vector
    // between the feature in f1 and its match in f2 for all inliers
    double u = 0;
    double v = 0;

    for (int i=0; i<inliers.size(); i++) {
        double xTrans, yTrans;

        // compute the translation implied by the ith inlier match
        // and store it in (xTrans,yTrans)
		//feature index
		int id1 = matches[inliers[i]].id1 - 1; 
		int id2 = matches[inliers[i]].id2 - 1;
		//compute distance
		xTrans = f2[id2].x - f1[id1].x;
		yTrans = f2[id2].y - f1[id1].y;	
        
        u += xTrans;
        v += yTrans;
    }

    u /= inliers.size();
    v /= inliers.size();

    M[0][0] = 1;
    M[0][1] = 0;
    M[0][2] = -u;
    M[1][0] = 0;
    M[1][1] = 1;
    M[1][2] = -v;
    M[2][0] = 0;
    M[2][1] = 0;
    M[2][2] = 1;

    return 0;
}

参考

Fischler M A , Bolles R C . Random sample consensus : a paradigm for model fitting with applications to image analysis and automated cartography[J]. Commun. Assoc. Comp. March, 1981, 24(6):381-395. ↩