## Copyright (C) 2024 David Legland
## All rights reserved.
## 
## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are met:
## 
##     1 Redistributions of source code must retain the above copyright notice,
##       this list of conditions and the following disclaimer.
##     2 Redistributions in binary form must reproduce the above copyright
##       notice, this list of conditions and the following disclaimer in the
##       documentation and/or other materials provided with the distribution.
## 
## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
## 
## The views and conclusions contained in the software and documentation are
## those of the authors and should not be interpreted as representing official
## policies, either expressed or implied, of the copyright holders.

function varargout = drawEllipsoid(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;
%
%   See also 
%   spheres, drawSphere, equivalentEllipsoid, ellipsoid, drawTorus, 
%   drawCuboid 
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2011-03-12, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2023 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;

drawAxes = false;
axesColor = 'k';
axesWidth = 2;


%% Extract input arguments

% extract handle of axis to draw on
[hAx, varargin] = parseAxisHandle(varargin{:});

% retrieve parameters of ellipsoid
elli = varargin{1};
varargin(1) = [];

% 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};
            
        case 'drawaxes'
            drawAxes = varargin{2};
            
        case 'axescolor'
            axesColor = varargin{2};
            
        case 'axeswidth'
            axesWidth = 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;


%% 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;


%% parametrisation of ellipsoid

% 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);

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


%% parametrisation of ellipses

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');

    % compute transformed ellipses
    [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

%% parametrisation of main axis edges

if drawAxes
    axesEndings = [-1 0 0; +1 0 0; 0 -1 0; 0 +1 0; 0 0 -1; 0 0 +1];
    axesEndings = transformPoint3d(axesEndings, trans);
end


%% Drawing 

ellipseOptions = {'color', ellipseColor, 'LineWidth', ellipseWidth};
axesOptions = {'color', axesColor, 'LineWidth', axesWidth};

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

    if drawEllipses
        plot3(hAx, xc1, yc1, zc1, ellipseOptions{:});
        plot3(hAx, xc2, yc2, zc2, ellipseOptions{:});
        plot3(hAx, xc3, yc3, zc3, ellipseOptions{:});
    end
    
    if drawAxes
        drawEdge3d(hAx, [axesEndings(1,:), axesEndings(2,:)], axesOptions{:});
        drawEdge3d(hAx, [axesEndings(3,:), axesEndings(4,:)], axesOptions{:});
        drawEdge3d(hAx, [axesEndings(5,:), axesEndings(6,:)], axesOptions{:});
    end
    
elseif nargout == 1
    % one output: draw the ellipsoid and return handle 
    varargout{1} = surf(x, y, z, options{:});
    if drawEllipses
        plot3(hAx, xc1, yc1, zc1, ellipseOptions{:});
        plot3(hAx, xc2, yc2, zc2, ellipseOptions{:});
        plot3(hAx, 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(hAx, xc1, yc1, zc1, ellipseOptions{:});
        plot3(hAx, xc2, yc2, zc2, ellipseOptions{:});
        plot3(hAx, xc3, yc3, zc3, ellipseOptions{:});
    end
    
elseif nargout == 4 && drawEllipses
    % Also returns handles to ellipses
    varargout{1} = surf(hAx, x, y, z, options{:});
    varargout{2} = plot3(hAx, xc1, yc1, zc1, ellipseOptions{:});
    varargout{3} = plot3(hAx, xc2, yc2, zc2, ellipseOptions{:});
    varargout{4} = plot3(hAx, xc3, yc3, zc3, ellipseOptions{:});
    
end