利用RANSAC方法找到二维与三维的最佳模型（仅线性）

        本文代码是对RANSAC介绍（Matlab版直线拟合+平面拟合）中代码的改进，RANSAC介绍（Matlab版直线拟合+平面拟合）中的代码没有完整地依据RANSAC的步骤进行，所以做了一些修改，加入了每次重新评估模型的部分。

        还是提一些RANSAC的原理和步骤，O(∩_∩)O~

        RANSAC 是在强干扰的情况下寻找以某种模型（如直线）出现的数据。为什么说是在强干扰下呢？？我们知道，最小二乘法可以用来建立模型拟合数据点，但是最小二乘法的目标是找到与所有点的整体误差最小的模型，也就是尽量去适应所有点。

        但这对于下面这种情形是不能使用的：（本图来源：我的数学之美（一）——RANSAC算法详解）



        很明显，我们用人眼可以看出，符合我们视觉的那条直线并不是通过最小二乘法得到的图中红色的那条。

       RANSAC就不同了，它可以获得所有点中最大的支持集合，用这个最大的支持集合来求出最佳的模型。

        RANSAC的流程如下：

        1、从数据集中随机选取一个子集，认为这其中的点都为局内点。

        2、找到一个模型适应于（怎么算适应呢？？用距离阈值或者其他方法）这些假设的局内点，也就是找到了例如y=k0x+b0

        3、用找到的这个模型测试数据集中所有的点，如果某个点适应于这个模型，即就暂且认为它是局内点。

        4、如果有足够多的点被归为这个模型，就认为这个模型合理。

        5、用当前归入该模型的所有局内点重新评估该模型（how？？比如通过这些局内点和模型的总体误差）。记录这个error。

        6、回到1，迭代固定次数，找到最小的error对应的模型。

        由此，得到最佳模型。 

        原博中的代码是没有第6步的，但我认为通过error评估模型是很重要的一步。仅仅根据归为模型的局内点的数目判断模型是否可靠并不妥当，因为对于一个模型的好坏，我们往往会从误差角度来判断一个模型对数据的拟合程度，拟合程度（对局内点）越好，才能说明这个模型越可靠。所以我在代码中加入了这一步。

        完整代码及对照注释如下：（原博对点的坐标似乎也搞错了，所以也做了改正）

%利用RANSAC拟合直线
clc;clear all;close all;

%%%二维直线拟合
%%%生成随机数据
%内点
mu=[0 0];  %均值
S=[1 2.5;2.5 8];  %协方差
data1=mvnrnd(mu,S,200);   %产生200个高斯分布数据
%外点
mu=[2 2];
S=[8 0;0 8];
data2=mvnrnd(mu,S,100);     %产生100个噪声数据
%合并数据
data=[data1',data2'];
iter = 100; %设置RANSAC固定的迭代次数，也就是随机选取多少次点的子集

%%% 绘制数据点
figure;plot(data(1,:),data(2,:),'o');hold on; % 显示数据点
number = size(data,2); % 总点数
bestParameter1=0; bestParameter2=0; % 初始化最佳匹配的参数
sigma = 1;
pretotal=0;     %初始化符合拟合模型的数据的个数
min_error = inf;     %初始化error
for i=1:iter
%%%随机选取一个子集，认为子集内的点都是局内点
idx = randperm(number,2);  %随机选取一个子集，即两个点
sample = data(:,idx);

%%%找到一个模型适应于这些假设的局内点（即刚才随机选取的两个点）
line = zeros(1,3);%拟合直线方程 y=kx+b
p1 = sample(:, 1);%第一个点
p2 = sample(:, 2);%第二个点

%      k=(y(1)-y(2))/(x(1)-x(2));      %直线斜率
%      b = y(1) - k*x(1);
%      line = [k -1 b]

k=(p2(2)-p1(2))/(p2(1)-p1(1));
b=p1(2) - k*p1(1);
line = [k -1 b];

%%%用找到的这个模型去测试其他所有的点，如果某个点适应于这个模型，就暂且认为它是局内点，这样得到依据此模型得到的局内点的数目
mask=abs(line*[data; ones(1,size(data,2))]);    %求每个数据到拟合直线的距离  对于每个点  相当于 | kxi+b-yi | 这种方式很巧妙
point_idx = zeros(1,size(data,2));
point_idx(mask<sigma) =1;
total=sum(point_idx);              %计算数据距离直线小于一定阈值的数据的个数

%     %%%如果有足够多的点被归类为这个模型，就认为这个模型合理。
%     %(本代码没有用当前归入的所有的局内点（利用错误率）去重新估计这个模型的过程)
%      if total>pretotal            %找到符合拟合直线数据最多的拟合直线
%          pretotal=total;
%          bestline=line;          %找到最好的拟合直线
%      end

%%%如果有足够多的点被归类为这个模型，就认为这个模型合理。
%(有重新评估模型的过程，即，用当前归入的所有的局内点（利用错误率）去重新估计这个模型)
if total>pretotal            %找到符合拟合直线数据最多的拟合直线
data_inliter = data(:,find(point_idx==1));
error =  sum(abs(line*[data_inliter ; ones(1,size(data_inliter ,2))])); %错误率
if error < min_error
pretotal=total;
bestline=line;          %找到最好的拟合直线
min_error = error;
end
end

end
%
%显示符合最佳拟合的数据
mask=abs(bestline*[data; ones(1,size(data,2))])<sigma;    %找到适应于这个最佳模型的点
hold on;
k=1;
for i=1:length(mask)
if mask(i)
inliers(1,k) = data(1,i);
k=k+1;
plot(data(1,i),data(2,i),'+');%适应于这个最佳模型的点被标为+
end
end

%%% 绘制最佳匹配曲线
bestParameter1 = -bestline(1)/bestline(2);  %最佳k
bestParameter2 = -bestline(3)/bestline(2);  %最佳b
xAxis = min(inliers(1,:)):max(inliers(1,:));
yAxis = bestParameter1*xAxis + bestParameter2;
plot(xAxis,yAxis,'r-','LineWidth',2);%画出最佳匹配曲线
title(['bestLine:  y =  ',num2str(bestParameter1),'x + ',num2str(bestParameter2)]);
fprintf('局内点的数目为：%d\n',pretotal);
fprintf('min_error为：%f\n',min_error);

        得到的结果如下：



        此时的局内点（误差率最小情况下，最大的支持集合）为：110个