RANSAC消除SIFT错配
RANSAC为RANdom SAmple Consensus的缩写，它是根据一组包含异常数据的样本数据集，计算出数据的数学模型参数，得到有效样本数据的算法。它于1981年由 Fischler和Bolles最先提出[1]。
RANSAC算法的基本假设是样本中包含正确数据(inliers，可以被模型描述的数据)，也包含异常数据(Outliers，偏离正常范 围很远、无法适应数学模型的数据)，即数据集中含有噪声。这些异常数据可能是由于错误的测量、错误的假设、错误的计算等产生的。同时RANSAC也假设， 给定一组正确的数据，存在可以计算出符合这些数据的模型参数的方法。
RANSAC基本思想描述如下：
①考虑一个最小抽样集的势为n的模型(n为初始化模型参数所需的最小样本数)和一个样本集P，集合P的样本数num(P)>n，从P中随机 抽取包含n个样本的P的子集S初始化模型M；
②余集SC=P\S中与模型M的误差小于某一设定阈值t的样本集以及S构成S*。S*认为是内点集，它们构成S的一致集(Consensus Set)；
③若#(S*)≥N，认为得到正确的模型参数，并利用集S*(内点inliers)采用最小二乘等方法重新计算新的模型M*；重新随机抽取新 的S，重复以上过程。
④在完成一定的抽样次数后，若为找到一致集则算法失败，否则选取抽样后得到的最大一致集判断内外点，算法结束。
由上可知存在两个可能的算法优化策略。①如果在选取子集S时可以根据某些已知的样本特性等采用特定的选取方案或有约束的随机选取来代替原来的 完全随机选取；②当通过一致集S*计算出模型M*后，可以将P中所有与模型M*的误差小于t的样本加入S*，然后重新计算M*。
RANSAC算法包括了3个输入的参数：①判断样本是否满足模型的误差容忍度t。t可以看作为对内点噪声均方差的假设，对于不同的输入数据需 要采用人工干预的方式预设合适的门限，且该参数对RANSAC性能有很大的影响；②随机抽取样本集S的次数。该参数直接影响SC中样本参与模型参数的检验 次数，从而影响算法的效率，因为大部分随机抽样都受到外点的影响；③表征得到正确模型时，一致集S*的大小N。为了确保得到表征数据集P的正确模型，一般 要求一致集足够大；另外，足够多的一致样本使得重新估计的模型参数更精确。
RANSAC算法经常用于计算机视觉中。例如，在立体视觉领域中同时解决一对相机的匹配点问题及基本矩阵的计算。
其伪代码如下
1 input:
2 data - a set of observations
3 model - a model that can be fitted to data
4 n - the minimum number of data required to fit the model
5 k - the number of iterations performed by the algorithm
6 t - a threshold value for determining when a datum fits a model
7 d - the number of close data values required to assert that a model fits well to data
8 output:
9 best_model - model parameters which best fit the data (or nil if no good model is found)
10 best_consensus_set - data point from which this model has been estimated
11 best_error - the error of this model relative to the data
12
13 iterations := 0
14 best_model := nil
15 best_consensus_set := nil
16 best_error := infinity
17 while iterations < k
18 maybe_inliers := n randomly selected values from data
19 maybe_model := model parameters fitted to maybe_inliers
20 consensus_set := maybe_inliers
21
22 for every point in data not in maybe_inliers
23 if point fits maybe_model with an error smaller than t
24 add point to consensus_set
25
26 if the number of elements in consensus_set is > d
27 (this implies that we may have found a good model,
28 now test how good it is)
29 this_model := model parameters fitted to all points in consensus_set
30 this_error := a measure of how well better_model fits these points
31 if this_error < best_error
32 (we have found a model which is better than any of the previous ones,
33 keep it until a better one is found)
34 best_model := this_model
35 best_consensus_set := consensus_set
36 best_error := this_error
37
38 increment iterations
39
40 return best_model, best_consensus_set, best_error
下面是C语言实现
View Code
1 void Ransac(IplImage *pload,IplImage *pnew,int *x_cord,int *y_cord,int length, double &c0,double &c1,double &c2,char *file_name_prefix2,int nums,int &min)
2 {
3 int *count = new int[RANSAC_TIMES]; //计数
4 memset(count,0x00,sizeof(int)*RANSAC_TIMES);
5 CvMat mat_x, mat_b, mat_c;
6 double p_mat_x[9], p_mat_b[3],p_mat_c[3];//要转化成数组
7 min = length; //求最小值用
8 int flag_linear =0;
9 int position = 0; //记录误差最小的位置
10
11 uchar *loadimageData = (uchar*)pload->imageData;
12 int step = pload->widthStep;
13 double a =0,b=0,c =0;
14 for(int i =0; i<RANSAC_TIMES; i++)
15 {
16 int randnums[3];
17 int n = 0;
18
19 for (int t=0;t<3;t++)
20 {
21 n = rand() % length;
22 randnums[t] = n;
23 p_mat_x[t*3] = 1.0;
24 p_mat_x[t*3+1] = x_cord[n];
25 p_mat_x[t*3+2] = x_cord[n] * x_cord[n];
26
27 p_mat_b[t] = y_cord[n];
28 p_mat_c[t] = 0;
29 }
30 cvInitMatHeader(&mat_x,3,3,CV_64FC1,p_mat_x);
31 cvInitMatHeader(&mat_b,3,1,CV_64FC1,p_mat_b);
32 cvInitMatHeader(&mat_c,3,1,CV_64FC1,p_mat_c);
33
34 flag_linear = cvSolve(&mat_x, &mat_b, &mat_c,CV_LU);
35 if ( flag_linear == 0)
36 {
37 continue;
38 }
39 double *temp = mat_c.data.db; //结果保存下来。
40 c = temp[0]; //常数项 //保留第i次结果
41 b = temp[1]; //一次项
42 a = temp[2]; //平方项
43 if (c2 < 0)
44 {
45 c2 = -c2;
46 }
47 double y_value = 0;
48 double x_square = 0;
49 int y_floor = 0;
50
51 double dis_error = 0;
52
53 for (int j=0;j<length;j++)
54 {
55 dis_error = fabs(y_cord[j] - ( a*(x_cord[j]*x_cord[j]) + b*x_cord[j] + c ));
56 if ( dis_error > RANSAC_THRESHOLD)
57 {
58 count[i]++;
59 }
60 }
61
62 if (min > count[i])
63 {
64 min = count[i];
65 position = i;
66 c0 = c;
67 c1 = b;
68 c2 = a;
69 }
70
71 } //50次循环结束。
72
73 double x_square = 0;
74 double y_value = 0;
75 int y_floor = 0;
76 int pixelposition;
77 for (int x=0; x<pnew->width; x++) //自变量取值范围
78 {
79 x_square = pow(double(x),2);
80 y_value = c2*x_square + c1*x + c0;
81 y_floor = cvFloor((y_value));
82
83 if ( y_floor >= 0 && y_floor < pnew->height )
84 {
85 ((uchar*)(pnew->imageData + pnew->widthStep*y_floor))[3*x]=255;
86 ((uchar*)(pnew->imageData + pnew->widthStep*y_floor))[3*x+1]=0;
87 ((uchar*)(pnew->imageData + pnew->widthStep*y_floor))[3*x+2]=0;
88 }
89 for (int y=0;y<pnew->height;y++)
90 {
91 pixelposition = y*step+x;
92 if ( 0 == loadimageData[pixelposition])
93 {
94 ((uchar*)(pnew->imageData + pnew->widthStep*y))[3*x]= 0;
95 ((uchar*)(pnew->imageData + pnew->widthStep*y))[3*x+1]=0;
96 ((uchar*)(pnew->imageData + pnew->widthStep*y))[3*x+2]=255;
97 }
98 }
99 }
100
101 delete []count;
102 cvSaveImage(file_name_prefix2,pnew);
103
104 }
