% taubincircle gives the 4 parameter vector of the apprximate Euclidean circle fit
% to a_1 X_i^2 + a_1 Y_i^2 + a_2 X_i + a_3 Y_i + a_4 = 0
% X & Y are lists of point coordinates and must be column vectors.
function a = taubincircle(X,Y)

   % normalize data
   mx = mean(X);
   my = mean(Y);
   sx = (max(X)-min(X))/2;
   sy = (max(Y)-min(Y))/2;
   s = max(sx,sy);
   x = (X-mx)/s;
   y = (Y-my)/s;

   % Build R matrix
   D = [ x.*x + y.*y  x  y  ones(size(x)) ];
   R = D'*D;

   % Build S matrix
   S = [[ 4*(sum(x.*x) + sum(y.*y)) 2*sum(x) 2*sum(y) 0]' ...
        [2*sum(x) 1 0 0]' [2*sum(y) 0 1 0]' [0 0 0 0]'];

   % Solve eigensystem
   [gevec, geval] = eig(R,S)%;

   % Find the smallest eigenvalue
   [MinV, MinC] = min(diag(geval));

   % Extract eigenvector corresponding to smallest eigenvalue
   A = gevec(:,MinC);

   % unnormalize
   a = [        ...
        A(1), ...
        A(2)*s - 2*A(1)*mx, ...
        A(3)*s - 2*A(1)*my, ...
        A(1)*(mx*mx+my*my) - A(2)*s*mx- A(3)*s*my - A(4)*s*s   ...
   ]';

