function ds=cosmo_surface_dataset(fn, varargin)
% Returns a dataset structure based on surface mesh data
%
% ds=cosmo_surface_dataset(filename, varargin)
%
% Inputs:
% filename filename of surface data to be loaded. Currently
% supported are '.niml.dset' (AFNI/SUMA NIML),
% 'smp' (BrainVoyager surface maps), and 'gii' (GIFTI).
% Also supported are structs as provided by
% afni_niml_readsimple, or objects provided by xff or
% the MATLAB/Octave GIfTI Library.
% 'targets', t Px1 targets for P samples; these will be stored in
% the output as ds.sa.targets
% 'chunks', c Px1 chunks for P samples; these will be stored in the
% the output as ds.sa.chunks
% Output:
% ds dataset struct
%
% Examples:
% % (this example requires the surfing toolbox)
% cosmo_skip_test_if_no_external('afni');
% %
% % construct AFNI NIML dataset struct
% cosmo_check_external('afni');
% niml=struct();
% niml.data=[1 2; 3 4; 5 6];
% niml.node_indices=[1 20 201];
% niml.stats={'Ttest(10)','Zscore()'};
% %
% % make surface dataset
% % (a filename of a NIML dataset in ASCII format is supported as well)
% ds=cosmo_surface_dataset(niml);
% cosmo_disp(ds)
% %|| .samples
% %|| [ 1 3 5
% %|| 2 4 6 ]
% %|| .sa
% %|| .stats
% %|| { 'Ttest(10)'
% %|| 'Zscore()' }
% %|| .fa
% %|| .node_indices
% %|| [ 1 2 3 ]
% %|| .a
% %|| .fdim
% %|| .labels
% %|| { 'node_indices' }
% %|| .values
% %|| { [ 2 21 202 ] }
%
%
% Notes:
% - this function is intended for datasets with surface data, i.e. with
% one or more values associated with each surface node. It does not
% support anatomical surface meshes that contain node coordinates and
% faces. To read and write such anatomical meshes, consider the surfing
% toolbox, github.com/nno/surfing
% - data can be mapped back to a surface file format using
% cosmo_map2surface
%
% Dependencies:
% - for Brainvoyager files (.smp), it requires the NeuroElf
% toolbox, available from: http://neuroelf.net
% - for AFNI/SUMA NIML files (.niml.dset) it requires the AFNI
% Matlab toolbox, available from: http://afni.nimh.nih.gov/afni/matlab/
% - for GIFTI files (.gii) it requires the MATLAB/Octave GIfTI Library,
% available from https://github.com/gllmflndn/gifti
%
% See also: cosmo_map2surface
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
defaults=struct();
defaults.targets=[];
defaults.chunks=[];
params = cosmo_structjoin(defaults, varargin);
ds=get_dataset(fn);
% set targets and chunks
ds=set_vec_sa(ds,'targets',params.targets);
ds=set_vec_sa(ds,'chunks',params.chunks);
% check consistency
cosmo_check_dataset(ds,'surface');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general helper functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ds=get_dataset(fn)
obj_to_clean=[];
% try to load from .mat file
obj_from_mat=try_dataset_from_mat(fn);
if isempty(obj_from_mat)
obj=fn;
else
obj=obj_from_mat;
obj_to_clean=obj;
end
% try to load from file
if ischar(obj)
obj=try_obj_from_file(obj);
obj_to_clean=obj;
end
[ds,fmt]=try_dataset_from_obj(obj);
if ~isempty(obj_to_clean)
apply_cleaner(obj_to_clean,fmt);
end
function obj=try_dataset_from_mat(fn)
obj=[];
if ischar(fn) && isempty(cosmo_strsplit(fn,'.mat',-1))
obj=importdata(fn);
end
function [obj,fmt]=try_obj_from_file(fn)
assert(ischar(fn));
[exts_cell,fmts]=get_format_attribute('exts');
endswith=@(s,e) isempty(cosmo_strsplit(s,e,-1));
endswith_any=@(s,ce) any(cellfun(@(x)endswith(s,x),ce));
match=cellfun(@(x)endswith_any(fn,x), exts_cell);
i=find_single_or_error(match, ...
@()sprintf('unknown extension for file %s',...
fn));
fmt=fmts{i};
reader=fmt.reader;
obj=reader(fn);
function [ds,fmt]=try_dataset_from_obj(obj)
if cosmo_check_dataset(obj,'surface',false)
% it is already a valid dataset
ds=obj;
fmt='';
return;
end
[matcher_cell,fmts]=get_format_attribute('matcher');
match=cellfun(@(f) f(obj),matcher_cell);
i=find_single_or_error(match, ...
@()sprintf('unknown object of type %s',...
class(obj)));
fmt=fmts{i};
ds=build_dataset(obj,fmt);
function ds=build_dataset(obj,fmt)
builder=fmt.builder;
[data,node_indices,sa]=builder(obj);
nfeatures=size(data,2);
if nfeatures~=numel(node_indices)
error(['The number of features (%d) does not match the number '...
'of node indices (%d)'],nfeatures,numel(node_indices));
end
ds=struct();
ds.samples=data;
% set sample attributes
ds.sa=sa;
% set feature attributes
ds.fa.node_indices=1:nfeatures;
% set dataset attributes
fdim=struct();
fdim.labels={'node_indices'};
fdim.values={node_indices(:)'};
ds.a.fdim=fdim;
function apply_cleaner(obj,fmt)
if ~isempty(fmt) && isfield(fmt,'cleaner')
cleaner=fmt.cleaner;
cleaner(obj);
end
function i=find_single_or_error(m, error_text_func)
i=find(m);
switch numel(i)
case 0
error(error_text_func());
case 1
% ok
otherwise
assert(false,'this should not happen');
end
function [elems,fmts]=get_format_attribute(fmt_key)
img_formats=get_img_formats();
keys=fieldnames(img_formats);
n_fmts=numel(keys);
elems=cell(n_fmts,1);
fmts=cell(n_fmts,1);
for k=1:n_fmts
key=keys{k};
fmt=img_formats.(key);
elems{k}=fmt.(fmt_key);
fmts{k}=fmt;
end
function ds=read_surface_dataset(fn)
[data,node_indices,sa]=read(fn);
nfeatures=size(data,2);
if nfeatures~=numel(node_indices)
error(['The number of features (%d) does not match the number '...
'of node indices (%d)'],nfeatures,numel(node_indices));
end
ds=struct();
ds.samples=data;
% set sample attributes
ds.sa=sa;
% set feature attributes
ds.fa.node_indices=1:nfeatures;
% set dataset attributes
fdim=struct();
fdim.labels={'node_indices'};
fdim.values={node_indices(:)'};
ds.a.fdim=fdim;
function ds=set_vec_sa(ds, label, values)
if isempty(values)
return;
end
if numel(values)==1
nsamples=size(ds.samples,1);
values=repmat(values,nsamples,1);
end
ds.sa.(label)=values(:);
function [data,node_indices,sa]=read(fn)
img_formats=get_img_formats();
endswith=@(s,e) isempty(cosmo_strsplit(s,e,-1));
names=fieldnames(img_formats);
n=numel(names);
for k=1:n
name=names{k};
img_format=img_formats.(name);
matcher=img_format.matcher;
builder=img_format.builder;
if matcher(fn)
[data,node_indices,sa]=builder(fn);
return
elseif ischar(fn)
exts=img_format.exts;
reader=img_format.reader;
for j=1:numel(exts)
if endswith(fn,exts{j})
cosmo_check_external(img_format.externals);
s=reader(fn);
if ~matcher(s)
error(['Could read file ''%s'', but the data '...
'is not matching the expected contents '...
'for a dataset of type %s'],...
fn, name);
end
[data,node_indices,sa]=builder(s);
return
end
end
end
end
if ischar(fn)
error('Unsupported extension in filename %s', fn);
else
error('Unsupported input of type %s', class(fn));
end
function img_formats=get_img_formats()
% define which formats are supports
% .exts indicates the extensions
% .matcher says whether a struct is of the type
% .reader should read a filaname and return a struct
% .externals are fed to cosmo_check_externals
img_formats=struct();
img_formats.niml_dset.exts={'.niml.dset'};
img_formats.niml_dset.matcher=@isa_niml_dset;
img_formats.niml_dset.reader=@read_niml_dset;
img_formats.niml_dset.builder=@build_niml_dset;
img_formats.niml_dset.externals={'afni'};
img_formats.bv_smp.exts={'.smp'};
img_formats.bv_smp.matcher=@isa_bv_smp;
img_formats.bv_smp.reader=@read_bv_smp;
img_formats.bv_smp.builder=@build_bv_smp;
img_formats.bv_smp.cleaner=@(x)x.ClearObject();
img_formats.bv_smp.externals={'neuroelf'};
img_formats.gii.exts={'.gii'};
img_formats.gii.matcher=@isa_gii;
img_formats.gii.reader=@read_gii;
img_formats.gii.builder=@build_gii;
img_formats.gii.externals={'gifti'};
img_formats.pymvpa.exts={};
img_formats.pymvpa.matcher=@isa_pymvpa;
img_formats.pymvpa.reader=@read_pymvpa;
img_formats.pymvpa.builder=@build_pymvpa;
img_formats.pymvpa.externals={};
function b=isa_pymvpa(x)
b=isstruct(x) && ...
isfield(x,'samples') && ...
isfield(x,'fa') && ...
any(cosmo_isfield(x.fa, pymvpa_get_node_indices_fields()));
function b=read_pymvpa(fn)
assert(false,'this function should not have been called');
function [data,node_indices,sa]=build_pymvpa(s)
data=s.samples;
index_fields=pymvpa_get_node_indices_fields();
common_fields=intersect(index_fields,fieldnames(s.fa));
assert(~isempty(common_fields));
first_field=common_fields{1};
node_indice_base0=s.fa.(first_field);
node_indices=double(node_indice_base0)+1;
sa=struct();
if isfield(s,'sa')
sa=convert_struct_with_3d_string_array_to_cellstr(s.sa,1);
end
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
function fields=pymvpa_get_node_indices_fields()
fields={'node_indices','center_ids'};
function b=isa_gii(x)
b=isa(x,'gifti') && isfield(x,'cdata');
function s=read_gii(fn)
s=gifti(fn);
function [data,node_indices,sa]=build_gii(s)
data=s.cdata';
if ~isa(data,'double')
data=double(data);
end
if isfield(s,'indices')
node_indices=double(s.indices);
else
nv=size(data,2);
node_indices=1:nv;
end
sa=struct();
function b=isa_niml_dset(x)
b=isstruct(x) && isfield(x,'data');
function s=read_niml_dset(fn)
s=afni_niml_readsimple(fn);
function [data,node_indices,sa]=build_niml_dset(s)
data=s.data';
if ~isa(data,'double')
data=double(data);
end
if isfield(s,'node_indices')
node_indices=s.node_indices+1; % base0 -> base1
else
nv=size(data,2);
node_indices=1:nv;
end
sa=struct();
if isfield(s,'stats')
sa.stats=cosmo_statcode(s.stats);
end
if isfield(s,'labels')
sa.labels=s.labels(:);
end
function b=isa_bv_smp(x)
b=isa(x,'xff') && isfield(x,'NrOfVertices') && isfield(x,'Map');
function s=read_bv_smp(fn)
s=xff(fn);
neuroelf_bless_wrapper(s);
function [data,node_indices,sa]=build_bv_smp(s)
nsamples=s.NrOfMaps;
nfeatures=s.NrOfVertices;
data=zeros(nsamples,nfeatures);
node_indices=1:nfeatures;
labels=cell(nsamples,1);
for k=1:nsamples
map=s.Map(k);
data(k,:)=map.SMPData(:);
labels{k}=map.Name;
end
sa=struct();
sa.labels=labels;
sa.stats=cosmo_statcode(s);
function result=neuroelf_bless_wrapper(arg)
% deals with recent neuroelf (>v1.1), where bless is deprecated
s=warning('off','neuroelf:xff:deprecated');
resetter=onCleanup(@()warning(s));
result=bless(arg);
