function Par = KarimakiOrig(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 is the original version, without corrective step % %-------------------------------------------------------------------------- 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); ZXY1 = [Z XY ones(n,1)]; M = ZXY1' * ZXY1; Cxx = det(M(2:2:4,2:2:4)); Cyy = det(M(3:4,3:4)); Cxy = det(M(2:2:4,3:4)); Czx = det(M(1:3:4,2:2:4)); Czy = det(M(1:3:4,3:4)); Czz = det(M(1:3:4,1:3:4)); frac = 2*(Czz*Cxy - Czx*Czy)/(Czz*(Cxx-Cyy)-Czx*Czx+Czy*Czy); phi = atan(frac)/2; if (Czz*Cxy - Czx*Czy)*sin(2*phi)<0 phi = phi+pi/2; end sp = sin(phi); cp = cos(phi); kappa = (Czx*sp - Czy*cp)/Czz; mu = (-M(1,4)*kappa + M(2,4)*sp - M(3,4)*cp)/n; rho = 2*kappa/sqrt(1-4*kappa*mu); d = 2*mu/(1+sqrt(1-4*kappa*mu)); Par = zeros(1,3); Par(1) = (d+1/rho)*sin(phi); Par(2) =-(d+1/rho)*cos(phi); Par(1:2) = Par(1:2) + pole; Par(3) = 1/abs(rho); end % KarimakiOrig