octave_packages/geometry-1.5.0/geom2d/cov2ellipse.m

   1 ## Copyright (c) 2012 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 %% -*- texinfo -*-
  17 %% @deftypefn {Function File} {@var{ellipse} = } cov2ellipse (@var{K})
  18 %% @deftypefnx {Function File} {[@var{ra} @var{rb} @var{theta}] = } cov2ellipse (@var{K})
  19 %% @deftypefnx {Function File} {@dots{} = } cov2ellipse (@dots{}, @samp{tol},@var{tol})
  20 %% Calculates ellipse parameters from covariance matrix.
  21 %%
  22 %% @var{K} must be symmetric positive (semi)definite. The optional argument
  23 %% @samp{tol} sets the tolerance for the verification of the
  24 %% positive-(semi)definiteness of the matrix @var{K} (see @command{isdefinite}).
  25 %%
  26 %% If only one output argument is supplied a vector defining a ellipse is returned
  27 %% as defined in @command{ellipses2d}. Otherwise the angle @var{theta} is given
  28 %% in radians.
  29 %%
  30 %% Run @code{demo cov2ellipse} to see an example.
  31 %%
  32 %% @seealso{ellipses2d, cov2ellipse, drawEllipse}
  33 %% @end deftypefn
  34
  35 function varargout = cov2ellipse (K, varargin);
  36
  37   parser = inputParser ();
  38   parser.FunctionName = "cov2ellipse";
  39   parser = addParamValue (parser,'Tol', 100*eps*norm (K, "fro"), @(x)x>0);
  40   parser = parse(parser,varargin{:});
  41
  42   if isdefinite (K,parser.Results.Tol) == -1
  43     print_usage
  44   end
  45   [R S W] = svd (K);
  46   theta = atan (R(1,1)/R(2,2));
  47   v     = sort (diag(S), 'ascend')';
  48
  49   if nargout == 1
  50     varargout{1} = [0 0 v theta*180/pi];
  51   elseif nargout == 3
  52     varargout{1} = v(1);
  53     varargout{2} = v(2);
  54     varargout{3} = theta;
  55   end
  56
  57 endfunction
  58
  59 %!demo
  60 %! K = [2 1; 1 2];
  61 %! L = chol(K,'lower');
  62 %! u = randn(1e3,2)*L';
  63 %!
  64 %! elli = cov2ellipse (K)
  65 %!
  66 %! figure(1)
  67 %! plot(u(:,1),u(:,2),'.r');
  68 %! hold on;
  69 %! drawEllipse(elli,'linewidth',2);
  70 %! hold off
  71 %! axis tight