最小二乘法直線擬合 java matlababla最小二乘法直線擬合
2023-06-13
最小二乘擬合
最小二乘擬合是一種數(shù)學(xué)上的相似性和優(yōu)化，利用已知數(shù)據(jù)獲得一條直線或曲線，使其在坐標(biāo)系上與已知數(shù)據(jù)之間的距離達(dá)到平方和最小。
2.RANSAC算法
3，直線擬合
一般方程AX在建模時(shí)使用直線。 BY C=0，隨機(jī)選擇兩點(diǎn)構(gòu)建直線模型，計(jì)算每個(gè)點(diǎn)到此直線的TLS。（Total Least Square），當(dāng)TLS低于一定閥值時(shí)，點(diǎn)是符合模型的點(diǎn)，當(dāng)點(diǎn)數(shù)最多時(shí)，模型是最好的直線模型。然后根據(jù)此時(shí)的直線參數(shù)繪制最終擬合直線。
4.橢圓擬合
橢圓定義方程用于建立模型：dist（P,A） dist（P,B）=DIST，P是橢圓的上一點(diǎn)，A和B是橢圓的兩個(gè)焦點(diǎn)。隨機(jī)選擇三點(diǎn)A，B,P構(gòu)建橢圓模型，計(jì)算從每個(gè)點(diǎn)到這兩個(gè)焦點(diǎn)的距離和與DIST的差異。當(dāng)差值低于一定閥值時(shí)，點(diǎn)為符合模型的點(diǎn)，點(diǎn)數(shù)最多時(shí)，模型為最佳橢圓模型，然后根據(jù)符合條件的點(diǎn)使用橢圓一般方程Ax2。 Bxy Cy2 Dx Ey F=0 根據(jù)函數(shù)式繪制最終擬合橢圓，并獲得符合點(diǎn)進(jìn)行指數(shù)擬合。
5.代碼matlab
(1)最小二乘擬合
View Code
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FILENAME LSF.m
% FUNCTION Input points with mouse,Least-squares fit of lines to
% 2D points
% DATE 2012-10-12
% AUTHOR zhangying
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
%% 鼠標(biāo)輸入點(diǎn)，結(jié)束enter鍵
axis([-10 10 -10 10]);
[x,y]=ginput; %讀出坐標(biāo)，直到按下回車，回到X，y坐標(biāo)點(diǎn)，
num=length(x); %計(jì)算點(diǎn)數(shù)量
%% 用最小二乘直接擬合
%找出數(shù)據(jù)的最佳函數(shù)，通過最小化誤差的平方來匹配
[p1,s1]=polyfit(x,y,1); %n=1為直線擬合 x,y是一個(gè)數(shù)據(jù)點(diǎn)，n是一個(gè)多項(xiàng)式級(jí)別，回到p是一個(gè)由高到低的多項(xiàng)式指數(shù)向量。
[p2,s2]=polyfit(x,y,num-2); %n>1為曲線擬合，找出頻率為n的多項(xiàng)式，對(duì)于數(shù)據(jù)點(diǎn)集（x,y）,平方和最小滿足差
[p3,s3]=polyfit(x,y,num-1); %x一定要無聊。在polyval中使用矩陣s來估計(jì)誤差。
xcurve=-10:2:10; %在這些點(diǎn)上尋求多項(xiàng)式數(shù)值。
p1curve=polyval(p1,xcurve); %多項(xiàng)式曲線求值，回到相應(yīng)的自變量xcurve給出指數(shù)P的多項(xiàng)式數(shù)值
p2curve=polyval(p2,xcurve);
p3curve=polyval(p3,xcurve);
%% 繪圖
plot(xcurve,p1curve,'g-',xcurve,p2curve,'b-',...
xcurve,p3curve,'r-',x,y,'*');
title(不同次數(shù)的最小二乘擬合');
legend('degree 1','degree num-2','degree num-1','points');基于RANSAC的直線擬合
View Code
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FILENAME RANSACF.m
% FUNCTION RANSAC fit of lines to 2D points
% DATE 2012-10-11
% AUTHOR zhangying
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear;
%% 隨機(jī)點(diǎn)生成
g_NumOfPoints = 500; % 點(diǎn)數(shù)
g_ErrPointPart = 0.4; % 噪音
g_NormDistrVar = 1; % 標(biāo)準(zhǔn)偏差
% 生成隨機(jī)點(diǎn)
theta = (rand(1) 1) * pi/6;
R = ( rand([1 g_NumOfPoints]) - 0.5) * 100;
DIST = randn([1 g_NumOfPoints]) * g_NormDistrVar;
Data = [cos(theta); sin(theta)] * R [-sin(theta); cos(theta)] * DIST;
Data(:, 1:floor(g_ErrPointPart * g_NumOfPoints)) = 2 * [max(abs(Data(1,:))), 0; 0, max(abs(Data(2,:)))] *...
(rand([2 floor(g_ErrPointPart * g_NumOfPoints)]) - 0.5);
plot(Data(1, :), Data(2, :), '.', 'Tag', 'DATA');
hold on;
%% 擬合RANSAC
% 擬合模型初始化
nSampLen = 2; %根據(jù)點(diǎn)數(shù)設(shè)置模型所依據(jù)的點(diǎn)數(shù)
nIter = 50; %最大循環(huán)次數(shù)
dThreshold = 2; %殘差閾值
nDataLen = size(Data, 2); %數(shù)據(jù)長(zhǎng)度
RANSAC_model = NaN; %跳過缺失模型
RANSAC_mask = zeros([1 nDataLen]); %全0矩陣，1表示符合模型，0表示不符合
nMaxInlyerCount = -1;
%% 主循環(huán)
for i = 1:nIter
% 取樣，選兩個(gè)不同的點(diǎn)
SampleMask = zeros([1 nDataLen]);
while sum( SampleMask ) ~= nSampLen
ind = ceil(nDataLen .* rand(1, nSampLen - sum(SampleMask))); %rand產(chǎn)生隨機(jī)數(shù)，ceil向離它最近的大整數(shù)圓整數(shù)。
SampleMask(ind) = 1;
end
Sample = find( SampleMask ); %找出非零元素的索引值，也就是建立模型的點(diǎn)
%% 建立模型,并且找出符合模型的點(diǎn)
ModelSet = feval（@TLS, Data(:, Sample)); %所有的計(jì)算都符合模型點(diǎn)最小二乘。
for iModel = 1:size(ModelSet, 3)
CurModel = ModelSet(:, :, iModel); %當(dāng)前模型對(duì)應(yīng)的直線參數(shù)
CurMask =( abs( CurModel * [Data; ones([1 size(Data, 2)])])< dThreshold);%直線距離低于閥值的點(diǎn)符合模型，標(biāo)記為1
nCurInlyerCount = sum(CurMask); %計(jì)算符合直線模型點(diǎn)數(shù)的數(shù)量。
%% 選最好的模型
if nCurInlyerCount > nMaxInlyerCount %最好的模型是符合模型點(diǎn)數(shù)最多的模型。
nMaxInlyerCount = nCurInlyerCount;
RANSAC_mask = CurMask;
RANSAC_model = CurModel;
end
end
end
%% 繪制最佳模型擬合結(jié)果
MinX=min(Data(1, :));
MaxX=max(Data(1, :));
MinX_Y=-(RANSAC_model(1).*MinX RANSAC_model(3))./RANSAC_model(2);
MaxX_Y=-(RANSAC_model(1).*MaxX RANSAC_model(3))./RANSAC_model(2);
plot([MinX MaxX],[MinX_Y MaxX_Y],'r-');
title(在噪聲條件下，ransac直線擬合');
%% 使用RANSAC方法擬合原理。
% RANSAC算法的輸入是一組可以解釋或適應(yīng)觀測(cè)數(shù)據(jù)的觀測(cè)數(shù)據(jù)的參數(shù)模型，以及一些可靠的參數(shù)值。
% RANSAC通過在數(shù)據(jù)中反復(fù)選擇一組隨機(jī)子集來實(shí)現(xiàn)目標(biāo)。
% RANSAC通過在數(shù)據(jù)中反復(fù)選擇一組隨機(jī)子集來實(shí)現(xiàn)目標(biāo)。選定的子集被假設(shè)為對(duì)局點(diǎn)，并通過以下方法進(jìn)行驗(yàn)證:
% 有一種模型適用于假設(shè)的對(duì)局點(diǎn)，即所有未知參數(shù)都可以從假設(shè)的對(duì)局點(diǎn)計(jì)算出來。
% 用1中獲得的模型對(duì)其他所有數(shù)據(jù)進(jìn)行測(cè)試，如果某一點(diǎn)適用于估計(jì)模型，則認(rèn)為它也是對(duì)局點(diǎn)。
% 如果有足夠多的點(diǎn)被歸類為假設(shè)的對(duì)局點(diǎn)，那么估計(jì)模型就足夠合理了。
% 接著，用所有假設(shè)的對(duì)局點(diǎn)再一次估計(jì)模型，因?yàn)樗皇潜怀跏技僭O(shè)的對(duì)局點(diǎn)估計(jì)過。
% 最后，通過估計(jì)對(duì)局點(diǎn)和模型之間的差錯(cuò)率來評(píng)估模型。
% 這個(gè)過程被反復(fù)執(zhí)行固定頻率，每次生成的模型要么因?yàn)橛螒螯c(diǎn)太少而被放棄，要么因?yàn)楸痊F(xiàn)在的模型更好而被選中。
%% 問題分析
% 關(guān)于畫直線：根據(jù)抽取的兩點(diǎn)畫出最后一條直線，結(jié)果不穩(wěn)定，每一條直線運(yùn)行后都會(huì)有一定的誤差。
% 直線參數(shù)值已在系統(tǒng)中計(jì)算，根據(jù)Ax計(jì)算。 By C=0，可選擇數(shù)據(jù)點(diǎn)中最左邊的點(diǎn)和最右邊的點(diǎn)來確定最后的直線，相對(duì)穩(wěn)定。
% 2.調(diào)用函數(shù)A時(shí)，如果有函數(shù)B作為參數(shù)，則在函數(shù)B之前添加@，函數(shù)B參數(shù)值另外傳輸，函數(shù)的執(zhí)行方式可以通過feval統(tǒng)一。
% 關(guān)于隨機(jī)點(diǎn)的生成：rand的應(yīng)用靈活多變。View Code
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FILENAME TLS.m
% FUNCTION Calculate Total Least Squares of input data
% DATE 2012-10-11
% AUTHOR unkown
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Total Least Squares TLS(Data)
%Return: [a,b,c] - line in a * x b * y c = 0 form
% where a ^ 2 b ^ 2 = 1
% TLS(X,Y) - approximates ALL points in array by one line
function line = TLS(Data)
if any( size(Data) == 0)
Line = [0, 0, 0];
return;
end
X = Data(1, :);
Y = Data(2, :);
len = length(X);
if size(X) ~= size(Y)
X = X';
end
M = [ mid(X .^ 2) - mid(X) ^ 2, sum(X .* Y) / len - mid(X) * mid(Y);...
sum(X .* Y) / len - mid(X) * mid(Y), mid(Y .^ 2) - mid(Y) ^ 2];
[ev,tmp] = eig(M);
ind = 2;
if ErrFunc(X, Y, ev(:, 1)) < ErrFunc(X, Y, ev(:, 2))
ind = 1;
end
line = [ev(1,ind), ev(2,ind),...
-ev(1,ind) * mid(X) - ev(2,ind) * mid(Y)];
return;
% Help function, calculates an error
function e = ErrFunc(X,Y,L)
maxE = 0;
e = 0;
c = -L(1) * mid(X) - L(2) * mid(Y);
for i = 1:length(X)
e = e ( L(1) * X(i) L(2) * Y(i) c )^2;
if (L(1) * X(i) L(2) * Y(i) c )^2 > maxE
maxE = (L(1) * X(i) L(2) * Y(i) c )^2;
end;
end
return;
% Middle value of vector X
function l = mid(X)
if length(X) > 0
l = sum(X) / length(X);
else
l = 0;
end;
return;以RANSAC為基礎(chǔ)的橢圓擬合
View Code
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FILENAME ellipsefit.m
% FUNCTIPN Least-squares fit of ellipse to 2D points
% DATE 2012-10-12
% AUTHOR zhangying
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
clc;
clear;
%% 生成 橢圓形帶噪音
% 參數(shù)初始化
g_NumOfPoints = 500; % 點(diǎn)數(shù)
g_ErrPointPart = 0.4; % 噪音
g_NormDistrVar = 3; % 標(biāo)準(zhǔn)偏差
a=10;b=20; %長(zhǎng)軸短軸
angle=60; %傾斜角
%% 橢圓生成
beta = angle * (pi / 180);
alpha = linspace(0, 360, g_NumOfPoints) .* (pi / 180);
X = (a * cos(alpha) * cos(beta)- b * sin(alpha) * sin(beta) ) wgn(1,length(alpha),g_NormDistrVar^2,'linear');
Y = (a * cos(alpha) * sin(beta) b * sin(alpha) * cos(beta) ) wgn(1,length(alpha),g_NormDistrVar^2,'linear');
Data=[X;Y];
plot(Data(1, :), Data(2, :), '.', 'Tag', 'DATA');
hold on;
%% 橢圓擬合結(jié)合RANSAC
%一般橢圓方程：Ax2 Bxy Cy2 Dx Ey F=0
%F=@(p,x)p(1)*x(:,1).^2 p(2)*x(:,1).*x(:,2) p(3)*x(:,2).^2 p(4)*x(:,1) p(5)
%% 參數(shù)初始化
nSampLen = 3; %根據(jù)點(diǎn)數(shù)設(shè)置模型所依據(jù)的點(diǎn)數(shù)
nDataLen = size(Data, 2); %數(shù)據(jù)長(zhǎng)度
nIter = 50; %最大循環(huán)次數(shù)
dThreshold = 2; %殘差閾值
nMaxInlyerCount=-1; %點(diǎn)數(shù)下限
A=zeros([2 1]);
B=zeros([2 1]);
P=zeros([2 1]);
%% 主循環(huán)
for i = 1:nIter
SampleMask = zeros([1 nDataLen]);
while sum( SampleMask ) ~= nSampLen
ind = ceil(nDataLen .* rand(1, nSampLen - sum(SampleMask))); %取樣，選擇nSampLen的不同點(diǎn)
SampleMask(ind) = 1;
end
Sample = find( SampleMask ); %找出索引值的非零元素，也就是建立模型的點(diǎn)
%% 建立模型，存儲(chǔ)建模所需的坐標(biāo)點(diǎn)、焦點(diǎn)和一個(gè)過橢圓點(diǎn)。
%橢圓定義方程：兩個(gè)定點(diǎn)之間的距離和常數(shù)
A(:,1)=Data(:,ind(1)); %焦點(diǎn)
B(:,1)=Data(:,ind(2)); %焦點(diǎn)
P(:,1)=Data(:,ind(3)); %橢圓上一點(diǎn)
DIST=sqrt((P(1,1)-A(1,1)).^2 (P(2,1)-A(2,1)).^2) sqrt((P(1,1)-B(1,1)).^2 (P(2,1)-B(2,1)).^2);
xx=[];
nCurInlyerCount=0; %初始化點(diǎn)數(shù)為0
%% 是否符合模型？
for k=1:g_NumOfPoints
CurModel=[A(1,1) A(2,1) B(1,1) B(2,1) DIST ];
pdist=sqrt((Data(1,k)-A(1,1)).^2 (Data(2,k)-A(2,1)).^2) sqrt((Data(1,k)-B(1,1)).^2 (Data(2,k)-B(2,1)).^2);
CurMask =(abs(DIST-pdist)< dThreshold); %直線距離低于閥值的點(diǎn)符合模型，標(biāo)記為1
nCurInlyerCount =nCurInlyerCount CurMask; %橢圓模型點(diǎn)的計(jì)算數(shù)量符合橢圓模型點(diǎn)數(shù)。
if(CurMask==1)
xx =[xx,Data(:,k)];
end
end
%% 選最好的模型
if nCurInlyerCount > nMaxInlyerCount %最好的模型是符合模型點(diǎn)數(shù)最多的模型。
nMaxInlyerCount = nCurInlyerCount;
Ellipse_mask = CurMask;
Ellipse_model = CurModel;
Ellipse_points = [A B P];
Ellipse_x =xx;
end
end
%% 橢圓由符合點(diǎn)擬合
%一般橢圓方程：Ax2 Bxy Cy2 Dx Ey F=0
F=@(p,x)p(1)*x(:,1).^2 p(2)*x(:,1).*x(:,2) p(3)*x(:,2).^2 p(4)*x(:,1) p(5)*x(:,2) p(6);
p0=[1 1 1 1 1 1];
x=Ellipse_x';
pr=nlinfit(x,zeros(size(x,1),1),F,p0); % 擬合指數(shù)，最小二乘法
xmin=min(x(:,1));
xmax=max(x(:,1));
ymin=min(x(:,2));
ymax=max(x(:,2));
%% 畫點(diǎn)作圖
plot(Ellipse_points(1,:),Ellipse_points(2,:),'r*');
hold on;
plot(Ellipse_x(1,:),Ellipse_x(2,:),'yo');
hold on;
ezplot(@(x,y)F(pr,[x,y]),[-1 xmin,1 xmax,-1 ymin,1 ymax]);
title(橢圓擬合'RANSAC);
legend('樣本點(diǎn)'，'抽取點(diǎn)'，'符合點(diǎn)'，'擬合曲線')
%% 問題分析
% 關(guān)于如何生成隨機(jī)點(diǎn)：基于標(biāo)準(zhǔn)橢圓，增加高斯白噪音--wgn();
% 如何建立橢圓模型：
% 方案一：橢圓通用方程：Ax2 Bxy Cy2 Dx Ey F=0.如果你認(rèn)為橢圓是由五個(gè)點(diǎn)決定的，那么橢圓擬合是由五個(gè)點(diǎn)帶入方程式進(jìn)行的，結(jié)果大部分都是
% 在繪制雙曲線的前提下，放棄；
% 方案二:使用橢圓定義:兩個(gè)定點(diǎn)之間的距離和常數(shù)為2a，選擇平面中的三個(gè)點(diǎn)、兩個(gè)焦點(diǎn)和一個(gè)過橢圓點(diǎn)來確定橢圓。
% 3.如何篩選符合條件的點(diǎn):此時(shí)，從計(jì)點(diǎn)到橢圓的距離過于復(fù)雜。如果定義，如果兩個(gè)焦點(diǎn)之間的距離和與2a之間的距離低于一定的閥值，則符合要求。
% 3.如何篩選符合條件的點(diǎn):此時(shí)，從計(jì)點(diǎn)到橢圓的距離過于復(fù)雜。如果定義，如果兩個(gè)焦點(diǎn)之間的距離和與2a之間的距離低于一定的閥值，則符合要求。
% 關(guān)于擬合函數(shù)：使用nlinfit，對(duì)輸入?yún)?shù)的維數(shù)有要求，需要x為N*P維，y為N*P維，n*一維，注意列向量。
% 5.關(guān)于如何畫橢圓:不像一般的繪圖指定X和Y，這個(gè)時(shí)候需要畫函數(shù)圖形。如果你在網(wǎng)上找到它，你應(yīng)該先建立函數(shù)F，然后利用它。
% ezplot(@(x,y)F(pr,[x,y])）可以畫出函數(shù)圖形。
% 數(shù)據(jù)是隨機(jī)生成的，每一個(gè)程序的運(yùn)行結(jié)果都會(huì)有所不同，在B中(:,1)=Data(:,ind(2));有時(shí)候數(shù)組越界會(huì)出錯(cuò)，但是單步調(diào)試沒有問題，
% 還沒有找到原因。
%%6.學(xué)習(xí)經(jīng)驗(yàn)
(1)不能過于依賴當(dāng)前的函數(shù)，需要多自己去思考算法，即使借用別人的函數(shù)，也要弄清楚原理和調(diào)用方法；
(2)matlab函數(shù)庫不熟悉，要多使用help。；
在編寫程序時(shí)，首先要對(duì)程序結(jié)構(gòu)進(jìn)行整體規(guī)劃，模塊化，組織清晰，自己要了解自己流程的每一步原理。
(4)理論的力量是無窮無盡的，代碼設(shè)計(jì)應(yīng)該基于理論的深刻理解。
本文為轉(zhuǎn)載內(nèi)容，我們尊重原作者對(duì)本文的作權(quán)。如有內(nèi)容錯(cuò)誤或侵權(quán)問題，歡迎原作者聯(lián)系我們更正或刪除內(nèi)容。
本文僅代表作者觀點(diǎn)，版權(quán)歸原創(chuàng)者所有，如需轉(zhuǎn)載請(qǐng)?jiān)谖闹凶⒚鱽碓醇白髡呙帧?/p>
免責(zé)聲明：本文系轉(zhuǎn)載編輯文章，僅作分享之用。如分享內(nèi)容、圖片侵犯到您的版權(quán)或非授權(quán)發(fā)布，請(qǐng)及時(shí)與我們聯(lián)系進(jìn)行審核處理或刪除，您可以發(fā)送材料至郵箱：service@tojoy.com