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. #
