function Par = Spath(XY,ParIni) %-------------------------------------------------------------------------- % % Geometric circle fit (minimizing orthogonal distances) by Spath algorithm % H. Spath, "Least-Squares Fitting By Circles", % Computing, Vol. 57, pages 179-185, (1996) % % Input: XY(n,2) is the array of coordinates of n points x(i)=XY(i,1), y(i)=XY(i,2) % ParIni = [a b R] is the initial guess (supplied by user) % % Output: Par = [a b R] is the fitting circle: % center (a,b) and radius R % %-------------------------------------------------------------------------- epsilon=0.00001; % tolerance (small threshold) IterMAX = 500; % maximal number of iterations; with a bad initial guess it may take >100 iterations centroid = mean(XY); % the centroid of the data set X = XY(:,1) - centroid(1); % centering data Y = XY(:,2) - centroid(2); % centering data ParNew = ParIni - [centroid , 0]; % centering the initial guess for iter=1:IterMAX % main loop, each run is one iteration ParOld = ParNew; Dx = X - ParOld(1); Dy = Y - ParOld(2); D = sqrt(Dx.*Dx + Dy.*Dy); Mu = mean(Dx./D); Mv = mean(Dy./D); Mr = mean(D); Radius = (Mu*ParOld(1) + Mv*ParOld(2) + Mr)/(1. - Mu*Mu - Mv*Mv); ParNew = [-Mu , -Mv , 1]*Radius; progress = (norm(ParNew-ParOld))/(norm(ParOld)+epsilon); if progress < epsilon, break, end % stopping rule end % the end of the main loop (over iterations) Par = ParOld + [centroid , 0]; % `decentering' the parameter vector for output end % Spath