matlab学习：最小二乘拟合&基于RANSAC的直线拟合&椭圆拟合

1.最小二乘拟合

最小二乘拟合是一种数学上的近似和优化，利用已知的数据得出一条直线或者曲线，使之在坐标系上与已知数据之间的距离的平方和最小。

2.RANSAC算法

参见王荣先老师的博文 /article/4959400.html

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        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('样本点','抽取点','符合点','拟合曲线')
%% 问题分析
%    1.关于如何生成随机点：在标准椭圆基础上，添加高斯白噪声--wgn();
%    2.关于如何建立椭圆模型：
%      方案一：椭圆一般方程：Ax2+Bxy+Cy2+Dx+Ey+F=0，认为由5个点确定一个椭圆，则用5个点代入方程式去做椭圆拟合，结果大多数
%      情况下画出双曲线，放弃；
%      方案二：利用椭圆定义：到两定点之间距离和为常数2a，选择平面内3个点，两个焦点，一个过椭圆的点，确定椭圆。
%    3.关于如何筛选符合条件的点：此时计点到椭圆距离过于复杂，用定义，到两焦点距离和与2a相差小于一定阈值，则符合。
%    4.关于拟合函数：使用nlinfit，对于输入参数的维数有要求，需要x为N*P维，y为n*1维，注意是列向量。
%    5.关于如何画椭圆：与一般的画图指定X和Y不同，此时要画的是函数图形，在网上查到，先建立函数F，再利用
%      ezplot(@(x,y)F(pr,[x,y])）可画出函数图形。
%    6.数据为随机产生，程序每次运行结果会不一样，在B(:,1)=Data(:,ind(2));时有时会数组越界出错，但单步调试时没问题，
%      原因还未找到。
%%

6.学习经验
（1）不能太依赖于现有函数，要多自己想算法，即便借用别人的函数，也要弄清楚原理及调用方式；

（2）matlab函数库不熟悉，要多用help；

（3）编写程序时要先总体规划好程序架构，模块化，条理清楚，自己要懂得自己程序的每一步原理。

（4）理论的力量是无穷的，要在理论深入理解的基础上进行代码优化。