function [varargout]=ellips(MI,SIGMA,R,N)
% [varargout]=ellips(MI,SIGMA,R,N)
%
% ELLIPS creates vectors or matrices for the ellipse plot by
%  MATLAB functions as PLOT, FILL in 2D or SURF, MESH in 3D.
%
%  An ellipsoid is described in terms of
%                  R^2=(X-MI)'*inv(SIGMA)*(X-MI)
%  where
%   MI [Mx1]    - is the center of ellipsoid. M is the dimension.
%   SIGMA [MxM] - determines on shape.
%   R [1x1]     - radius of ellipsoid.
%  The ellipsoid is interpolated by N lines (N=20 is default).
%
%  Order of the matrix SIGMA (as well as length MI) determines
%  the number of returned variables:
%  For 2D (M=2)
%    [X,Y]=ellips(MI,SIGMA,R,N)
%  and for 3D (M=3) it is
%    [X,Y,Z]=ellips(MI,SIGMA,R,N)
%
%  Example:
%    [x,y]=ellipse([1;1],[1 0.5;0.5 1],1);
%    plot(x,y);
%
% See also PLOT, FILL, MESH, SURFACE.
%
% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz
% Written Vojtech Franc (diploma thesis) 27.02.2000, 4.4.2000
% Modifications
% 24. 6.00 V. Hlavac, comments polished.

MI=MI(:);

% input arguments processing
if nargin < 3,
   help('ellips');
   return;
elseif nargin < 4
   N=20;
end

% dimension
DIM=size(SIGMA,1);
INVSIGMA=inv(SIGMA);

if DIM == 3,
   % 3D ellipsoid, creates the matrix X,Y,Z appropriate for using in
   % the functions mesh,surf, etc.

   [X,Y,Z]=sphere(N);

   X=R*X;
   Y=R*Y;
   Z=R*Z;

   [A,p]=chol(INVSIGMA);
   if p ~= 0,
      error('Covariance matrix must be positive definite');
   end
   A=inv(A);

   [ROWS,COLUMNS]=size(X);

   % transforms the sphere to the ellipse according to the matrix INVSIGMA
   for i=1:ROWS,
      P=[X(i,:);Y(i,:);Z(i,:)];
      Q=A*P;

      % if the translation is given then translate points
      Q=Q+repmat(MI(:),1,COLUMNS);

      X(i,:)=Q(1,:);
      Y(i,:)=Q(2,:);
      Z(i,:)=Q(3,:);
   end % for i=1:ROWS

   % return variables
   varargout{1}=X;
   varargout{2}=Y;
   varargout{3}=Z;

elseif DIM == 2,

   from=0;
   to=2*pi;
   step=(to-from)/N;
   X=cos([from:step:to]);
   Y=sin([from:step:to]);

   X=R*X;
   Y=R*Y;

   [A,p]=chol(INVSIGMA);
   if p ~= 0,
      error('Covariance matrix must be positive definite');
   end
   A=inv(A);

   P=[X;Y];
   Q=A*P;

   Q=Q+repmat(MI(:),1,N+1);

   X=Q(1,:);
   Y=Q(2,:);

   % return variables
   varargout{1}=X;
   varargout{2}=Y;

else
   help('ellips');
end;