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 (type 'help angles3d'
   for more details).

   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, equivalentEllipsoid, ellipsoid, angles3d
   drawSphere, drawCuboid, drawTorus, drawCylinder

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