# 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
%
%
% Input:
%
% 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
%  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 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.

## 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
%       .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');
%                                             '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 ]
%     %||     [ 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.
%