qian1.m

一个拟合椭圆的程序

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

dy=y1;

[m,n]=size(dx);
dx1=dx;
dx2=dx.^2;
dx3=dx.^3;
dx4=dx.^4;

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

% find P in QX=P
P=-[dx2*dy2';dx1*dy3';dx1*dy2';sum(dy3);sum(dy2)];

% find Q in QX=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 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,1)
plot(x,y1,x,y2)

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

%fitted the curve




%x=lyxleastsquare(A,b);
%D=[0,2;0,-2;1,0;-1,0;0.5,1.732;-0.5,1.732;0.5,-1.732;-0.5,-1.732]';
%lyxellipse1(D)