function [mri] = ft_read_mri(filename, varargin)

filename = fetch_url(filename);

dataformat = ft_getopt(varargin, 'dataformat');

outputfield = ft_getopt(varargin, 'outputfield', 'anatomy');

spmversion = ft_getopt(varargin, 'spmversion');

readbids = ft_getopt(varargin, 'readbids', 'ifmakessense');

if isempty(spmversion)

if ~ft_hastoolbox('spm')

spmversion = 'spm12';

else

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

end

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

if islogical(readbids)

% it should be either yes/no/ifmakessense

if readbids

readbids = 'yes';

else

readbids = 'no';

end

if ~exist(filename, 'file')

ft_error('file ''%s'' does not exist', filename);

if isempty(dataformat)

% only do the autodetection if the format was not specified

dataformat = ft_filetype(filename);

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

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;

inflated = false;

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;

if exist('seg', 'var')

% store the segmented data

mri.seg = seg;

if exist('hdr', 'var')

% store the header with all file format specific details

mri.hdr = hdr;

if exist('transform', 'var')

% store the homogeneous transformation matrix if present

mri.transform = transform;

if exist('unit', 'var')

% determine the geometrical units in which it is expressed

mri.unit = unit;

% estimate the units from the data

mri = ft_determine_units(mri);

if exist('coordsys', 'var')

% add a descriptive label for the coordinate system

mri.coordsys = coordsys;

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;

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;

if inflated

% compressed file has been unzipped on the fly, clean up

delete(filename);

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

filecontent = load(filename, varname);

value = filecontent.(varname); % read the variable named according to the input specification

clear filecontent

function bool = isrighthanded(orient)

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