【请教】拟合椭圆参数方程/椭圆拟合

sunyj09

各位大侠，大家好：       我已经提取了椭圆的坐标值，如何应用MATLAB拟合出参数方程啊
       望各位不吝赐教！在线等！！！先谢过了

qibbxxt

这个说来也不难，1.线性回归，得到5个参数；2.根据数学公式得出中心，长短轴和角度
给你一个有用的m程序，应该可以解决你的问题

1. function report=ellipsefit(XY)
   
2. %ELLIPSEFIT - form 2D ellipse fit to given x,y data
   
3. %
   
4. % report=ellipsefit(XY)
   
5. %
   
6. %in:
   
7. %
   
8. % XY: Input matrix of 2D coordinates to be fit. Each column XY(:,i) is [xi;yi]
   
9. %
   
10. %out: Finds the ellipse fitting the input data parametrized both as
    
11. % A*x^2+B*x*y C*y^2+D*x+E*y=1 and [x-x0,y-y0]*Q*[x-x0;y-y0]=1
    
12. %
    
13. % report: a structure output with the following fields
    
14. %
    
15. % report.Q: the matrix Q
    
16. % report.d: the vector [x0,y0]
    
17. % report.ABCDE: the vector [A,B,C,D,E]
    
18. % report.AxesDiams: The minor and major ellipse diameters
    
19. % report.theta: The counter-clockwise rotation of the ellipse.
    
20. %
    
21. %NOTE: The code will give errors if the data fit traces out a non-elliptic or
    
22. % degenerate conic section.
    
23. 
    
24. 
    
25. X=XY(1,:).';
    
26. Y=XY(2,:).';
    
27. 
    
28. M= [X.^2, X.*Y, Y.^2, X, Y, -ones(size(X,1),1)];
    
29. 
    
30. ABCDEF=null(M);
    
31. 
    
32. if size(ABCDEF,2)>1
    
33.    
    
34.     error 'Data cannot be fit with unique ellipse'
    
35. else
    
36. 
    
37.     ABCDEF=num2cell(ABCDEF);
    
38. end
    
39. 
    
40. [A,B,C,D,E,F]=deal(ABCDEF{:});
    
41. 
    
42. 
    
43. Q=[A, B/2;B/2 C];
    
44. x0=-Q\[D;E]/2;
    
45. dd=F+x0'*Q*x0;
    
46. 
    
47. Q=Q/dd;
    
48. 
    
49. [R,eigen]=eig(Q);
    
50. eigen=eigen([1,4]);
    
51. 
    
52. if ~all(eigen>=0), error 'Fit produced a non-elliptic conic section'; end
    
53. 
    
54. idx=eigen>0;
    
55. eigen(idx)=1./eigen(idx);
    
56. AxesDiams = 2*sqrt(eigen);
    
57. 
    
58. theta=atand(tand(-atan2(R(1),R(2))*180/pi));
    
59. 
    
60. 
    
61. report.Q=Q;
    
62. report.d=x0(:).';
    
63. report.ABCDE=[A, B, C, D, E]/F;
    
64. report.AxesDiams=sort(AxesDiams(:)).';
    
65. report.theta=theta;

复制代码

sunyj09

嗯 谢谢了！我试试

hankezhou

解决问题的人。。很专业啊。。从今天开始~天天关注各版最新进展~

千里马1898

这位大侠程序是个好东西 可是怎么运行呢？我是个新手现在记者用这个程序拟合散点，但是我在matlab对话框中输入XY[5 6 4 0 -3 -4 0;0 3 7 5 0 -4 -3]后 提示出错啊！怎么才能像上面一样得出拟合的椭圆来呢？ 2# qibbxxt

千里马1898

具体怎么运行还请指点

千里马1898

运行出现如下提示 ：还请达人指点
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,:).';

bainhome

用点坐标数据和“1”的移项构造标准形式椭圆系数矩阵，用null命令做正交基，思想还是最小二乘，但程序代码简单的多，不错。
构思很漂亮！作者对矩阵思想的理解很深啊，自叹不如。

千里马1898

但是运行起来好像有点问题啊！楼上的高手能不能指点下哦！

bainhome

1. p1=[92.8 228.3];p2=[251.5 233.8];p3=[129.2 157.8];p4=[220.5 160.5];p5=[180.6 252.2];XY=[p1;p2;p3;p4;p5]';

复制代码

1. 
   
2. report=ellipsefit2(XY)
   
3. 
   
4. report =
   
5.             Q: [2x2 double]
   
6.             d: [173.5491 202.7305]
   
7.         ABCDE: [-5.8005e-006 9.6955e-007 -2.2233e-005 0.0018 0.0088]
   
8.     AxesDiams: [98.8620 193.8521]
   
9.         theta: -88.3117

复制代码

你的问题？还是程序运行的问题？
这实在是个问题。

千里马1898

是我对这个程序理解上有点误区，我以为这个是可以把很多点拟合成一个近似椭圆的，原来它值支持5个点的拟合，但是上面的说明语句中貌似有吧很多点拟合成一个近似椭圆的意思！
看来我还得再找找能把很多点拟合成椭圆的程序，多谢楼上的斑竹指点！

bainhome

支持至少5个点的最小二乘椭圆拟合程序：

1. function [varargout]=ellipsefit(x,y)
   
2. %ELLIPSEFIT Stable Direct Least Squares Ellipse Fit to Data.
   
3. % [Xc,Yc,A,B,Phi,P]=ELLIPSEFIT(X,Y) finds the least squares ellipse that
   
4. % best fits the data in X and Y. X and Y must have at least 5 data points.
   
5. % Xc and Yc are the x- and y-axis center of the ellipse respectively.
   
6. % A and B are the major and minor axis of the ellipse respectively.
   
7. % Phi is the radian angle of the major axis with respect to the x-axis.
   
8. % P is a vector containing the general conic parameters of the ellipse.
   
9. % The conic representation of the ellipse is given by:
   
10. %
    
11. % P(1)*x^2 + P(2)*x*y + P(3)*y^2 + P(4)*x + P(5)*y + P(6) = 0
    
12. %
    
13. % S=ELLIPSEFIT(X,Y) returns the output data in a structure with field names
    
14. % equal to the variable names given above, e.g., S.Xc, S.Yc, S.A, S.B,
    
15. % S.Phi and S.P
    
16. %
    
17. % Reference: R. Halif and J. Flusser, "Numerically Stable Direct Least
    
18. % Squares FItting of Ellipses," Department of Software Engineering, Charles
    
19. % University, Czech Republic, 2000.
    
20. 
    
21. % Conversion from conic to conventional ellipse equation inspired by
    
22. % fit_ellipse.m on MATLAB Central
    
23. 
    
24. % D.C. Hanselman, University of Maine, Orono, ME 04469
    
25. % Mastering MATLAB 7
    
26. % 2005-02-28
    
27. % Rotation angle fixed 2005-08-09
    
28. 
    
29. %--------------------------------------------------------------------------
    
30. x=x(:); % convert data to column vectors
    
31. y=y(:);
    
32. if numel(x)~=numel(y) || numel(x)<5
    
33.    error('X and Y Must be the Same Length and Contain at Least 5 Values.')
    
34. end
    
35. 
    
36. D1=[x.*x x.*y y.*y]; % quadratic terms
    
37. D2=[x y ones(size(x))]; % linear terms
    
38. S1=D1'*D1;
    
39. S2=D1'*D2;
    
40. 
    
41. [Q2,R2]=qr(D2,0);
    
42. if condest(R2)>1.0e10
    
43.    warning('ellipsefit',...
    
44.       'Data is Poorly Conditioned and May Not Represent an Ellipse.')
    
45. end
    
46. T=-R2\(R2'\S2'); % -inv(S3) * S2'
    
47. 
    
48. M=S1+S2*T;
    
49. CinvM=[M(3,:)/2; -M(2,:); M(1,:)/2];
    
50. [V,na]=eig(CinvM);
    
51. c=4*V(1,:).*V(3,:) - V(2,:).^2;
    
52. A1=V(:,c>0);
    
53. P=[A1; T*A1];
    
54. 
    
55. % correct signs if needed
    
56. P=sign(P(1))*P;
    
57. 
    
58. Phi=atan(P(2)/(P(3)-P(1)))/2;
    
59. c=cos(Phi);
    
60. s=sin(Phi);
    
61. 
    
62. % rotate the ellipse parallel to x-axis
    
63. Pr=zeros(6,1);
    
64. Pr(1)=P(1)*c*c - P(2)*c*s + P(3)*s*s;
    
65. Pr(2)=2*(P(1)-P(3))*c*s + (c^2-s^2)*P(2);
    
66. Pr(3)=P(1)*s*s + P(2)*s*c + P(3)*c*c;
    
67. Pr(4)=P(4)*c - P(5)*s;
    
68. Pr(5)=P(4)*s + P(5)*c;
    
69. Pr(6)=P(6);
    
70. 
    
71. % extract other data
    
72. XcYc=[c s;-s c]*[-Pr(4)/(2*Pr(1));-Pr(5)/(2*Pr(3))];
    
73. Xc=XcYc(1);
    
74. Yc=XcYc(2);
    
75. F=-Pr(6) + Pr(4)^2/(4*Pr(1)) + Pr(5)^2/(4*Pr(3));
    
76. AB=sqrt(F./Pr(1:2:3));
    
77. A=AB(1);
    
78. B=AB(2);
    
79. Phi=-Phi;
    
80. if A<B % x-axis not major axis, so rotate it pi/2
    
81.    Phi=Phi-sign(Phi)*pi/2;
    
82.    A=AB(2);
    
83.    B=AB(1);
    
84. end
    
85. S.Xc=Xc;
    
86. S.Yc=Yc;
    
87. S.A=A;
    
88. S.B=B;
    
89. S.Phi=Phi;
    
90. S.P=P;
    
91. if nargout==1
    
92.    varargout{1}=S;
    
93. else
    
94.    outcell=struct2cell(S);
    
95.    varargout=outcell(1:nargout);
    
96. end

复制代码

qibbxxt

11# 千里马1898
我在前面说的很清楚，先是线性回归得到系数，然后用数学公式得到参数，你只需要将上面的程序中的

1. null

复制代码

部分用

1. regress

复制代码

代替，就可以了，希望你真正的弄明白程序的含义，学习，改进，达到自己想要的效果，不要想着在网上直接找到源程序，这样对你学习没有什么帮助的

stats01

千里马1898 发表于 2010-9-16 20:59
运行出现如下提示 ：还请达人指点
XY =

clear,clc
x=[5     6     4     0    -3    -4     0]';
y=[0     3     7     5     0    -4    -3]';
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=[0.879099539880687,-0.940582944470011,-0.592917179136204,-1.36681691508122,-13.9710207325040];
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')