Example of a Surface-based searchlight
import AnatSearchlight.searchlight as sl
import numpy as np
import nibabel as nib
To define a searchlight we need a:
Pial Surface (individual)
Gray Matter Mask (individual)
Functional mask image (mask.nii from GLM) telling us which voxels are availabe for searchlight analysis
If we want to calculate the searchlight only for a part of a the surface, we can use a roi_specification (i.e. excluding the medial wall).
The searchlight can be specified to have a specific radius in mm (on the surface) or a specific number of voxels.
# Give it an appropriate structure name: Used for generating cifti-files for output
mySearchlight = sl.SearchlightSurface('left_cortex')
# These are the inndividual surfaces
surf = ['example_data/sub-02_space-32k_hemi-L_pial.surf.gii',
'example_data/sub-02_space-32k_hemi-L_white.surf.gii']
# This is the mask image for the functional data space
voxel_mask = 'example_data/sub-02_ses-s1_mask.nii'
# Optionally, you can define a roi-mask for the searchlight region on the surface. Centers outside of the mask will be ignored.
# Instead of a gifti-image you can also use a numpy array with the indices of the vertices that should be included in the searchlight.
roi_mask = 'example_data/tpl-fs32k_hemi-L_mask.label.gii'
# Here we defining
mySearchlight.define(surf,voxel_mask,roi=roi_mask,maxradius=20,maxvoxels=np.inf)
mySearchlight.save('example_data/searchlight_surf.h5')
Processing center 0 of 29759
Processing center 1000 of 29759
Processing center 2000 of 29759
Processing center 3000 of 29759
Processing center 4000 of 29759
Processing center 5000 of 29759
Processing center 6000 of 29759
Processing center 7000 of 29759
Processing center 8000 of 29759
Processing center 9000 of 29759
Processing center 10000 of 29759
Processing center 11000 of 29759
Processing center 12000 of 29759
Processing center 13000 of 29759
Processing center 14000 of 29759
Processing center 15000 of 29759
Processing center 16000 of 29759
Processing center 17000 of 29759
Processing center 18000 of 29759
Processing center 19000 of 29759
Processing center 20000 of 29759
Processing center 21000 of 29759
Processing center 22000 of 29759
Processing center 23000 of 29759
Processing center 24000 of 29759
Processing center 25000 of 29759
Processing center 26000 of 29759
Processing center 27000 of 29759
Processing center 28000 of 29759
Processing center 29000 of 29759
# You can now use the pre-computed searchlight.
mySearchlight = sl.load('example_data/searchlight_surf.h5')
print(f'Average size of the searchlight: {mySearchlight.radius.mean():.1f} mm')
print(f'Average number of voxels: {mySearchlight.nvoxels.mean():.1f}')
Average size of the searchlight: 20.0 mm
Average number of voxels: 325.5
# Define the MVPA function that will be used in the searchlight
# The functions take a n_images x n_voxels array as input.
# and return either a scalar or vector as output.
# Here we use simple example that calculate the mean activity in the searchlight region.
def mvpa_mean_function(data):
# Example of MVPA function that returns a scalar as an output argument
return np.mean(data)
def mvpa_multi_function(data):
# Example of MVPA function that returns a vector as an output argument
return np.mean(data,axis=1)
## Run the searchlight with the defined scalar MVPA function
datafiles = [f"example_data/sub-02_ses-s1_run-{r:02d}_reg-{s:02d}_beta.nii" for r in range(1,3) for s in range(2)]
results = mySearchlight.run(datafiles, mvpa_mean_function,verbose=False)
mySearchlight.data_to_cifti(results,outfilename='example_data/output.dscalar.nii')
<nibabel.cifti2.cifti2.Cifti2Image at 0x2071d7f99d0>
## Run the searchlight with the defined scalar vector function
datafiles = [f"example_data/sub-02_ses-s1_run-{r:02d}_reg-{s:02d}_beta.nii" for r in range(1,3) for s in range(2)]
results = mySearchlight.run(datafiles, mvpa_multi_function,verbose=False)
mySearchlight.data_to_cifti(results,outfilename='example_data/output_vector.dscalar.nii')
<nibabel.cifti2.cifti2.Cifti2Image at 0x2071f9a6d30>