function [nbrhood,vo,fo,out2in]=cosmo_surficial_neighborhood(ds, surfs, varargin)
% neighborhood definition for surface-based searchlight
%
% [nbrhood,vo,fo,out2in]=cosmo_surficial_neighborhood(ds, surfs, ...)
%
% Inputs:
% ds fMRI-like volumetric or surface-based dataset
% surfs cell with specification of surfaces. The following
% formats are valid options:
% # Volumetric datasets:
% * Single surface (Caret, BrainVoyager)
% - {fn, [start,stop]}
% - {fn, [start,stop], niter}
% - {v,f,[start,stop], niter}
% where
% + fn is the filename of a single surface
% + start and stop are distances towards and away
% from the center of a hemisphere. For example,
% if surface coordinates are in mm (which is
% typical), start=-3 and stop=2 means selecting
% voxels that are 3 mm or closer to the surface
% on the white-matter side, up to voxels that are
% 2 mm from the surface on the pial matter side
% + niter is the number of subsampling (mesh
% decimation) iterations to use for the output
% surface.
% + v are Px3 coordinates for P nodes
% + f are Qx3 node indices for Q faces
% * Twin surface (FreeSurfer)
% - {fn1, fn2}
% - {fn1, fn2, fn_o}
% - {fn1, fn2, niter}
% - {v1,v2,f}
% - {v1,v2,f,v_o,f_o}
% - {v1,v2,f,niter}
% where
% + fn{1,2} are filenames for pial and white
% surfaces
% + fn_o is the filename of an intermediate
% (node-wise average of a pial and white
% surface) output surface, and can be of a lower
% resolution than the pial and white surfaces.
% This can be achieved by using MapIcosahedron
% to generate meshes of different densities, e.g.
% ld=64 for the pial and white surface and ld=16
% for the intermediate surface. It is required
% that for every node in the intermediate surface
% there is a corresponding node (at the same
% spatial location) on the node-wise average of
% the pial and white surface.
% + v{1,2} are Px3 coordinates for P nodes of the
% pial and white surfaces
% + f are Qx3 node indices for Q faces on the pial
% and white surfaces
% + v_o and f_o are the nodes and vertices of an
% intermediate output surface. It is required
% that the vertices in v_o are a subset of the
% pair-wise averages of vertices in v1 and v2. A
% typical use case is using standard topologies
% from AFNI's MapIcosahedron, where the v1 and v2
% are high-res with X linear divisions, v_o is
% low-res with Y linear divisions, and Y<X and X
% is an integer multiple of Y.
% + niter is the number of subsampling (mesh
% decimation) iterations to use for the output
% surface.
% # Surface-based datasets:
% - {fn}
% - {fn, niter}
% - {v, f, niter}
% select neighbors:
% 'radius', r } either within radius r (usually in mm), grow
% 'count', c } the radius to get the nearest c neighbors,
% 'area', a } grow the radius to get neighbors whose area size
% } is less than or equal area a (usually in mm^2)
% 'direct', true } or get direct neighbors only (nodes who share an
% } edge)
% These four options are mutually exclusive
% 'metric',metric distance metric along cortical surface. One of
% 'geodesic' (default), 'dijkstra' or 'euclidean'.
% 'line_def',line_def definition of lines from inner to outer surface
% See surfing_nodeidxs2coords
% 'progress', p Show progress every p centers (default: 100)
% 'vol_def', v Volume definition with fields .mat (4x4 affine
% voxel-to-world matrix) and .dim (1x3 number of
% voxels in each dimension). If omittied, voldef is
% based on ds
% 'subsample_min_ratio', r When niter is provided and positive, it
% defines the minimum ratio between surfaces
% of old and new faces. This should help in
% preventing the generation of triangles with
% angles close to 180 degrees (default: 0.2).
% 'center_ids', c Center ids in output surface of the nodes for
% which the neighborhood should be computed. Empty
% means all nodes are used (default: [])
%
% Outputs:
% nbrhood neighborhood struct with fields
% .neighbors Mx1 cell with .neighbors{k} the
% indices of features (relative to ds)
% in the neighborhood of node k
% .a.fdim.labels set to {'node_indices'}
% .a.fdim.values set to {center_ids} with center_ids
% Mx1 ids of each neighborhood
% .fa.node_indices identical to center_ids
% .fa.radius } if ds is a surface-based dataset,
% .fa.area } these fields are present and contain
% } the radius and total surface area of
% the selected nodes
% v_o NVx3 coordinates for N nodes of the output surface
% f_o NFx3 node indices for Q faces of the output surface
% out2in Node mapping from output to input surface.
%
% Example:
% % this example uses BrainVoyager Surfaces
% anat_fn='SUB02_MPRAGE_ISO_IIHC_TAL.vmr';
% surf_fn='SUB02_MPRAGE_ISO_IIHC_TAL_LH_RECOSM.srf';
%
% % read dataset and surface
% ds=cosmo_fmri_dataset(anat_fn,'mask',anat_fn);
% fprintf('Dataset has %d samples and %d features\n',size(ds.samples));
%
% [v,f]=surfing_read(surf_fn);
% fprintf('Input surface has %d nodes and %d faces\n',size(v,1),size(f,1))
%
% % define neighborhood parameters
% count=100; % 100 voxels per searchlight
% niter=10; % 10 mesh decimation algorithms
%
% % surface definition; voxels are selected that are 3 mm or closer to the
% % surface on the white-matter side, up to voxels that are 1 mm from the
% % surface on the pial matter side.
%
% % Note: for FreeSurfer surfaces one could use
% % {surf_white_fn, surf_pial_fn [,niter]}
% surfs={surf_fn,[-3 1],niter};
% [nbrhood,vo,fo]=cosmo_surficial_neighborhood(ds,surfs,'radius',radius);
% fprintf('Neighborhood has %d elements\n', numel(nbrhood.neighbors))
% fprintf('Output surface has %d nodes and %d faces\n',size(vo,1),size(fo,1))
%
% % write decimated (smaller) output surface in ASCII format
% surfing_write(['small_surf.asc'],vo,fo);
%
% % define a measure that counts the number of features in each searchlight
% % (in more practical applications the measure could be, for example,
% % cosmo_correlation_measure or cosmo_crossvalidation_measure)
% voxel_counter=@(x,opt) cosmo_structjoin('samples',size(x.samples,2));
%
% % run a searchlight
% res=cosmo_searchlight(ds,voxel_counter,'nbrhood',nbrhood);
% fprintf('Searchlights have mean %.3f +/- %.3f features\n',...
% mean(res.samples), std(res.samples));
%
%
% % store surface dataset (in AFNI/SUMA NIML format) with the number of
% % features (voxels) for each center node of the output surface
% s=struct();
% s.node_indices=res.a.fdim.values{1}(res.fa.node_indices);
% s.data=res.samples';
% surfing_write('voxcount.niml.dset',s);
%
% % store volume dataset (in NIFTI format) with the number of times each
% % voxel was selected by a searchlight
% vol=ds;
% vol.samples(:)=0;
% for k=1:numel(nbrhood.neighbors);
% nbrs=nbrhood.neighbors{k};
% vol.samples(nbrs)=vol.samples(nbrs)+1;
% end
% cosmo_map2fmri(vol,'voxcount.nii');
% cosmo_map2fmri(ds,'anat.nii');
%
% Notes
% - Higher values of niter, or using a less dense mesh for fn_o or
% (v_o, f_o) means that the output surface has fewer nodes, which will
% decrease execution time for searchlights.
% - This function requires the surfing toolbox, surfing.sourceforge.net or
% github.com/nno/surfing
%
% See also: surfing, surfing_nodeidxs2coords, cosmo_searchlight
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
cosmo_check_external('surfing');
check_input(surfs, varargin{:});
defaults=struct();
defaults.metric='geodesic';
defaults.line_def=[10 0 1];
defaults.progress=100;
defaults.vol_def=[];
defaults.subsample_min_ratio=.2;
defaults.center_ids=[];
opt=cosmo_structjoin(defaults,varargin{:});
ds_type=get_ds_type(ds);
if strcmp(opt.metric,'geodesic') && ~isfield(opt,'direct')
cosmo_check_external('fast_marching');
end
% get surfaces
[v1,v2,f,vo,fo]=parse_surfs(surfs, ds_type);
ds_is_surface=strcmp(ds_type, 'surface');
one_surface=numel(v2)==2 || (ds_is_surface && isempty(v2));
% define intermediate high-res surface
if one_surface
vi=v1;
else
vi=(v1+v2)*.5;
end
switch numel(vo)
case 0
% no low-res surface, so use intermediate surface as source
vo=vi;
fo=f;
case 1
% subsample the surface
[vo,fo]=surfing_subsample_surface(vi,f,vo,...
opt.subsample_min_ratio, opt.progress);
end
% find mapping from low to high res surface
% if they are identical then low2high is the identity
out2in=surfing_maplow2hires(vo', vi');
% set center ids
center_ids=opt.center_ids(:)'; % ensure row vector
if isempty(opt.center_ids)
center_ids=1:size(vo,1);
else
out2in=out2in(center_ids);
end
% either surface or volumetric input dataset
ds_is_surface=cosmo_check_dataset(ds,'surface',false);
ds_is_volume=cosmo_check_dataset(ds,'fmri',false);
if sum([ds_is_surface, ds_is_volume])~=1
error('Unsupported dataset: expected surface or fmri dataset');
end
if ds_is_surface
% todo: use out2in to allow for subset of node indices as center
% (the current implementation computes neighborhoods for all
% nodes)
if ~isequal(center_ids,1:size(vo,1))
error('Unsupported center_ids option');
end
if ~isequal(out2in',1:numel(out2in))
%error('Unsupported other output surface than input surface');
end
nbrhood=surface_to_surface_neighborhood(ds,vi,f,opt);
elseif ds_is_volume
% get circle definition
circle_def=get_circle_def(opt);
% set volume definition
if isempty(opt.vol_def)
vol_def=ds.a.vol;
else
vol_def=opt.vol_def;
end
% define voxel mask based on features in dataset
ds_one=cosmo_slice(ds,1); % arbitrary sample
ds_one.samples(:)=1; % all set to one
ds_unflat=cosmo_unflatten(ds_one); % missing features have value of zero
vol_def.mask=ds_unflat~=0; % define mask
n2v=surfing_voxelselection(v1',v2',f',circle_def,vol_def,...
out2in,opt.line_def,...
opt.metric,0+opt.progress);
ncenters=numel(center_ids);
assert(ncenters==numel(n2v));
% determine unique center ids used
[unq_center_ids,unused,center_ids_mapping]=unique(center_ids);
% find mapping from features
nf_full=numel(vol_def.mask);
nf_mask=sum(vol_def.mask(:));
all2msk_indices=zeros(1,nf_full);
all2msk_indices(vol_def.mask)=1:nf_mask;
% store neighborhood information
nbrhood=struct();
nbrhood.a.fdim.labels={'node_indices'};
nbrhood.a.fdim.values={unq_center_ids};
nbrhood.fa.node_indices=center_ids_mapping(:)';
% set neighbors while considering the masked nature of the input dataset
neighbors=cell(ncenters,1);
for k=1:ncenters
msk_indices=all2msk_indices(double(n2v{k}));
assert(all(msk_indices>0));
neighbors{k}=msk_indices;
end
nbrhood.neighbors=neighbors;
else
assert(false,'this should never happen');
end
% compatibility wrapper for the surfing toolbox
nbrhood=ensure_neighbors_row_vectors(nbrhood);
origin=struct();
origin.a=ds.a;
origin.fa=ds.fa;
nbrhood.origin=origin;
cosmo_check_neighborhood(nbrhood,ds);
function nbrhood=surface_to_surface_neighborhood(ds,vertices,faces,opt)
circle_def=get_circle_def(opt);
dim_label='node_indices';
[two,fdim_index,attr_name,dim_name]=cosmo_dim_find(ds,...
dim_label,true);
fdim_nodes=ds.a.(dim_name).values{fdim_index};
fa_indices=ds.(attr_name).(dim_label);
if ~isequal(sort(fdim_nodes), unique(fdim_nodes))
error('values in .a.%s.values{%d} are not all unique',...
dim_name, fdim_index);
end
% - unq_nodes contains the node index associated with each unique
% node_indices feature in the dataset
[fa_idxs,unq_fa_indices]=cosmo_index_unique(fa_indices');
unq_nodes=fdim_nodes(unq_fa_indices);
nvertices=size(vertices,1);
too_large_index=find(unq_nodes>nvertices,1);
if any(too_large_index)
error(['surface has %d vertices, but maximum .fa.%s '...
'is %d'],...
nvertices,dim_label,fa_node_ids(too_large_index));
end
ignore_vertices_msk=true(nvertices,1);
ignore_vertices_msk(unq_nodes)=false;
vertices(ignore_vertices_msk,:)=NaN;
% set mapping from nodes to feature ids
ncenters=numel(unq_fa_indices);
node2feature_ids=cell(1,nvertices);
for k=1:ncenters
node=unq_nodes(k);
if ignore_vertices_msk(node)
feature_ids=cell(1,0);
else
feature_ids=fa_idxs{k};
end
node2feature_ids{node}=feature_ids;
end
% run node selection
[n2ns,radii]=surfing_nodeselection(vertices',faces',circle_def,...
opt.metric,opt.progress);
% set output
nbrhood=struct();
nbrhood.a.fdim.labels={'node_indices'};
nbrhood.a.fdim.values={fdim_nodes};
nbrhood.neighbors=cell(ncenters,1);
nbrhood.fa.radius=zeros(1,ncenters);
nbrhood.fa.node_indices=zeros(1,ncenters);
% add area fieldname
node2area=surfing_surfacearea(vertices,faces);
nbrhood.fa.area=zeros(1,ncenters);
for k=1:ncenters
center_node=unq_nodes(k);
nbrhood.fa.node_indices(k)=unq_fa_indices(k);
nbrhood.fa.radius(k)=radii(center_node);
nbrhood.fa.area(k)=sum(node2area(n2ns{center_node}));
if ignore_vertices_msk(center_node)
around_feature_ids=zeros(1,0);
else
around_nodes=n2ns{center_node};
around_feature_ids=cat(1,node2feature_ids{around_nodes});
end
nbrhood.neighbors{k}=around_feature_ids(:)';
end
function nbrhood=ensure_neighbors_row_vectors(nbrhood)
% compatibility wrapper for the surfing toolbox
for k=1:numel(nbrhood.neighbors)
if size(nbrhood.neighbors{k},1)~=1
nbrhood.neighbors{k}=nbrhood.neighbors{k}';
end
end
function [v1,v2,f,vo,fo]=parse_surfs(surfs, ds_type)
% helper function to get surfaces
ds_is_surface=strcmp(ds_type,'surface');
[v1,v2,f,vo,fo]=parse_surfs_arguments(ds_is_surface, surfs);
if numel(v2)~=2 && (numel(f)==2 || isempty(f))
% swap position
temp=f;
f=v2;
v2=temp;
end
check_surf_arguments(ds_is_surface,v1,v2,f,vo,fo);
function check_surf_arguments(ds_is_surface,v1,v2,f,vo,fo)
if isempty(v1) || (isempty(v2) && ~ds_is_surface) || isempty(f)
error('Not enough arguments for surfaces and topology');
end
% c2 can be vector with two elements indicating the size of the curved
% cylinder with a circle as basis; if not, it should be a surface with
% the same size as c1
one_surface=numel(v2)==2;
if ~(one_surface || ds_is_surface || isequal(size(v1),size(v2)))
error('Size mismatch between surfaces: %dx%d != %dx%d',...
size(v1,1), size(v1,2),...
size(v2,1), size(v2,2));
end
surfing_check_surface(v1,f);
if ~(one_surface || ds_is_surface)
surfing_check_surface(v2,f);
end
if isempty(fo)
if numel(vo)>1
% if not a scaler (niter) then throw an error
error('Topology missing for output surface');
end
else
surfing_check_surface(vo,fo);
end
function [v1,v2,f,vo,fo]=parse_surfs_arguments(ds_is_surface,surfs)
if ~iscell(surfs)
error('surfs argument must be a cell');
end
n=numel(surfs);
% space for output
v1=[];
v2=[];
f=[];
vo=[];
fo=[];
for k=1:n
s=surfs{k};
if ischar(s)
% filename; read the surface
[c,f_]=surfing_read(s);
if isempty(v1)
v1=c;
f=f_;
elseif isempty(v2) && ~ds_is_surface
if ~isequal(f_,f)
error('Topology mismatch between the two surfaces');
end
v2=c;
elseif isempty(vo)
vo=c;
fo=f_;
else
error('Superfluous argument at position %d', k);
end
elseif isnumeric(s)
% numeric array; coordinates or faces
if isempty(v1)
v1=s;
elseif isempty(v2)
v2=s;
elseif isempty(f)
f=s;
elseif isempty(vo)
vo=s;
elseif isempty(fo)
fo=s;
else
error('Superfluous argument at position %d', k);
end
else
error('Expected surface filename or array at position %d', k);
end
end
function ds_type=get_ds_type(ds)
supported_ds_types={'surface','fmri'};
n=numel(supported_ds_types);
for k=1:n
supported_ds_type=supported_ds_types{k};
if cosmo_check_dataset(ds,supported_ds_type,false)
ds_type=supported_ds_type;
return
end
end
% maybe it is not a dataset at all, check this possibility
cosmo_check_dataset(ds);
% it is a valid dataset but not fmri or surface;
% try to give an appropriate error message
for k=1:n
supported_ds_type=supported_ds_types{k};
if is_ds_type_like(ds, supported_ds_type)
% let it throw an error
cosmo_check_dataset(ds,supported_ds_type);
end
end
error('Unknown dataset type, supported are: %s.', ...
cosmo_strjoin(supported_ds_types,', '));
function tf=is_ds_type_like(ds, type)
% helper function that uses heuristics to find type of dataset
tf=false;
switch type
case 'surface'
if cosmo_isfield(ds,'fa.node_indices')
tf=true;
return;
end
if cosmo_isfield(ds,'a.fdim.values') && ...
isequal({'node_indices'},ds.a.fdim.values)
tf=true;
return;
end
case 'fmri'
if any(cosmo_isfield(ds,{'fa.i','fa.j','fa.k'}))
tf=true;
return;
end
if cosmo_isfield(ds,'a.fdim.values') && ...
any(cosmo_match({'i','j','k'},ds.a.fdim.values))
tf=true;
return
end
otherwise
error('unsupported type %s', type);
end
function check_input(surfs, varargin)
% give deprecation notice
if isnumeric(surfs)
error(['Second argument must be a cell with a surface '...
'definition (as of Jan 2015, the syntax for this '...
'function has changed)']);
end
function circle_def=get_circle_def(opt)
metric2circle_def_func=struct();
metric2circle_def_func.radius=@(x)x;
metric2circle_def_func.count=@(x)[0 x]; % fixed initial radius
metric2circle_def_func.area=@(x)[5 Inf x]; % default initial radius
metric2circle_def_func.direct=@get_direct_circle_def;
% options are mutually exclusive
metrics=fieldnames(metric2circle_def_func);
metric_msk=cosmo_isfield(opt,metrics);
if sum(metric_msk)~=1
error('Use one of these arguments to define neighbors: %s',...
cosmo_strjoin(metrics, ', '));
end
metric=metrics{metric_msk};
func=metric2circle_def_func.(metric);
param=opt.(metric);
if ~isscalar(param) || param<0
error('value for ''%s'' must be a scalar not less than zero',...
metric);
end
circle_def=func(opt.(metric));
function circle_def=get_direct_circle_def(radius)
% false or 0 give a zero radius, otherwise NaN
if isnan(radius) || radius
circle_def=NaN;
else
circle_def=[5 1];
end
