最小二乘法直线拟合 java matlababla最小二乘法直线拟合
2023-06-13
最小二乘拟合
最小二乘拟合是一种数学上的相似性和优化，利用已知数据获得一条直线或曲线，使其在坐标系上与已知数据之间的距离达到平方和最小。
2.RANSAC算法
3，直线拟合
一般方程AX在建模时使用直线。 BY C=0，随机选择两点构建直线模型，计算每个点到此直线的TLS。（Total Least Square），当TLS低于一定阀值时，点是符合模型的点，当点数最多时，模型是最好的直线模型。然后根据此时的直线参数绘制最终拟合直线。
4.椭圆拟合
椭圆定义方程用于建立模型：dist（P,A） dist（P,B）=DIST，P是椭圆的上一点，A和B是椭圆的两个焦点。随机选择三点A，B,P构建椭圆模型，计算从每个点到这两个焦点的距离和与DIST的差异。当差值低于一定阀值时，点为符合模型的点，点数最多时，模型为最佳椭圆模型，然后根据符合条件的点使用椭圆一般方程Ax2。 Bxy Cy2 Dx Ey F=0 根据函数式绘制最终拟合椭圆，并获得符合点进行指数拟合。
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;
%% 鼠标输入点，结束enter键
axis([-10 10 -10 10]);
[x,y]=ginput; %读出坐标，直到按下回车，回到X，y坐标点，
num=length(x); %计算点数量
%% 用最小二乘直接拟合
%找出数据的最佳函数，通过最小化误差的平方来匹配
[p1,s1]=polyfit(x,y,1); %n=1为直线拟合 x,y是一个数据点，n是一个多项式级别，回到p是一个由高到低的多项式指数向量。
[p2,s2]=polyfit(x,y,num-2); %n>1为曲线拟合，找出频率为n的多项式，对于数据点集（x,y）,平方和最小满足差
[p3,s3]=polyfit(x,y,num-1); %x一定要无聊。在polyval中使用矩阵s来估计误差。
xcurve=-10:2:10; %在这些点上寻求多项式数值。
p1curve=polyval(p1,xcurve); %多项式曲线求值，回到相应的自变量xcurve给出指数P的多项式数值
p2curve=polyval(p2,xcurve);
p3curve=polyval(p3,xcurve);
%% 绘图
plot(xcurve,p1curve,'g-',xcurve,p2curve,'b-',...
xcurve,p3curve,'r-',x,y,'*');
title(不同次数的最小二乘拟合');
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;
%% 随机点生成
g_NumOfPoints = 500; % 点数
g_ErrPointPart = 0.4; % 噪音
g_NormDistrVar = 1; % 标准偏差
% 生成随机点
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; %根据点数设置模型所依据的点数
nIter = 50; %最大循环次数
dThreshold = 2; %残差阈值
nDataLen = size(Data, 2); %数据长度
RANSAC_model = NaN; %跳过缺失模型
RANSAC_mask = zeros([1 nDataLen]); %全0矩阵，1表示符合模型，0表示不符合
nMaxInlyerCount = -1;
%% 主循环
for i = 1:nIter
% 取样，选两个不同的点
SampleMask = zeros([1 nDataLen]);
while sum( SampleMask ) ~= nSampLen
ind = ceil(nDataLen .* rand(1, nSampLen - sum(SampleMask))); %rand产生随机数，ceil向离它最近的大整数圆整数。
SampleMask(ind) = 1;
end
Sample = find( SampleMask ); %找出非零元素的索引值，也就是建立模型的点
%% 建立模型,并且找出符合模型的点
ModelSet = feval（@TLS, Data(:, Sample)); %所有的计算都符合模型点最小二乘。
for iModel = 1:size(ModelSet, 3)
CurModel = ModelSet(:, :, iModel); %当前模型对应的直线参数
CurMask =( abs( CurModel * [Data; ones([1 size(Data, 2)])])< dThreshold);%直线距离低于阀值的点符合模型，标记为1
nCurInlyerCount = sum(CurMask); %计算符合直线模型点数的数量。
%% 选最好的模型
if nCurInlyerCount > nMaxInlyerCount %最好的模型是符合模型点数最多的模型。
nMaxInlyerCount = nCurInlyerCount;
RANSAC_mask = CurMask;
RANSAC_model = CurModel;
end
end
end
%% 绘制最佳模型拟合结果
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算法的输入是一组可以解释或适应观测数据的观测数据的参数模型，以及一些可靠的参数值。
% RANSAC通过在数据中反复选择一组随机子集来实现目标。
% RANSAC通过在数据中反复选择一组随机子集来实现目标。选定的子集被假设为对局点，并通过以下方法进行验证:
% 有一种模型适用于假设的对局点，即所有未知参数都可以从假设的对局点计算出来。
% 用1中获得的模型对其他所有数据进行测试，如果某一点适用于估计模型，则认为它也是对局点。
% 如果有足够多的点被归类为假设的对局点，那么估计模型就足够合理了。
% 接着，用所有假设的对局点再一次估计模型，因为它只是被初始假设的对局点估计过。
% 最后，通过估计对局点和模型之间的差错率来评估模型。
% 这个过程被反复执行固定频率，每次生成的模型要么因为游戏点太少而被放弃，要么因为比现在的模型更好而被选中。
%% 问题分析
% 关于画直线：根据抽取的两点画出最后一条直线，结果不稳定，每一条直线运行后都会有一定的误差。
% 直线参数值已在系统中计算，根据Ax计算。 By C=0，可选择数据点中最左边的点和最右边的点来确定最后的直线，相对稳定。
% 2.调用函数A时，如果有函数B作为参数，则在函数B之前添加@，函数B参数值另外传输，函数的执行方式可以通过feval统一。
% 关于随机点的生成：rand的应用灵活多变。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为基础的椭圆拟合
View Code
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FILENAME ellipsefit.m
% FUNCTIPN Least-squares fit of ellipse to 2D points
% DATE 2012-10-12
% AUTHOR zhangying
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
clc;
clear;
%% 生成 椭圆形带噪音
% 参数初始化
g_NumOfPoints = 500; % 点数
g_ErrPointPart = 0.4; % 噪音
g_NormDistrVar = 3; % 标准偏差
a=10;b=20; %长轴短轴
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;
%% 椭圆拟合结合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)
%% 参数初始化
nSampLen = 3; %根据点数设置模型所依据的点数
nDataLen = size(Data, 2); %数据长度
nIter = 50; %最大循环次数
dThreshold = 2; %残差阈值
nMaxInlyerCount=-1; %点数下限
A=zeros([2 1]);
B=zeros([2 1]);
P=zeros([2 1]);
%% 主循环
for i = 1:nIter
SampleMask = zeros([1 nDataLen]);
while sum( SampleMask ) ~= nSampLen
ind = ceil(nDataLen .* rand(1, nSampLen - sum(SampleMask))); %取样，选择nSampLen的不同点
SampleMask(ind) = 1;
end
Sample = find( SampleMask ); %找出索引值的非零元素，也就是建立模型的点
%% 建立模型，存储建模所需的坐标点、焦点和一个过椭圆点。
%椭圆定义方程：两个定点之间的距离和常数
A(:,1)=Data(:,ind(1)); %焦点
B(:,1)=Data(:,ind(2)); %焦点
P(:,1)=Data(:,ind(3)); %椭圆上一点
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; %初始化点数为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); %直线距离低于阀值的点符合模型，标记为1
nCurInlyerCount =nCurInlyerCount CurMask; %椭圆模型点的计算数量符合椭圆模型点数。
if(CurMask==1)
xx =[xx,Data(:,k)];
end
end
%% 选最好的模型
if nCurInlyerCount > nMaxInlyerCount %最好的模型是符合模型点数最多的模型。
nMaxInlyerCount = nCurInlyerCount;
Ellipse_mask = CurMask;
Ellipse_model = CurModel;
Ellipse_points = [A B P];
Ellipse_x =xx;
end
end
%% 椭圆由符合点拟合
%一般椭圆方程：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); % 拟合指数，最小二乘法
xmin=min(x(:,1));
xmax=max(x(:,1));
ymin=min(x(:,2));
ymax=max(x(:,2));
%% 画点作图
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('样本点'，'抽取点'，'符合点'，'拟合曲线')
%% 问题分析
% 关于如何生成随机点：基于标准椭圆，增加高斯白噪音--wgn();
% 如何建立椭圆模型：
% 方案一：椭圆通用方程：Ax2 Bxy Cy2 Dx Ey F=0.如果你认为椭圆是由五个点决定的，那么椭圆拟合是由五个点带入方程式进行的，结果大部分都是
% 在绘制双曲线的前提下，放弃；
% 方案二:使用椭圆定义:两个定点之间的距离和常数为2a，选择平面中的三个点、两个焦点和一个过椭圆点来确定椭圆。
% 3.如何筛选符合条件的点:此时，从计点到椭圆的距离过于复杂。如果定义，如果两个焦点之间的距离和与2a之间的距离低于一定的阀值，则符合要求。
% 3.如何筛选符合条件的点:此时，从计点到椭圆的距离过于复杂。如果定义，如果两个焦点之间的距离和与2a之间的距离低于一定的阀值，则符合要求。
% 关于拟合函数：使用nlinfit，对输入参数的维数有要求，需要x为N*P维，y为N*P维，n*一维，注意列向量。
% 5.关于如何画椭圆:不像一般的绘图指定X和Y，这个时候需要画函数图形。如果你在网上找到它，你应该先建立函数F，然后利用它。
% ezplot(@(x,y)F(pr,[x,y])）可以画出函数图形。
% 数据是随机生成的，每一个程序的运行结果都会有所不同，在B中(:,1)=Data(:,ind(2));有时候数组越界会出错，但是单步调试没有问题，
% 还没有找到原因。
%%6.学习经验
(1)不能过于依赖当前的函数，需要多自己去思考算法，即使借用别人的函数，也要弄清楚原理和调用方法；
(2)matlab函数库不熟悉，要多使用help。；
在编写程序时，首先要对程序结构进行整体规划，模块化，组织清晰，自己要了解自己流程的每一步原理。
(4)理论的力量是无穷无尽的，代码设计应该基于理论的深刻理解。
本文为转载内容，我们尊重原作者对本文的作权。如有内容错误或侵权问题，欢迎原作者联系我们更正或删除内容。
本文仅代表作者观点，版权归原创者所有，如需转载请在文中注明来源及作者名字。
免责声明：本文系转载编辑文章，仅作分享之用。如分享内容、图片侵犯到您的版权或非授权发布，请及时与我们联系进行审核处理或删除，您可以发送材料至邮箱：service@tojoy.com