lyle.m

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

function [z, a, b, alpha, step] = lyle (X, show);
%LYLE   Iterative algebraic ellipse fit
%
% [z, a, b, alpha, step] = lyle (X, show{0});
% fit ellipse with algebraic method using curvature weights.
%
% X: points given <X(i,1), X(i,2)>
% show: if (show == 1), test output
%
% z, a, b, alpha: ellipse found
% step: nof iterations

  delta = 1;
  omega = 2.0;
  myeps = 1e-3;

  new = [0;0;0;0;0];
  old = new;
  W = ones(size(X,1), 1);
  step = 0;
  m = size(X, 1);

  while (delta > myeps),
    step  = step + 1;
    if (step > 20),
      disp ('warning: number of steps exceeded limit');
      break;
    end
    myeps = omega*myeps;
    [z, a, b, alpha, err] = algellipse (X, W, show);
    if (err),
      disp ('warning: found non-ellipse');
      break;
    end
    new = [z; a; b; alpha];
    if (step == 1),
      delta = 1;
    else
      delta = norm (new - old);
    end
    old = new;
    
    c = cos(alpha); s = sin(alpha);
    Q = [c -s; s c];
    % compute initial approximations for phi_i
    du  = Q'*[ X(:,1)-z(1) X(:,2)-z(2)]';
    phi = angle(du(1,:)/a  + sqrt(-1)*du(2,:)/b)';
  
    %% weights is inverse of ellipse radius
    C = cos(phi); S = sin(phi);
    W = (a*b) ./ (a^2*S.^2 + b^2*C.^2).^(3/2);

  end % while

  z  = new(1:2);
  a  = new(3);
  b  = new(4);
  alpha = new(5);

end % lyle