Code covered by the BSD License

# geom3d

### David Legland (view profile)

19 Jun 2009 (Updated )

Library to handle 3D geometric primitives: create, intersect, display, and make basic computations

### Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

drawEllipsoid(elli, varargin)
```function varargout = drawEllipsoid(elli, varargin)
%DRAWELLIPSOID Draw a 3D ellipsoid
%
%   drawEllipsoid(ELLI)
%   Displays a 3D ellipsoid on current axis. ELLI is given by:
%   [XC YC ZC A B C PHI THETA PSI],
%   where (XC, YC, ZC) is the ellipsoid center, A, B and C are the half
%   lengths of the ellipsoid main axes, and PHI THETA PSI are Euler angles
%   representing ellipsoid orientation, in degrees.
%
%   drawEllipsoid(..., 'drawEllipses', true)
%   Also displays the main 3D ellipses corresponding to XY, XZ and YZ
%   planes.
%
%
%   Example
%     figure; hold on;
%     drawEllipsoid([10 20 30   50 30 10   5 10 0]);
%     axis equal;
%
%     figure; hold on;
%     elli = [[10 20 30   50 30 10   5 10 0];
%     drawEllipsoid(elli, 'FaceColor', 'r', ...
%         'drawEllipses', true, 'EllipseColor', 'b', 'EllipseWidth', 3);
%     axis equal;
%
%   spheres, drawSphere, inertiaEllipsoid, ellipsoid, drawTorus, drawCuboid
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-03-12,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

%% Default values

% number of meridians
nPhi    = 32;

% number of parallels
nTheta  = 16;

% settings for drawing ellipses
drawEllipses = false;
ellipseColor = 'k';
ellipseWidth = 1;

%% Extract input arguments

% Parse the input (try to extract center coordinates and radius)
if nargin == 0
% no input: assumes ellipsoid with default shape
elli = [0 0 0 5 4 3 0 0 0];
end

% default set of options for drawing meshes
options = {'FaceColor', 'g', 'linestyle', 'none'};

while length(varargin) > 1
switch lower(varargin{1})
case 'nphi'
nPhi = varargin{2};

case 'ntheta'
nTheta = varargin{2};

case 'drawellipses'
drawEllipses = varargin{2};

case 'ellipsecolor'
ellipseColor = varargin{2};

case 'ellipsewidth'
ellipseWidth = varargin{2};

otherwise
% assumes this is drawing option
options = [options varargin(1:2)]; %#ok<AGROW>
end

varargin(1:2) = [];
end

%% Parse numerical inputs

% Extract ellipsoid parameters
xc  = elli(:,1);
yc  = elli(:,2);
zc  = elli(:,3);
a   = elli(:,4);
b   = elli(:,5);
c   = elli(:,6);
k   = pi / 180;
ellPhi   = elli(:,7) * k;
ellTheta = elli(:,8) * k;
ellPsi   = elli(:,9) * k;

%% parametrisation

% spherical coordinates
theta   = linspace(0, pi, nTheta+1);
phi     = linspace(0, 2*pi, nPhi+1);

% convert to cartesian coordinates
sintheta = sin(theta);
x = cos(phi') * sintheta;
y = sin(phi') * sintheta;
z = ones(length(phi),1) * cos(theta);

if drawEllipses
% parametrisation for ellipses
nVertices = 120;
t = linspace(0, 2*pi, nVertices+1);

% XY circle
xc1 = cos(t');
yc1 = sin(t');
zc1 = zeros(size(t'));

% XZ circle
xc2 = cos(t');
yc2 = zeros(size(t'));
zc2 = sin(t');

% YZ circle
xc3 = zeros(size(t'));
yc3 = cos(t');
zc3 = sin(t');
end

%% Coordinates computation

% convert unit basis to ellipsoid basis
sca     = createScaling3d(a, b, c);
rotZ    = createRotationOz(ellPhi);
rotY    = createRotationOy(ellTheta);
rotX    = createRotationOx(ellPsi);
tra     = createTranslation3d([xc yc zc]);

% concatenate transforms
trans   = tra * rotZ * rotY * rotX * sca;

% transform mesh vertices
[x, y, z] = transformPoint3d(x, y, z, trans);

% transform the 3D polygons corresponding to ellipse cross sections
if drawEllipses
[xc1, yc1, zc1] = transformPoint3d(xc1, yc1, zc1, trans);
[xc2, yc2, zc2] = transformPoint3d(xc2, yc2, zc2, trans);
[xc3, yc3, zc3] = transformPoint3d(xc3, yc3, zc3, trans);
end

%% Drawing

ellipseOptions = {'color', ellipseColor, 'LineWidth', ellipseWidth};

% Process output
if nargout == 0
% no output: draw the ellipsoid
surf(x, y, z, options{:});

if drawEllipses
plot3(xc1, yc1, zc1, ellipseOptions{:});
plot3(xc2, yc2, zc2, ellipseOptions{:});
plot3(xc3, yc3, zc3, ellipseOptions{:});
end

elseif nargout == 1
% one output: draw the ellipsoid and return handle
varargout{1} = surf(x, y, z, options{:});
if drawEllipses
plot3(xc1, yc1, zc1, ellipseOptions{:});
plot3(xc2, yc2, zc2, ellipseOptions{:});
plot3(xc3, yc3, zc3, ellipseOptions{:});
end

elseif nargout == 3
% 3 outputs: return computed coordinates
varargout{1} = x;
varargout{2} = y;
varargout{3} = z;
if drawEllipses
plot3(xc1, yc1, zc1, ellipseOptions{:});
plot3(xc2, yc2, zc2, ellipseOptions{:});
plot3(xc3, yc3, zc3, ellipseOptions{:});
end

elseif nargout == 4 && drawEllipses
% Also returns handles to ellipses
varargout{1} = surf(x, y, z, options{:});
varargout{2} = plot3(xc1, yc1, zc1, ellipseOptions{:});
varargout{3} = plot3(xc2, yc2, zc2, ellipseOptions{:});
varargout{4} = plot3(xc3, yc3, zc3, ellipseOptions{:});

end

```