function ds=cosmo_fmri_dataset(filename, varargin)
% load an fmri volumetric dataset
%
% ds = cosmo_fmri_dataset(filename, [,'mask',mask],...
% ['targets',targets],...
% ['chunks',chunks])
%
% Inputs:
% filename One of:
% * filename of fMRI dataset, it should end with one of:
% .nii, .nii.gz NIfTI
% .hdr, .img ANALYZE
% +{orig,tlrc}.{HEAD,BRIK[.gz]} AFNI
% .vmr, .vmp, .vtc, .glm, .msk BrainVoyager
% .mat SPM (SPM.mat)
% .mat:beta SPM beta
% .mat:con SPM contrast
% .mat:spm SPM stats
% .mat Matlab file with
% CoSMoMVPA or PyMVPA
% dataset.
% * xff structure (from neuroelf's xff)
% * nifti structure (from load_untouch_nii)
% * FieldTrip source MEEG structure
% * SPM structure
% * CoSMoMVPA fMRI or MEEG source dataset structure
% * PyMVPA fMRI dataset structure (exported using PyMVPA's
% cosmo.map2cosmo function)
% 'mask', m Any input as for filename (in which case the output must
% contain a single volume), or one of:
% '-all' exclude features where all values are
% zero or NaN
% '-any' exclude features where any value is
% zero or NaN
% '-auto' require that '-all' and '-any' exclude the
% same features, and exclude the
% corresponding features
% true equivalent to '-auto'
% false do not apply a mask
% The mask must have voxels at the same coordinates as the
% data indicated by filename, although it may
% have a different orientation (e.g. RAI, LPI, AIR).
% Only voxels that are non-zero and not NaN are selected
% from the data indicated by filename.
% If 'mask' is not given, then no mask is applied and a
% warning message (suggesting to use a mask) is printed if
% at least 5% of the values are non{zero,finite}.
% 'targets', t optional Tx1 numeric labels of experimental
% conditions (where T is the number of samples (volumes)
% in the dataset)
% 'chunks, c optional Tx1 numeric labels of chunks, typically indices
% of runs of data acquisition
% 'volumes', v optional vector with indices of volumes to load. If
% empty or not provided, then all volumes are loaded.
% 'block_size', b optional block size by which data is read (if
% supported by the format; currently NIfTI, ANALYZE, AFNI
% and SPM. If this option is provided *and* a mask is
% provided, then data is loaded in chunks (subsets of
% volumes) that contain at most block_size elements each;
% only data that survives the mask is then selected before
% the next block is loaded.
% The default value is 20,000,000, corresponding to ~160
% megabytes of memory required for a block (using
% numbers with double (64 bit) precsision).
% The rationale for this option is to reduce memory
% requirements, at the expensive of a possible increase of
% duration of disk reading operations.
%
%
% Returns:
% ds dataset struct with the following fields:
% .samples NxM matrix containing the data loaded from filename, for
% N samples (observations, volumes) and M features (spatial
% locations, voxels).
% If the original nifti file contained data with X,Y,Z,T
% elements in the three spatial and one temporal dimension
% and no mask was applied, then .samples will have
% dimensions N x M, where N = T and M = X*Y*Z. If a mask
% was applied then M (M<=X*Y*Z) is the number of non-zero
% voxels in the mask input dataset.
% .a struct with dataset-relevent data.
% .a.fdim.labels dimension labels, set to {'i','j','k'}
% .a.fdim.values dimension values, set to {1:X, 1:Y, 1:Z}
% .a.vol.dim 1x3 vector indicating the number of voxels in the 3
% spatial dimensions.
% .a.vol.mat 4x4 voxel-to-world transformation matrix (base-1).
% .a.vol.dim 1x3 number of voxels in each spatial dimension
% .sa struct for holding sample attributes (e.g., sa.targets,
% sa.chunks)
% .fa struct for holding feature attributes
% .fa.{i,j,k} indices of voxels (in voxel space).
%
% Notes:
% - Most MVPA applications require that .sa.targets (experimental
% condition of each sample) and .sa.chunks (partitioning of the samples
% in independent sets) are set, either by using this function or
% manually afterwards.
% - Data can be mapped to the volume using cosmo_map2fmri
% - SPM data can also be specified as filename:format, where format
% can be 'beta', 'con' or 'spm' (e.g. 'SPM.mat:beta', 'SPM.mat:con', or
% 'SPM.mat:spm') to load beta, contrast, or statistic images,
% respectively. When using 'beta', estimates for motion parameters and
% intercepts (which in most cases are estimates of no interest) are
% not returned. If format is omitted it is set to 'beta'.
% - If SPM data contains a field .Sess (session) then .sa.chunks are set
% according to its contents
% - If a mask is supplied, then all features that are in the mask are
% returned, even if some voxels contain NaN. To remove such features,
% consider applying cosmo_remove_useless_data to the output of this
% function.
%
% Dependencies:
% - for NIfTI, analyze (.hdr/.img) and SPM.mat files, it requires the
% NIfTI toolbox by Jimmy Shen
% (note that his toolbox is included in CoSMoMVPA in /externals)
% - for Brainvoyager files (.vmp, .vtc, .msk, .glm), it requires the
% NeuroElf toolbox, available from: http://neuroelf.net
% - for AFNI files (+{orig,tlrc}.{HEAD,BRIK[.gz]}) it requires the AFNI
% Matlab toolbox, available from: https://github.com/afni/AFNI
%
% Examples:
% % load nifti file
% ds=cosmo_fmri_dataset('mydata.nii');
%
% % load gzipped nifti file
% ds=cosmo_fmri_dataset('mydata.nii.gz');
%
% % load ANALYZE file and apply brain mask
% ds=cosmo_fmri_dataset('mydata.hdr','mask','brain_mask.hdr');
%
% % load AFNI file with 6 'bricks' (values per voxel, e.g. beta
% % values); set chunks (e.g. runs) and targets (experimental
% % conditions); use a mask
% ds=cosmo_fmri_dataset('mydata+tlrc', 'chunks', [1 1 1 2 2 2]', ...
% 'targets', [1 2 3 1 2 3]', ...
% 'mask', 'masks/brain_mask+tlrc);
%
% % load BrainVoyager VMR file in directory 'mydata', and apply an
% % automask that removes all features (voxels) that are zero or
% % non-finite for all samples
% ds=cosmo_fmri_dataset('mydata/mydata.vmr', 'mask', true);
%
% % load two datasets, one for odd runs, the other for even runs, and
% % combine them into one dataset. Note that the chunks are set here,
% % but the targets are not - for subsequent analyses this may have to
% % be done manually
% ds_even=cosmo_fmri_dataset('data_even_runs.glm','chunks',1);
% ds_odd=cosmo_fmri_dataset('data_odd_runs.glm','chunks',2);
% ds=cosmo_stack({ds_even,ds_odd});
%
% % load beta values from SPM GLM analysis stored
% % in a file SPM.mat.
% % If SPM.mat contains a field .Sess (sessions) then .sa.chunks
% % is set according to the contents of .Sess.
% ds=cosmo_fmri_dataset('path/to/SPM.mat');
%
% % as above, and apply an automask to remove voxels that
% % are zero or non-finite in all samples.
% ds=cosmo_fmri_dataset('path/to/SPM.mat','mask',true);
%
% % load contrast beta values from SPM GLM file SPM.mat
% ds=cosmo_fmri_dataset('path/to/SPM.mat:con');
%
% % load contrast statistic values from SPM GLM file SPM.mat
% ds=cosmo_fmri_dataset('path/to/SPM.mat:spm');
%
% % load PyMVPA dataset 'pymvpa_ds' that was exported using
% % PyMVPA's cosmo.map2cosmo function in Python with:
% %
% % >>> from mvpa2.datasets import cosmo
% % >>> cosmo.map2cosmo(pymvpa_ds,'data.mat')
% %
% cosmo_ds=cosmo_fmri_dataset('data.mat')
%
% See also: cosmo_map2fmri
%
% part of the NIfTI code is based on code by Robert W Cox, 2003,
% dedicated to the public domain.
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
% set defaults
defaults.mask=[];
defaults.targets=[];
defaults.chunks=[];
defaults.volumes=[];
defaults.block_size=2e7;
% set parameters
params = cosmo_structjoin(defaults, varargin{:});
% create dataset from filename
ds=convert_to_dataset(filename, params);
% set chunks and targets
ds=set_sa_vec(ds,params,'targets');
ds=set_sa_vec(ds,params,'chunks');
if param_specifies_auto_mask(params)
auto_mask=compute_auto_mask(ds.samples, params.mask);
ds=cosmo_slice(ds,auto_mask,2);
else
warn_if_nan_present_in_dataset(ds)
end
function ds_all=convert_to_dataset(fn, params)
% main data reader function
if string_endswith(fn,'.mat')
data=fast_import_data(fn);
else
data=fn;
end
img_formats_collection=get_img_formats;
label=find_img_format(data, img_formats_collection);
img_format=img_formats_collection.(label);
% make sure the required externals exists
externals=img_format.externals;
cosmo_check_external(externals);
% get the helper functions for this format
data_converter=img_format.data_converter;
has_data_reader=isfield(img_format,'data_reader');
% verify the input
input_is_filename=~img_format.matcher(data);
if input_is_filename && ~ischar(data)
error('illegal input of type ''%s''', class(data));
end
if input_is_filename
% read header from file
raw_header=img_format.header_reader(data);
elseif isfield(img_format, 'struct_header_reader')
raw_header=img_format.struct_header_reader(fn,data);
else
% input is a struct or object, no reading required
raw_header=data;
end
% get dataset header
[ds_hdr_full,nsamples]=img_format.header_converter(raw_header,...
params);
% now we have enough information about the dataset to load the mask
% and compare its orientation and position of voxels. If there is
% a mismatch, an error is raised
mask_ds=get_user_ds_mask(params, ds_hdr_full);
% only read in multiple blocks if it is supported for this format, and
% if there is a mask; otherwise what is the point?
has_mask=~isempty(mask_ds);
load_multiple_blocks=has_data_reader && has_mask;
nvoxels=prod(ds_hdr_full.a.vol.dim(1:3));
volumes_cell=partition_volumes(nsamples, nvoxels,...
load_multiple_blocks, params);
n_blocks=numel(volumes_cell);
% allocate space for output
ds_cell=cell(n_blocks,1);
% each block contains a subset of volumes; when they are combined
% the contain all volumes that have to be loaded.
% To reduce memory usage when a mask is supplied *and* the file
% format reader allows selecting a subset of volumes:
% - for each block, read all data from the corresponding volumes
% - apply the mask to the data
% - only store the result (which is much smaller)
% After this has been done for each block, results are stacked to
% form a single dataset
for block=1:n_blocks
volumes=volumes_cell{block};
% if all volumes are read, call the data reader or converter
% with empty input, so that it knows that all volumes can be read
volumes_or_empty=volumes;
is_all_volumes=isequal(volumes, 1:nsamples);
if is_all_volumes
volumes_or_empty=[];
end
if input_is_filename && has_data_reader
% selecting subset of volumes is done by the data_reader,
% so that only part of the whole file has to be read
data_reader=img_format.data_reader;
% make a copy of the filename or SPM struct
is_first_block=block==1;
if is_first_block
original_data=data;
end
data=data_converter(data_reader(original_data, raw_header, ...
volumes_or_empty), []);
else
% all data is probably already in memory, so select subset
% of volumes through the converter
data=data_converter(raw_header, volumes_or_empty);
end
if isfield(img_format,'convert_volume') && ...
~img_format.convert_volume
ds=data;
else
ds=flatten_data_array(data, ds_hdr_full.a.vol);
end
clear data; % reduce memory usage
if has_mask
if block==1
% only get the mask once, as it is the same for all blocks
ds_ids_mask=get_binary_dataset_mask(mask_ds,ds);
end
ds=cosmo_slice(ds, ds_ids_mask, 2);
end
ds_hdr=ds_hdr_full;
% get sample attributes for these volumes from the header
if isfield(ds_hdr,'sa')
ds_hdr.sa=cosmo_slice(ds_hdr.sa, volumes,1,'struct');
end
% update from header
ds=cosmo_structjoin(ds_hdr, ds);
% store block
ds_cell{block}=ds;
end
ds_all=cosmo_stack(ds_cell,[],[],false);
cosmo_check_dataset(ds_all,'fmri');
function ids_mask=get_binary_dataset_mask(ds_mask, ds)
% return a binary mask that can be used to slice ds
% to contain only features indexed by ds_mask
%
% This function also works if ds_mask and ds do not have the same
% features (voxels), and/or if the same location is indexed by multiple
% features.
% this should always be fine
check_datasets_in_same_space(ds_mask, ds);
assert(islogical(ds_mask.samples));
[mask_lin_ids,dim]=get_linear_feature_ids(ds_mask);
[lin_ids, dim_]=get_linear_feature_ids(ds);
assert(isequal(dim,dim_));
% allow duplicate feature ids in either ds_mask and/or ds
n_ids_mask=prod(dim);
mask=false(1,n_ids_mask);
mask(mask_lin_ids)=samples_to_binary_mask(ds_mask.samples);
ids_mask=mask(lin_ids);
function mask=samples_to_binary_mask(samples)
mask=samples~=0 & ~isnan(samples);
function warn_if_nan_present_in_dataset(ds)
nsamples=size(ds.samples,1);
for k=1:nsamples
if any(isnan(ds.samples(k,:)))
cosmo_warning(['The input dataset has NaN (not a number) '...
'values, which may cause the output '...
'of subsequent analyses to contain NaNs as '...
'well. For many use cases, NaNs are not '...
'desirabe. To remove features (voxels) '...
'with NaN values, consider using:\n\n'...
' ds_clean=cosmo_remove_useless_data(ds)'...
'\n\nwhere ds is the output from this '...
'function (%s)'],mfilename());
return
end
end
function [lin_ids, vol_dim]=get_linear_feature_ids(ds)
% get linear ids for each feature
keys={'i','j','k'};
n_keys=numel(keys);
sub_ids=cell(1,n_keys);
vol_dim=ds.a.vol.dim(1:3);
for k=1:n_keys;
key=keys{k};
[dim, unused, unused, unused, values]=cosmo_dim_find(ds,key,true);
if dim~=2
error('Unexpected key ''%s'' in sample dimension',key);
end
assert(numel(values)==vol_dim(k));
sub_ids{k}=values(ds.fa.(key));
end
lin_ids=sub2ind(vol_dim,sub_ids{:});
function tf=param_specifies_auto_mask(params)
% return true if an auto mask is specified
mask_param=params.mask;
tf=(islogical(mask_param) && ...
mask_param) || ...
(ischar(mask_param) && ...
numel(mask_param)>0 && ...
mask_param(1)=='-');
function tf=param_specifies_user_mask(params)
% return true if a user mask is specified
mask_param=params.mask;
tf=~(isempty(mask_param) || ...
islogical(mask_param) || ...
param_specifies_auto_mask(params));
function ds_mask=get_user_ds_mask(params, ds_hdr)
% returns a dataset containing the user mask,
% or [] if no user mask was supplied
mask_param=params.mask;
if ~param_specifies_user_mask(params);
ds_mask=[];
return;
end
% get the mask in dataset form
opt=struct();
opt.mask=false;
ds_mask=convert_to_dataset(mask_param,opt);
% only support single volume
nsamples_mask=size(ds_mask.samples,1);
if nsamples_mask~=1
error('mask must have a single volume, found %d',...
nsamples_mask);
end
% ensure they are in the same space
ds_mask=align_mask_to_ds_space(ds_mask, ds_hdr);
% set the samples to a boolean array
ds_mask.samples=samples_to_binary_mask(ds_mask.samples);
function ds_mask=align_mask_to_ds_space(ds_mask, ds_hdr)
% throw an error if the mask is in a different space as ds_hdr
% based on the header, make a minimal dataset (with no samples)
% so that its orientation can be obtained for ds_hdr
vol=ds_hdr.a.vol;
data=zeros([vol.dim(1:3) 0]);
ds_hdr=flatten_data_array(data, vol);
% if in another orientation, try to match to orientation of the mask
ds_orient=cosmo_fmri_orientation(ds_hdr);
if ~isequal(ds_orient, cosmo_fmri_orientation(ds_mask))
ds_mask=cosmo_fmri_reorient(ds_mask, ds_orient);
end
check_datasets_in_same_space(ds_mask, ds_hdr);
function check_datasets_in_same_space(ds_mask, ds_hdr)
% ensure the mask is compatible with the dataset
if ~isequal(ds_mask.a.fdim,ds_hdr.a.fdim)
error('.a.fdim mismatch between data and mask');
end
% check voxel-to-world mapping
max_delta=1e-4; % allow for minor tolerance
delta=max(abs(ds_mask.a.vol.mat(:)-ds_hdr.a.vol.mat(:)));
if delta>max_delta
error(['voxel dimension mismatch between data and mask:'...
'max difference is %.5f > %.5f'],...
delta,max_delta);
end
function ds_data=flatten_data_array(data, vol)
if all(cosmo_isfield(data,{'samples','a.vol.dim'}))
ds_data=data;
return;
end
% see how many dimensions there are, and their size
data_size = size(data);
ndim = numel(data_size);
if ndim<4
% simple reshape operation
data=reshape(data,[1 data_size]);
elseif ndim==4
data=shiftdim(data,3);
else
error('need 3 or 4 dimensions, found %d', ndim);
end
% number of values in 3 spatial dimensions
full_size=[size(data) 1 1];
ni=full_size(2);
nj=full_size(3);
nk=full_size(4);
% make a dataset
ds_data=cosmo_flatten(data,{'i';'j';'k'},{1:ni;1:nj;1:nk});
ds_data.a.vol=vol;
function volumes_cell=partition_volumes(n_volumes_total, n_voxels, ...
do_partial, params)
% return a cell with indices of volumes to load as specified
% by params.volumes.
% If params.volumes is empty, all volume indices are returned.
% If do_partial=true then each element in volume_cell has a
% limited number of indices so that the volumes in each cell
% element correspond to at most params.block_size elements (number
% of voxels times number of volumes).
% If do_partial=false, a cell with a single element with all volumes
% is returned.
if ~isfield(params,'volumes') || isempty(params.volumes)
volumes=1:n_volumes_total;
else
volumes=params.volumes;
end
if do_partial
block_size=params.block_size;
n_volumes_per_block=floor(block_size / n_voxels);
if n_volumes_per_block<1
n_volumes_per_block=1;
end
n_volumes_to_load=numel(volumes);
n_blocks=ceil(n_volumes_to_load / n_volumes_per_block);
if n_blocks<1
n_blocks=n_volumes_to_load;
n_volumes_per_block=1;
end
volumes_cell=cell(1, n_blocks);
first_index=1;
for block=1:n_blocks
last_index=min(n_volumes_to_load, ...
first_index+n_volumes_per_block-1);
volumes_cell{block}=volumes(first_index:last_index);
first_index=last_index+1;
end
assert(all(cellfun(@numel,volumes_cell)>=1));
else
volumes_cell={volumes};
end
function result=fast_import_data(fn)
x=load(fn);
keys=fieldnames(x);
if isfield(x,'samples') && isnumeric(x.samples)
result=x;
elseif numel(keys)==1
result=x.(keys{1});
else
error('Cannot load .mat file %s with multiple variables: %s',...
fn, cosmo_strjoin(keys,', '));
end
function ds=set_sa_vec(ds,p,fieldname)
% helper: sets a sample attribute as a vector
% throws an error if it has the wrong size
v=p.(fieldname);
if isequal(size(v),[0 0])
% ignore '[]', but not zeros(10,0)
return;
end
nsamples=size(ds.samples,1);
n=numel(v);
if n==1
% singleton element - repeat nsamples times.
v=repmat(v,nsamples,1);
n=nsamples;
end
if ~(n==0 || n==nsamples)
error('size mismatch for %s: expected %d values, found %d', ...
fieldname, nsamples, n);
end
ds.sa.(fieldname)=v(:);
function tf=string_endswith(s, tail)
tf=ischar(s) && ~isempty(s) && isempty(cosmo_strsplit(s, tail, -1));
function img_format=find_img_format(filename, img_formats)
% helper: find image format of filename fn
fns=fieldnames(img_formats);
n=numel(fns);
for k=1:n
fn=fns{k};
if ischar(filename)
exts=img_formats.(fn).exts;
m=numel(exts);
for j=1:m
ext=exts{j};
if string_endswith(filename,ext)
img_format=fn;
return
end
end
else
% it could be a struct - try that
matcher=img_formats.(fn).matcher;
if matcher(filename)
img_format=fn;
return
end
end
end
if ischar(filename)
desc=sprintf('file ''%s''',filename);
else
desc=sprintf('<%s> input',class(filename));
end
error('Could not find image format for %s',desc);
function auto_mask=compute_auto_mask(data, mask_type)
% mask_type can be 'any', 'all', 'auto', or ''
% When using 'auto', 'any' and 'all' should give the same mask
% When using '', a warning is shown when the percentage of
% non{zero,finite} features exceeds pct_thrshold
if isequal(mask_type,true)
mask_type='-auto';
end
pct_threshold=5;
to_remove=~samples_to_binary_mask(data);
% take as a mask anywhere where any feature is nonzero.
if cosmo_match({mask_type},{'-any','-auto',''})
to_remove_any=any(to_remove,1);
end
if cosmo_match({mask_type},{'-all','-auto',''})
to_remove_all=all(to_remove,1);
end
switch mask_type
case {'-auto',''}
%
any_equals_all=isequal(to_remove_any, to_remove_all);
n=numel(to_remove_any);
n_any=sum(to_remove_any(:));
n_all=sum(to_remove_all(:));
pct_any=100*n_any/n;
pct_all=100*n_all/n;
do_mask_suggestion=pct_all>pct_threshold && ...
strcmp(mask_type,'');
if any_equals_all
if do_mask_suggestion
me_name=mfilename();
msg=sprintf(['%d (%.1f%%) features are non'...
'{zero,finite} in all samples (and no '...
'features have non-{zero,finite} '...
'values in some samples and not '...
'in others)\n'...
'To use a mask excluding these '...
'features: %s(...,''mask'',-auto'')\n'],...
n_all,pct_all,me_name);
cosmo_warning(msg);
to_remove=[];
else
to_remove=to_remove_any;
end
else
% give error or warning
me_name=mfilename();
msg=sprintf(['%d (%.1f%%) features are non{zero,'...
'finite} in all samples\n'...
'%d (%.1f%%) features are non{zero,'...
'finite} in at least one sample\n'...
'To use a mask excluding '...
'features:\n'...
'- where *all* values are non{zero,finite}:'...
' %s(...,''mask'',-all'')\n'...
'- where *any* value is non{zero,finite}:'...
' %s(...,''mask'',-any'')\n'],...
n_all,pct_all,n_any,pct_any,me_name,me_name);
if strcmp(mask_type,'-auto');
error('automatic mask failed:\n%s',msg);
else
% give a warning
cosmo_warning(msg);
% set mask to empty; a mask will not be applied
to_remove=[];
end
end
case '-any'
to_remove=to_remove_any;
case '-all'
to_remove=to_remove_all;
otherwise
error('illegal mask specification ''-%s''', mask_type);
end
auto_mask=~to_remove(:)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Definition of supported data formats
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function img_formats=get_img_formats()
% define which formats are supported.
%
% The intention is to support a wide variety of formats through
% different external toolboxes, yet support these in a uniform way. If
% the toolbox supports loading a subset of the data
% (NIfTI, AFNI, SPM), then a .data_reader can be defined which
% only loads the requested data; otherwise (BrainVoyager, non-SPM
% .mat files) all data is read and the requested data selected
% afterwards
%
% In the definitions below:
% - native header: contains at least minimal header information in
% the format's representation
% - minimal dataset: contains at least minimal header information in
% dataset format (in .sa.samples, .a.vol.dim and
% .a.vol.mat)
% - nvolumes number of volumes available in the input
% - native header } native header information and data stored in
% and data: } native representation
% - volumes indices of volumes to read
% - 4D data volume data in X x Y x Z x T
%
% Each format is defined by these fields:
% .exts file name extensions
% .externals required externals
% .header_reader file name -> native header
% [.struct_header_reader] (file name, native header) -> native header
% (used for SPM.mat, to avoid re-loading the file)
% .header_converter (native header, params) -> (minimal dataset,
% nvolumes)
% [.data_reader] (file name,
% native header, volumes) -> (native header and
% data)
% .data_converter (native header and data,
% volumes) -> 4D data array
% [.convert_volume] if false, then the native header format must be
% an fmri dataset
%
% Notes:
% - if .data_reader is absent, then .header_reader must return a
% header that contains the data and which can be converter by
% .data_converter.
% - if .convert_volume is not present, then the output frmo
% .data_converter is assumed to return a COSMoMVPA dataset struct
% (instead of 4D data)
img_formats=struct();
img_formats.nii.exts={'.nii.gz','.nii','.hdr','.img'};
img_formats.nii.externals={'nifti'};
img_formats.nii.matcher=@isa_nii;
img_formats.nii.header_reader=@read_nii_header;
img_formats.nii.header_converter=@convert_nii_header;
img_formats.nii.data_reader=@read_nii_data;
img_formats.nii.data_converter=@convert_nii_data;
img_formats.bv_glm.exts={'.glm'};
img_formats.bv_glm.externals={'neuroelf'};
img_formats.bv_glm.matcher=@isa_bv_glm;
img_formats.bv_glm.header_reader=@read_bv_glm_header;
img_formats.bv_glm.header_converter=@convert_bv_glm_header;
img_formats.bv_glm.data_converter=@convert_bv_glm_data;
img_formats.bv_msk.exts={'.msk'};
img_formats.bv_msk.externals={'neuroelf'};
img_formats.bv_msk.matcher=@isa_bv_msk;
img_formats.bv_msk.header_reader=@read_bv_msk_header;
img_formats.bv_msk.header_converter=@convert_bv_msk_header;
img_formats.bv_msk.data_converter=@convert_bv_msk_data;
img_formats.bv_vtc.exts={'.vtc'};
img_formats.bv_vtc.externals={'neuroelf'};
img_formats.bv_vtc.matcher=@isa_bv_vtc;
img_formats.bv_vtc.header_reader=@read_bv_vtc_header;
img_formats.bv_vtc.header_converter=@convert_bv_vtc_header;
img_formats.bv_vtc.data_converter=@convert_bv_vtc_data;
img_formats.bv_vmp.exts={'.vmp'};
img_formats.bv_vmp.externals={'neuroelf'};
img_formats.bv_vmp.matcher=@isa_bv_vmp;
img_formats.bv_vmp.header_reader=@read_bv_vmp_header;
img_formats.bv_vmp.header_converter=@convert_bv_vmp_header;
img_formats.bv_vmp.data_converter=@convert_bv_vmp_data;
img_formats.bv_vmr.exts={'.vmr'};
img_formats.bv_vmr.externals={'neuroelf'};
img_formats.bv_vmr.matcher=@isa_bv_vmr;
img_formats.bv_vmr.header_reader=@read_bv_vmr_header;
img_formats.bv_vmr.header_converter=@convert_bv_vmr_header;
img_formats.bv_vmr.data_converter=@convert_bv_vmr_data;
img_formats.spm.exts={'mat:con','mat:beta','mat:spm'};
img_formats.spm.externals=img_formats.nii.externals;
img_formats.spm.matcher=@isa_spm;
img_formats.spm.struct_header_reader=@read_spm_struct_header;
img_formats.spm.header_reader=@read_spm_header;
img_formats.spm.header_converter=@convert_spm_header;
img_formats.spm.data_reader=@read_spm_data;
img_formats.spm.data_converter=@convert_spm_data;
img_formats.afni.exts={'+orig','+orig.HEAD','+orig.BRIK',...
'+orig.BRIK.gz','+tlrc','+tlrc.HEAD',...
'+tlrc.BRIK','+tlrc.BRIK.gz'};
img_formats.afni.externals={'afni'};
img_formats.afni.matcher=@isa_afni;
img_formats.afni.header_reader=@read_afni_header;
img_formats.afni.header_converter=@convert_afni_header;
img_formats.afni.data_reader=@read_afni_data;
img_formats.afni.data_converter=@convert_afni_data;
img_formats.ft_source.exts=cell(0);
img_formats.ft_source.externals=cell(0);
img_formats.ft_source.matcher=@isa_ft_source;
img_formats.ft_source.header_reader=@read_ft_source_header;
img_formats.ft_source.header_converter=@convert_ft_source_header;
img_formats.ft_source.data_converter=@convert_ft_source_data;
img_formats.ft_source.convert_volume=false; % already dataset
% fMRI volumetric datasets exported with PyMVPA's cosmo module
img_formats.pymvpa_fmri_ds.exts=cell(0);
img_formats.pymvpa_fmri_ds.matcher=@isa_pymvpa_fmri;
img_formats.pymvpa_fmri_ds.externals=cell(0);
img_formats.pymvpa_fmri_ds.header_reader=@read_pymvpa_ds_header;
img_formats.pymvpa_fmri_ds.header_converter=@convert_pymvpa_ds_header;
img_formats.pymvpa_fmri_ds.data_converter=@convert_pymvpa_ds_data;
img_formats.pymvpa_fmri_ds.convert_volume=false; % already dataset
img_formats.cosmo_fmri_ds.exts=cell(0);
img_formats.cosmo_fmri_ds.matcher=@(x)isa_cosmo_fmri(x) && ...
~isa_pymvpa_fmri(x);
img_formats.cosmo_fmri_ds.externals=cell(0);
img_formats.cosmo_fmri_ds.header_reader=@read_cosmo_ds_header;
img_formats.cosmo_fmri_ds.header_converter=@convert_cosmo_ds_header;
img_formats.cosmo_fmri_ds.data_converter=@convert_cosmo_ds_data;
img_formats.cosmo_fmri_ds.convert_volume=false; % already dataset
%%%%%%%%%%%%%%%%%%%%%%%%
% General
function data=slice_4d(data, volumes)
data_size=size(data);
ndim=numel(data_size);
if ndim>4
% Could be the AFNI NIFTI conversion syndrome, where the 4th
% dimension is singleton and the fifth one contains the data.
% Such data is accepted and treated as if the fifth dimension is
% the fourth one.
time_size=data_size(4:end);
if sum(time_size>1)>1
error(['More than one singleton dimension found in '...
'time dimension; this is currently not '...
'supported. If you want to be able '...
'to load such data, please get in touch '...
'with the CoSMoMVPA developers']);
end
ntime=prod(time_size);
data=reshape(data,[data_size(1:3),ntime]);
end
if ~isempty(volumes)
data=data(:,:,:,volumes);
end
if ~isa(data,'double')
data=double(data);
end
function ds=slice_dataset_volumes(ds, volumes)
if ~isempty(volumes)
ds=cosmo_slice(ds,volumes,1);
end
function require_singleton_volume(volumes)
if ~isequal(volumes,[]) && ~isequal(volumes,1)
error('Only a single volume is supported for this data format');
end
function hdr=get_and_check_data(hdr, loader_func, check_func, varargin)
% is hdr is a char, load it using loader with optional arguments
% from varargin
% For other input, the input is returned.
%
% in any case the output is checked using check_func.
if ischar(hdr)
hdr=loader_func(hdr, varargin{:});
end
if ~check_func(hdr)
error('Illegal input of type %s - failed to pass %s',...
class(hdr), func2str(check_func));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% format-specific helper functions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
% NIfTI
%%%%%%%%%%%%%%%%%%%%%%%%
% helpers
function [mx, xform]=get_nifti_transform(hdr, varargin)
% Get LPI affine transformation from NIfTI file
%
% Input:
% fn nifti filename
%
% Output:
% mx 4x4 affine transformation matrix from voxel to world
% coordinates. voxel indices as base0, not base1 as in
% CoSMoMVPA
% xform string with xform based on sform or qform code
%
% Notes:
% - this function is experimental
% - initial testing suggests agreement with MRIcron (thanks to Chris
% Rorden for providing this software)
% - functionality in the subfunctions are based on nftii1_io.h in
% AFNI, written Robert W Cox (2003), public domain dedication;
% http://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1_io.c
% - to convert voxel coordinates (i,j,k) to (x,y,z), compute
% lpi_mx*[x y z 1]'. For the reverse, compute inv(lpi_mx)*[i j k 1]'
%
% NNO Dec 2014
defaults.nifti_form=[];
opt=cosmo_structjoin(defaults,varargin);
if isempty(opt.nifti_form)
[mx, nifti_form]=nifti_matrix_from_auto(hdr);
else
[mx, nifti_form]=nifti_matrix_using_method(hdr,opt.nifti_form);
end
% prioritize sformch
switch nifti_form
case 'sform'
key=hdr.hist.sform_code;
case 'qform'
key=hdr.hist.qform_code;
otherwise
key=0;
end
xform=cosmo_fmri_convert_xform('nii',key);
function [mx, nifti_form]=nifti_matrix_from_auto(hdr)
% get matrix automatically, assumes that qform and sform (if present)
% are identical
max_delta_s_and_q=1e-3; % maximum allowed difference between s and q
has_sform=hdr.hist.sform_code>0;
has_qform=hdr.hist.qform_code>0;
if has_sform;
nifti_form='sform';
elseif has_qform
nifti_form='qform';
else
nifti_form='pixdim';
end
mx=nifti_matrix_using_method(hdr,nifti_form);
if has_sform && has_qform
mx_q=nifti_matrix_using_method(hdr,'qform');
max_diff=max(abs(mx(:)-mx_q(:)));
if max_diff>max_delta_s_and_q
str_mx=matrix2string(mx);
str_mx_s=matrix2string(mx_q);
url=['http://nifti.nimh.nih.gov/nifti-1/documentation/'...
'nifti1fields/nifti1fields_pages/qsform.html'];
error(['the affine matrices mapping voxel-to-world '...
'coordinates according to the sform and qform '...
'in the NIfTI header differ '...
'by %d, exceeding the treshold %d.\n\n'...
'The sform matrix is:\n\n%s\n\n',...
'The qform matrix is:\n\n%s\n\n'...
'To resolve this, set the ''nifti_form'' '...
'option to either:\n'...
' ''pixdim'' (method 1), or\n'...
' ''qform'' (method 2), or\n'...
' ''sform'' (method 3).\n\n'...
'For more information, see:\n %s\n\n',...
'If you have absolutely no idea what to use, '...
'try ''sform''\n\n'],...
max_diff, max_delta_s_and_q,...
str_mx_s, str_mx, url);
end
end
function s=matrix2string(mx)
float_pat='%6.3f';
line_pat=repmat([float_pat ' '],1,4);
line_pat(end)=sprintf('\n');
mx_pat=repmat(line_pat,1,3);
mx_3x4=mx(1:3,1:4);
s=sprintf(mx_pat,mx_3x4');
function [mx, method]=nifti_matrix_using_method(hdr, method)
switch method
case 'pixdim'
mx=nifti_matrix_from_pixdim(hdr);
case 'qform'
mx=nifti_matrix_from_qform(hdr);
case 'sform'
mx=nifti_matrix_from_sform(hdr);
otherwise
error('illegal method %s', method);
end
assert(isequal(size(mx),[4 4]));
function mx=nifti_matrix_from_pixdim(hdr)
mx=[[diag(hdr.dime.pixdim(2:4)); 0 0 0] [0;0;0;1]];
function mx=nifti_matrix_from_qform(hdr)
% convert quaternion to affine matrix
%
% based on "quatern_to_mat44" in nifti1_io.c by Robert W. Cox, 2003,
% which he dedicated to the public domain.
%
% http://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1_io.c
qfac=hdr.dime.pixdim(1);
dx=hdr.dime.pixdim(2);
dy=hdr.dime.pixdim(3);
dz=hdr.dime.pixdim(4);
qb=hdr.hist.quatern_b;
qc=hdr.hist.quatern_c;
qd=hdr.hist.quatern_d;
qx=hdr.hist.qoffset_x;
qy=hdr.hist.qoffset_y;
qz=hdr.hist.qoffset_z;
% ported from MRIcron
b=qb;
c=qc;
d=qd;
a=1-(b^2+c^2+d^2);
if a<1e-7
a=1/sqrt(b^2+c^2+d^2);
b=b*a;
c=c*a;
d=d*a;
a=0;
else
a=sqrt(a);
end
if dx>0
xd=dx;
else
xd=1;
end
if dy>0
yd=dy;
else
yd=1;
end
if dz>0
zd=dz;
else
zd=1;
end
if qfac<0
zd=-zd;
end
% construct affine matrix
mx=zeros(4,4);
mx(1,1)= (a*a+b*b-c*c-d*d) * xd ;
mx(1,2)= 2 * (b*c-a*d ) * yd ;
mx(1,3)= 2 * (b*d+a*c ) * zd ;
mx(2,1)= 2 * (b*c+a*d ) * xd ;
mx(2,2)= (a*a+c*c-b*b-d*d) * yd ;
mx(2,3)= 2 * (c*d-a*b ) * zd ;
mx(3,1)= 2 * (b*d-a*c ) * xd ;
mx(3,2)= 2 * (c*d+a*b ) * yd ;
mx(3,3)= (a*a+d*d-c*c-b*b) * zd ;
mx(1,4)=qx;
mx(2,4)=qy;
mx(3,4)=qz;
mx(4,:)=[0 0 0 1];
function mx=nifti_matrix_from_sform(hdr)
% set the srow values
mx=[hdr.hist.srow_x;...
hdr.hist.srow_y;...
hdr.hist.srow_z;...
[0 0 0 1]];
function scaling=nifti_get_scaling_factor(hdr)
% get scaling factor, if present for this dataset
% scaling=[intercept slope] if present, otherwise []
is_datatype=any(hdr.dime.datatype==[2,4,8,16,64,256,512,768]);
is_nonidentity=hdr.dime.scl_inter~=0 || hdr.dime.scl_slope~=1;
is_nonzero=hdr.dime.scl_slope~=0;
if is_datatype && is_nonidentity && is_nonzero
scaling=[hdr.dime.scl_inter hdr.dime.scl_slope];
else
scaling=[];
end
% NIfTI input
% -----------
function b=isa_nii(hdr)
b=isstruct(hdr) && isfield(hdr,'img') && isnumeric(hdr.img) && ...
isfield(hdr,'hdr') && isfield(hdr.hdr,'dime') && ...
isfield(hdr.hdr.dime,'dim') && isnumeric(hdr.hdr.dime.dim);
function nii=read_nii_header(fn)
nii.hdr=load_untouch_header_only(fn);
function [hdr_ds,nsamples]=convert_nii_header(nii, params)
hdr=nii.hdr;
dim=hdr.dime.dim(2:4);
% get original affine matrix
[mat, xform]=get_nifti_transform(hdr, params);
% make matrix base1 friendly
mat(1:3,4)=mat(1:3,4)+mat(1:3,1:3)*[-1 -1 -1]';
vol=struct();
vol.mat=mat;
vol.xform=xform;
vol.dim=dim;
sa=struct();
nsamples=hdr.dime.dim(5);
hdr_ds=struct();
hdr_ds.a.vol=vol;
hdr_ds.sa=sa;
function nii=read_nii_data(fn, hdr_from_nii, volumes)
% get original header
hdr=hdr_from_nii.hdr;
% load the data
nii=load_untouch_nii(fn, volumes);
assert(isequal(hdr.dime.dim(1:4),nii.hdr.dime.dim(1:4)));
hdr.dime.dim=nii.hdr.dime.dim;
hdr.dime.glmax=nii.hdr.dime.glmax;
hdr.dime.glmin=nii.hdr.dime.glmin;
assert(isequal(hdr,nii.hdr));
function data=convert_nii_data(nii, volumes)
% get scaling factor
scaling=nifti_get_scaling_factor(nii.hdr);
data=slice_4d(nii.img, volumes);
% apply scaling
if ~isempty(scaling)
data(:)=scaling(1)+scaling(2)*data;
end
%%%%%%%%%%%%%%%%%%%%%%%%
% AFNI
%%%%%%%%%%%%%%%%%%%%%%%%
% helpers
function vol=get_vol_afni(hdr)
% afni volume info
orient='LPI'; % always return LPI-based matrix
% origin and basis vectors in world space
k=[0 0 0;eye(3)];
[unused,i]=AFNI_Index2XYZcontinuous(k,hdr,orient);
% basis vectors in voxel space
e1=i(2,:)-i(1,:);
e2=i(3,:)-i(1,:);
e3=i(4,:)-i(1,:);
% change from base0 (afni) to base1 (SPM/Matlab)
o=i(1,:)-(e1+e2+e3);
% create matrix
mat=[e1;e2;e3;o]';
% set 4th row
mat(4,:)=[0 0 0 1];
vol=struct();
vol.mat=mat;
vol.dim=hdr.DATASET_DIMENSIONS(1:3);
vol.xform=cosmo_fmri_convert_xform('afni',hdr.SCENE_DATA(1));
% AFNI HEAD & BRIK input
% ----------------------
function b=isa_afni(hdr)
b=isstruct(hdr) && isfield(hdr,'DATASET_DIMENSIONS') && ...
isfield(hdr,'DATASET_RANK');
function head=read_afni_header(fn)
[err,head]=BrikInfo(fn);
if err
error('Could not read %s', fn);
end
function [hdr_ds,nsamples]=convert_afni_header(head, params)
if numel(head.DATASET_RANK)<2
error('illegal AFNI header: DATASET_RANK');
end
nsamples=head.DATASET_RANK(2);
% set sample attributes
sa=struct();
if isfield(head,'BRICK_LABS') && ~isempty(head.BRICK_LABS);
% if present, get labels
labels=cosmo_strsplit(head.BRICK_LABS,'~');
if numel(labels)==nsamples+1 && isempty(labels{end})
labels=labels(1:(end-1));
end
sa.labels=labels(:);
end
if isfield(head,'BRICK_STATAUX') && ~isempty(head.BRICK_STATAUX);
% if present, get stat codes
sa.stats=cosmo_statcode(head);
end
hdr_ds.sa=sa;
% get volume info
hdr_ds.a.vol=get_vol_afni(head);
function head=read_afni_data(fn, head, volumes)
opt=struct();
opt.Frames=volumes;
[err,data,head_data,err_msg]=BrikLoad(fn,opt);
if err
error('Error reading afni file: %s', err_msg);
end
assert(isequal(head_data,head));
head.img=data;
function data=convert_afni_data(head, volumes)
data=slice_4d(head.img, volumes);
%%%%%%%%%%%%%%%%%%%%%%%%
% BrainVoyager
%%%%%%%%%%%%%%%%%%%%%%%%
% helpers
function z=xff_struct(x)
% helper function: applies xff and returns a struct
% this avoids clearing of the object (which xff seems like doing)
y=xff(x);
z=getcont(y);
% BoundingBox is a method; copy its output to the struct
z.BoundingBox=y.BoundingBox;
% Clear object
y.ClearObject();
function vol=get_vol_bv(hdr)
% bv vol info
bbox=hdr.BoundingBox;
mat=neuroelf_bvcoordconv_wrapper([],'bvx2tal',bbox);
dim=bbox.DimXYZ;
vol=struct();
vol.mat=mat;
vol.dim=dim;
vol.xform=cosmo_fmri_convert_xform('bv',NaN);
function mat=neuroelf_bvcoordconv_wrapper(varargin)
% converts BV bounding box to affine transformation matrix
% helper functions that deals with both new neuroelf (version 1.0)
% and older versions.
% the old version provides a 'bvcoordconv' .m file
% the new version privides this function in the neuroelf class
has_bvcoordconv=~isempty(which('bvcoordconv'));
% set function handle
if has_bvcoordconv
f=@bvcoordconv;
else
n=neuroelf();
f=@n.bvcoordconv;
end
mat=f(varargin{:});
% BV GLM input
% ------------
function b=isa_bv_glm(hdr)
b=(isa(hdr,'xff') || isstruct(hdr)) && isfield(hdr,'Predictor') && ...
isfield(hdr,'GLMData') && isfield(hdr,'DesignMatrix');
function hdr=read_bv_glm_header(fn)
hdr=get_and_check_data(fn, @xff_struct, @isa_bv_glm);
function [hdr_ds,nsamples]=convert_bv_glm_header(hdr, params)
nsamples=hdr.NrOfPredictors;
% get sample attributes
name1=cell(nsamples,1);
name2=cell(nsamples,1);
rgb=zeros(nsamples,3);
for k=1:nsamples
p=hdr.Predictor(k);
name1{k}=p.Name1;
name2{k}=p.Name2;
rgb(k,:)=p.RGB(1,:);
end
sa=struct();
sa.Name1=name1;
sa.Name2=name2;
sa.RGB=rgb;
vol=get_vol_bv(hdr);
hdr_ds=struct();
hdr_ds.sa=sa;
hdr_ds.a.vol=vol;
function data=convert_bv_glm_data(hdr, volumes)
% ignore filename, because all data is already in hdr
data=slice_4d(hdr.GLMData.BetaMaps, volumes);
% BV volumetric map input
% -----------------------
function b=isa_bv_vmp(hdr)
b=(isa(hdr,'xff') || isstruct(hdr)) && ...
isfield(hdr,'Map') && isstruct(hdr.Map) && ...
isfield(hdr,'VMRDimX') && isfield(hdr,'NrOfMaps');
function hdr=read_bv_vmp_header(fn)
hdr=get_and_check_data(fn, @xff_struct, @isa_bv_vmp);
function [hdr_ds,nsamples]=convert_bv_vmp_header(hdr, params)
nsamples=hdr.NrOfMaps;
sa=struct();
labels=cell(nsamples,1);
for k=1:nsamples
map=hdr.Map(k);
labels{k}=map.Name;
end
sa.labels=labels;
sa.stats=cosmo_statcode(hdr);
hdr_ds=struct();
hdr_ds.sa=sa;
hdr_ds.a.vol=get_vol_bv(hdr);
function data=convert_bv_vmp_data(hdr, volumes)
nsamples=numel(hdr.Map);
if isempty(volumes)
volumes=1:nsamples;
end
nvolumes=numel(volumes);
data_cell=cell(nvolumes,1);
for k=1:nvolumes
map=hdr.Map(volumes(k)).VMPData;
data_cell{k}=map;
end
data=cat(4,data_cell{:});
% BV mask input
% -------------
function b=isa_bv_msk(hdr)
b=(isa(hdr,'xff') || isstruct(hdr)) && isfield(hdr, 'Mask');
function hdr=read_bv_msk_header(fn)
hdr=get_and_check_data(fn, @xff_struct, @isa_bv_msk);
function [hdr_ds,nsamples]=convert_bv_msk_header(hdr, params)
nsamples=1;
hdr_ds=struct();
hdr_ds.sa=struct();
hdr_ds.a.vol=get_vol_bv(hdr);
function data=convert_bv_msk_data(hdr, volumes)
require_singleton_volume(volumes);
data=double(hdr.Mask);
% BV volume time course input
% ---------------------------
function b=isa_bv_vtc(hdr)
b=(isa(hdr,'xff') || isstruct(hdr)) && isfield(hdr, 'VTCData');
function hdr=read_bv_vtc_header(fn)
hdr=get_and_check_data(fn, @xff_struct, @isa_bv_vtc);
function [hdr_ds,nsamples]=convert_bv_vtc_header(hdr, params)
nsamples=size(hdr.VTCData,1);
hdr_ds=struct();
hdr_ds.sa=struct();
hdr_ds.a.vol=get_vol_bv(hdr);
function data=convert_bv_vtc_data(hdr, volumes)
data=slice_4d(shiftdim(hdr.VTCData,1), volumes);
% BV volumetric MR (anatomy) input
% --------------------------------
function b=isa_bv_vmr(hdr)
b=(isa(hdr,'xff') || isstruct(hdr)) && isfield(hdr,'VMRData');
function hdr=read_bv_vmr_header(fn)
hdr=get_and_check_data(fn, @xff_struct, @isa_bv_vmr);
function [hdr_ds,nsamples]=convert_bv_vmr_header(hdr, params)
nsamples=1;
hdr_ds=struct();
hdr_ds.sa=struct();
hdr_ds.a.vol=get_vol_bv(hdr);
function data=convert_bv_vmr_data(hdr, volumes)
data=slice_4d(hdr.VMRData,volumes);
%%%%%%%%%%%%%%%%%%%%%%%%
% SPM structure
%%%%%%%%%%%%%%%%%%%%%%%%
% SPM input
% ---------
function b=isa_spm(hdr)
b=isstruct(hdr) && isfield(hdr,'xX') && isfield(hdr.xX,'X') && ...
isnumeric(hdr.xX.X) && isfield(hdr,'SPMid');
function hdr=read_spm_header(fn)
hdr=read_spm_struct_header(fn, []);
function hdr=read_spm_struct_header(fn_with_input_type, spm_header)
% input can be 'SPM.mat', 'SPM.mat:beta', 'SPM.mat:con', or
% 'SPM.mat:spm'. Output is the SPM struct with extra fields 'path' and
% 'input_type' added, so that convert_spm_header can get the correct
% fields
if isa_spm(fn_with_input_type)
fn_with_input_type='SPM.mat';
end
[fn,input_type]=get_spm_input_type(fn_with_input_type);
if isempty(spm_header)
% .mat file must be loaded still
spm_struct=load(fn);
if ~isstruct(spm_struct) || ...
~isequal(fieldnames(spm_struct),{'SPM'})
error('expected data with struct ''SPM'' in file ''%s''',fn);
end
spm_header=spm_struct.SPM;
end
hdr=spm_header;
hdr.path=fileparts(fn);
hdr.input_type=input_type;
get_and_check_data(hdr, [], @isa_spm);
function [fn,input_type]=get_spm_input_type(fn_with_input_type)
sep=':';
input_type=cosmo_strsplit(fn_with_input_type,sep,-1);
fn=fn_with_input_type;
switch input_type
case {'beta','con','spm'}
fn=fn(1:(end-numel(input_type)-numel(sep)));
otherwise
input_type='beta'; % the default; function will crash
% if fn is not a proper filename
end
function [hdr_ds,nsamples]=convert_spm_header(spm_struct, params)
if isfield(spm_struct,'input_type')
input_type=spm_struct.input_type;
else
input_type='beta';
end
if isfield(spm_struct,'path')
path=spm_struct.path;
else
path='';
end
% get data of interest
switch input_type
case 'beta'
input_vols=spm_struct.Vbeta;
input_labels=spm_struct.xX.name';
case 'con'
input_vols=[spm_struct.xCon.Vcon];
input_labels={spm_struct.xCon.name}';
case 'spm'
input_vols=[spm_struct.xCon.Vspm];
input_labels={spm_struct.xCon.name}';
otherwise
error('illegal data type %s', input_type);
end
n_input=numel(input_vols);
assert(numel(input_labels)==n_input);
sa=struct();
if isfield(spm_struct,'Sess') && strcmp(input_type,'beta')
% single subject GLM with betas; will use only betas of interest
% and set chunks based on runs
nruns=numel(spm_struct.Sess);
nbeta=numel(spm_struct.Vbeta);
sessions=zeros(nbeta,1);
beta_index=zeros(nbeta,1);
for k=1:nruns
sess=spm_struct.Sess(k);
sess_idxs=[sess.Fc.i];
row_idxs=sess.col(sess_idxs);
sessions(row_idxs)=k;
beta_index(row_idxs)=row_idxs;
end
keep_vol_msk=sessions>0;
sa.chunks=sessions(keep_vol_msk);
sa.beta_index=beta_index(keep_vol_msk);
else
% anything else: use all volumes
keep_vol_msk=true(n_input,1);
end
sa.labels=input_labels(keep_vol_msk);
sa.fname=cellfun(@(fn)fullfile(path,fn),...
{input_vols(keep_vol_msk).fname}',...
'UniformOutput',false);
nsamples=sum(keep_vol_msk);
if nsamples==0
error('Illegal: empty input');
end
% get volume info
first_vol=input_vols(1);
vol=struct();
vol.dim=first_vol.dim;
vol.mat=first_vol.mat;
hdr_ds=struct();
hdr_ds.a.vol=vol;
hdr_ds.sa=sa;
function hdr=read_spm_data(fn, hdr, volumes)
% store volumes, so that when convert_spm_data is called, this
% information is used
hdr.volumes=volumes;
function data=convert_spm_data(hdr, volumes)
% get SPM info
hdr_ds=convert_spm_header(hdr);
% see which files to load
file_names=hdr_ds.sa.fname;
n_files=numel(file_names);
if isfield(hdr, 'volumes')
% call after read_spm_data; volumes are already selected
% this must be done by the calling function
assert(isempty(volumes));
volumes=hdr.volumes;
end
if isempty(volumes)
volumes=1:n_files;
end
n_volumes=numel(volumes);
% allocate space for output
dim=hdr_ds.a.vol.dim;
data=zeros([dim n_volumes]);
%
for k=1:n_volumes
file_name=file_names{volumes(k)};
nii=load_untouch_nii(file_name);
data(:,:,:,k)=nii.img;
end
%%%%%%%%%%%%%%%%%%%%%%%%
% FieldTrip source
% FieldTrip source input
% ----------------------
function b=isa_ft_source(hdr)
b=isstruct(hdr) && ((isfield(hdr,'inside') && ...
isfield(hdr,'pos')) || ...
(cosmo_check_dataset(hdr,false) && ...
cosmo_isfield(hdr,'fa.pos')));
function hdr=read_ft_source_header(fn)
hdr=get_and_check_data(fn, @fast_import_data, @isa_ft_source);
function [ds,nsamples]=convert_ft_source_header(hdr, params)
if isfield(hdr,'inside') && isfield(hdr,'pos')
% fieldtrip struct
ds_meeg=cosmo_meeg_dataset(hdr);
else
% must be dataset struct with field .fa.pos
cosmo_isfield(hdr,'fa.pos',true);
ds_meeg=hdr;
end
cosmo_check_dataset(ds_meeg,'meeg');
ds=cosmo_vol_grid_convert(ds_meeg,'tovol');
nsamples=size(ds.samples,1);
function ds=convert_ft_source_data(ds_meeg, volumes)
ds=convert_ft_source_header(ds_meeg, struct());
ds=slice_dataset_volumes(ds, volumes);
%%%%%%%%%%%%%%%%%%%%%%%%
% PyMVPA fMRI datasset
function tf=isa_pymvpa_fmri(ds)
tf=isstruct(ds) && isfield(ds,'samples') && ...
all(cosmo_isfield(ds,{'a.imgaffine',...
'a.voxel_dim',...
'a.voxel_eldim',...
'fa.voxel_indices'}));
function ds=read_pymvpa_ds_header(ds)
ds=get_and_check_data(ds, [], @isa_pymvpa_fmri);
function [ds, nsamples]=convert_pymvpa_ds_header(pymvpa_ds, params)
% set volume information
vol=struct();
mat_base0=pymvpa_ds.a.imgaffine;
offset=mat_base0(1:3,1:3)*ones(3,1);
mat_base1=mat_base0;
mat_base1(1:3,4)=mat_base1(1:3,4)-offset;
vol.mat=mat_base1;
dim=pymvpa_ds.a.voxel_dim;
vol.dim=double(dim(:)');
vol.xform='unknown';
% set feature dimensions
dim_labels={'i';'j';'k'};
n_dim_labels=numel(dim_labels);
fdim=struct();
fdim.labels=dim_labels;
fdim.values=arrayfun(@(x)1:x,vol.dim(:),'UniformOutput',false);
% copy over data from PyMPVA struct
ds=pymvpa_ds;
% ensure samples are double
if ~isa(ds.samples,'double')
ds.samples=double(ds.samples);
end
% remove PyMVPA volume-specific fields
ds.a=rmfield(ds.a,{'imgaffine','voxel_dim','voxel_eldim'});
if isfield(ds.a,'mapper')
ds.a=rmfield(ds.a,'mapper');
end
% set vol and fdim
ds.a.vol=vol;
ds.a.fdim=fdim;
% set voxel indices
ijk_base0=double(ds.fa.voxel_indices);
ijk_base1=ijk_base0+1;
for k=1:n_dim_labels
ds.fa.(dim_labels{k})=ijk_base1(k,:);
end
ds.fa=rmfield(ds.fa,'voxel_indices');
% deal with scipy's 3d character arrays
ds.fa=convert_struct_with_3d_string_array_to_cellstr(ds.fa,2);
if isfield(ds,'sa')
ds.sa=convert_struct_with_3d_string_array_to_cellstr(ds.sa,1);
end
% set number of samples
nsamples=size(ds.samples,1);
function ds=convert_pymvpa_ds_data(ds, volumes)
ds=convert_pymvpa_ds_header(ds, struct());
ds=slice_dataset_volumes(ds,volumes);
function c=convert_3d_string_array_to_cellstr(arr,dim)
sz=size(arr);
assert(numel(sz)==3);
if dim==1
new_sz_idxs=[1 3];
else
new_sz_idxs=[2 3];
end
arr_2d=reshape(arr,sz(new_sz_idxs));
c=cellstr(arr_2d);
if dim==2
c=c';
end
function s=convert_struct_with_3d_string_array_to_cellstr(s,dim)
% deal with scipy's character arrays that represent 2d cell strings
fns=fieldnames(s);
for k=1:numel(fns)
fn=fns{k};
value=s.(fn);
if ischar(value) && numel(size(value))==3
s.(fn)=convert_3d_string_array_to_cellstr(value,dim);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%
% CoSMoMVPA datasset
% CoSMoMVPA dataset input
% -----------------------
function tf=isa_cosmo_fmri(ds)
tf=cosmo_check_dataset(ds,'fmri',false);
function ds=read_cosmo_ds_header(ds)
ds=get_and_check_data(ds, [], @isa_cosmo_fmri);
function [ds, nsamples]=convert_cosmo_ds_header(ds, params)
ds=get_and_check_data(ds, [], @isa_cosmo_fmri);
nsamples=size(ds.samples,1);
function ds=convert_cosmo_ds_data(ds, volumes)
ds=slice_dataset_volumes(ds,volumes);
