cosmo surficial neighborhood hdrΒΆ

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}
%    'radius', r         } select neighbors either within radius r, grow
%    'count', c          } the radius to get neighbors are c locations,
%    'direct', true      } or get direct neighbors only (nodes who share an
%                        } edge)
%                        } These three 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
%     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.           #