function Par = KMvH(XY)

%--------------------------------------------------------------------------
%  
%     Consistent circle fit (Kukush-Markovsky-van Huffel)
%      A. Kukush, I. Markovsky, S. Van Huffel,
%      "Consistent estimation in an implicit quadratic measurement error model",
%       Comput. Statist. Data Anal., Vol. 47, pp. 123-147, (2004)
%
%     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*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);

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));

end   %  KMvH