function [data] = ft_determine_coordsys(data, varargin)

% FT_DETERMINE_COORDSYS plots a geometrical object, allowing you to perform

% a visual check on the coordinatesystem, the units and on the anatomical

% labels for the coordinate system axes.

%

% Use as

% [dataout] = ft_determine_coordsys(datain, ...)

% where the input data structure can be either

% - an anatomical MRI

% - a cortical or head surface mesh

% - an electrode, gradiometer or optode definition

% - a volume conduction model of the head

% or most other FieldTrip structures that represent geometrical information.

%

% Additional optional input arguments should be specified as key-value pairs

% and can include

% interactive = string, 'yes' or 'no' (default = 'yes')

% axisscale = scaling factor for the reference axes and sphere (default = 1)

% clim = lower and upper anatomical MRI limits (default = [0 1])

%

% This function will pop up a figure that allows you to check whether the

% alignment of the object relative to the coordinate system axes is correct

% and what the anatomical labels of the coordinate system axes are. You

% should switch on the 3D rotation option in the figure panel to rotate and

% see the figure from all angles. To change the anatomical labels of the

% coordinate system, you should press the corresponding keyboard button.

%

% Recognized and supported coordinate systems are 'ctf', 'bti', '4d', 'yokogawa',

% 'eeglab', 'neuromag', 'itab', 'acpc', 'spm', 'mni', 'fsaverage', 'tal', 'scanras',

% 'scanlps', 'dicom'.

%

% Furthermore, supported coordinate systems that do not specify the origin are 'ras',

% 'als', 'lps', etc. See https://www.fieldtriptoolbox.org/faq/coordsys for more

% details.

%

% See also FT_CONVERT_COORDSYS, FT_DETERMINE_UNITS, FT_CONVERT_UNITS, FT_PLOT_AXES, FT_PLOT_XXX

% Copyright (C) 2015-2021, Jan-Mathijs Schoffelen

%

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

dointeractive = ft_getopt(varargin, 'interactive', 'yes');

axisscale = ft_getopt(varargin, 'axisscale', 1); % this is used to scale the axmax and rbol

clim = ft_getopt(varargin, 'clim', [0 1]); % this is used to scale the orthoplot

fontsize = ft_getopt(varargin, 'fontsize'); % this is passed to ft_plot_axes

data = ft_checkdata(data, 'hasunit', 'yes');

dtype = ft_datatype(data);

% the high-level data structures are detected with ft_datatype, but there are

% also some low-level data structures that need to be supported here

if strcmp(dtype, 'unknown')

if isfield(data, 'fid') || (isfield(data, 'tri') && isfield(data, 'pos'))

dtype = 'headshape';

elseif isfield(data, 'hex') && isfield(data, 'pos')

dtype = 'mesh';

elseif isfield(data, 'tet') && isfield(data, 'pos')

dtype = 'mesh';

elseif ~strcmp(ft_headmodeltype(data), 'unknown')

dtype = 'headmodel';

elseif ~strcmp(ft_senstype(data), 'unknown')

dtype = 'sens';

end

elseif strcmp(dtype, 'volume+label')

% we don't care about the labels here

dtype = 'volume';

elseif strcmp(dtype, 'mesh+label')

% we don't care about the labels here

dtype = 'mesh';

end

if isfield(data, 'coordsys') && ~isempty(data.coordsys)

% ensure that it is in lower case

data.coordsys = lower(data.coordsys);

% print the interpretation of the coordinate system

[labelx, labely, labelz] = coordsys2label(data.coordsys, 2, 0);

fprintf('The positive x-axis is pointing towards %s\n', labelx);

fprintf('The positive y-axis is pointing towards %s\n', labely);

fprintf('The positive z-axis is pointing towards %s\n', labelz);

end

% plot the geometrical object

% the plotting style depends on the data content

figure;

switch dtype

case 'volume'

funparam = [];

if isfield(data, 'anatomy')

funparam = data.anatomy;

elseif isfield(data, 'gray')

funparam = data.gray;

elseif isfield(data, 'white')

funparam = data.white;

elseif isfield(data, 'brick0')

funparam = data.brick0; % used for an atlas

elseif isfield(data, 'brick1')

funparam = data.brick1; % used for an atlas

else

% try to determine it automatically

fn = fieldnames(data);

for i=1:length(fn)

if isequal(size(data.(fn{i})), data.dim)

funparam = data.(fn{i});

break;

end

end

end

if isempty(funparam)

ft_error('don''t know which volumetric parameter to plot');

end

% the volumetric data needs to be interpolated onto three orthogonal planes

% determine a resolution that is close to, or identical to the original resolution

[corner_vox, corner_head] = cornerpoints(data.dim, data.transform);

diagonal_head = norm(range(corner_head));

diagonal_vox = norm(range(corner_vox));

resolution = (diagonal_head+eps)/(diagonal_vox+eps); % this is in units of "data.unit"

% scale funparam between 0 and 1

if ~isa(funparam, 'double') % avoid integer datatypes to allow for scaling

funparam = double(funparam);

end

dmin = min(funparam(:));

dmax = max(funparam(:));

funparam = (funparam-dmin)./(dmax-dmin);

clear ft_plot_slice

ft_plot_ortho(funparam, 'transform', data.transform, 'unit', data.unit, 'resolution', resolution, 'style', 'intersect', 'clim', clim);

axis vis3d

view([110 36]);

case 'source'

if isfield(data, 'inside') && ~isfield(data, 'tri')

% only plot the source locations that are inside the volume conduction model

ft_plot_mesh(data.pos(data.inside, :));

else

ft_plot_mesh(data, 'edgecolor','none', 'facecolor', [0.6 0.8 0.6], 'facealpha', 0.6);

end

camlight;

case 'dip'

ft_plot_mesh(data, 'edgecolor','none', 'facecolor', 'none');

camlight;

case 'headshape'

ft_plot_headshape(data);

camlight;

case {'mesh', 'source+mesh'}

ft_plot_mesh(data);

camlight;

case 'headmodel'

ft_plot_headmodel(data);

camlight;

case {'grad' 'elec' 'opto' 'sens'}

ft_plot_sens(data, 'label', 'label');

camlight;

case {'raw', 'timelock', 'freq', 'mvar', 'freqmvar', 'comp'}

% the data may contain a gradiometer or electrode definition

if isfield(data, 'grad')

ft_plot_sens(data.grad);

elseif isfield(data, 'elec')

ft_plot_sens(data.elec, 'label', 'label');

elseif isfield(data, 'opto')

ft_plot_sens(data.opto, 'label', 'label');

end

case 'unknown'

end % switch dtype{k}

if isfield(data, 'tri')

% this makes the 3-D object easier to understand

camlight

lighting gouraud

end

% plot the 3-D axes, labels, and sphere at the origin

ft_plot_axes(data, 'axisscale', axisscale, 'fontsize', fontsize);

if istrue(dointeractive)

% ensure the figure is updated prior to asking the question

% this was needed for FT_ELECTRODEPLACEMENT in combination with MATLAB 2022a

drawnow

if ~isfield(data, 'coordsys') || isempty(data.coordsys)

% default is yes

fprintf('The coordinate system is not specified.\n')

value = smartinput('Do you want to change the anatomical labels for the axes [Y, n]? ', 'y');

else

% default is no

fprintf('The coordinate system is specified as "%s".\n', data.coordsys)

value = smartinput('Do you want to change the anatomical labels for the axes [y, N]? ', 'n');

end

if strcmpi(value, 'n')

return

end

% interactively determine orientation

orientation = ' ';

while ~any(strcmp(orientation(1), {'r', 'l', 'a', 'p', 's', 'i'}))

orientation(1) = smartinput('What is the anatomical label for the positive X-axis [r, l, a, p, s, i]? ', '');

end

while ~any(strcmp(orientation(2), {'r', 'l', 'a', 'p', 's', 'i'}))

orientation(2) = smartinput('What is the anatomical label for the positive Y-axis [r, l, a, p, s, i]? ', '');

end

while ~any(strcmp(orientation(3), {'r', 'l', 'a', 'p', 's', 'i'}))

orientation(3) = smartinput('What is the anatomical label for the positive Z-axis [r, l, a, p, s, i]? ', '');

end

% interactively determine origin

origin = ' ';

while ~any(strcmp(origin, {'a', 'i', 's', 'n'}))

origin = input('Is the origin of the coordinate system at the a(nterior commissure), i(nterauricular), s(scanner origin), n(ot a landmark)? ', 's');

end

% some coordinate systems are identical or very similar, see https://www.fieldtriptoolbox.org/faq/coordsys

if origin=='a' && strcmp(orientation, 'ras')

coordsys = 'acpc'; % also used for spm, mni, tal

elseif origin=='i' && strcmp(orientation, 'als')

coordsys = 'ctf'; % also used for 4d, bti, eeglab

elseif origin=='i' && strcmp(orientation, 'ras')

coordsys = 'neuromag'; % also used for itab

elseif origin=='s' && strcmp(orientation, 'ras')

coordsys = 'scanras'; % also used for nifti

elseif origin=='s' && strcmp(orientation, 'lps')

coordsys = 'scanlps'; % also used for dicom

else

% only use the orientation, not the origin

coordsys = orientation;

end

data.coordsys = coordsys;

end % if interactive