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