## Copyright (C) 2001, James B. Rawlings and John W. Eaton
##
## This program is free software; you can redistribute it and/or
## modify it under the terms of the GNU General Public License as
## published by the Free Software Foundation; either version 3, or (at
## your option) any later version.
##
## This program is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; see the file COPYING.  If not, write to
## the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
## MA 02111-1307, USA.

## -*- texinfo -*-
## @deftypefn {Function File} ellipse (@var{a}, @var{level})
## @deftypefnx{Function File} ellipse (@var{a}, @var{level}, @var{shift})
## @deftypefnx{Function File} ellipse (@var{a}, @var{level}, @var{shift}, @var{n})
## @deftypefnx{Function File} ellipse (@var{a}, @var{level}, @var{shift}, @var{n}, @dots{})
## @deftypefnx{Function File} address@hidden =} ellipse (@dots{})
## @deftypefnx{Function File} address@hidden, @var{y}, @var{major}, @var{minor}, @
## @var{bbox}] =} ellipse (@dots{})
##
## Plot an ellipse with principal axes given the Eigen vectors of a 2 by 2 matrix.
##
## The Eigen vectors of @var{a} multiplied by the corresponding Eigen values and
## @var{level} are used as the principal axes of the plotted ellipse. By default
## @var{level} is 1. If it is a vector, several ellipses are plotted.
##
## The optional argument @var{shift} controls the center of the ellipse. By default
## this is a two-vector of zeros.
##
## The optional argument @var{n} controls the number of points used to plot the
## ellipse. By default 100 points are used.
##
## The rest of the arguments are passed to @code{plot}, and can be used to control
## the visual style of the plot.
##
## Besides the usual style arguments that can be passed to plot, @code{ellipse}
## also accepts the following strings that adds further graphics to the figure.
##
## @table @t
## @item "majoraxis"
## Also draw the major axis of the ellipse.
## @item "minoraxis"
## Also draw the minor axis of the ellipse.
## @item "boundingbox"
## Also draw the bounding box of the ellipse.
## @end table
##
## @noindent
## Input arguments following any of these strings will be used to control the
## style of the extra graphics elements.
##
## If one output argument is requested, a graphics handle of the plot is returned.
##
## If more than one output argument is requested, no plot is generated. Instead
## the data used to generate the plot is returned. @var{x} and @var{y} are
## @var{n}-vectors containing the points on the ellipse. @var{major} and @var{minor}
## contains the axes of the ellipse, and @var{bbox} contains the bounding box of
## the ellipse.
##
## The following example generates non-isotropic normally distributed data, and
## plots the ellipses with Mahalanobis' distance of 1 and 3 to the center.
##
## @example
## ## Generate data
## X = randn (100, 2);
## scale = 2;
## X (:, 1) *= scale;
## theta = 0.4;
## R = [cos(theta), -sin(theta); sin(theta), cos(theta)];
## X = X * R;
##
## ## Estimate mean and covariance matrix
## mu = mean (X);
## C = cov (X);
##
## ## Plot data and ellipses
## figure
## plot (X (:, 1), X (:, 2), 'ro');
## hold on, ellipse (inv (C), [1, 3], mu, :, "k"); hold off
## axis equal
## @end example
##
## @noindent
## Note the use of @code{:} to denote the default value of @var{n}. To add a
## green major axis of width 2, and a blue bounding box to the figure, replace
## the call to @code{ellipse} with
##
## @example
## ellipse (inv (C), [1, 3], mu, :, "k", "majoraxis", "g", ...
##          "linewidth", 2, "boundingbox", "b");
## @end example
## @seealso{plot}
## @end deftypefn

function varargout = ellipse (semiaxes, level = 1, shift = [0, 0], n = 100, varargin)
  ## Check input
  if (nargin < 1)
    error ("ellipse: not enough input arguments");
  endif
  
  if (!isnumeric (semiaxes) || !isequal (size (semiaxes), [2, 2]))
    error ("ellipse: first input argument must be a 2x2 matrix");
  endif
  
  if (!isvector (level))
    error ("ellipse: second input argument must be a vector");
  endif
  
  if (!isnumeric (shift) || numel (shift) != 2)
    error ("ellipse: third input argument must be a two-vector");
  endif
  shift = shift (:).';
  
  if (!isscalar (n) || n < 0 || n != round (n))
    error ("ellipse: fourth input argument must be a positive integer");
  endif
  
  ## Possibly override default options
  available_options = {"minoraxis", "majoraxis", "boundingbox"};
  num_options = length (available_options);
  option_idx = zeros (1, num_options);
  for idx = 1:length (varargin)
    for k = 1:num_options
      opt = available_options {k};
      if (strcmpi (opt, varargin {idx}))
        option_idx (k) = idx;
      endif
    endfor
  endfor
  
  [option_idx, idx] = sort (option_idx);
  available_options = available_options (idx);
  option_idx = [option_idx, length(varargin)+1];
  
  ellipse_opt = minoraxis_opt = majoraxis_opt = boundingbox_opt = {};

  last_idx = option_idx (find (option_idx > 0, 1));
  ellipse_opt = varargin (1:last_idx-1);
  for k = 1:num_options
    opt = available_options {k};
    if (option_idx (k) > 0)
      val = varargin (last_idx+1:option_idx (k+1)-1);
      last_idx = option_idx (k+1);
      eval (sprintf ("%s_opt = val;", opt));
    endif
  endfor
  
  ## We're ready to do the actual work
  [v, l] = eig (semiaxes);
  
  dl = diag (l);
  if (any (imag (dl)) || any (dl <= 0))
    error ("ellipse: first input argument must be positive definite");
  endif

  ## Generate contour data.
  a = 1 / sqrt (dl (1));
  b = 1 / sqrt (dl (2));

  a = a * level (:);
  b = b * level (:);
    
  t = linspace (0, 2*pi, n);

  xt = a * cos (t);  xt = xt.';
  yt = b * sin (t);  yt = yt.';

  ## Rotate the contours.
  ra = atan2 (v (2, 1), v (1, 1));

  cos_ra = cos (ra);
  sin_ra = sin (ra);

  x = xt * cos_ra - yt * sin_ra + shift (1);
  y = xt * sin_ra + yt * cos_ra + shift (2);

  ## Endpoints of the major and minor axes.
  [mx, ii] = max (level);
    
  minor = (v * diag ([a(ii), b(ii)])).';
  major = minor;

  major (2, :) = -major (1,:);
  minor (1, :) = -minor (2,:);

  t = [1; 1] * shift;

  major = major + t;
  minor = minor + t;

  ## Bounding box for the ellipse using magic formula.
  ainv = inv (semiaxes);
  xbox = level (ii) * sqrt (ainv (1,1));
  ybox = level(ii) * sqrt (ainv (2,2));

  bbox = [xbox ybox; xbox -ybox; -xbox -ybox; -xbox ybox; xbox ybox];

  t = [1; 1; 1; 1; 1] * shift;
  bbox = bbox + t;
  
  ## Should we plot the ellipse, or return it?
  if (nargout <= 1)
    plot_args = {x, y, ellipse_opt{:}};

    if (!isempty (minoraxis_opt))
      plot_args = {plot_args{:}, minor(:,1), minor(:,2), minoraxis_opt{:}};
    endif

    if (!isempty (majoraxis_opt))
      plot_args = {plot_args{:}, major(:,1), major(:,2), majoraxis_opt{:}};
    endif

    if (!isempty (boundingbox_opt))
      plot_args = {plot_args{:}, bbox(:,1), bbox(:,2), boundingbox_opt{:}};
    endif
    handle = plot (plot_args {:});

    if (nargout == 1)
      varargout {1} = handle;
    endif
  else
    varargout {1} = x;
    varargout {2} = y;
    varargout {3} = major;
    varargout {4} = minor;
    varargout {5} = bbox;
  endif
endfunction

%!demo
%! ## Generate data
%! X = randn (100, 2);
%! scale = 2;
%! X (:, 1) *= scale;
%! theta = 0.4;
%! R = [cos(theta), -sin(theta); sin(theta), cos(theta)];
%! X = X * R;
%!
%! ## Estimate mean and covariance matrix
%! mu = mean (X);
%! C = cov (X);
%!
%! ## Plot data and ellipses
%! figure
%! plot (X (:, 1), X (:, 2), 'ro', 'linewidth', 3);
%! hold on, ellipse (inv (C), [1, 3], mu, :, "k"); hold off
%! axis equal

%!demo
%! ## Generate data
%! X = randn (100, 2);
%! scale = 2;
%! X (:, 1) *= scale;
%! theta = 0.4;
%! R = [cos(theta), -sin(theta); sin(theta), cos(theta)];
%! X = X * R;
%!
%! ## Estimate mean and covariance matrix
%! mu = mean (X);
%! C = cov (X);
%!
%! ## Plot data and ellipses. This time also with a major axis and a bounding box.
%! figure
%! plot (X (:, 1), X (:, 2), 'ro', 'linewidth', 3);
%! hold on
%! ellipse (inv (C), [1, 3], mu, :, "k", "majoraxis", "g", ...
%!          "linewidth", 2, "boundingbox", "b");
%! hold off
%! axis equal