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 % % 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