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
