function [Par,Var] = KMvD(XY)

%--------------------------------------------------------------------------
%  
%     Consistent version of simple algebraic circle fit (Kasa method)
%
%     Input:  XY(n,2) is array of coordinates of n points <XY(i,1),XY(i,2)>
%
%     Output: Par = [a b R] is the fitting circle:
%                           center (a,b) and radius R
%
%--------------------------------------------------------------------------

n = size(XY,1);      % number of data points
Z = XY(:,1).*XY(:,1) + XY(:,2).*XY(:,2);

ZXY1 = [Z XY ones(n,1)];
M0 = ZXY1' * ZXY1;
M1 = [8*M0(1,4) 4*M0(2,4) 4*M0(3,4) 2*n; 4*M0(2,4) n 0 0; 4*M0(3,4) 0 n 0; 2*n 0 0 0];
M2 = [8*n 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 0];

Vmin = 0;
XYcent = XY - repmat(mean(XY),n,1);
scatter = XYcent'*XYcent;
[E,D] = eig(XYcent'*XYcent);
Vmax = min(diag(D));
epsilon = 0.00001*Vmax;

while (Vmax - Vmin > epsilon)
       V = (Vmin + Vmax)/2;
       M = M0 - V*M1 + V*V*M2;
       Eval   = eig(M);
       if min(Eval)>0
          Vmin = V;
       else
          Vmax = V;
       end
end

Var = V;

%[Evec,Eval] = eig(M);
%[Evalmin,im] = min(diag(Eval));
%Evecmin = Evec(:,im);
%P = Evecmin(2:4,1)/Evecmin(1,1);

XY1 = [XY ones(n,1)];
K = XY1' * XY1 - V * [n 0 0; 0 n 0; 0 0 0];
N = XY1' * (Z - 4*V) + [0 0 2*V*n]';

P = K \ N;

Par = zeros(1,3);
Par(1) = P(1)/2;
Par(2) = P(2)/2;
Par(3) = sqrt(Par(1)*Par(1)+Par(2)*Par(2)+P(3));