Skip to content
Permalink
Browse files

EXC: add surface-based searchlight exercise

  • Loading branch information...
nno committed Jul 25, 2019
1 parent f8e5914 commit ad315c930df27789e78158d3301cc767038902bd
@@ -34,7 +34,7 @@ Detailed
contents_demo
get_started
download
labman2017
nmsm2019
faq
contact
contribute
@@ -52,7 +52,9 @@ Full documentation
cimec2014
rhul2016
cimec2016
prni2016
labman2017
nmsm2019
faq
contact
contribute
@@ -0,0 +1,45 @@
.. # For CoSMoMVPA's license terms and conditions, see #
# the COPYING file distributed with CoSMoMVPA #
.. _`ex_surface_searchlight`:

Surface-based fMRI searchlight
==============================

Reading material
----------------

- :cite:`OWDD11`: One of the early papers using a surface-based searchlight, with comparison between volumetric and surface-based searchlight..

Required toolboxes
------------------
- surfing toolbox: https://github.com/nno/surfing
- AFNI Matlab toolbox: https://sscc.nimh.nih.gov/pub/dist/tgz/afni_matlab.tgz (or part of the full AFNI distribution in the ``src/matlab`` directory, see https://github.com/afni/afni)

Background
++++++++++
In this exercise, data is analyzed from one participant who pressed buttons with either the index or middle finger in blocks. We try to infer where in the brain
Surface-models were reconstructed using FreeSurfer's ``recon-all``; further processing was done using AFNI and the script ``prep_afni_surf.py`` that is part of PyMVPA. Surfaces representing the left and right hemispheres were merged as described in the FAQ.


Exercise
++++++++

For this exercise use the `digit` dataset.

Part 1 (cortical thickness)
---------------------------
Load the anatomical surface models for the outer (pial) and inner (white) surfaces that separate the grey matter from non-gray matter. Compute, for each node on the surface, the cortical thickness, and then plot the thickness on a 3D surface model.

Part 2 (Classificiation analysis)
---------------------------------
Load the anatomical surface models and the functional data.
Define a surface-based neighborhood with approximately 100 voxels per searchlight center. Then run the searchlight with a classifier to distinguish between the different digit presses and visualize the results.

Hint: :ref:`run_surface_searchlight_skl`

Solution: :ref:`run_surface_searchlight` / :pb:`surface_searchlight`


.. include:: links.txt

@@ -26,6 +26,7 @@ Contents:
ex_meeg_searchlight
ex_meeg_time_generalization
ex_rsa_tutorial
ex_surface_searchlight
ex_multiple_comparisons


@@ -201,7 +201,7 @@ Friday
-------------- ---------------------------------------------------------------------------------------------------
12:30 Lunch break
-------------- ---------------------------------------------------------------------------------------------------
14:00 Surface-based searchlight. (Exercise to be written)
14:00 Surface-based searchlight. :doc:`ex_surface_searchlight`
-------------- ---------------------------------------------------------------------------------------------------
15:30 Coffee break
-------------- ---------------------------------------------------------------------------------------------------
@@ -0,0 +1,279 @@
%% Demo: fMRI surface-based searchlights with LDA classifier
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% + 'digit'
% A participant made finger pressed with the index and middle finger of
% the right hand during 4 runs in an fMRI study. Each run was divided in
% 4 blocks with presses of each finger and analyzed with the GLM,
% resulting in 2*4*4=32 t-values
%
% A searchlight is run with a 100 voxel searchlight, using a
% disc for which the metric radius varies from node to node.
%
% This example requires the surfing toolbox, github.com/nno/surfing
%
% This example may take quite some time to run. For faster execution but
% lower spatial precision, set ld=16 below; for slower execution use ld=64.
%
% If you use this code for a publication, please cite:
% Oosterhof, N.N., Wiestler, T, Downing, P.E., & Diedrichsen, J. (2011)
% A comparison of volume-based and surface-based information mapping.
% Neuroimage. DOI:10.1016/j.neuroimage.2010.04.270
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #

%% Check externals
cosmo_check_external('surfing');
cosmo_check_external('afni');

%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config=cosmo_config();

digit_study_path=fullfile(config.tutorial_data_path,'digit');
readme_fn=fullfile(digit_study_path,'README');
cosmo_type(readme_fn);

output_path=config.output_data_path;

%%

% resolution parameter for input surfaces
% 64 is for high-quality results; use 16 for fast execution
surface_ld=16;

% Define twin surface filenames (FreeSurfer)
pial_fn=fullfile(digit_study_path,...
sprintf('ico%d_mh.pial_al.asc', surface_ld));
white_fn=fullfile(digit_study_path,...
sprintf('ico%d_mh.smoothwm_al.asc', surface_ld));


% read the surface in pial_fn using surfing_read, and assign the vertices
% and faces to variables pial_v and pial_f, respectively
% >@@>
[pial_v, pial_f] = surfing_read(pial_fn);
% <@@<

fprintf('The pial surface has %d vertices, %d faces\n',...
size(pial_v,1), size(pial_f,1))

% do the same for the white_fn, assign the faces and vertices to
% white_v and white_f
% >@@>
[white_v, white_f] = surfing_read(white_fn);
% <@@<

fprintf('The white surface has %d vertices, %d faces\n',...
size(pial_v,1), size(pial_f,1))

% verify that the face information in pial_f and white_f are the same
assert(isequal(white_f, pial_f));

% show the content of the surfaces
fprintf('pial_v\n');
cosmo_disp(pial_v)

fprintf('pial_f\n');
cosmo_disp(pial_f)

fprintf('white_v\n');
cosmo_disp(white_v)

fprintf('white_f\n');
cosmo_disp(white_f)


%% Part 1: compute thickness of the cortex

% compute the element-wise difference in coordinates between pial_v
% and white_v, and assign to a variable delta
% >@@>
delta = pial_v - white_v;
% <@@<

% square the differences element-wise, and assign to delta_squared.
% hint: use ".^2"
% >@@>
delta_squared = delta .^ 2;
% <@@<

% compute the thickness squared, by summing the elements in
% delta_squared along the second dimension. Assign to thickness_squared
% >@@>
thickness_squared = sum(delta_squared,2);
% <@@<

% finally compute the thickness by taking the square root, assign
% the result to thickness
% >@@>
thickness = sqrt(thickness_squared);
% <@@<
%%
% plot a histogram of the thickness values
% >@@>
hist(thickness,100);
xlabel('thickness (mm)');
% <@@<

%% For visualization purposed, read inflated surface
inflated_fn=fullfile(digit_study_path,...
sprintf('ico%d_mh.inflated_alCoMmedial.asc', surface_ld));
[infl_v,infl_f]=surfing_read(inflated_fn);
fprintf('The inflated surface has %d vertices, %d faces\n',...
size(infl_v,1), size(infl_f,1))

%% visualize surface in Matlab using AFNI Matlab toolbox
nvertices=size(infl_v,1);

min_thickness=1;
max_thickness=4;
show_edge=false;

opt=struct();
opt.ShowEdge=show_edge;
opt.Zlim=[min_thickness,max_thickness]; % this does not seem to work
opt.Dim='3D';

if show_edge
t='with edges';
else
t='without edges';
end

desc='thickness';
header=strrep([desc ' ' t],'_',' ');

range_adj_thickness=max(min_thickness,...
min(thickness,max_thickness));

DispIVSurf(infl_v,infl_f,1:nvertices,...
range_adj_thickness,0,opt);
title(sprintf('Inflated %s', header));


%%
ds_thickness=struct();
ds_thickness.fa.node_indices=1:nvertices;
ds_thickness.samples=thickness(:)';
ds_thickness.a.fdim.labels={'node_indices'};
ds_thickness.a.fdim.values={(1:nvertices)'};

% July 2019: strangely enough these gifti files cannot be read by AFNI SUMA.
% https://github.com/CoSMoMVPA/CoSMoMVPA/issues/186
if cosmo_check_external('gifti',false)
output_fn=fullfile(config.output_data_path,'thickness.gii');
cosmo_map2surface(ds_thickness,output_fn,'encoding','ASCII');
end

% AFNI output
output_fn=fullfile(config.output_data_path,'thickness.niml.dset');
cosmo_map2surface(ds_thickness,output_fn);

%% Part 3: run surface-based searchlight

% Load volumetric functional data
data_path=digit_study_path;
data_fn=fullfile(data_path,'glm_T_stats_perblock+orig');

% set targets
targets=repmat(1:2,1,16)'; % class labels: 1 2 1 2 1 2 1 2 1 2 ... 1 2
chunks=floor(((1:32)-1)/8)+1; % run labels: 1 1 1 1 1 1 1 1 2 2 ... 4 4

% load functional data
ds = cosmo_fmri_dataset(data_fn,'targets',targets,'chunks',chunks);

% remove zero elements
zero_msk=all(ds.samples==0,1);
ds = cosmo_slice(ds, ~zero_msk, 2);

fprintf('Dataset has %d samples and %d features\n', size(ds.samples));

% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds);

%% set measure arguments
% Assign to measure a function handle to cosmo_cross_validation_measure
% >@@>
measure = @cosmo_crossvalidation_measure;
% <@@<

% use as arguments for the measure:
% - classifier: cosmo_classify_naive_bayes
% - partitions: odd-even partitions
% assign these to a struct measure_args
% >@@>
measure_args = struct();
measure_args.classifier = @cosmo_classify_naive_bayes;
measure_args.partitions = cosmo_oddeven_partitioner(ds);
% <@@<


%%
% Set neighborhood parameters
% Make a cell with the outer surface vertices (pial_v),
% the inner surface vertices (white_v), and the face
% indices (pial_f or white_f). Assign to a variable surface_def
% >@@>
surface_def={pial_v,white_v,pial_f};
% <@@<

% Define a surface-based neighborhood (using cosmo_surficial_neighborhood)
% with approximately 100 voxels per searchlight. Assign the result
% to nbrhood.

% >@@>
nbrhood=cosmo_surficial_neighborhood(ds,surface_def,'count',100);
% <@@<

% visualize the neighborhood using cosmo_disp.
% What is the feature attribute of the neighborhood?

%%
% run surface-based searchlight using the variables define above
% Assign the result to ds_sl
% >@@>
ds_sl=cosmo_searchlight(ds,nbrhood,measure,measure_args);
% <@@<

%% plot surface in Matlab using AFNI Matlab toolbox

nvertices=size(infl_v,1);
show_edge=false;

opt=struct();
opt.ShowEdge=show_edge;
opt.Dim='3D';

if show_edge
t='with edges';
else
t='without edges';
end

desc='classification accuracy';
header=strrep([desc ' ' t],'_',' ');

DispIVSurf(infl_v,infl_f,1:nvertices,...
ds_sl.samples',0,opt);
title(sprintf('Inflated %s', header));

% July 2019: strangely enough these gifti files cannot be read by AFNI SUMA.
% https://github.com/CoSMoMVPA/CoSMoMVPA/issues/186
if cosmo_check_external('gifti',false)
output_fn=fullfile(config.output_data_path,'digit_accuracy.gii');
cosmo_map2surface(ds_sl,output_fn,'encoding','ASCII');
end

% AFNI output
output_fn=fullfile(config.output_data_path,'digit_accuracy.niml.dset');
cosmo_map2surface(ds_sl,output_fn);

% Show citation information
cosmo_check_external('-cite');

0 comments on commit ad315c9

Please sign in to comment.
You can’t perform that action at this time.