parcircle.m

椭圆拟合的相关介绍与数学运算方法

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
  end

end % parcircle