function Par = RTKD(XY,pole)

%--------------------------------------------------------------------------
%  
%     RTKD (Rusu-Tico-Kuosmanen-Delp) circle fit
%   C. Rusu, M. Tico, P. Kuosmanen, and E. J. Delp,
%   "Classical geometrical approach to circle fitting—review and new developments",
%    J. Electron. Imaging, Vol. 12, pp. 179-193, (2003)
%
%     Input:  XY(n,2) is array of coordinates of n points <XY(i,1),XY(i,2)>
%
%     Output: Par = [a b R] is the fitting circle:
%                           center (a,b) and radius R
%
%     This is a non-iterative version
%
%--------------------------------------------------------------------------

n = size(XY,1);      % number of data points
XY = XY - repmat(pole,n,1);   % transfer coordinates to the pole
Z = XY(:,1).*XY(:,1) + XY(:,2).*XY(:,2);

ZXY = [Z XY];
M = ZXY' * ZXY;
C = [det(M(1:2,1:2)) det(M(1:2,1:2:3)); det(M(1:2,1:2:3)) det(M(1:2:3,1:2:3))];
[E,D] = eig(C);
[Dsort,ID] = sort(diag(D));
if (Dsort(1)<0) 
    disp('Error in RTKD: the smallest e-value is negative...')
end
BC = E(:,ID(1))';
A = -BC*M(2:3,1)/M(1,1);

Par = zeros(1,3);
Par(1:2) = -BC/2/A;
Par(3) = norm(BC)/2/abs(A);
Par(1:2) = Par(1:2) + pole;

end   %  RTKD