function ellipsoid = makeEllipsoid(grid_size, position, radii, plot_ellipsoid, binary)
%MAKEELLIPSOID Create a binary map of a filled ellipsoid within a 3D grid.
%
% DESCRIPTION:
%     makeEllipsoid creates a binary map of a filled ellipsoid within a
%     three-dimensional grid (the ellipsoid position is denoted by 1's in the
%     matrix with 0's elsewhere). 
%
% USAGE:
%     ellipsoid = makeEllipsoid(grid_size, position, radii)
%     ellipsoid = makeEllipsoid(grid_size, position, radii, plot_ellipsoid)
%     ellipsoid = makeEllipsoid(grid_size, position, radii, [], binary)
%
% INPUTS:
%     grid_size       - size of the 3D grid given as a three element
%                       vector [Nx, Ny, Nz] [grid points] 
%     position        - centre of the ellipsoid given as a three element
%                       vector [cx, cy, cz] [grid points]  
%     radii           - radii of ellipsoid in each direction (semi-axis
%                       lengths) given as a three element vector [rx, ry,
%                       rz] [grid points]  
%
% OPTIONAL INPUTS:
%     plot_ellipsoid  - Boolean controlling whether the ellipsoid is
%                       plotted using voxelPlot (default = false)
%     binary          - Boolean controlling whether the ellipsoid map is
%                       returned as a double precision matrix (false) or
%                       a logical matrix (true) (default = false)
%
% OUTPUTS:
%     ellipsoid       - 3D binary map of a filled ellipsoid
%
% ABOUT:
%     author          - Bradley Treeby
%     date            - 13th February 2019
%     last update     - 13th February 2019
%       
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2019 Bradley Treeby
%
% See also makeBall, makeSphere 

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
% 
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details. 
% 
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% check for plot input
if nargin < 4 || isempty(plot_ellipsoid)
    plot_ellipsoid = false;
end

% check for binary input
if nargin < 5 || isempty(binary)
    binary = false;
end

% force integer grid values
Nx = round(grid_size(1));
Ny = round(grid_size(2));
Nz = round(grid_size(3));

% assign center and radius
cx = position(1);
cy = position(2);
cz = position(3);
rx = radii(1);
ry = radii(2);
rz = radii(3);

% check for zero values for center position, and set to middle of grid
if cx == 0
    cx = floor(Nx / 2) + 1;
end
if cy == 0
    cy = floor(Ny / 2) + 1;
end
if cz == 0
    cz = floor(Nz / 2) + 1;
end

% create grid axes
[x_index, y_index, z_index] = ndgrid(1:Nx, 1:Ny, 1:Nz) ;

% create ellipsoid based on code from https://stackoverflow.com/questions/36420023/generating-a-3d-binary-mask-of-geometric-shapes-in-matlab
ellipsoid = ((x_index - cx).^2 / (rx.^2)) + ((y_index - cy).^2 / (ry.^2)) + ((z_index - cz).^2 / (rz.^2)) <= 1;

% convert to double precision if not binary
if ~binary
    ellipsoid = double(ellipsoid);
end

% plot results
if plot_ellipsoid
    voxelPlot(double(ellipsoid));
end