function [z,a,b,alpha]=AlgebraicEllipseFit(X);
% ALGEBRAICELLIPSEFIT ellipse fit, minimizing the algebraic distance
%   [z,a,b,alpha]=AlgebraicEllipseFit(X) fits an ellipse by minimizing
%   the algebraic distance to given points P_i=[X(i,1), X(i,2)] in the
%   least squares sense x'A x + bb'x + c=0. z is the center, a,b are
%   the main axes and alpha the angle between a and the x-axis.

[U S V]=svd([X(:,1).^2  X(:,1).*X(:,2) X(:,2).^2 ...
	     X(:,1) X(:,2) ones(size(X(:,1)))]);
u=V(:,6); A=[u(1) u(2)/2; u(2)/2 u(3)]; 
bb=[u(4); u(5)]; c=u(6);
[Q,D]=eig(A);
alpha=atan2(Q(2,1),Q(1,1));
bs=Q'*bb; zs=-(2*D)\bs; z=Q*zs;
h=-bs'*zs/2-c; a=sqrt(h/D(1,1)); b=sqrt(h/D(2,2));