# HG changeset patch
# User Soren Hauberg <address@hidden>
# Date 1233428405 -3600
# Node ID 8b338522286d8208e5e086effb27375ea5dd6853
# Parent  4385bb503467d6cbd834378dd4023b1f5052b858
scripts/plot/ellipse.m: new function

diff -r 4385bb503467 -r 8b338522286d scripts/ChangeLog
--- a/scripts/ChangeLog	Thu Jan 29 18:13:06 2009 -0500
+++ b/scripts/ChangeLog	Sat Jan 31 20:00:05 2009 +0100
@@ -1,3 +1,7 @@
+2009-01-31  Soren Hauberg  <address@hidden>
+
+	* plot/ellipse.m: New function.
+
 2009-01-29  John W. Eaton  <address@hidden>
 
 	* miscellaneous/fileparts.m: Match all possible file separators.
diff -r 4385bb503467 -r 8b338522286d scripts/plot/Makefile.in
--- a/scripts/plot/Makefile.in	Thu Jan 29 18:13:06 2009 -0500
+++ b/scripts/plot/Makefile.in	Sat Jan 31 20:00:05 2009 +0100
@@ -98,6 +98,7 @@
   cylinder.m \
   diffuse.m \
   gnuplot_drawnow.m \
+  ellipse.m \
   ellipsoid.m \
   errorbar.m \
   ezcontourf.m \
diff -r 4385bb503467 -r 8b338522286d scripts/plot/ellipse.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/ellipse.m	Sat Jan 31 20:00:05 2009 +0100
@@ -0,0 +1,276 @@
+## 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{n})
+## @deftypefnx{Function File} ellipse (@var{a}, @var{level}, @var{n}, @var{shift})
+## @deftypefnx{Function File} ellipse (@var{a}, @var{level}, @var{n}, @var{shift}, @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 eigenvectors of a 2 by 2 matrix.
+##
+## The eigenvectors of @var{a} multiplied by the corresponding eigenvalues 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, n = 100, shift = [0, 0], 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
+