lyxellipse.m

一个拟合椭圆的程序

clear
N=300;
x=-1:2/400:1;
y1=sqrt(4*(1-x.*x));
y2=-sqrt(4*(1-x.*x));

dx=x;
dy=y1;

[m,n]=size(dx);

dx1=dx;%dx1^1
dx2=dx.^2;%dx^2
dx3=dx.^3;%dx3
dx4=dx.^4;%dx4

dy1=dy;%dy^1
dy2=dy.^2;
dy3=dy.^3;
dy4=dy.^4;%dy^4

% find P in Qcoe=P,coe means the coefficient of the quadratic curve
P=-[dx2*dy2';dx1*dy3';dx1*dy2';sum(dy3);sum(dy2)];

% find Q in Qcoe=P
Q(1,1)=sum(dx4);
Q(1,2)=dx3*dy1';
Q(1,3)=sum(dx3);
Q(1,4)=dx2*dy1';
Q(1,5)=sum(dx2);

Q(2,1)=dx3*dy1';;
Q(2,2)=dx2*dy2';
Q(2,3)=dx2*dy1';
Q(2,4)=dx1*dy2';
Q(2,5)=dx1*dy1';

Q(3,1)=sum(dx3);
Q(3,2)=dx2*dy1';
Q(3,3)=sum(dx2);
Q(3,4)=dx1*dy1';
Q(3,5)=sum(dx1);

Q(4,1)=dx2*dy1';
Q(4,2)=dx1*dy2';
Q(4,3)=dx1*dy1';
Q(4,4)=sum(dy2);
Q(4,5)=sum(dy1);

Q(5,1)=sum(dx2);
Q(5,2)=dx1*dy1';
Q(5,3)=sum(dx1);
Q(5,4)=sum(dy1);
Q(5,5)=n;

Deta=det(Q);

C=inv(Q)*P;
Cxx=C(1)
Cyy=1;
Cxy=C(2)
Cx=C(3)
Cy=C(4)
C0=C(5)

%curve the data given
subplot(2,1,1)
plot(x,y1,x,y2)

%curve equation Cxx*X^2+Y^2+Cxy*XY+Cx*X+Cy*Y+C0=0

xx=-1:2/N:1;
yy1=-(Cxy*xx+Cy);
yy1=yy1+sqrt((Cxy*xx+Cy).*(Cxy*xx+Cy)-4*(Cxx*xx.*xx+Cx*xx+C0));
yy1=yy1/2;

yy2=-(Cxy*xx+Cy);
yy2=yy2-sqrt((Cxy*xx+Cy).*(Cxy*xx+Cy)-4*(Cxx*xx.*xx+Cx*xx+C0));
yy2=yy2/2;

subplot(2,1,2)
plot(xx,yy1,xx,yy2)