function Par = PrattSVD(XY)

%--------------------------------------------------------------------------
%  
%     Algebraic circle fit by Pratt
%      V. Pratt, "Direct least-squares fitting of algebraic surfaces",
%      Computer Graphics, Vol. 21, pages 145-152 (1987)
%
%     Input:  XY(n,2) is the array of coordinates of n points x(i)=XY(i,1), y(i)=XY(i,2)
%
%     Output: Par = [a b R] is the fitting circle:
%                           center (a,b) and radius R
%
%     Note: this is a version optimized for stability, not for speed
%
%--------------------------------------------------------------------------

centroid = mean(XY);   % the centroid of the data set

[U,S,V]=svd([(XY(:,1)-centroid(1)).^2+(XY(:,2)-centroid(2)).^2,...
        XY(:,1)-centroid(1), XY(:,2)-centroid(2), ones(size(XY,1),1)],0);

if (S(4,4)/S(1,1) < 1e-12)   %  singular case
    A = V(:,4);
    disp('Pratt singular case');
else                         %  regular case
    W=V*S;
    Binv = [0 0 0 -0.5; 0 1 0 0; 0 0 1 0; -0.5 0 0 0];
    [E,D] = eig(W'*Binv*W);
    [Dsort,ID] = sort(diag(D));
    A = E(:,ID(2));
    for i=1:4
        S(i,i)=1/S(i,i); 
    end
    A = V*S*A;
end

Par = -(A(2:3))'/A(1)/2 + centroid;
Par = [Par , sqrt(A(2)^2+A(3)^2-4*A(1)*A(4))/abs(A(1))/2];

end   %  PrattSVD