bookstein.m

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

function [z, a, b, alpha] = bookstein (X, show);
%BOOKSTEIN       Algebraic ellipse fit 
%
% [z, a, b, alpha] = bookstein (X, show{0});
%
% Approximate ellipse to points <X(i,1),X(i,2)>.
% Invariant under euclidian transformation, see
% BOOKSTEIN, "Fitting conic secions to scattered data", in
% Computer graphics & image processing 9, 56-71 (1979).
%
% X: given points <X(i,1),X(i,2)>
% show: if (show == 1), test output
%
% z, a, b, alpha: parameters for ellipse found

  if (nargin < 2), show = 0; end;

  m      = size(X,1);
  A      = [X(:,1).^2  X(:,1).*X(:,2) X(:,2).^2 ...
            X(:,1) X(:,2) ones(size(X(:,1)))];
  
  S = A'*A;
  T = S(1:3,1:3) - S(1:3,4:6)*(S(4:6,4:6)'\S(4:6,1:3));
  T = diag([1,2,1])*T;
  [V, D] = eig(T);
    
  emin = 0;
  kmin = 0;
  for k = 1:3,
    A = V(1,k); B = V(2,k); C = V(3,k);
    I0 = (A + C);
    I1 = (A*C - B^2/4);
    if (I1 <= 0),
      % this is not an ellipse !
    else
      val   = (I0^2 - 4*I1)/(I0^2 - 2*I1);
      if (emin == 0 | val < emin),
        emin = val;
        kmin = k;
      end
    end
  end
  if (kmin == 0), kmin = 1; end; % not an ellipse 
  y1 = V(:,kmin);
  y2 = -(S(4:6,4:6)')\(S(1:3,4:6)'*y1);
  u  = [y1; y2];
  
  [z, a, b, alpha, err] = ellipse_params (u, show);

end % bookstein