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