Daniel@0: % PLOTCOV3 - Plots a covariance ellipsoid with axes for a trivariate
Daniel@0: %            Gaussian distribution.
Daniel@0: %
Daniel@0: % Usage:
Daniel@0: %   [h, s] = plotcov3(mu, Sigma[, OPTIONS]);
Daniel@0: % 
Daniel@0: % Inputs:
Daniel@0: %   mu    - a 3 x 1 vector giving the mean of the distribution.
Daniel@0: %   Sigma - a 3 x 3 symmetric positive semi-definite matrix giving
Daniel@0: %           the covariance of the distribution (or the zero matrix).
Daniel@0: %
Daniel@0: % Options:
Daniel@0: %   'conf'      - a scalar between 0 and 1 giving the confidence
Daniel@0: %                 interval (i.e., the fraction of probability mass to
Daniel@0: %                 be enclosed by the ellipse); default is 0.9.
Daniel@0: %   'num-pts'   - if the value supplied is n, then (n + 1)^2 points
Daniel@0: %                 to be used to plot the ellipse; default is 20.
Daniel@0: %   'plot-opts' - a cell vector of arguments to be handed to PLOT3
Daniel@0: %                 to contol the appearance of the axes, e.g., 
Daniel@0: %                 {'Color', 'g', 'LineWidth', 1}; the default is {}
Daniel@0: %   'surf-opts' - a cell vector of arguments to be handed to SURF
Daniel@0: %                 to contol the appearance of the ellipsoid
Daniel@0: %                 surface; a nice possibility that yields
Daniel@0: %                 transparency is: {'EdgeAlpha', 0, 'FaceAlpha',
Daniel@0: %                 0.1, 'FaceColor', 'g'}; the default is {}
Daniel@0: % 
Daniel@0: % Outputs:
Daniel@0: %   h     - a vector of handles on the axis lines
Daniel@0: %   s     - a handle on the ellipsoid surface object
Daniel@0: %
Daniel@0: % See also: PLOTCOV2
Daniel@0: 
Daniel@0: % Copyright (C) 2002 Mark A. Paskin
Daniel@0: %
Daniel@0: % This program is free software; you can redistribute it and/or modify
Daniel@0: % it under the terms of the GNU General Public License as published by
Daniel@0: % the Free Software Foundation; either version 2 of the License, or
Daniel@0: % (at your option) any later version.
Daniel@0: %
Daniel@0: % This program is distributed in the hope that it will be useful, but
Daniel@0: % WITHOUT ANY WARRANTY; without even the implied warranty of
Daniel@0: % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Daniel@0: % General Public License for more details.
Daniel@0: %
Daniel@0: % You should have received a copy of the GNU General Public License
Daniel@0: % along with this program; if not, write to the Free Software
Daniel@0: % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
Daniel@0: % USA.
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: 
Daniel@0: function [h, s] = plotcov3(mu, Sigma, varargin)
Daniel@0: 
Daniel@0: if size(Sigma) ~= [3 3], error('Sigma must be a 3 by 3 matrix'); end
Daniel@0: if length(mu) ~= 3, error('mu must be a 3 by 1 vector'); end
Daniel@0: 
Daniel@0: [p, ...
Daniel@0:  n, ...
Daniel@0:  plot_opts, ...
Daniel@0:  surf_opts] = process_options(varargin, 'conf', 0.9, ...
Daniel@0: 					'num-pts', 20, ...
Daniel@0: 			                'plot-opts', {}, ...
Daniel@0: 			                'surf-opts', {});
Daniel@0: h = [];
Daniel@0: holding = ishold;
Daniel@0: if (Sigma == zeros(3, 3))
Daniel@0:   z = mu;
Daniel@0: else
Daniel@0:   % Compute the Mahalanobis radius of the ellipsoid that encloses
Daniel@0:   % the desired probability mass.
Daniel@0:   k = conf2mahal(p, 3);
Daniel@0:   % The axes of the covariance ellipse are given by the eigenvectors of
Daniel@0:   % the covariance matrix.  Their lengths (for the ellipse with unit
Daniel@0:   % Mahalanobis radius) are given by the square roots of the
Daniel@0:   % corresponding eigenvalues.
Daniel@0:   if (issparse(Sigma))
Daniel@0:     [V, D] = eigs(Sigma);
Daniel@0:   else
Daniel@0:     [V, D] = eig(Sigma);
Daniel@0:   end
Daniel@0:   if (any(diag(D) < 0))
Daniel@0:     error('Invalid covariance matrix: not positive semi-definite.');
Daniel@0:   end
Daniel@0:   % Compute the points on the surface of the ellipsoid.
Daniel@0:   t = linspace(0, 2*pi, n);
Daniel@0:   [X, Y, Z] = sphere(n);
Daniel@0:   u = [X(:)'; Y(:)'; Z(:)'];
Daniel@0:   w = (k * V * sqrt(D)) * u;
Daniel@0:   z = repmat(mu(:), [1 (n + 1)^2]) + w;
Daniel@0: 
Daniel@0:   % Plot the axes.
Daniel@0:   L = k * sqrt(diag(D));
Daniel@0:   h = plot3([mu(1); mu(1) + L(1) * V(1, 1)], ...
Daniel@0: 	    [mu(2); mu(2) + L(1) * V(2, 1)], ...
Daniel@0: 	    [mu(3); mu(3) + L(1) * V(3, 1)], plot_opts{:});
Daniel@0:   hold on;
Daniel@0:   h = [h; plot3([mu(1); mu(1) + L(2) * V(1, 2)], ...
Daniel@0: 		[mu(2); mu(2) + L(2) * V(2, 2)], ...
Daniel@0: 		[mu(3); mu(3) + L(2) * V(3, 2)], plot_opts{:})];
Daniel@0:   h = [h; plot3([mu(1); mu(1) + L(3) * V(1, 3)], ...
Daniel@0: 		[mu(2); mu(2) + L(3) * V(2, 3)], ...
Daniel@0: 		[mu(3); mu(3) + L(3) * V(3, 3)], plot_opts{:})];
Daniel@0: end
Daniel@0: 
Daniel@0: s = surf(reshape(z(1, :), [(n + 1) (n + 1)]), ...
Daniel@0: 	 reshape(z(2, :), [(n + 1) (n + 1)]), ...
Daniel@0: 	 reshape(z(3, :), [(n + 1) (n + 1)]), ...
Daniel@0: 	 surf_opts{:});
Daniel@0: 
Daniel@0: if (~holding) hold off; end