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