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