轉載：RANSAC的演算法 （ 原貼有少量錯誤，訂正連結下面）
阿新 • • 發佈：2018-11-26
由於文章是轉載的，無法在原文編輯。更正後的程式碼連結： https://blog.csdn.net/qq_35752161/article/details/84502542
（2）通過這兩個點，我們可以計算出這兩個點所表示的模型方程y=ax+b。
（3）我們將所有的資料點套到這個模型中計算誤差。
（4）我們找到所有滿足誤差閾值的點，第4幅圖中可以看到，只有很少的點支援這個模型。
（5）然後我們再重複（1）~（4）這個過程，直到達到一定迭代次數後，我們選出那個被支援的最多的模型，作為問題的解。如下圖所示：
RANSAC介紹（Matlab版直線擬合+平面擬合）
本人郵箱：[email protected]，歡迎交流討論，
歡迎轉載，轉載請註明網址http://blog.csdn.net/u010128736/
一、RANSAC介紹
隨機抽樣一致演算法（RANdom SAmple Consensus,RANSAC）,採用迭代的方式從一組包含離群的被觀測資料中估算出數學模型的引數。 維基百科中的RANSAC該演算法最早由Fischler和Bolles於1981年提出。RANSAC演算法假設資料中包含正確資料和異常資料（或稱為噪聲）。正確資料記為內點（inliers），異常資料記為外點（outliers）。同時RANSAC也假設，給定一組正確的資料，存在可以計算出符合這些資料的模型引數的方法。該演算法核心思想就是隨機性和假設性，隨機性是根據正確資料出現概率去隨機選取抽樣資料，根據大數定律，隨機性模擬可以近似得到正確結果。假設性是假設選取出的抽樣資料都是正確資料，然後用這些正確資料通過問題滿足的模型，去計算其他點，然後對這次結果進行一個評分。
RANSAC演算法被廣泛應用在計算機視覺領域和數學領域，例如直線擬合、平面擬合、計算影象或點雲間的變換矩陣、計算基礎矩陣
二、演算法基本思想
如下圖所示，存在很多離散的點，而我們認為這些點構成了一條直線。當然，人眼能很清晰地擬合出這條直線，找到外點。但要讓計算機找到這條直線，在很久之前是很難的，RACSAC的出現是通過數學之美解決這一難題的重要發明。
通過例項對演算法基本思想進行描述：
（1）首先，我們知道要得到一個直線模型，我們需要兩個點唯一確定一個直線方程。所以第一步我們隨機
（2）通過這兩個點，我們可以計算出這兩個點所表示的模型方程y=ax+b。
（3）我們將所有的資料點套到這個模型中計算誤差。
（4）我們找到所有滿足誤差閾值的點，第4幅圖中可以看到，只有很少的點支援這個模型。
（5）然後我們再重複（1）~（4）這個過程，直到達到一定迭代次數後，我們選出那個被支援的最多的模型，作為問題的解。如下圖所示：
以上是對RANSAC演算法的基本思想的介紹，我們可以發現，雖然這個資料集中外點和內點的比例幾乎相等，但是RANSAC演算法還是能找到最合適的解。試想一下，這個問題如果使用最小二乘法進行優化，由於噪聲資料的干擾，我們得到的結果肯定是一個錯誤的結果，如下圖所示。這是由於最小二乘法是一個將外點參與討論的代價優化問題，而RANSAC是一個使用內點進行優化的問題。經實驗驗證，對於包含80%誤差的資料集，RANSAC的效果遠優於直接的最小二乘法。
三、RANSAC的數學推導
我們假設內點在整個資料集中的概率為t，即：
四、利用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)]);
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
結果如圖所示：
五、利用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)]);
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
結果如圖所示：
