function Par = KarimakiAlg(XY,pole) %-------------------------------------------------------------------------- % % Karimaki circle fit % V. Karimaki, % "Effective circle fitting for particle trajectories", % Nucl. Instr. Meth. Phys. Res. A., Vol. 305, pp. 187--191, (1991) % % Input: XY(n,2) is array of coordinates of n points % % Output: Par = [a b R] is the fitting circle: % center (a,b) and radius R % % This version is modified from the original one; it uses algebraic parameters % No corective steps % %-------------------------------------------------------------------------- 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]; means = mean(ZXY); ZXY0 = ZXY - repmat(means,n,1); M = ZXY0' * ZXY0; 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 Karimaki: the smallest e-value is negative...') end BC = E(:,ID(1))'; A = -BC*M(2:3,1)/M(1,1); Dp = -A*means(1) - BC*means(2:3)'; Par = zeros(1,3); Par(1:2) = -BC/2/A; Par(3) = sqrt(norm(BC)^2-4*A*Dp)/2/abs(A); Par(1:2) = Par(1:2) + pole; end % KarimakiAlg