drawEllipsoid

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;

   See also
   spheres, drawSphere, equivalentEllipsoid, ellipsoid, drawTorus,
   drawCuboid

0001 function varargout = drawEllipsoid(elli, varargin)
0002 %DRAWELLIPSOID Draw a 3D ellipsoid.
0003 %
0004 %   drawEllipsoid(ELLI)
0005 %   Displays a 3D ellipsoid on current axis. ELLI is given by:
0006 %   [XC YC ZC A B C PHI THETA PSI],
0007 %   where (XC, YC, ZC) is the ellipsoid center, A, B and C are the half
0008 %   lengths of the ellipsoid main axes, and PHI THETA PSI are Euler angles
0009 %   representing ellipsoid orientation, in degrees.
0010 %
0011 %   drawEllipsoid(..., 'drawEllipses', true)
0012 %   Also displays the main 3D ellipses corresponding to XY, XZ and YZ
0013 %   planes.
0014 %
0015 %
0016 %   Example
0017 %     figure; hold on;
0018 %     drawEllipsoid([10 20 30   50 30 10   5 10 0]);
0019 %     axis equal;
0020 %
0021 %     figure; hold on;
0022 %     elli = [10 20 30   50 30 10   5 10 0];
0023 %     drawEllipsoid(elli, 'FaceColor', 'r', ...
0024 %         'drawEllipses', true, 'EllipseColor', 'b', 'EllipseWidth', 3);
0025 %     axis equal;
0026 %
0027 %   See also
0028 %   spheres, drawSphere, equivalentEllipsoid, ellipsoid, drawTorus,
0029 %   drawCuboid
0030 %
0031 
0032 % ------
0033 % Author: David Legland
0034 % e-mail: david.legland@inrae.fr
0035 % Created: 2011-03-12,    using Matlab 7.9.0.529 (R2009b)
0036 % Copyright 2011 INRA - Cepia Software Platform.
0037 
0038 
0039 %% Default values
0040 
0041 % number of meridians
0042 nPhi    = 32;
0043 
0044 % number of parallels
0045 nTheta  = 16;
0046 
0047 % settings for drawing ellipses
0048 drawEllipses = false;
0049 ellipseColor = 'k';
0050 ellipseWidth = 1;
0051 
0052 drawAxes = false;
0053 axesColor = 'k';
0054 axesWidth = 2;
0055 
0056 
0057 %% Extract input arguments
0058 
0059 % Parse the input (try to extract center coordinates and radius)
0060 if nargin == 0
0061     % no input: assumes ellipsoid with default shape
0062     elli = [0 0 0 5 4 3 0 0 0];
0063 end
0064 
0065 % default set of options for drawing meshes
0066 options = {'FaceColor', 'g', 'linestyle', 'none'};
0067 
0068 while length(varargin) > 1
0069     switch lower(varargin{1})
0070         case 'nphi'
0071             nPhi = varargin{2};
0072             
0073         case 'ntheta'
0074             nTheta = varargin{2};
0075 
0076         case 'drawellipses'
0077             drawEllipses = varargin{2};
0078             
0079         case 'ellipsecolor'
0080             ellipseColor = varargin{2};
0081             
0082         case 'ellipsewidth'
0083             ellipseWidth = varargin{2};
0084             
0085         case 'drawaxes'
0086             drawAxes = varargin{2};
0087             
0088         case 'axescolor'
0089             axesColor = varargin{2};
0090             
0091         case 'axeswidth'
0092             axesWidth = varargin{2};
0093             
0094         otherwise
0095             % assumes this is drawing option
0096             options = [options varargin(1:2)]; %#ok<AGROW>
0097     end
0098 
0099     varargin(1:2) = [];
0100 end
0101 
0102 
0103 %% Parse numerical inputs
0104 
0105 % Extract ellipsoid parameters
0106 xc  = elli(:,1);
0107 yc  = elli(:,2);
0108 zc  = elli(:,3);
0109 a   = elli(:,4);
0110 b   = elli(:,5);
0111 c   = elli(:,6);
0112 k   = pi / 180;
0113 ellPhi   = elli(:,7) * k;
0114 ellTheta = elli(:,8) * k;
0115 ellPsi   = elli(:,9) * k;
0116 
0117 
0118 %% Coordinates computation
0119 
0120 % convert unit basis to ellipsoid basis
0121 sca     = createScaling3d(a, b, c);
0122 rotZ    = createRotationOz(ellPhi);
0123 rotY    = createRotationOy(ellTheta);
0124 rotX    = createRotationOx(ellPsi);
0125 tra     = createTranslation3d([xc yc zc]);
0126 
0127 % concatenate transforms
0128 trans   = tra * rotZ * rotY * rotX * sca;
0129 
0130 
0131 %% parametrisation of ellipsoid
0132 
0133 % spherical coordinates
0134 theta   = linspace(0, pi, nTheta+1);
0135 phi     = linspace(0, 2*pi, nPhi+1);
0136 
0137 % convert to cartesian coordinates
0138 sintheta = sin(theta);
0139 x = cos(phi') * sintheta;
0140 y = sin(phi') * sintheta;
0141 z = ones(length(phi),1) * cos(theta);
0142 
0143 % transform mesh vertices
0144 [x, y, z] = transformPoint3d(x, y, z, trans);
0145 
0146 
0147 %% parametrisation of ellipses
0148 
0149 if drawEllipses
0150     % parametrisation for ellipses
0151     nVertices = 120;
0152     t = linspace(0, 2*pi, nVertices+1);
0153 
0154     % XY circle
0155     xc1 = cos(t');
0156     yc1 = sin(t');
0157     zc1 = zeros(size(t'));
0158 
0159     % XZ circle
0160     xc2 = cos(t');
0161     yc2 = zeros(size(t'));
0162     zc2 = sin(t');
0163 
0164     % YZ circle
0165     xc3 = zeros(size(t'));
0166     yc3 = cos(t');
0167     zc3 = sin(t');
0168 
0169     % compute transformed ellipses
0170     [xc1, yc1, zc1] = transformPoint3d(xc1, yc1, zc1, trans);
0171     [xc2, yc2, zc2] = transformPoint3d(xc2, yc2, zc2, trans);
0172     [xc3, yc3, zc3] = transformPoint3d(xc3, yc3, zc3, trans);
0173 end
0174 
0175 %% parametrisation of main axis edges
0176 
0177 if drawAxes
0178     axesEndings = [-1 0 0; +1 0 0; 0 -1 0; 0 +1 0; 0 0 -1; 0 0 +1];
0179     axesEndings = transformPoint3d(axesEndings, trans);
0180 end
0181 
0182 
0183 %% Drawing
0184 
0185 ellipseOptions = {'color', ellipseColor, 'LineWidth', ellipseWidth};
0186 axesOptions = {'color', axesColor, 'LineWidth', axesWidth};
0187 
0188 % Process output
0189 if nargout == 0
0190     % no output: draw the ellipsoid
0191     surf(x, y, z, options{:});
0192 
0193     if drawEllipses
0194         plot3(xc1, yc1, zc1, ellipseOptions{:});
0195         plot3(xc2, yc2, zc2, ellipseOptions{:});
0196         plot3(xc3, yc3, zc3, ellipseOptions{:});
0197     end
0198     
0199     if drawAxes
0200         drawEdge3d([axesEndings(1,:), axesEndings(2,:)], axesOptions{:});
0201         drawEdge3d([axesEndings(3,:), axesEndings(4,:)], axesOptions{:});
0202         drawEdge3d([axesEndings(5,:), axesEndings(6,:)], axesOptions{:});
0203     end
0204     
0205 elseif nargout == 1
0206     % one output: draw the ellipsoid and return handle
0207     varargout{1} = surf(x, y, z, options{:});
0208     if drawEllipses
0209         plot3(xc1, yc1, zc1, ellipseOptions{:});
0210         plot3(xc2, yc2, zc2, ellipseOptions{:});
0211         plot3(xc3, yc3, zc3, ellipseOptions{:});
0212     end
0213     
0214 elseif nargout == 3
0215     % 3 outputs: return computed coordinates
0216     varargout{1} = x; 
0217     varargout{2} = y; 
0218     varargout{3} = z; 
0219     if drawEllipses
0220         plot3(xc1, yc1, zc1, ellipseOptions{:});
0221         plot3(xc2, yc2, zc2, ellipseOptions{:});
0222         plot3(xc3, yc3, zc3, ellipseOptions{:});
0223     end
0224     
0225 elseif nargout == 4 && drawEllipses
0226     % Also returns handles to ellipses
0227     varargout{1} = surf(x, y, z, options{:});
0228     varargout{2} = plot3(xc1, yc1, zc1, ellipseOptions{:});
0229     varargout{3} = plot3(xc2, yc2, zc2, ellipseOptions{:});
0230     varargout{4} = plot3(xc3, yc3, zc3, ellipseOptions{:});
0231     
0232 end
0233