Define voxel selection for a searchlight¶
In a searchlight (Kriegeskorte et al 2006) one computes a (spatially) locally contrained measure for each feature in a dataset.
An important building block for a searchlight is the voxel selection part: for each feature, decide which features are in its neighborhood. Typically the neighborhood is defined as a sphere of a certain radius.
Compute voxel offsets in a sphere¶
In turn, computing the relative offsets from a voxel location is a building block for feature selection part. In this exercise, write a function that returns the voxel indices relative to the origin. That is, given a radius r, it returns a Px3 array where every row (i,j,k) is unique and it holds that (i,j,k) is at most at distance r from the origin (0,0,0). This function has the following signature:
function [offsets,distances]=cosmo_sphere_offsets(radius)
% computes sub index offsets for voxels in a sphere
%
% [offsets,distances]=cosmo_sphere_offsets(radius)
%
% Input:
% radius radius of the sphere (in voxel units)
%
% Output
% offsets Px3 sub indices relative to the origin (0,0,0).
% offsets(p,:)=[i, j, k] means that the euclidian distance
% between points at (i,j,k) and the origin is less than or
% or equal to radius
% distances Px1 distances from the origin (in voxel units).
%
% Example:
% % compute offsets for voxels that share a side or edge, but not those
% % that only share a corner (because sqrt(2) < 1.5 < sqrt(3)).
% [o, d]=cosmo_sphere_offsets(1.5);
% cosmo_disp(o);
% %|| [ 0 0 0
% %|| -1 0 0
% %|| 0 -1 0
% %|| : : :
% %|| 1 0 -1
% %|| 1 0 1
% %|| 1 1 0 ]@19x3
% cosmo_disp(d)
% %|| [ 0
% %|| 1
% %|| 1
% %|| :
% %|| 1.41
% %|| 1.41
% %|| 1.41 ]@19x1
%
% Notes:
% - this function computes distances in voxel space, not in world space.
% If voxels have the same size sof mm along all dimensions (i.e. are
% of size s x s x s mm^2 for some value of s; this property is known as
% "isotropic") then using radius=d/s will select voxels within
% distance d mm. This function is less suitable for spherical offsets
% in world space when voxels are non-isotropic.
% - offsets and distances are sorted by distance.
%
% See also sub2ind, ind2sub, cosmo_spherical_neighborhood
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
Then, for each radius from 2 to 6 in steps of .5, plot the indices and show how many voxels are in a sphere of that raidus.
Hint: cosmo sphere offsets skl / run sphere offsets skl
Solution: cosmo sphere offsets / run sphere offsets / Matlab output: run_sphere_offsets
Using the sphere offsets for voxel selection¶
If you are up for quite an advanced exercise: write a function that performs voxel selection (if not, don’t spend your time on this and just look at the answer). It should have the following signature:
function nbrhood=cosmo_spherical_neighborhood(ds, varargin)
% computes neighbors for a spherical searchlight
%
% nbrhood=cosmo_spherical_neighborhood(ds, opt)
%
% Inputs
% ds a dataset struct, either:
% - in fmri form (from cosmo_fmri_dataset), when
% ds.fa has the fields .i, .j and .k
% - in meeg source form (from cosmo_meeg_dataset),
% when ds.fa has the field .pos. In this case, the
% features must have positions that can be
% converted to a grid.
% 'radius', r } either use a radius of r, or select
% 'count', c } approximately c voxels per searchlight
% Notes:
% - These two options are mutually exclusive
% - When using this option for an fmri dataset, the
% radius r is expressed in voxel units; for an meeg
% source dataset, the radius r is in whatever units
% the source dataset uses for the positions
% 'progress', p show progress every p features (default: 1000)
%
% Outputs
% nbrhood dataset-like struct without .sa or .samples, with:
% .a dataset attributes, from dataset.a
% .fa feature attributes with the same fields as fs.fa,
% and in addition the fields:
% .nvoxels 1xP number of voxels in each searchlight
% .radius 1xP radius in voxel units
% .center_ids 1xP feature center id
% .neighbors Px1 cell so that center2neighbors{k}==nbrs contains
% the feature ids of the neighbors of feature k
% If the dataset has a field ds.fa.inside, then
% features that are not inside are not included as
% neighbors in the output
% .origin Has fields .a and .fa from input dataset
%
%
% Example:
% ds=cosmo_synthetic_dataset('type','fmri');
% radius=1; % radius=3 is typical for 'real-world' searchlights
% nbrhood=cosmo_spherical_neighborhood(ds,'radius',radius,...
% 'progress',false);
% cosmo_disp(nbrhood)
% %|| .a
% %|| .fdim
% %|| .labels
% %|| { 'i' 'j' 'k' }
% %|| .values
% %|| { [ 1 2 3 ] [ 1 2 ] [ 1 ] }
% %|| .vol
% %|| .mat
% %|| [ 2 0 0 -3
% %|| 0 2 0 -3
% %|| 0 0 2 -3
% %|| 0 0 0 1 ]
% %|| .dim
% %|| [ 3 2 1 ]
% %|| .xform
% %|| 'scanner_anat'
% %|| .fa
% %|| .nvoxels
% %|| [ 3 4 3 3 4 3 ]
% %|| .radius
% %|| [ 1 1 1 1 1 1 ]
% %|| .center_ids
% %|| [ 1 2 3 4 5 6 ]
% %|| .i
% %|| [ 1 2 3 1 2 3 ]
% %|| .j
% %|| [ 1 1 1 2 2 2 ]
% %|| .k
% %|| [ 1 1 1 1 1 1 ]
% %|| .neighbors
% %|| { [ 1 4 2 ]
% %|| [ 2 1 5 3 ]
% %|| [ 3 2 6 ]
% %|| [ 4 1 5 ]
% %|| [ 5 4 2 6 ]
% %|| [ 6 5 3 ] }
% %|| .origin
% %|| .a
% %|| .fdim
% %|| .labels
% %|| { 'i'
% %|| 'j'
% %|| 'k' }
% %|| .values
% %|| { [ 1 2 3 ]
% %|| [ 1 2 ]
% %|| [ 1 ] }
% %|| .vol
% %|| .mat
% %|| [ 2 0 0 -3
% %|| 0 2 0 -3
% %|| 0 0 2 -3
% %|| 0 0 0 1 ]
% %|| .dim
% %|| [ 3 2 1 ]
% %|| .xform
% %|| 'scanner_anat'
% %|| .fa
% %|| .i
% %|| [ 1 2 3 1 2 3 ]
% %|| .j
% %|| [ 1 1 1 2 2 2 ]
% %|| .k
% %|| [ 1 1 1 1 1 1 ]
%
%
% Notes:
% - this function can return neighborhoods with either a fixed number of
% features, or a fixed radius. When used with a searchlight, the
% former has the advantage that the number of features is less
% variable (especially near edges of the brain, in an fmri dataset),
% which can make it easier to compare result in different regions as
% the number of features can affect
% pattern discriminablity. The latter has the advantage that the
% smoothness of the output maps under the null hypothesis can be more
% uniformly smooth.
%
% See also: cosmo_fmri_dataset, cosmo_meeg_dataset, cosmo_searchlight
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
As a start you can ignore the requirement that a negative value for radius should select a certain number of voxels; in other words just focus on positive values for radius.
A more advanced exercise: modify this function so that it also supports negative values for radius, which return neighborhoods of the same size (i.e. same number of voxels) at each location.
Solution: cosmo spherical neighborhood
