pause(2); % wait 2 seconds
[xi,yi] = xinput([-2 2 -2 2]);
axis('equal');

%% fit general ellipse
F = [xi.*yi yi.^2  xi yi ones(size(xi))];
p = LinearRegression(F,-xi.^2);

m  = [2 p(1);p(1) 2*p(2)];
x0 = -m\[p(3);p(4)]
a11= 1/(x0(1)^2+p(1)*x0(1)*x0(2)+p(2)*x0(2)^2-p(5));
[V,la] = eig(a11*m/2);

alpha = atan2(V(2,1)/V(1,1),1)*180/pi
a = 1/sqrt(la(1,1)) 
b = 1/sqrt(la(2,2))

np = 37;
phi = linspace(0,2*pi,np);
xx = V*([a*cos(phi); b*sin(phi)])+x0*ones(1,np);

xg = xx(1,:);yg = xx(2,:);

%% constraint regression
FC = [xi.^2  yi.^2  xi.*yi xi yi ones(size(xi))];
p = RegressionConstraint(FC,1);

xy0 = -[2*p(1), p(3);p(3), 2*p(2)]\[p(4);p(5)];
x0 = xy0(1); y0 = xy0(2);

lambda = p(1)*x0^2 + p(2)*y0^2 + p(3)*x0*y0-1;
p = p/lambda;

[V,la] = eig([p(1), p(3)/2;p(3)/2,p(2)]);
alpha  = atan2(V(2,1)/V(1,1),1)*180/pi
a = 1/sqrt(la(1,1)) 
b = 1/sqrt(la(2,2))

np = 37; phi = linspace(0,2*pi,np);

xx = V*([a*cos(phi); b*sin(phi)]) + xy0*ones(1,np);
xC  = xx(1,:); yC = xx(2,:);

plot(xi,yi,'*r',xg,yg,'b',xC,yC,'g');
