最小二乘椭圆拟合ovalfit m sci 代码

 //m sci 代码

function [xc,yc,Ra,Rb,theta,res] = ovalfit(x,y)

//椭圆方程：x^2+res(1) xy+res(2) y^2+res(3) x+res(4) y+res(5)=0;

//长半轴Ra，短半轴Rb，中心(xc1,yc1); 长轴与x轴夹角theta.

N=length(x);
g=zeros(5,5);
b=zeros(5);
g(1,1)=sum(x.*x.*y.*y);
g(1,2)=sum(x.*y.*y.*y);
g(1,3)=sum(x.*x.*y);
g(1,4)=sum(x.*y.*y);
g(1,5)=sum(x.*y);

g(2,1)=sum(x.*y.*y.*y);
g(2,2)=sum(y.*y.*y.*y);
g(2,3)=sum(x.*y.*y);
g(2,4)=sum(y.*y.*y);
g(2,5)=sum(y.*y);

g(3,1)=sum(x.*x.*y);
g(3,2)=sum(x.*y.*y);
g(3,3)=sum(x.*x);
g(3,4)=sum(x.*y);
g(3,5)=sum(x);

g(4,1)=sum(x.*y.*y);
g(4,2)=sum(y.*y.*y);
g(4.3)=sum(x.*y);
g(4,4)=sum(y.*y);
g(4,5)=sum(y);

g(5,1)=sum(x.*y);
g(5,2)=sum(y.*y);
g(5,3)=sum(x);
g(5,4)=sum(y);
g(5,5)=N;

b(1)=-sum(x.*x.*x.*y);
b(2)=-sum(x.*x.*y.*y);
b(3)=-sum(x.*x.*x);
b(4)=-sum(x.*x.*y);
b(5)=-sum(x.*x);

res=zeros(5);
res=b'/g;
A=res(1);
B=res(2);
C=res(3);
D=res(4);
E=res(5);

Ra=sqrt(2*(A*C*D-B*C*C-D*D+4*B*E)/(A*A-4*B)/(B-sqrt(A*A+(1-B)*(1-B))+1));
Rb=sqrt(2*(A*C*D-B*C*C-D*D+4*B*E)/(A*A-4*B)/(B+sqrt(A*A+(1-B)*(1-B))+1));
theta=atan(sqrt((ra*ra-rb*rb*B)/(ra*ra*B-rb*rb)));
xc=(2*B*C-A*D)/(A*A-4*B);
yc=(2*B*D-A*D)/(A*A-4*B);

endfunction

 

调用示例：

N=6;
t=0.001:1/N:1.0;
PI=3.14159265358979;
x=0.2*sin(t*2*PI);
y=1.7*cos(t*2*PI);

 

[xc1,yc1,Ra,Rb,theta,res]=ovalfit(x,y)

 

椭圆方程：x^2+res(1) xy+res(2) y^2+res(3) x+res(4) y+res(5)=0;

长半轴Ra，短半轴Rb，中心(xc1,yc1); 长轴与x轴夹角theta.