ft_read_mri.m

function [mri] = ft_read_mri(filename, varargin)

% FT_READ_MRI reads anatomical and functional MRI data from different file formats.
% The output data is structured in such a way that it is compatible with
% FT_DATATYPE_VOLUME.
%
% Use as
%   [mri] = ft_read_mri(filename, ...)
%
% Additional options should be specified in key-value pairs and can include
%   'dataformat'  = string specifying the file format, determining the low-level
%                   reading routine to be used. If no explicit format is given,
%                   it is determined automatically from the filename.
%   'volumes'     = vector with the volume indices to read from a 4D nifti (only for 'nifti_spm')
%   'outputfield' = string specifying the name of the field in the structure in which the
%                   numeric data is stored (default = 'anatomy')
%   'fixel2voxel' = string, the operation to apply to the fixels belonging to the
%                   same voxel, can be 'max', 'min', 'mean' (only for 'mrtrix_mif', default = 'max')
%   'indexfile'   = string, pointing to a fixel index file, if not present in the same directory
%                   as the functional data (only for 'mrtrix_mif')
%   'spmversion'  = string, version of SPM to be used (default = 'spm12')
%   'readbids'    = string, 'yes', no', or 'ifmakessense', whether to read information from
%                   the BIDS sidecar files (default = 'ifmakessense')
%
% The supported dataformats are
%   'afni_head'/'afni_brik'      uses AFNI code
%   'analyze_img'/'analyze_hdr'  uses SPM code
%   'analyze_old'                uses Darren Webber's code
%   'asa_mri'
%   'ctf_mri'
%   'ctf_mri4'
%   'ctf_svl'
%   'dicom'                      uses SPM code, same as dicom_spm
%   'dicom_spm'                  uses SPM code, also allows for multiframe dicoms (one volume in a file)     
%   'dicom_fs'                   uses FreeSurfer code, which relies on MATLAB image processing toolbox
%   'dicom_old'                  uses MATLAB image processing toolbox code
%   'freesurfer_mgh'             uses FreeSurfer code
%   'freesurfer_mgz'             uses FreeSurfer code
%   'jnifti_jnii'
%   'jnifti_bnii'
%   'matlab'                     assumes a MATLAB *.mat file containing a struct
%   'minc'                       uses SPM, this requires SPM5 or older
%   'mrtrix_mif'                 uses mrtrix code
%   'neuromag_fif'               uses MNE toolbox code
%   'neuromag_fif_old'           uses meg-pd toolbox
%   'nifti'                      uses FreeSurfer code
%   'nifti_spm'                  uses SPM code
%   'seg3d_mat'                  MATLAB file from Seg3D with a scirunnrrd structure
%   'yokogawa_mri'
%
% The following MRI file formats are supported
%   AFNI (*.head, *.brik)
%   ANT - Advanced Neuro Technology (*.mri)
%   Analyze (*.img, *.hdr)
%   CTF (*.svl, *.mri version 4 and 5)
%   DICOM (*.dcm, *.ima)
%   FreeSurfer (*.mgz, *.mgh)
%   MATLAB (*.mat)
%   MINC (*.mnc)
%   Mrtrix image format (*.mif)
%   NIFTi (*.nii) and zipped NIFTi (*.nii.gz)
%   Neuromag/Elekta/Megin (*.fif)
%   Yokogawa (*.mrk, incomplete)
%
% If you have a series of DICOM files, please provide the name of any of the files in
% the series (e.g. the first one). The files corresponding to the whole volume will
% be found automatically.
%
% The output MRI may have a homogenous transformation matrix that converts the
% coordinates of each voxel (in xgrid/ygrid/zgrid) into head coordinates.
%
% If the input file is a 4D nifti, and you wish to load in just a subset of the
% volumes, for example due to memory constraints, you should use as dataformat 'nifti_spm',
% which uses the optional key-value pair 'volumes' = vector, with the indices of the
% to-be-read volumes. The order of the indices is ignored, and the volumes will be
% sorted according to the original numeric indices, i.e., 1:10 yields the same as 10:-1:1
%
% See also FT_DATATYPE_VOLUME, FT_WRITE_MRI, FT_READ_DATA, FT_READ_HEADER, FT_READ_EVENT

% Copyright (C) 2008-2024, 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$

% optionally get the data from the URL and make a temporary local copy
filename = fetch_url(filename);

% get the options
dataformat  = ft_getopt(varargin, 'dataformat');
outputfield = ft_getopt(varargin, 'outputfield', 'anatomy');
spmversion  = ft_getopt(varargin, 'spmversion');
readbids    = ft_getopt(varargin, 'readbids', 'ifmakessense');

% 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

% the following is added for backward compatibility of using 'format' rather than 'dataformat'
format    = ft_getopt(varargin, 'format');
if ~isempty(format)
  ft_warning('the option ''format'' will be deprecated soon, please use ''dataformat'' instead');
  if isempty(dataformat)
    dataformat  = format;
  end
end

% for backward compatibility with https://github.com/fieldtrip/fieldtrip/issues/1585
if islogical(readbids)
  % it should be either yes/no/ifmakessense
  if readbids
    readbids = 'yes';
  else
    readbids = 'no';
  end
end

% test whether the file exists
if ~exist(filename, 'file')
  ft_error('file ''%s'' does not exist', filename);
end

if isempty(dataformat)
  % only do the autodetection if the format was not specified
  dataformat = ft_filetype(filename);
end

% deal with data that is organized according to BIDS
% this has to happen prior to the unzipping of the file
if strcmp(readbids, 'yes') || strcmp(readbids, 'ifmakessense')
  [p, f, x] = fileparts(filename);
  % check whether it is a BIDS dataset with a json sidecar file
  isbids = startsWith(f, 'sub-');
  if isbids
    % try to read the metadata from the BIDS sidecar files
    sidecar = bids_sidecar(filename);
    if ~isempty(sidecar)
      mri_json = ft_read_json(sidecar);
    end
  end
end

if strcmp(dataformat, 'compressed') || (strcmp(dataformat, 'freesurfer_mgz') && ispc) || any(filetype_check_extension(filename, {'gz', 'zip', 'tar', 'tgz'}))
  % the file is compressed, unzip on the fly,
  % -freesurfer mgz files get special treatment only on a pc
  % -compressed AFNI BRIKS also need the HEAD copied over to the temp dir

  filename_old = filename;
  filename    = inflate_file(filename_old);
  if strcmp(dataformat, 'freesurfer_mgz')
    filename_old = filename;
    filename     = [filename '.mgh'];
    movefile(filename_old, filename);
  end
  if ~strcmp(dataformat, 'nifti_spm')
    % replace it with the filetype's default format, but don't overwrite in
    % case dataformat was nifti_spm
    dataformat = ft_filetype(filename);
  end
  if strcmp(dataformat, 'afni_brik')
    [p, f, e] =fileparts(filename_old);
    copyfile(fullfile(p, strrep(f, 'BRIK', 'HEAD')), strrep(filename, 'BRIK', 'HEAD'));
  end
  inflated = true;
else
  inflated = false;
end

% test for the presence of some external functions from other toolboxes
hasspm2  = ft_hastoolbox('spm2');    % see http://www.fil.ion.ucl.ac.uk/spm/
hasspm5  = ft_hastoolbox('spm5');    % see http://www.fil.ion.ucl.ac.uk/spm/
hasspm8  = ft_hastoolbox('spm8');    % see http://www.fil.ion.ucl.ac.uk/spm/
hasspm12 = ft_hastoolbox('spm12');   % see http://www.fil.ion.ucl.ac.uk/spm/

switch dataformat
  case 'ctf_mri'
    [img, hdr] = read_ctf_mri(filename);
    transform = hdr.transformMRI2Head;
    coordsys  = 'ctf';

    if issubfield(hdr, 'fiducial.head')
      fid.label = fieldnames(hdr.fiducial.head);
      for i=1:length(fid.label)
        fid.pos(i,:) = hdr.fiducial.head.(fid.label{i});
      end
    elseif issubfield(hdr, 'fiducial.mri')
      fid.label = fieldnames(hdr.fiducial.mri);
      for i=1:length(fid.label)
        fid.pos(i,:) = hdr.fiducial.mri.(fid.label{i});
      end
      % convert from voxel to headcoordinates
      fid.pos = ft_warp_apply(transform, fid.pos);
    end

  case 'ctf_mri4'
    [img, hdr] = read_ctf_mri4(filename);
    transform = hdr.transformMRI2Head;
    coordsys  = 'ctf';

    if issubfield(hdr, 'fiducial.head')
      fid.label = fieldnames(hdr.fiducial.head);
      for i=1:length(fid.label)
        fid.pos(i,:) = hdr.fiducial.head.(fid.label{i});
      end
    elseif issubfield(hdr, 'fiducial.mri')
      fid.label = fieldnames(hdr.fiducial.mri);
      for i=1:length(fid.label)
        fid.pos(i,:) = hdr.fiducial.mri.(fid.label{i});
      end
      % convert from voxel to headcoordinates
      fid.pos = ft_warp_apply(transform, fid.pos);
    end

  case 'ctf_svl'
    [img, hdr] = read_ctf_svl(filename);
    transform = hdr.transform;

  case 'asa_mri'
    [img, seg, hdr] = read_asa_mri(filename);
    transform = hdr.transformMRI2Head;
    if isempty(seg)
      % in case seg exists, it will be added to the output
      clear seg
    end

  case 'minc'
    if ~(hasspm2 || hasspm5)
      fprintf('the SPM2 or SPM5 toolbox is required to read *.mnc files\n');
      ft_hastoolbox('spm2', 1);
    end
    % use the functions from SPM
    hdr = spm_vol_minc(filename);
    img = spm_read_vols(hdr);
    transform = hdr.mat;

  case 'nifti_spm'
    if ~(hasspm5 || hasspm8 || hasspm12)
      fprintf('the SPM5 or later toolbox is required to read *.nii files\n');
      ft_hastoolbox(spmversion, 1);
    end
    volumes = ft_getopt(varargin, 'volumes', []);

    % use the functions from SPM
    hdr = spm_vol_nifti(filename);
    if isempty(volumes)
      img = double(hdr.private.dat);
    else
      volumes = sort(intersect(volumes, 1:size(hdr.private.dat,4)));
      img = double(hdr.private.dat(:,:,:,volumes));
    end
    transform = hdr.mat;
    unit = 'mm';

    try
      % The nifti header allows three methods for specifying the coordinates of the
      % voxels. For two of them (sform and qform), the header contains a numerical
      % code from 0 to 4 that specifies the coordinate system. SPM8 and SPM12 use a
      % field in the private header of the nifti object to code this.
      switch hdr.private.mat_intent
        case 'UNKNOWN'
          coordsys = 'unknown';
        case 'Scanner'
          coordsys = 'scanras';
        case 'Aligned'
          coordsys = 'aligned';
        case 'Talairach'
          coordsys = 'tal';
        case 'MNI152'
          coordsys = 'mni152';
      end

      % We cannot trust it yet, see https://github.com/fieldtrip/fieldtrip/issues/1879
      ft_notice('the coordinate system appears to be ''%s''\n', coordsys);
      clear coordsys
    end


  case {'analyze_img' 'analyze_hdr'}
    if ~(hasspm8 || hasspm12)
      fprintf('the SPM8 or newer toolbox is required to read analyze files\n');
      ft_hastoolbox(spmversion, 1);
    end

    % use the image file instead of the header
    filename((end-2):end) = 'img';
    % use the functions from SPM to read the Analyze MRI
    hdr = spm_vol(filename);
    img = spm_read_vols(hdr);
    transform = hdr.mat;

  case 'analyze_old'
    % use the functions from Darren Weber's mri_toolbox to read the Analyze MRI
    ft_hastoolbox('mri', 1);     % from Darren Weber, see http://eeg.sourceforge.net/

    avw = avw_img_read(filename, 0); % returned volume is LAS*
    img = avw.img;
    hdr = avw.hdr;

    % The default Analyze orientation is axial unflipped (LAS*), which means
    % that the resulting volume is according to the radiological convention.
    % Most other fMRI and EEG/MEG software (except Mayo/Analyze) uses
    % neurological conventions and a right-handed coordinate system, hence
    % the first axis of the 3D volume (right-left) should be flipped to make
    % the coordinate system comparable to SPM
    ft_warning('flipping 1st dimension (L-R) to obtain volume in neurological convention');
    img = flip(img, 1);

    transform      = diag(hdr.dime.pixdim(2:4));
    transform(4,4) = 1;

  case {'afni_brik' 'afni_head'}
    % needs afni
    ft_hastoolbox('afni', 1);    % see http://afni.nimh.nih.gov/

    [err, hdr] = BrikInfo(filename);

    % check the precision of the data, and if scaling is required. If the precision is other than float,
    %and no scaling is required, then return the data in its native precision, let the low level code take
    % care of that
    if any(hdr.BRICK_FLOAT_FACS~=0)
      opts.OutPrecision = '';
    else
      opts.OutPrecision = '*';
      opts.Scale = 0;
    end
    [err, img, hdr, ErrMessage] = BrikLoad(filename, opts);
    if err
      ft_error('could not read AFNI file');
    end

    if isfield(hdr, 'ORIENT_SPECIFIC')
      [err, orient, flipvec] = AFNI_OrientCode(hdr.ORIENT_SPECIFIC);
      % FIXME, I don't understand why the orient vector needs to be like
      % this: it seems the opposite of what is reflected in the coordsys
      % (see below), but it seems to yield internally consistent results
    else
      % afni volume info
      orient = 'LPI'; % hope for the best
    end

    % origin and basis vectors in world space
    [unused, ix] = AFNI_Index2XYZcontinuous([0 0 0; eye(3)], hdr, orient);

    % basis vectors in voxel space
    e1 = ix(2,:) - ix(1,:);
    e2 = ix(3,:) - ix(1,:);
    e3 = ix(4,:) - ix(1,:);

    % change from base0 (afni) to base1 (SPM/Matlab)
    o = ix(1,:) - (e1+e2+e3);

    % create matrix
    transform = [e1;e2;e3;o]';
    transform = cat(1, transform, [0 0 0 1]);

    coordsys = lower(hdr.Orientation(:,2)');
    if contains(filename, 'TTatlas') || (isfield(hdr, 'TEMPLATE_SPACE') && ~isempty(hdr.TEMPLATE_SPACE))
      if isfield(hdr, 'TEMPLATE_SPACE') && ~isempty(hdr.TEMPLATE_SPACE)
        space = hdr.TEMPLATE_SPACE;
      else
        space = 'tal'; % accommodate the case when this is not specified in the hdr, make assumption
      end
      if startsWith(space, 'tt_') || startsWith(space, 'TT_')
        space = 'tal';
      elseif startsWith(space, 'mni')
        space = 'mni';
      elseif startsWith(space, 'tlrc')
        % according to the documentation tlrc is rather generic as a
        % specification of the space, but originally it meant tal.
        ft_warning('space ''tlrc'' might be ambiguous, here assuming the coordsys to be ''tal''');
        space = 'tal';
      end
      if ismember(space, {'tal' 'mni'})
        % xyz orientation should be RAS
        if ~strcmp(coordsys, 'ras')
          ft_warning('the template space suggests that the image is in %s coordinates, but the xyz orientation %s does not match this', space, coordsys);
          xxx2ras = true;
        else
          coordsys = space;
        end
      end
    end

  case 'neuromag_fif'
    % needs mne toolbox
    ft_hastoolbox('mne', 1);

    % use the mne functions to read the Neuromag MRI
    hdr = fiff_read_mri(filename);
    img_t = cat(3, hdr.slices.data);
    img = permute(img_t,[2 1 3]);
    hdr.slices = rmfield(hdr.slices, 'data'); % remove the image data to save memory

    % information below is from MNE - fiff_define_constants.m
    % coordinate system 4 - is the MEG head coordinate system (fiducials)
    % coordinate system 5 - is the MRI coordinate system
    % coordinate system 2001 - MRI voxel coordinates
    % coordinate system 2002 - Surface RAS coordinates (is mainly vertical shift, no rotation to 2001)
    % MEG sensor positions come in system 4
    % MRI comes in system 2001

    transform = eye(4);
    if isfield(hdr, 'trans') && issubfield(hdr.trans, 'trans')
      if (hdr.trans.from == 4) && (hdr.trans.to == 5)
        transform = hdr.trans.trans;
      else
        ft_warning('W: trans does not transform from 4 to 5.');
        ft_warning('W: Please check the MRI fif-file');
      end
    else
      ft_warning('W: trans structure is not defined.');
      ft_warning('W: Maybe coregistration is missing?');
    end
    if isfield(hdr, 'voxel_trans') && issubfield(hdr.voxel_trans, 'trans')
      % centers the coordinate system and switches from mm to m
      if (hdr.voxel_trans.from == 2001) && (hdr.voxel_trans.to == 5)
        % matlab_shift compensates for the different index conventions between C and MATLAB

        % the lines below is old code (prior to Jan 3, 2013) and only works with 1 mm resolution MRIs
        %   matlab_shift = [ 0 0 0 0.001; 0 0 0 -0.001; 0 0 0 0.001; 0 0 0 0];
        % transform transforms from 2001 to 5 and further to 4
        %   transform = transform\(hdr.voxel_trans.trans+matlab_shift);

        % the lines below should work with arbitrary resolution
        matlab_shift = eye(4);
        matlab_shift(1:3,4) = [-1,-1,-1];
        transform = transform\(hdr.voxel_trans.trans * matlab_shift);

        coordsys  = 'neuromag';
        mri.unit  = 'm';
      else
        ft_warning('W: voxel_trans does not transform from 2001 to 5.');
        ft_warning('W: Please check the MRI fif-file');
      end
    else
      ft_warning('W: voxel_trans structure is not defined.');
      ft_warning('W: Please check the MRI fif-file');
    end

  case 'neuromag_fif_old'
    % needs meg_pd toolbox
    ft_hastoolbox('meg-pd', 1);

    % use the meg_pd functions to read the Neuromag MRI
    [img,coords] = loadmri(filename);
    dev = loadtrans(filename,'MRI','HEAD');
    transform  = dev*coords;
    hdr.coords = coords;
    hdr.dev    = dev;

  case {'dicom' 'dicom_spm'}
    % requires spm12
    ft_hastoolbox('spm12', 1);
    dicomhdr = spm_dicom_headers(filename); % read the header for a single file
    if isfield(dicomhdr{1}, 'SharedFunctionalGroupSequence')
      % this seems to be a multiframe image, we don't need to read in the separate headers for each of the slices
    else
      % this is a 'traditional' image, one file per slice
      [f, p, e] = fileparts(filename);
      d         = dir(fullfile(f,sprintf('*%s',e)));
      filenames = {d.name};
      for k = 1:numel(filenames)
        filenames{k} = fullfile(f, filenames{k});
      end
      hdr1 = spm_dicom_headers(filename);  % read the header for the predefined slice
      hdrs = spm_dicom_headers(filenames); % read the headers for all slices in the folder 
      for k = 1:numel(hdrs)
        seriesnumber(k,1) = hdrs{k}.SeriesNumber;
      end
      sel      = seriesnumber==hdr1{1}.SeriesNumber;
      dicomhdr = hdrs(sel);
    end
    tdir  = tempdir;
    fname = spm_dicom_convert(dicomhdr, 'all', 'flat', 'nii', tdir, 1);
    
    % convert  the cell(-array) dicomhdr as generated by SPM into a
    % struct(-array)
    for k = 1:numel(dicomhdr)
      dicomhdr_struct(k) = dicomhdr{k};
    end
    clear dicomhdr;
    
    mri   = ft_read_mri(fname.files{1});
    mri.coordsys = 'scanras';
    mri.unit     = 'mm';
    mri.hdr      = dicomhdr_struct;
    delete(fname.files{1});
    
  case 'dicom_fs'
    % this returns a right-handed volume with the transformation matrix stored in the file headers
    % see https://github.com/fieldtrip/website/pull/444

    % needs the freesurfer toolbox
    ft_hastoolbox('freesurfer', 1);
    [dcmdir, junk1, junk2] = fileparts(filename);
    if isempty(dcmdir)
      dcmdir = '.';
    end
    [img, transform,hdr, mr_params] = load_dicom_series(dcmdir,dcmdir,filename);
    transform = vox2ras_0to1(transform);
    coordsys  = 'scanras';
    unit      = 'mm';

  case 'dicom_old'
    % this uses the Image processing toolbox
    % the DICOM files represent a stack of slices, and possibly even multiple volumes
    orig = dicominfo(filename);
    dim(1) = orig.Rows;
    dim(2) = orig.Columns;

    % this works for the Siemens scanners at the FCDC
    [p, f] = fileparts(filename);
    tok = tokenize(f, '.');
    for i=5:length(tok)
      tok{i} = '*';
    end
    filename = sprintf('%s.', tok{:});  % reconstruct the filename with wildcards and '.' between the segments
    filename = filename(1:end-1);       % remove the last '.'
    dirlist  = dir(fullfile(p, filename));
    dirlist  = {dirlist.name};

    if isempty(dirlist)
      % this is for the Philips data acquired at KI
      ft_warning('could not determine list of dicom files, trying with *.dcm');
      dirlist  = dir(fullfile(p, '*.dcm'));
      dirlist  = {dirlist.name};
    end

    if isempty(dirlist)
      ft_warning('could not determine list of dicom files, trying with *.ima');
      dirlist  = dir(fullfile(p, '*.ima'));
      dirlist  = {dirlist.name};
    end

    if length(dirlist)==1
      % try something else to get a list of all the slices
      dirlist = dir(fullfile(p, '*'));
      dirlist = {dirlist(~[dirlist.isdir]).name};
    end

    keep = false(1, length(dirlist));
    for i=1:length(dirlist)
      filename = char(fullfile(p, dirlist{i}));
      if ~strcmp(dataformat, 'dicom_old')
        keep(i) = false;
        fprintf('skipping ''%s'' because of incorrect filetype\n', filename);
      end
      % read the header information
      info     = dicominfo(filename);
      if info.SeriesNumber~=orig.SeriesNumber
        keep(i) = false;
        fprintf('skipping ''%s'' because of different SeriesNumber\n', filename);
      else
        keep(i) = true;
        hdr(i)  = info;
      end
    end
    % remove the files that were skipped
    hdr     = hdr(keep);
    dirlist = dirlist(keep);

    % pre-allocate enough space for the subsequent slices
    dim(3) = length(dirlist);
    img    = zeros(dim(1), dim(2), dim(3));
    for i=1:length(dirlist)
      filename = char(fullfile(p, dirlist{i}));
      ft_info('reading image data from ''%s''\n', filename);
      img(:,:,i) = dicomread(hdr(i));
    end

    % reorder and concatenate the slices
    [z, indx]   = sort(cell2mat({hdr.SliceLocation}));
    hdr = hdr(indx);
    img = img(:,:,indx);

    % construct a homgeneous transformation matrix that performs the scaling from voxels to mm
    transform = dicom2transform(hdr);
    coordsys  = 'dicom'; % identical to scanlps, see https://www.fieldtriptoolbox.org/faq/coordsys/#details-of-the-dicom-coordinate-system
    unit      = 'mm';

    % this makes the mapping of voxels to patient coordinates consistent with Horos
    img = permute(img, [2, 1, 3]);

  case {'nifti', 'nifti_gz', 'freesurfer_mgz', 'freesurfer_mgh'}
    ft_hastoolbox('freesurfer', 1);
    tmp = MRIread(filename);
    ndims = numel(size(tmp.vol));
    if ndims==3
      img = permute(tmp.vol, [2 1 3]);
      % FIXME although this is probably correct
      % see the help of MRIread, anecdotally columns and rows seem to need a swap
      % in order to match the transform matrix (alternatively a row switch of the
      % latter can be done)
    elseif ndims==4
      img = permute(tmp.vol, [2 1 3 4]);
    end
    hdr = rmfield(tmp, 'vol');
    transform = hdr.vox2ras1;
    unit = 'mm';

    if isfield(hdr, 'niftihdr')
      % The nifti header allows three methods for specifying the coordinates of the
      % voxels. For two of them (sform and qform), the header contains a numerical
      % code that specifies the coordinate system.
      coordsys = {'unknown', 'scanras', 'aligned', 'tal', 'mni152'}; % corresponding to 0, 1, 2, 3, 4
      if isequal(hdr.vox2ras0, hdr.niftihdr.sform)
        coordsys = coordsys{hdr.niftihdr.sform_code + 1};
      elseif isequal(hdr.vox2ras0, hdr.niftihdr.qform)
        coordsys = coordsys{hdr.niftihdr.qform_code + 1};
      else
        coordsys = 'unknown';
      end

      % We cannot trust it yet, see https://github.com/fieldtrip/fieldtrip/issues/1879
      ft_notice('the coordinate system appears to be ''%s''\n', coordsys);
      clear coordsys
    end

  case 'yokogawa_mri'
    ft_hastoolbox('yokogawa', 1);
    fid = fopen_or_error(filename, 'rb');
    mri_info = GetMeg160MriInfoM(fid);
    patient_info = GetMeg160PatientInfoFromMriFileM(fid);
    [data_style, model, marker, image_parameter, normalize, besa_fiducial_point] = GetMeg160MriFileHeaderInfoM(fid);
    fclose(fid);

    % gather all meta-information
    hdr.mri_info = mri_info;
    hdr.patient_info = patient_info;
    hdr.data_style = data_style;
    hdr.model = model;
    hdr.marker = marker;
    hdr.image_parameter = image_parameter;
    hdr.normalize = normalize;
    hdr.besa_fiducial_point = besa_fiducial_point;

    ft_error('FIXME yokogawa_mri implementation is incomplete');

  case 'seg3d_mat'
    scirunnrrd = loadvar(filename, 'scirunnrrd');
    mri.anatomy = scirunnrrd.data;
    mri.dim = [scirunnrrd.axis.size];
    % construct the homogenous transformation matrix
    mri.transform = eye(4);
    mri.transform(1,1) = scirunnrrd.axis(1).spacing;
    mri.transform(2,2) = scirunnrrd.axis(2).spacing;
    mri.transform(3,3) = scirunnrrd.axis(3).spacing;
    mri.transform(1,4) = scirunnrrd.axis(1).min - scirunnrrd.axis(1).spacing;
    mri.transform(2,4) = scirunnrrd.axis(2).min - scirunnrrd.axis(2).spacing;
    mri.transform(3,4) = scirunnrrd.axis(3).min - scirunnrrd.axis(3).spacing;

  case 'matlab'
    mri = loadvar(filename, 'mri');

  case {'mif' 'mrtrix_mif'}
    ft_hastoolbox('mrtrix', 1);
    tmp = read_mrtrix(filename);

    % check if it's sparse fixeldata
    isfixel = numel(tmp.dim==3) && tmp.dim(3)==1;

    if ~isfixel
      mri.hdr     = removefields(tmp, {'data'});
      mri.(outputfield) = tmp.data;
      mri.dim     = tmp.dim(1:length(size(tmp.data)));
      mri.transform = tmp.transform;
      mri.transform(1:3,1:3) = diag(tmp.vox(1:3))*mri.transform(1:3,1:3);
    else
      fix2vox_fun = ft_getopt(varargin, 'fixel2voxel', 'max');
      indexfile   = ft_getopt(varargin, 'indexfile');
      if isempty(indexfile)
        % assume the index file to be in the same directory as the data file
        [p,f,e]   = fileparts(filename);
        indexfile = fullfile(p,'index.mif');
      end
      index     = read_mrtrix(indexfile);
      tmpdata   = reshape(index.data, [], 2);

      vox_index = find(tmpdata(:,1)>0);
      num_index = tmpdata(vox_index,1);
      fix_index = tmpdata(vox_index,2)+1;

      % create a mapping matrix of fixel2voxel -> currently this only works
      % for scalar fixel data.
      tmpdata = zeros(numel(num_index), max(num_index));
      for k = 1:numel(num_index)
        tmpdata(k, 1:num_index(k)) = fix_index(k)-1 + (1:num_index(k));
      end
      tmpdata(tmpdata==0) = nan;
      tmpdata             = tmpdata.'; % transpose is intended
      tmpdata(isfinite(tmpdata)) = tmp.data(tmpdata(isfinite(tmpdata)));

      switch fix2vox_fun
        case 'magmax'
          tmpx    = nanmin(tmpdata,[],1).';
          tmpdata = nanmax(tmpdata,[],1).';
          tmpdata(abs(tmpx)>tmpdata) = tmpx(abs(tmpx)>tmpdata);
        case 'max'
          tmpdata = nanmax(tmpdata,[],1).';
        case 'min'
          tmpdata = nanmin(tmpdata,[],1).';
        case 'mean'
          tmpdata = nanmean(tmpdata,1).';
        case 'none'
          tmpdata = tmpdata.';
        otherwise
          ft_error('unsupported fixel2voxel operation requested');
      end

      mri.hdr  = removefields(tmp, {'data'});
      mri.(outputfield) = zeros([prod(index.dim(1:3)) size(tmpdata,2)]);
      mri.(outputfield)(vox_index,:) = tmpdata;
      mri.(outputfield)   = reshape(mri.(outputfield), [index.dim(1:3) size(tmpdata,2)]);
      mri.dim       = size(mri.(outputfield)) ;%index.dim(1:length(size(tmp.data)));
      mri.transform = tmp.transform;
      mri.transform(1:3,1:3) = diag(tmp.vox(1:3))*mri.transform(1:3,1:3);
    end

  case {'neurojson_jnii' 'neurojson_bnii'}
    % this depends on two external toolboxes
    ft_hastoolbox('jsonlab', 1);
    ft_hastoolbox('jnifti', 1);

    jnii = loadjnifti(filename);

    mri.hdr     = jnii.NIFTIHeader;
    mri.(outputfield) = jnii.NIFTIData;
    mri.dim     = jnii.NIFTIHeader.Dim;
    mri.unit    = jnii.NIFTIHeader.Unit.L; % units of length

    % see https://brainder.org/2012/09/23/the-nifti-file-format/
    coordsys = {'scanras', 'aligned', 'tal', 'mni152'};

    if jnii.NIFTIHeader.SForm>0
      mri.coordsys  = coordsys{jnii.NIFTIHeader.SForm};
      mri.transform = [
        jnii.NIFTIHeader.Affine
        0 0 0 1
        ];

    elseif jnii.NIFTIHeader.QForm>0
      mri.coordsys = coordsys{jnii.NIFTIHeader.QForm};

      % this is adapted from freesurfer/load_nifti_hdr.m
      b = jnii.NIFTIHeader.Quatern.b;
      c = jnii.NIFTIHeader.Quatern.c;
      d = jnii.NIFTIHeader.Quatern.d;
      x = jnii.NIFTIHeader.QuaternOffset.x;
      y = jnii.NIFTIHeader.QuaternOffset.y;
      z = jnii.NIFTIHeader.QuaternOffset.z;

      a = 1.0 - (b*b + c*c + d*d);
      if(abs(a) < 1.0e-7)
        a = 1.0 / sqrt(b*b + c*c + d*d);
        b = b*a;
        c = c*a;
        d = d*a;
        a = 0.0;
      else
        a = sqrt(a);
      end

      r11 = a*a + b*b - c*c - d*d;
      r12 = 2.0*b*c - 2.0*a*d;
      r13 = 2.0*b*d + 2.0*a*c;
      r21 = 2.0*b*c + 2.0*a*d;
      r22 = a*a + c*c - b*b - d*d;
      r23 = 2.0*c*d - 2.0*a*b;
      r31 = 2.0*b*d - 2*a*c;
      r32 = 2.0*c*d + 2*a*b;
      r33 = a*a + d*d - c*c - b*b;

      if(jnii.NIFTIHeader.VoxelSize(1) < 0.0)
        r13 = -r13;
        r23 = -r23;
        r33 = -r33;
      end

      R = [r11 r12 r13; r21 r22 r23; r31 r32 r33];
      S = diag(jnii.NIFTIHeader.VoxelSize(2:4));
      T = [x y z]';
      mri.transform = [R*S T; 0 0 0 1];

    else
      mri.coordsys = 'unknown';
      % Method 1
      mri.transform      = diag(jnii.NIFTIHeader.VoxelSize(2:4));
      mri.transform(4,4) = 1;
    end

    % ensure that this is double precision and not uint8
    mri.transform = double(mri.transform);

    % ensure that this is double precision and not uint8
    mri.dim = double(mri.dim);

    % these are already part of the output structure and should not be reassigned
    clear coordsys transform unit

  otherwise
    ft_error('unrecognized filetype ''%s'' for ''%s''', dataformat, filename);
end % switch dataformat

if exist('img', 'var')
  % determine the size of the volume in voxels
  mri.dim = size(img);
  if numel(mri.dim)>3
    mri.dim = mri.dim(1:3); % dim should be the spatial 3D dim
  end
  % store the anatomical data
  mri.(outputfield) = img;
end

if exist('seg', 'var')
  % store the segmented data
  mri.seg = seg;
end

if exist('hdr', 'var')
  % store the header with all file format specific details
  mri.hdr = hdr;
end

if exist('transform', 'var')
  % store the homogeneous transformation matrix if present
  mri.transform = transform;
end

if exist('unit', 'var')
  % determine the geometrical units in which it is expressed
  mri.unit = unit;
else
  % estimate the units from the data
  mri = ft_determine_units(mri);
end

if exist('coordsys', 'var')
  % add a descriptive label for the coordinate system
  mri.coordsys = coordsys;
end

if exist('xxx2ras', 'var') && xxx2ras==true
  % this is needed for AFNI formatted data, where the created voxels-to-world
  % mapping matrix is diagonal for the 3x3 rotation part (i.e. ijk should
  % be ras, in order for the tal/mni coordsys to make sense
  mri = ft_convert_coordsys(mri, 'ras', 0);
  mri.coordsys = space;
end

if (strcmp(readbids, 'yes') || strcmp(readbids, 'ifmakessense')) && isbids
  % the BIDS sidecar files extend/overrule the information that is present in the file header itself
  if exist('mri_json', 'var')
    fid =[];
    if isfield(mri_json, 'AnatomicalLandmarkCoordinates')
      fid.label = fieldnames(mri_json.AnatomicalLandmarkCoordinates);
      for i=1:length(fid.label)
        fid.pos(i,:) = mri_json.AnatomicalLandmarkCoordinates.(fid.label{i});
      end
      if isfield(mri_json, 'AnatomicalLandmarkCoordinateSystem') && isfield(mri_json, 'AnatomicalLandmarkCoordinateUnits')
        fid.coordsys = mri_json.AnatomicalLandmarkCoordinateSystem;
        fid.unit     = mri_json.AnatomicalLandmarkCoordinateUnits;
      else
        % assume that it is according to https://bids-specification.readthedocs.io/en/stable/glossary.html#anatomicallandmarkcoordinates-sense-2-metadata
        fid.pos = fid.pos+1;                             % convert from 0 to 1 offset
        fid.pos = ft_warp_apply(mri.transform, fid.pos); % convert from voxel to head coordinates
      end
    end % if fiducials present
  end % if mri_json
end % if readbids

if exist('fid', 'var') && ~isempty(fid)
  % store the fiducial details
  mri.fid = fid;
end

if inflated
  % compressed file has been unzipped on the fly, clean up
  delete(filename);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function value = loadvar(filename, varname)
var = whos('-file', filename);
if length(var)==1
  filecontent = load(filename); % read the one variable in the file, regardless of how it is called
  value       = filecontent.(var.name);
  clear filecontent
else
  filecontent = load(filename, varname);
  value       = filecontent.(varname);  % read the variable named according to the input specification
  clear filecontent
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bool = isrighthanded(orient)
bool = ismember(orient, {'ALS' 'RAS' 'PRS' 'LPS' 'SAL' 'SRA' 'SPR' 'SLP'});