RANSAC算法的基本思想和Matlab实现

Posted by Kriz on 2017-04-13

写cv作业的时候涉及到了RANSAC。时间过去得有点久忘了不少,网上的资料也大多比较抽象,所以在这边稍微记录一下自己的理解。


RANSAC是干什么的?

假设分到了个活儿,要求根据输入的两个或更多点坐标,拟合出一条直线来。

两个点的情况当然简单,中学生表示操纵各种直线方程分分钟得解,顺便如果直线上存在更多点,使用相同的方式进行验证也是洒洒水(考虑到计算机性能因素,一般使用向量进行判断)。

更加接近现实情况的当然是多个不在同一条直线上的点集,在噪声不明显的情况下,最小二乘法也足够用了。事实上,对于等式个数大于未知数个数的情况,最小二乘法基本是第一选择,毕竟是就算曲线也拟合给你看的存在。

然而涉猎广泛的最小二乘法也是会翻车的,一个很明显的前提条件便是,它对误差的要求稍显严格。如果在噪声偏多的情况下进行应用,茫茫然不知所措的程序八成会返回一个充满迷幻气息的结果,就像这样:

oops

而RANSAC(Random Sample Consensus,随机抽样一致)就是专门用于在高噪声影响环境下完成此类线性回归的算法。在大多数情况下,它的适配结果都是非常精准且稳定的。


RANSAC是怎么干的?

首先来看一下权威解释,如果有醍醐灌顶茅塞顿开之感,那么之后的文章就不用看了:

RANSAC基本思想描述如下:

①考虑一个最小抽样集的势为n的模型(n为初始化模型参数所需的最小样本数)和一个样本集P,集合P的样本数#(P)>n,从P中随机抽取包含n个样本的P的子集S初始化模型M;

②余集SC=P\S中与模型M的误差小于某一设定阈值t的样本集以及S构成S。S认为是内点集,它们构成S的一致集(Consensus Set);

③若#(S*)≥N,认为得到正确的模型参数,并利用集S*(内点inliers)采用最小二乘等方法重新计算新的模型M*;重新随机抽取新的S,重复以上过程。

④在完成一定的抽样次数后,若未找到一致集则算法失败,否则选取抽样后得到的最大一致集判断内外点,算法结束。

总之我个人一直觉得这个解释的毒性很大,虽然费劲地多读几遍大概是能理解的,然而总有些云里雾里。毕竟不给例子讲道理都是耍流氓,下面的部分我会结合一个示范来进行RANSAC步骤的基本描述。

先让我们初始化一个噪声巨大无比的观测数据点集,大概像酱:

ransac

基本版的RANSAC是个非常暴躁的算法,它会在这一片点里随机取定两个(选取用来拟合模型的项可以自定义,一般取两个):

ransac

两点确定一条直线,于是我们得到了这样的结果。

ransac

好的拟合完成了。

……当然不是这样。随机取到的两个点基本无法成功拟合,但是否能在这样的操作中看出些端倪了?

没错,只要把所有点对都进行一次计算,比较一下它们返回的某个可以用于衡量拟合结果好坏的量,取最优即可。

这个量怎样取比较合适?容易想到,如果大多数点都离我们由两点获得的直线很近,那么这样的拟合就比较令人信服。因此,先计算每个点到这条直线的距离:

ransac

紧接着,我们借由设定好的一个合适的阈值,来判断有多少点是我们能够接受的。

ransac

可以看出,只有3个点到直线的距离是我们能够容忍的。这个数目意味着什么?我们目前只进行了一次这个过程,还没有参考。接下来,另外一个随机点对会被选取,并重复相同的过程:

ransac

当所有的点对都被测试完成,我们选取其中可接受距离阈值内包含点最多的拟合模型,当做最终的结果模型。可以看到,虽然算法本身很狂躁,但效果还是很不错的。

ransac


Matlab实现RANSAC

稍微翻了翻好像没发现现成的轮子,好像很多都集成到Computer Vision System Toolbox里了。那么手撸一个好了。

需求里给了少量离散点,所以我在这里就直接把它们包含进来了。如果要进行更严谨的测试,可以使用mvnrnd函数进行有效数据和噪声的模拟,这里不细说。

1
2
3
4
5
6
7
8
% Input coordinates
data = [-2 0; 0 0.9; 2 2.0; 3 6.5; 4 2.9; 5 8.8; 6 3.95; 8 5.03; 10 5.97; 12 7.1; 13 1.2; 14 8.2; 16 8.5; 18 10.1];
figure;
plot(data(:,1), data(:,2), 'ro');
axis([-2 18 0 12])
axis equal
hold on

初始化:

1
2
3
4
5
% Main process
max = 100; % Number of loops
thresh = 1; % Bearable distance to line
pretotal = 0; % Number of points which are inside the threshold

循环体的处理,最优点对数据存在两个perfect里:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
for i = 1:max
rand_index = round(randi([1 14],1,2));
if rand_index(1) == rand_index(2)
continue;
end
pA = [data(rand_index(1),1) data(rand_index(1),2)];
pB = [data(rand_index(2),1) data(rand_index(2),2)];
get_line = generate_line([pA;pB]);
dist=abs(get_line*[data ones(size(data,1),1)]');
total = sum(dist<thresh);
if total > pretotal
pretotal = total;
perfect_pA = pA;
perfect_pB = pB;
end
end

最后绘制直线。讲道理绘制过程不应该这么麻烦的,然而line函数不知道怎么回事一直调教失败,就只能臃肿而不优雅地直接plot线段了:

1
2
3
4
5
6
% Draw line
k = (perfect_pB(2)-perfect_pA(2))/(perfect_pB(1)-perfect_pA(1));
p_start = [data(1,1), perfect_pA(2)-k*(perfect_pA(1)-data(1,1))];
p_end = [data(end,1), perfect_pA(2)-k*(perfect_pA(1)-data(end,1))];
plot([p_start(1),p_end(1)],[p_start(2),p_end(2)]);

generate_line函数的定义,用于输出直线一般式的参数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function get_line=generate_line(data)
x = data(:,1);
y = data(:,2);
k=(y(1)-y(2))/(x(1)-x(2));
a=sqrt(1-1/(1+k^2));
b=sqrt(1-a^2);
if k>0
b=-b;
end
c=-a*x(1)-b*y(1);
get_line=[a b c];
end

输出为成功拟合的直线,显示正常。


参考文献

http://www.cnblogs.com/yin52133/archive/2012/07/21/2602562.html

http://blog.sina.com.cn/s/blog_875c3b2f0106huux.html

http://www.cnblogs.com/tiandsp/archive/2013/06/03/3115743.html