一、RANSAC介紹
随機抽樣一緻算法（RANdom SAmple Consensus,RANSAC）,采用疊代的方式從一組包含離群的被觀測資料中估算出數學模型的參數。維基百科中的RANSAC該算法最早由Fischler和Bolles于1981年提出。RANSAC算法假設資料中包含正确資料和異常資料（或稱為噪聲）。正确資料記為内點（inliers），異常資料記為外點（outliers）。同時RANSAC也假設，給定一組正确的資料，存在可以計算出符合這些資料的模型參數的方法。該算法核心思想就是随機性和假設性，随機性是根據正确資料出現機率去随機選取抽樣資料，根據大數定律，随機性模拟可以近似得到正确結果。假設性是假設選取出的抽樣資料都是正确資料，然後用這些正确資料通過問題滿足的模型，去計算其他點，然後對這次結果進行一個評分。
RANSAC算法被廣泛應用在計算機視覺領域和數學領域，例如直線拟合、平面拟合、計算圖像或點雲間的變換矩陣、計算基礎矩陣等方面，使用的非常多。本文将在對RANSAC介紹完後，附兩段直線拟合以及平面拟合的matlab代碼。關于計算機視覺中基于RANSAC架構的矩陣求解問題，在OpenCV中都有對應的函數接口，如果以後有機會再把這個整理出來。
二、算法基本思想
如下圖所示，存在很多離散的點，而我們認為這些點構成了一條直線。當然，人眼能很清晰地拟合出這條直線，找到外點。但要讓計算機找到這條直線，在很久之前是很難的，RACSAC的出現是通過數學之美解決這一難題的重要發明。
通過執行個體對算法基本思想進行描述：
（1）首先，我們知道要得到一個直線模型，我們需要兩個點唯一确定一個直線方程。是以第一步我們随機選擇兩個點。
（2）通過這兩個點，我們可以計算出這兩個點所表示的模型方程y=ax+b。
（3）我們将所有的資料點套到這個模型中計算誤差。
（4）我們找到所有滿足誤差門檻值的點，第4幅圖中可以看到，隻有很少的點支援這個模型。
（5）然後我們再重複（1）~（4）這個過程，直到達到一定疊代次數後，我們選出那個被支援的最多的模型，作為問題的解。如下圖所示：
以上是對RANSAC算法的基本思想的介紹，我們可以發現，雖然這個資料集中外點和内點的比例幾乎相等，但是RANSAC算法還是能找到最合适的解。試想一下，這個問題如果使用最小二乘法進行優化，由于噪聲資料的幹擾，我們得到的結果肯定是一個錯誤的結果，如下圖所示。這是由于最小二乘法是一個将外點參與讨論的代價優化問題，而RANSAC是一個使用内點進行優化的問題。經實驗驗證，對于包含80%誤差的資料集，RANSAC的效果遠優于直接的最小二乘法。
三、RANSAC的數學推導
我們假設内點在整個資料集中的機率為t，即：
确定該問題的模型需要n個點，這個n是根據問題而定義的，例如拟合直線時n為2，平面拟合時n=3，求解點雲之間的剛性變換矩陣時n=3，求解圖像之間的射影變換矩陣是n=4，等等。k表示疊代次數，即随機選取n個點計算模型的次數。P為在這些參數下得到正确解的機率。為友善表示，我們可以得到，n個點都是内點的機率為tntn，則n個點中至少有一個是外點的機率為
表示k次随機抽樣中都沒有找到一次全是内點的情況，這個時候得到的是錯誤解，也就是：
内點機率t是一個先驗值，可以給出一些魯棒的值。同時也可以看出，即使t給的過于樂觀，也可以通過增加疊代次數k，來保證正确解的機率P。同樣的，我們可以通過上面式子計算出來疊代次數k，即我們假設需要正确機率為P（例如我們需要99%的機率取到正确解），則：
四、利用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;
%%% 繪制資料點
figure;plot(data(1,:),data(2,:),'o');hold on; % 顯示資料點
number = size(data,2); % 總點數
bestParameter1=0; bestParameter2=0; % 最佳比對的參數
sigma = 1;
pretotal=0; %符合拟合模型的資料的個數
for i=1:iter
%%% 随機選擇兩個點
idx = randperm(number,2);
sample = data(:,idx);
%%%拟合直線方程 y=kx+b
line = zeros(1,3);
x = sample(:, 1);
y = sample(:, 2);
k=(y(1)-y(2))/(x(1)-x(2)); %直線斜率
b = y(1) - k*x(1);
line = [k -1 b]
mask=abs(line*[data; ones(1,size(data,2))]); %求每個資料到拟合直線的距離
total=sum(mask<sigma); %計算資料距離直線小于一定門檻值的資料的個數
if total>pretotal %找到符合拟合直線資料最多的拟合直線
pretotal=total;
bestline=line; %找到最好的拟合直線
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);
bestParameter2 = -bestline(3)/bestline(2);
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)]);
結果如圖所示：
五、利用RANSAC拟合平面
clc;clear all;close all;
%%%三維平面拟合
%%%生成随機資料
%内點
mu=[0 0 0]; %均值
S=[2 0 4;0 4 0;4 0 8]; %協方差
data1=mvnrnd(mu,S,300); %産生200個高斯分布資料
%外點
mu=[2 2 2];
S=[8 1 4;1 8 2;4 2 8]; %協方差
data2=mvnrnd(mu,S,100); %産生100個噪聲資料
%合并資料
data=[data1',data2'];
iter = 1000;
%%% 繪制資料點
figure;plot3(data(1,:),data(2,:),data(3,:),'o');hold on; % 顯示資料點
number = size(data,2); % 總點數
bestParameter1=0; bestParameter2=0; bestParameter3=0; % 最佳比對的參數
sigma = 1;
pretotal=0; %符合拟合模型的資料的個數
for i=1:iter
%%% 随機選擇三個點
idx = randperm(number,3);
sample = data(:,idx);
%%%拟合直線方程 z=ax+by+c
plane = zeros(1,3);
x = sample(1,:);
y = sample(2,:);
z = sample(3,:);
a = ((z(1)-z(2))*(y(1)-y(3)) - (z(1)-z(3))*(y(1)-y(2)))/((x(1)-x(2))*(y(1)-y(3)) - (x(1)-x(3))*(y(1)-y(2)));
b = ((z(1) - z(3)) - a * (x(1) - x(3)))/(y(1)-y(3));
c = z(1) - a * x(1) - b * y(1);
plane = [a b -1 c]
mask=abs(plane*[data; ones(1,size(data,2))]); %求每個資料到拟合平面的距離
total=sum(mask<sigma); %計算資料距離平面小于一定門檻值的資料的個數
if total>pretotal %找到符合拟合平面資料最多的拟合平面
pretotal=total;
bestplane=plane; %找到最好的拟合平面
end
end
%顯示符合最佳拟合的資料
mask=abs(bestplane*[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);
inliers(2,k) = data(2,i);
plot3(data(1,i),data(2,i),data(3,i),'r+');
k = k+1;
end
end
%%% 繪制最佳比對平面
bestParameter1 = bestplane(1);
bestParameter2 = bestplane(2);
bestParameter3 = bestplane(4);
xAxis = min(inliers(1,:)):max(inliers(1,:));
yAxis = min(inliers(2,:)):max(inliers(2,:));
[x,y] = meshgrid(xAxis, yAxis);
z = bestParameter1 * x + bestParameter2 * y + bestParameter3;
surf(x, y, z);
title(['bestPlane: z = ',num2str(bestParameter1),'x + ',num2str(bestParameter2),'y + ',num2str(bestParameter3)]);
結果如圖所示：
cankao:
https://blog.csdn.net/u010128736/article/details/53422070
http://www.cnblogs.com/xrwang/archive/2011/03/09/ransac-1.html