【请教】拟合椭圆参数方程/椭圆拟合
本帖最后由 messenger 于 2011-1-19 17:05 编辑各位大侠，大家好： 我已经提取了椭圆的坐标值，如何应用MATLAB拟合出参数方程啊
望各位不吝赐教！在线等！！！先谢过了 本帖最后由 qibbxxt 于 2010-9-9 14:26 编辑
这个说来也不难，1.线性回归，得到5个参数；2.根据数学公式得出中心，长短轴和角度
给你一个有用的m程序，应该可以解决你的问题
function report=ellipsefit(XY)
%ELLIPSEFIT - form 2D ellipse fit to given x,y data
%
% report=ellipsefit(XY)
%
%in:
%
% XY: Input matrix of 2D coordinates to be fit. Each column XY(:,i) is
%
%out: Finds the ellipse fitting the input data parametrized both as
% A*x^2+B*x*y C*y^2+D*x+E*y=1 and *Q*=1
%
% report: a structure output with the following fields
%
% report.Q: the matrix Q
% report.d: the vector
% report.ABCDE: the vector
% report.AxesDiams: The minor and major ellipse diameters
% report.theta: The counter-clockwise rotation of the ellipse.
%
%NOTE: The code will give errors if the data fit traces out a non-elliptic or
% degenerate conic section.
X=XY(1,:).';
Y=XY(2,:).';
M= ;
ABCDEF=null(M);
if size(ABCDEF,2)>1
error 'Data cannot be fit with unique ellipse'
else
ABCDEF=num2cell(ABCDEF);
end
=deal(ABCDEF{:});
Q=;
x0=-Q\/2;
dd=F+x0'*Q*x0;
Q=Q/dd;
=eig(Q);
eigen=eigen();
if ~all(eigen>=0), error 'Fit produced a non-elliptic conic section'; end
idx=eigen>0;
eigen(idx)=1./eigen(idx);
AxesDiams = 2*sqrt(eigen);
theta=atand(tand(-atan2(R(1),R(2))*180/pi));
report.Q=Q;
report.d=x0(:).';
report.ABCDE=/F;
report.AxesDiams=sort(AxesDiams(:)).';
report.theta=theta; 嗯 谢谢了！我试试 解决问题的人。。很专业啊。。从今天开始~天天关注各版最新进展~ 这位大侠程序是个好东西 可是怎么运行呢？我是个新手现在记者用这个程序拟合散点，但是我在matlab对话框中输入XY后 提示出错啊！怎么才能像上面一样得出拟合的椭圆来呢？ 2# qibbxxt 具体怎么运行还请指点 运行出现如下提示 ：还请达人指点
XY =
5 6 4 0 -3 -4 0
0 3 7 5 0 -4 -3
??? Input argument "XY" is undefined.
Error in ==> ellipsefit at 24
X=XY(1,:).'; 用点坐标数据和“1”的移项构造标准形式椭圆系数矩阵，用null命令做正交基，思想还是最小二乘，但程序代码简单的多，不错。
构思很漂亮！作者对矩阵思想的理解很深啊，自叹不如。 但是运行起来好像有点问题啊！楼上的高手能不能指点下哦！ p1=;p2=;p3=;p4=;p5=;XY=';
report=ellipsefit2(XY)
report =
Q:
d:
ABCDE: [-5.8005e-006 9.6955e-007 -2.2233e-005 0.0018 0.0088]
AxesDiams:
theta: -88.3117你的问题？还是程序运行的问题？
这实在是个问题。 是我对这个程序理解上有点误区，我以为这个是可以把很多点拟合成一个近似椭圆的，原来它值支持5个点的拟合，但是上面的说明语句中貌似有吧很多点拟合成一个近似椭圆的意思！
看来我还得再找找能把很多点拟合成椭圆的程序，多谢楼上的斑竹指点！ 支持至少5个点的最小二乘椭圆拟合程序：function =ellipsefit(x,y)
%ELLIPSEFIT Stable Direct Least Squares Ellipse Fit to Data.
% =ELLIPSEFIT(X,Y) finds the least squares ellipse that
% best fits the data in X and Y. X and Y must have at least 5 data points.
% Xc and Yc are the x- and y-axis center of the ellipse respectively.
% A and B are the major and minor axis of the ellipse respectively.
% Phi is the radian angle of the major axis with respect to the x-axis.
% P is a vector containing the general conic parameters of the ellipse.
% The conic representation of the ellipse is given by:
%
% P(1)*x^2 + P(2)*x*y + P(3)*y^2 + P(4)*x + P(5)*y + P(6) = 0
%
% S=ELLIPSEFIT(X,Y) returns the output data in a structure with field names
% equal to the variable names given above, e.g., S.Xc, S.Yc, S.A, S.B,
% S.Phi and S.P
%
% Reference: R. Halif and J. Flusser, "Numerically Stable Direct Least
% Squares FItting of Ellipses," Department of Software Engineering, Charles
% University, Czech Republic, 2000.
% Conversion from conic to conventional ellipse equation inspired by
% fit_ellipse.m on MATLAB Central
% D.C. Hanselman, University of Maine, Orono, ME 04469
% Mastering MATLAB 7
% 2005-02-28
% Rotation angle fixed 2005-08-09
%--------------------------------------------------------------------------
x=x(:); % convert data to column vectors
y=y(:);
if numel(x)~=numel(y) || numel(x)<5
error('X and Y Must be the Same Length and Contain at Least 5 Values.')
end
D1=; % quadratic terms
D2=; % linear terms
S1=D1'*D1;
S2=D1'*D2;
=qr(D2,0);
if condest(R2)>1.0e10
warning('ellipsefit',...
'Data is Poorly Conditioned and May Not Represent an Ellipse.')
end
T=-R2\(R2'\S2'); % -inv(S3) * S2'
M=S1+S2*T;
CinvM=;
=eig(CinvM);
c=4*V(1,:).*V(3,:) - V(2,:).^2;
A1=V(:,c>0);
P=;
% correct signs if needed
P=sign(P(1))*P;
Phi=atan(P(2)/(P(3)-P(1)))/2;
c=cos(Phi);
s=sin(Phi);
% rotate the ellipse parallel to x-axis
Pr=zeros(6,1);
Pr(1)=P(1)*c*c - P(2)*c*s + P(3)*s*s;
Pr(2)=2*(P(1)-P(3))*c*s + (c^2-s^2)*P(2);
Pr(3)=P(1)*s*s + P(2)*s*c + P(3)*c*c;
Pr(4)=P(4)*c - P(5)*s;
Pr(5)=P(4)*s + P(5)*c;
Pr(6)=P(6);
% extract other data
XcYc=*[-Pr(4)/(2*Pr(1));-Pr(5)/(2*Pr(3))];
Xc=XcYc(1);
Yc=XcYc(2);
F=-Pr(6) + Pr(4)^2/(4*Pr(1)) + Pr(5)^2/(4*Pr(3));
AB=sqrt(F./Pr(1:2:3));
A=AB(1);
B=AB(2);
Phi=-Phi;
if A<B % x-axis not major axis, so rotate it pi/2
Phi=Phi-sign(Phi)*pi/2;
A=AB(2);
B=AB(1);
end
S.Xc=Xc;
S.Yc=Yc;
S.A=A;
S.B=B;
S.Phi=Phi;
S.P=P;
if nargout==1
varargout{1}=S;
else
outcell=struct2cell(S);
varargout=outcell(1:nargout);
end 11# 千里马1898
我在前面说的很清楚，先是线性回归得到系数，然后用数学公式得到参数，你只需要将上面的程序中的null部分用regress代替，就可以了，希望你真正的弄明白程序的含义，学习，改进，达到自己想要的效果，不要想着在网上直接找到源程序，这样对你学习没有什么帮助的 千里马1898 发表于 2010-9-16 20:59 static/image/common/back.gif
运行出现如下提示 ：还请达人指点
XY =
clear,clc
x=';
y=';
fx1=@(b,x)-(b(3) + b(4)*x + (b(3)^2 + 2*b(3)*b(4)*x + b(4)^2*x.^2 - 4*b(1)*x.^2 - 4*b(1)*b(2)*x - 4*b(1)*b(5)).^(1/2))/(2*b(1));
fx2=@(b,x)-(b(3) + b(4)*x - (b(3)^2 + 2*b(3)*b(4)*x + b(4)^2*x.^2 - 4*b(1)*x.^2 - 4*b(1)*b(2)*x - 4*b(1)*b(5)).^(1/2))/(2*b(1));
b=;
n=length(y);
figure(1),clf
plot(x,y,'o','markerfacecolor','k','markersize',10)
hold on
rg=range(x);
x1=min(x)-rg/40:rg/250:max(x)+rg/8;
y1=fx1(b,x1);
y2=fx2(b,x1);
plot(x1,y1,'r-',x1,y2,'r-','linewidth',3)
legend('data','fit','location','best')
页:
[1]