wate2.m

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

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

  delta = 1;
  omega = 0;
  myeps = 1e-6;

  m     = size(X, 1);
  sumX  = sum(X)/m;
  X     = X - ones(m,1)*sumX;
  normX = norm(X, inf);
  X     = X/normX;

  oldu = zeros(6,1);
  W1   = [ones(size(X,1), 1)];
  A1   = [X(:,1).^2  X(:,1).*X(:,2) X(:,2).^2 ...
          X(:,1) X(:,2) ones(size(X(:,1)))];
  A2   = [1 0 -1 0 0 0; ...
          0 1  0 0 0 0];
  step = 0;

  while (1), %% breaks

    [U S V] = svd([diag(W1) * A1; omega * norm(W1,inf) * A2]);
    u = V(:,6);
    if (norm (oldu - u) <= myeps), break; end;

    [z, a, b, alpha, err] = ellipse_params (u);
    if (err),
      %% the result is not an ellipse, adjust weight !
      disp ('warning: weighted');
      if (omega == 0), omega = 1;
      else             omega = 4*omega;
      end
    else  
      oldu = u;
    
      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)';
      du  = du';

      phi_geom = ellipse_phi (du, a, b, phi);

      %% weights is "real distance / algebraic weight"
      C      = cos(phi);      S      = sin(phi);
      C_geom = cos(phi_geom); S_geom = sin(phi_geom);
     
      D_geom = sqrt ((du(:,1) - a*C_geom).^2 + (du(:,2) - b*S_geom).^2);
      D      = (du(:,1).^2 + du(:,2).^2) ./ (a^2*C.^2 + b^2*S.^2) - ...
               ones (size(du,1), 1);
      W1     = D_geom./D;

    end % if (is ellipse ?)

    step = step + 1;
    if (step > 20),
      disp ('warning: number of steps exceeeded limit');
      break;
    end
  end % while

  z = z*normX;
  a = a*normX;
  b = b*normX;
  z = z + sumX';
  
  if (show == 3),
    M = abs((diag(W1)*(A1*u))./D_geom);
    M = M/norm(M,-inf);
    disp ('output of M = wQ/g, M = M/norm(M), norm(M - 1)');
    norm (M - 1)
  end

end % wate2