function [hx, hy, hz] = ft_plot_ortho(dat, varargin)

% FT_PLOT_ORTHO plots three orthographic slices through a 3-D volume and interpolates
% the data if needed.
%
% Use as
%   ft_plot_ortho(dat, ...)
% or
%   ft_plot_ortho(dat, mask, ...)
% where dat and mask are equal-sized 3-D arrays.
%
% Additional options should be specified in key-value pairs and can be
%   'style'        = string, 'subplot' or 'intersect' (default = 'subplot')
%   'orientation'  = 3x3 matrix specifying the directions orthogonal through the planes which will be plotted
%   'parents'      = (optional) 3-element vector containing the handles of the axes for the subplots (when style = 'subplot')
%   'surfhandle'   = (optional) 3-element vector containing the handles of the surfaces for each of the sublots (when style = 'subplot'). Parents and surfhandle are mutually exclusive
%   'update'       = (optional) 3-element boolean vector with the axes that should be updated (default = [true true true])
%   'coordsys'     = string, assume the data to be in the specified coordinate system (default = 'unknown')
%
% The following options are supported and passed on to FT_PLOT_SLICE
%   'clim'                = [min max], lower and upper color limits
%   'facealpha'           = transparency when no mask is specified, between 0 and 1 (default = 1)
%   'transform'           = 4x4 homogeneous transformation matrix specifying the mapping from voxel space to the coordinate system in which the data are plotted
%   'location'            = 1x3 vector specifying the intersection point at which the three slices will be plotted. The coordinates should be expressed in the coordinate system of the data.
%   'datmask'             = 3D-matrix with the same size as the matrix dat, serving as opacitymap if the second input argument to the function contains a matrix, this will be used as the mask
%   'maskstyle'           = string, 'opacity' or 'colormix', defines the rendering
%   'background'          = needed when maskstyle is 'colormix', 3D-matrix with the same size as the data matrix, serving as grayscale image that provides the background
%   'interpmethod'        = string specifying the method for the interpolation, see INTERPN (default = 'nearest')
%   'colormap'            = string, see COLORMAP
%   'unit'                = string, can be 'm', 'cm' or 'mm' (default is automatic)
%   'intersectmesh'       = triangulated mesh, see FT_PREPARE_MESH
%   'intersectcolor'      = string, color specification
%   'intersectlinestyle'  = string, line specification
%   'intersectlinewidth'  = number
%
% See also FT_PLOT_SLICE, FT_PLOT_MONTAGE, FT_PLOT_MESH, FT_SOURCEPLOT

% Copyrights (C) 2010, Jan-Mathijs Schoffelen
% Copyrights (C) 2022-2024, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip 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 General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if isstruct(dat)
  if isfield(dat, 'transform') && isfield(dat, 'dim')
    % the input is an MRI structure, call this function recursively
    varargin = ft_setopt(varargin, 'transform', dat.transform);
    if isfield(dat, 'coordsys')
      varargin = ft_setopt(varargin, 'coordsys', dat.coordsys);
    end
    if isfield(dat, 'unit')
      varargin = ft_setopt(varargin, 'unit', dat.unit);
    end
    fn = fieldnames(dat);
    for i=1:numel(fn)
      if isequal(size(dat.(fn{i})), dat.dim)
        ft_info('plotting %s', fn{i});
        [hx, hy, hz] = ft_plot_ortho(dat.(fn{i}), varargin{:});
      end % if
    end % for
    return
  else
    ft_error('unsupported input structure');
  end
end

% parse first input argument(s), it is either
% (dat, varargin)
% (dat, msk, varargin)
% (dat, [], varargin)
% this is done in ft_plot_slice

sellist = 1:numel(varargin);
if ~isempty(sellist)
  if isempty(varargin{1}) || isnumeric(varargin{1})
    sellist(1) = [];
  end
end

% get the optional input arguments
% other options such as location and transform are passed along to FT_PLOT_SLICE
style     = ft_getopt(varargin(sellist), 'style', 'subplot');
ori       = ft_getopt(varargin(sellist), 'orientation', eye(3));
coordsys  = ft_getopt(varargin, 'coordsys');

if strcmp(style, 'subplot')
  parents     = ft_getopt(varargin(sellist), 'parents');
  surfhandle  = ft_getopt(varargin(sellist), 'surfhandle');
  patchhandle = ft_getopt(varargin(sellist), 'patchhandle');
  update      = ft_getopt(varargin(sellist), 'update', [true true true]);
  if ~isempty(surfhandle) && ~isempty(parents)
    ft_error('if specifying handles, you should either specify handles to the axes or to the surface objects, not both');
  end
end

if ~isa(dat, 'double')
  dat = cast(dat, 'double');
end

% determine the orientation key-value pair
keys = varargin(sellist(1:2:end));

sel  = find(strcmp('orientation', keys));
if isempty(sel)
  % add orientation key-value pair if it does not exist
  sel             = numel(keys)+1;
  varargin{2*sel-1} = 'orientation';
  varargin{2*sel} = [];
end

switch style
  case 'subplot'
    if isempty(parents) && isempty(surfhandle)
      Hx = subplot(2,2,1);
      Hy = subplot(2,2,2);
      Hz = subplot(2,2,4);
    elseif ~isempty(parents) && isempty(surfhandle)
      Hx = parents(1);
      Hy = parents(2);
      Hz = parents(3);
    elseif isempty(parents) && ~isempty(surfhandle)
      % determine the parents from the surface handle and use the
      % surfhandle for efficient visualization (overwriting existing data)
      if update(1), Hx = get(surfhandle(1), 'parent'); else Hx = []; end
      if update(2), Hy = get(surfhandle(2), 'parent'); else Hy = []; end
      if update(3), Hz = get(surfhandle(3), 'parent'); else Hz = []; end
    end
    
    if ~isempty(Hx)
      if ~isempty(surfhandle) && update(1)
        varargin(sellist) = ft_setopt(varargin(sellist), 'surfhandle', surfhandle(1));
      end
      if ~isempty(patchhandle) && update(1)
        varargin(sellist) = ft_setopt(varargin(sellist), 'patchhandle', patchhandle(1));
      end
      % swap the first 2 dimensions because of meshgrid vs ndgrid issues
      varargin{2*sel} = ori(2,:);
      set(gcf,'currentaxes',Hx);
      hx = ft_plot_slice(dat, varargin{:}, 'tag', 'x');
      set(Hx, 'view', [0 0]); %, 'xlim', [0.5 size(dat,1)-0.5], 'zlim', [0.5 size(dat,3)-0.5]);
      if isempty(parents)
        % only change axis behavior if no parents are specified
        axis off
      end
    end
    
    if ~isempty(Hy)
      if ~isempty(surfhandle) && update(2)
        varargin(sellist) = ft_setopt(varargin(sellist), 'surfhandle', surfhandle(2));
      end
      if ~isempty(patchhandle) && update(2)
        varargin(sellist) = ft_setopt(varargin(sellist), 'patchhandle', patchhandle(2));
      end
      varargin{2*sel} = ori(1,:);
      set(gcf,'currentaxes',Hy);
      hy = ft_plot_slice(dat, varargin{:}, 'tag', 'y');
      set(Hy, 'view', [90 0]); %, 'ylim', [0.5 size(dat,2)-0.5], 'zlim', [0.5 size(dat,3)-0.5]);
      if isempty(parents)
        % only change axis behavior if no parents are specified
        axis off
      end
    end
    
    if ~isempty(Hz)
      if ~isempty(surfhandle) && update(3)
        varargin(sellist) = ft_setopt(varargin(sellist), 'surfhandle', surfhandle(3));
      end
      if ~isempty(patchhandle) && update(3)
        varargin(sellist) = ft_setopt(varargin(sellist), 'patchhandle', patchhandle(3));
      end
      varargin{2*sel} = ori(3,:);
      set(gcf,'currentaxes',Hz);
      hz = ft_plot_slice(dat, varargin{:}, 'tag', 'z');
      set(Hz, 'view', [0 90]); %, 'xlim', [0.5 size(dat,1)-0.5], 'ylim', [0.5 size(dat,2)-0.5]);
      if isempty(parents)
        % only change axis behavior if no parents are specified
        axis off
      end
    end
    
  case 'intersect'
    holdflag = ishold;
    if ~holdflag
      hold on
    end
    
    varargin{2*sel} = ori(1,:);
    hx = ft_plot_slice(dat, varargin{:}, 'tag', 'x');
    
    varargin{2*sel} = ori(2,:);
    hy = ft_plot_slice(dat, varargin{:}, 'tag', 'y');
    
    varargin{2*sel} = ori(3,:);
    hz = ft_plot_slice(dat, varargin{:}, 'tag', 'z');

    axis equal; axis tight; axis vis3d; axis off
    view(3);
    
    if ~holdflag
      hold off
    end

    if ~isempty(coordsys) && ~strcmp(coordsys, 'unknown')
      % add a context sensitive menu to change the 3d viewpoint to top|bottom|left|right|front|back
      menu_viewpoint(gca, coordsys)
    end
    
  otherwise
    ft_error('unsupported style %s', style);
    
end % switch style

% if strcmp(interactive, 'yes')
%   flag = 1;
%   while flag
%
%   end
% end
