function [z, r, phi, step] = parcircle (X, z, r, show);%PARCIRCLE      Geometric circle fit%% [z, r, phi, step] = parcircle (X, z, r, show{0});% computes the best fit circle in parameterform % x = z(1) + r  cos(phi), y = z(2) + r sin(phi)%% X: given points <X(i,1), X(i,2)>% z, r: starting values for circle% show: if (show == 1), test output%% z, r: circle found% phi: phi(i) angle to nearest point on circle% step: nof iterations  m = size(X,1);  % compute initial approximations for phi_i  phi  = angle(X(:,1)-z(1)+ i*(X(:,2)-z(2)));  step = 0;  h    = 1;   while (norm(h) > 1e-5),    step  = step+1;    if (step > 100),      disp ('warning: number of iterations exceeded limit');      break;    end  % form Jacobian     S =  diag(sin(phi));    C =  diag(cos(phi));    A = [ -ones(size(phi)) zeros(size(phi)) -cos(phi)];    B = [ zeros(size(phi)) -ones(size(phi)) -sin(phi)];    J = [r*S A; -r*C B];  % QR decomposition of Jacobian     Q = [S C; -C S];    [U R] = qr(C*A+S*B);    RR =  R(1:3, 1:3);    RRR = [r*eye(m) S*A-C*B; zeros(size(A')) RR];  % transform right hand  side    Y = [X(:,1)-z(1)-r*cos(phi); X(:,2)-z(2)-r*sin(phi)];    Y = [eye(m) zeros(m); zeros(m) U]'*Q'*Y;  % solve triangular system    h = -RRR\Y(1:m+3);  % update solution     phi = phi + h(1:m);    z = z +  h(m+1:m+2);    r = r+h(m+3);    if (show == 1),      drawcircle (z,r);      [z' r phi']    end  endend % parcircle