function [object] = ft_convert_coordsys(object, target, varargin)

% FT_CONVERT_COORDSYS changes the coordinate system of the input object to the

% specified coordinate system. The coordinate system of the input object is

% determined from the 'coordsys' field in the input data, or needs to be determined

% and specified interactively by the user.

%

% Use as

% [output] = ft_convert_coordsys(input, target)

% [output] = ft_convert_coordsys(input, target, method)

% [output] = ft_convert_coordsys(input, target, method, template)

% to determine and convert the coordinate system.

%

% With the optional method input argument you can determine whether to use SPM for an

% affine or non-linear transformation.

% method = 0: only an approximate coregistration (default for non-MRI data)

% method = 1: an approximate coregistration, followed by spm_affreg

% method = 2: an approximate coregistration, followed by spm_normalise (default for MRI data)

%

% The following input data structures are supported

% electrode or gradiometer array, see FT_DATATYPE_SENS

% volume conduction model, see FT_DATATYPE_HEADMODEL

% source model, see FT_DATATYPE_SOURCE and FT_PREPARE_SOURCEMODEL

% anatomical mri, see FT_DATATYPE_VOLUME

% segmented mri, see FT_DATATYPE_SEGMENTATION

% anatomical or functional atlas, see FT_READ_ATLAS

%

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

%

% Note that the conversion will be an automatic and approximate conversion, not

% taking into account differences in individual anatomies/differences in conventions

% where to put the fiducials.

%

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

% Undocumented options

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

% Copyright (C) 2005-2023, Robert Oostenveld & 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$

if nargin>3 && isnumeric(varargin{1})

% old-style with 2 extra input arguments

tmp = {'method', varargin{1}, 'template', varargin{2}};

varargin = tmp;

elseif nargin>2 && isnumeric(varargin{1})

% old-style with 1 extra input argument

tmp = {'method', varargin{1}};

varargin = tmp;

end

method = ft_getopt(varargin, 'method'); % default is handled below

templatefile = ft_getopt(varargin, 'template'); % default is handled in the SPM section

feedback = ft_getopt(varargin, 'feedback', 'no');

spmversion = ft_getopt(varargin, 'spmversion');

% use the version that is on the path, or default to spm12

if isempty(spmversion)

if ~ft_hastoolbox('spm')

spmversion = 'spm12';

else

spmversion = lower(spm('ver'));

end

end

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

% determine the coordinate system of the input object

object = ft_determine_coordsys(object, 'interactive', 'yes');

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

% the call to ft_determine_coordsys should have taken care of this, but

% it is possible that the user aborted the coordinate system

% determination. See http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2526

ft_error('the coordinate system of the geometrical object is not specified');

end

ft_notice('the coordinate system is determined as ''%s''', object.coordsys);

end

if isequal(object.coordsys, target)

% no further conversion is needed

return

end

if isempty(method)

if isfield(object, 'transform') && isfield(object, 'anatomy')

% the default for an anatomical MRI is to start with an approximate alignment,

% followed by a call to spm_normalise for a better quality alignment

method = 2;

else

% the default for all other objects is to do only an approximate alignment

method = 0;

end

end

if ~ismember(object.coordsys, {'spm', 'mni', 'fsaverage', 'tal'}) && ismember(target, {'spm', 'mni', 'fsaverage', 'tal'})

% the input appears to be an individual subject MRI which has not been rescaled

% the target is a template coordinate system

% see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3304 and http://bit.ly/2sw7eC4

ft_warning('Not applying any scaling, using ''acpc'' instead of ''%s''. See http://bit.ly/2sw7eC4', target);

target = 'acpc';

end

if ~isfield(object, 'unit') || isempty(object.unit)

% determine the units of the input object

object = ft_determine_units(object);

ft_notice('the units are determined as ''%s''', object.unit);

end

% all of the internal logic inside this function requires that the units are in millimeter

% it will be converted back at the end of this function

original = object;

object = ft_convert_units(object, 'mm');

%--------------------------------------------------------------------------

% Do an initial affine registration (rigid body) alignment to the target coordinate

% system. This deals with 90-degree rotations between RAS and ALS, or between CTF and

% NEUROMAG, but also with approximate alignment between external NAS/LPA/RPA

% landmark-based coordinates to internal ACPCP landmark-based coordinate systems.

% these are the 48 generic axis orientation triplets, these specify the axes but no origin

% a = anterior

% p = posterior

% l = left

% r = right

% s = superior

% i = inferior

generic = {

'als'; 'ali'; 'ars'; 'ari';...

'pls'; 'pli'; 'prs'; 'pri';...

'las'; 'lai'; 'ras'; 'rai';...

'lps'; 'lpi'; 'rps'; 'rpi';...

'asl'; 'ail'; 'asr'; 'air';...

'psl'; 'pil'; 'psr'; 'pir';...

'sal'; 'ial'; 'sar'; 'iar';...

'spl'; 'ipl'; 'spr'; 'ipr';...

'sla'; 'ila'; 'sra'; 'ira';...

'slp'; 'ilp'; 'srp'; 'irp';...

'lsa'; 'lia'; 'rsa'; 'ria';...

'lsp'; 'lip'; 'rsp'; 'rip'}';

if ismember(object.coordsys, generic) && strcmp(target, 'acpc') && method>0

% converting anatomical MRI data from a generic coordinate system to ACPC with SPM should work

% do an initial alignment of the anatomical MRI data to RAS to bring it closer to ACPC

initial = ft_affinecoordinates(object.coordsys, 'ras');

else

initial = ft_affinecoordinates(object.coordsys, target);

end

object = ft_transform_geometry(initial, object);

object.coordsys = target;

%--------------------------------------------------------------------------

% Do a second round of affine registration (rigid body) to get improved

% alignment with ACPC coordinate system. This is needed because there may be

% different conventions defining LPA and RPA. The affine registration may

% fail however, e.g. if the initial alignment is not close enough. In that

% case SPM will throw an error.

%

% We expect the template MRIs to be expressed in millimeter and to be

% approximately aligned with ACPC.

if method>0

if isdeployed && isempty(templatefile)

ft_error('you need to specify a template filename for the coregistration');

end

if ~isfield(object, 'transform') || ~isfield(object, 'anatomy')

ft_error('affine or non-linear transformation are only supported for anatomical MRIs');

end

if ~isfield(object, 'unit') || ~strcmp(object.unit, 'mm')

ft_error('affine or non-linear transformation require the anatomial MRI to be expressed in mm');

end

% the following requires SPM to be on the path

ft_hastoolbox(spmversion, 1);

end

if method==1

% use spm_affreg

if isempty(templatefile)

switch lower(spm('ver'))

case 'spm2'

if isdeployed

if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using method==1'); end

else

templatefile = fullfile(spm('Dir'),'templates','T1.mnc');

end

case 'spm8'

if isdeployed

if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using method==1'); end

elseif isempty(templatefile)

templatefile = fullfile(spm('Dir'),'templates','T1.nii');

end

case 'spm12'

if isdeployed

if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using method==1'); end

else

templatefile = fullfile(spm('Dir'),'toolbox','OldNorm','T1.nii');

if ~exist('spm_affreg', 'file')

addpath(fullfile(spm('Dir'),'toolbox','OldNorm'));

end

end

fprintf('using ''OldNorm'' affine registration\n');

otherwise

ft_error('unsupported SPM version');

end

end

template = ft_read_mri(templatefile);

tname1 = [tempname, '.img'];

tname2 = [tempname, '.img'];

V1 = ft_write_mri(tname1, object.anatomy, 'transform', object.transform, 'spmversion', spm('ver'), 'dataformat', 'nifti_spm');

V2 = ft_write_mri(tname2, template.anatomy, 'transform', template.transform, 'spmversion', spm('ver'), 'dataformat', 'nifti_spm');

% the below, using just spm_affreg does not work robustly enough in some cases

flags.regtype = 'rigid';

[M, scale] = spm_affreg(V1,V2,flags);

% some juggling around with the transformation matrices

vox2head = object.transform;

vox2acpc2 = M \ V1.mat;

acpc2head2 = vox2head / vox2acpc2;

% update the transformation matrix

object.transform = vox2acpc2;

% delete the temporary files

delete(tname1); delete(strrep(tname1, 'img', 'hdr'));

delete(tname2); delete(strrep(tname2, 'img', 'hdr'));

elseif method==2

% use spm_normalise

if isempty(templatefile)

switch lower(spm('ver'))

case 'spm2'

if isdeployed

if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using method==2'); end

else

templatefile = fullfile(spm('Dir'),'templates','T1.mnc');

end

case 'spm8'

if isdeployed

if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using method==2'); end

else

templatefile = fullfile(spm('Dir'),'templates','T1.nii');

end

case 'spm12'

% this uses the 'OldNorm' functionality, so the path needs to be added, can only be done if non-deployed.

if isdeployed

if nargin<3, ft_error('you need to specify a template filename when in deployed mode and using method==2'); end

else

templatefile = fullfile(spm('Dir'),'toolbox','OldNorm','T1.nii');

if ~exist('spm_normalise', 'file')

addpath(fullfile(spm('Dir'),'toolbox','OldNorm'));

end

end

fprintf('using ''OldNorm'' normalisation\n');

otherwise

ft_error('unsupported SPM version');

end

end

template = ft_read_mri(templatefile);

tname1 = [tempname, '.img'];

tname2 = [tempname, '.img'];

V1 = ft_write_mri(tname1, object.anatomy, 'transform', object.transform, 'unit', object.unit, 'spmversion', spm('ver'), 'dataformat', 'nifti_spm');

V2 = ft_write_mri(tname2, template.anatomy, 'transform', template.transform, 'unit', template.unit, 'spmversion', spm('ver'), 'dataformat', 'nifti_spm');

flags.nits = 0; % set number of non-linear iterations to zero

flags.regtype = 'rigid';

params = spm_normalise(V2,V1,[],[],[],flags);

% some juggling around with the transformation matrices

vox2head = object.transform;

acpc2head2 = V1.mat*params.Affine/V2.mat;

vox2acpc2 = acpc2head2\vox2head;

% update the transformation matrix

object.transform = vox2acpc2;

% delete the temporary files

delete(tname1); delete(strrep(tname1, 'img', 'hdr'));

delete(tname2); delete(strrep(tname2, 'img', 'hdr'));

end

if istrue(feedback)

% give some graphical feedback

ft_determine_coordsys(object, 'interactive', 'no', 'fontsize', 15);

% also add the original axes

transform = object.transform * inv(original.transform);

ft_plot_axes([], 'transform', transform, 'unit', 'mm', 'coordsys', original.coordsys, 'fontsize', 15);

end

% all of the internal logic inside this function required that the units were in millimeter

% convert back to the original units

object = ft_convert_units(object, original.unit);