椭圆拟合算法总结

版权声明：本文为博主原创文章，未经博主允许不得转载。

椭圆拟合算法

椭圆拟合在医疗图像、工程设计、交通识别、生物识别、人工智能、林火蔓延等领域具有广泛应用，本文总结前人研究综合出基于MATLAB的三种椭圆拟合的算法：直接计算方法，二次多项式拟合方法，标准方程拟合方法。

 

一，二次多项式拟合方法

利用圆锥曲线二次函数拟合，求算五个参数。关键在于估计的参数单位. p0=[1 1 1 1 1 1]为系数。

 

[plain]  view plain  copy

1. % 拟合椭圆型曲线段  
2. close all;  
3. clear; clc; close all;  
4. % 设出圆锥曲线方程  
5. 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);%椭圆一般方程  
6. % 离散数据点  
7. Up=xlsread('E:\\Up.xls');%excel文件路径  
8. UpX=Up(:,1);  
9. UpY=Up(:,2);  
10.    
11. % p0系数初值  
12. p0=[1 1 1 1 1 1];  
13. warning off  
14. % 拟合系数，最小二乘方法  
15. p =nlinfit(Up,zeros(size(Up,1),1),F,p0);  
16. plot(UpX,UpY,'r.');  
17. hold on;  
18.    
19. UpMinx=min(UpX);  
20. UpMaxx=max(UpX);  
21. UpMiny=min(UpY);  
22. UpMaxy=max(UpY);  
23.    
24. % 作图  
25. ezplot(@(x,y)F(p,[x,y]),[-1+UpMinx,1+UpMaxx,-1+UpMiny,1+UpMaxy]);  
26. title('曲线拟合');  
27. legend('样本点','拟合曲线');  

 

二，直接计算方法

 

[plain]  view plain  copy

1. S=xlsread('E:\\Up.xls');  
2. r=S(:,1);  
3. 　　t=S(:,2);  
4. 　　x=r.*cos((90-t)*pi/180);%转换成直角坐标  
5. 　    y=r.*sin((90-t)*pi/180);%转换成直角坐标  
6.    
7. XY=[x,y];  
8. A = EllipseDirectFit(XY);  
9. % A = [a b c d e f]  
10. %  ax^2 + bxy + cy^2 + dx + ey + f = 0  
11.    
12.    
13. function A = EllipseDirectFit(XY);  
14. centroid = mean(XY);   % thecentroid of the data set  
15.    
16. D1 = [(XY(:,1)-centroid(1)).^2,(XY(:,1)-centroid(1)).*(XY(:,2)-centroid(2)),...  
17.       (XY(:,2)-centroid(2)).^2];  
18. D2 = [XY(:,1)-centroid(1),XY(:,2)-centroid(2), ones(size(XY,1),1)];  
19. S1 = D1'*D1;  
20. S2 = D1'*D2;  
21. S3 = D2'*D2;  
22. T = -inv(S3)*S2';  
23. M = S1 + S2*T;  
24. M = [M(3,:)./2; -M(2,:);M(1,:)./2];  
25. [evec,eval] = eig(M);  
26. cond =4*evec(1,:).*evec(3,:)-evec(2,:).^2;  
27. A1 = evec(:,find(cond>0));  
28. A = [A1; T*A1];  
29. A4 =A(4)-2*A(1)*centroid(1)-A(2)*centroid(2);  
30. A5 =A(5)-2*A(3)*centroid(2)-A(2)*centroid(1);  
31. A6 =A(6)+A(1)*centroid(1)^2+A(3)*centroid(2)^2+...  
32.      A(2)*centroid(1)*centroid(2)-A(4)*centroid(1)-A(5)*centroid(2);  
33. A(4) = A4;  A(5) = A5; A(6) = A6;  
34. A = A/norm(A);  
35.    
36. end  %  EllipseDirectFit  

 

三、标准方程拟合方法

设椭圆拟合方程：f(x,y)=(x-h)²/a²+(y-k)²/b²-1，即椭圆中心为（h,k），长半轴为a，短半轴为b

[plain]  view plain  copy

1. X=[xi yi]  
2. n=length(X(:,1));  
3. y=ones(n,1);  
4. b=[rand(1) rand(1)  rand(1)];  
5. fun=inline('(X(:,1)-a(1)).^2./a(3).^2+(X(:,2)-a(2).^2./a(4).^2-1','a','X');  
6. [a,r,J] = nlinfit(X,y,fun,b);  
7. vpa(a,10)  
8. h=a(1),k=a(2)  
9. b=a(3),a=a(4)