Merge pull request #26 from wmvanvliet/rois

Add support for ROIs

Author: wmvanvliet

examples/plot_rsa_roi.py

@@ -0,0 +1,117 @@
+#!/usr/bin/env python
+# coding: utf-8
+"""
+Source-level RSA using ROI's
+============================
+
+In this example, we use anatomical labels as Regions Of Interest (ROIs). Rather
+than using a searchlight, we compute DSMs for each ROI and then compute RSA
+with a single model DSM.
+
+The dataset will be the MNE-sample dataset: a collection of 288 epochs in which
+the participant was presented with an auditory beep or visual stimulus to
+either the left or right ear or visual field.
+"""
+# sphinx_gallery_thumbnail_number=2
+# Import required packages
+import mne
+import mne_rsa
+
+mne.set_log_level(True)  # Be less verbose
+mne.viz.set_3d_backend('pyvista')
+
+###############################################################################
+# We'll be using the data from the MNE-sample set. To speed up computations in
+# this example, we're going to use one of the sparse source spaces from the
+# testing set.
+sample_root = mne.datasets.sample.data_path(verbose=True)
+testing_root = mne.datasets.testing.data_path(verbose=True)
+sample_path = sample_root / 'MEG' / 'sample'
+testing_path = testing_root / 'MEG' / 'sample'
+subjects_dir = sample_root / 'subjects'
+
+###############################################################################
+# Creating epochs from the continuous (raw) data. We downsample to 100 Hz to
+# speed up the RSA computations later on.
+raw = mne.io.read_raw_fif(sample_path / 'sample_audvis_filt-0-40_raw.fif')
+events = mne.read_events(sample_path / 'sample_audvis_filt-0-40_raw-eve.fif')
+event_id = {'audio/left': 1,
+            'audio/right': 2,
+            'visual/left': 3,
+            'visual/right': 4}
+epochs = mne.Epochs(raw, events, event_id, preload=True)
+epochs.resample(100)
+
+###############################################################################
+# It's important that the model DSM and the epochs are in the same order, so
+# that each row in the model DSM will correspond to an epoch. The model DSM
+# will be easier to interpret visually if the data is ordered such that all
+# epochs belonging to the same experimental condition are right next to
+# each-other, so patterns jump out. This can be achieved by first splitting the
+# epochs by experimental condition and then concatenating them together again.
+epoch_splits = [epochs[cl] for cl in ['audio/left', 'audio/right',
+                                      'visual/left', 'visual/right']]
+epochs = mne.concatenate_epochs(epoch_splits)
+
+###############################################################################
+# Now that the epochs are in the proper order, we can create a DSM based on the
+# experimental conditions. This type of DSM is referred to as a "sensitivity
+# DSM". Let's create a sensitivity DSM that will pick up the left auditory
+# response when RSA-ed against the MEG data. Since we want to capture areas
+# where left beeps generate a large signal, we specify that left beeps should
+# be similar to other left beeps. Since we do not want areas where visual
+# stimuli generate a large signal, we specify that beeps must be different from
+# visual stimuli. Furthermore, since in areas where visual stimuli generate
+# only a small signal, random noise will dominate, we also specify that visual
+# stimuli are different from other visual stimuli. Finally left and right
+# auditory beeps will be somewhat similar.
+
+
+def sensitivity_metric(event_id_1, event_id_2):
+    """Determine similarity between two epochs, given their event ids."""
+    if event_id_1 == 1 and event_id_2 == 1:
+        return 0  # Completely similar
+    if event_id_1 == 2 and event_id_2 == 2:
+        return 0.5  # Somewhat similar
+    elif event_id_1 == 1 and event_id_2 == 2:
+        return 0.5  # Somewhat similar
+    elif event_id_1 == 2 and event_id_1 == 1:
+        return 0.5  # Somewhat similar
+    else:
+        return 1  # Not similar at all
+
+
+model_dsm = mne_rsa.compute_dsm(epochs.events[:, 2], metric=sensitivity_metric)
+mne_rsa.plot_dsms(model_dsm, title='Model DSM')
+
+###############################################################################
+# This example is going to be on source-level, so let's load the inverse
+# operator and apply it to obtain a cortical surface source estimate for each
+# epoch. To speed up the computation, we going to load an inverse operator from
+# the testing dataset that was created using a sparse source space with not too
+# many vertices.
+inv = mne.minimum_norm.read_inverse_operator(
+    f'{testing_path}/sample_audvis_trunc-meg-eeg-oct-4-meg-inv.fif')
+epochs_stc = mne.minimum_norm.apply_inverse_epochs(epochs, inv, lambda2=0.1111)
+
+###############################################################################
+# ROIs need to be defined as ``mne.Label`` objects. Here, we load the APARC
+# parcellation generated by FreeSurfer and treat each parcel as an ROI.
+rois = mne.read_labels_from_annot(parc='aparc', subject='sample',
+                                  subjects_dir=subjects_dir)
+
+###############################################################################
+# Performing the RSA. To save time, we don't use a searchlight over time, just
+# over the ROIs. The results are returned not only as a NumPy `ndarray`, but
+# also as an `mne.SourceEstimate` object, where each vertex beloning to the
+# same ROI has the same value.
+rsa_vals, stc = mne_rsa.rsa_stcs_rois(epochs_stc, model_dsm, inv['src'], rois,
+                                      temporal_radius=None, n_jobs=1,
+                                      verbose=False)
+
+###############################################################################
+# To plot the RSA values on a brain, we can use one of MNE-RSA's own
+# visualization functions.
+brain = mne_rsa.plot_roi_map(rsa_vals, rois, subject='sample',
+                             subjects_dir=subjects_dir)
+brain.show_view('lateral', distance=600)

mne_rsa/__init__.py

@@ -1,10 +1,11 @@
 __version__ = '0.8dev'
-from .source_level import rsa_stcs, dsm_stcs, rsa_nifti, dsm_nifti
+from .source_level import (rsa_stcs, dsm_stcs, rsa_stcs_rois, rsa_nifti,
+                           dsm_nifti)
 from .sensor_level import rsa_evokeds, rsa_epochs, dsm_evokeds, dsm_epochs
 from .searchlight import searchlight
 from .rsa import rsa, rsa_gen, rsa_array
 from .dsm import compute_dsm, compute_dsm_cv, dsm_array
-from .viz import plot_dsms, plot_dsms_topo
+from .viz import plot_dsms, plot_dsms_topo, plot_roi_map
 from .folds import create_folds
 
 # This function is useful to have nearby

mne_rsa/searchlight.py

@@ -1,3 +1,4 @@
+"""Classes and functions having to do with creating searchlights."""
 import numpy as np
 from mne.utils import logger
 
@@ -42,13 +43,27 @@ class searchlight:
         your distance computations.
 
         Defaults to ``None``.
-    spatial_radius : floats | None
-        The spatial radius of the searchlight patch in meters. All source
-        points within this radius will belong to the searchlight patch. Set to
-        None to only perform the searchlight over time. When this parameter is
-        set, the ``dist`` parameter must also be specified. Defaults to
-        ``None``.
-    temporal_radius : float | None
+    spatial_radius : float | list of list of int | None
+        This controls how spatial patches will be created. There are several
+        ways to do this:
+
+        The first way is to specify a spatial radius in meters. In this case,
+        the ``dist`` parameter must also be specified. This will create a
+        searchlight where each patch contains all source points within this
+        radius.
+
+        The second way is to specify a list of predefined patches. In this
+        case, each element of the list should itself be a list of integer
+        indexes along the spatial dimension of the data array. Each element of
+        this list will become a separate patch using the data at the specified
+        indices.
+
+        The third way is to set this to ``None``, which will disable the making
+        of spatial patches and only perform the searchlight over time. This can
+        be thought of as pooling everything into a single spatial patch.
+
+        Defaults to``None``.
+    temporal_radius : int | None
         The temporal radius of the searchlight patch in samples. Set to
         ``None`` to only perform the searchlight over sensors/source points.
         Defaults to ``None``.
@@ -126,21 +141,29 @@ def __init__(self, shape, dist=None, spatial_radius=None,
 
         # Will we be creating spatial searchlight patches?
         if self.spatial_radius is not None:
-            if self.dist is None:
-                raise ValueError('A spatial radius was requested, but no '
-                                 'distance information was specified '
-                                 '(=dist parameter).')
             if self.series_dim is None:
                 raise ValueError('Cannot create spatial searchlight patches: '
                                  f'the provided data shape ({shape}) has no '
                                  'spatial dimension.')
-            if self.sel_series is None:
-                self.sel_series = np.arange(shape[self.series_dim])
-
-            # Compressed Sparse Row format is optimal for our computations
-            from scipy.sparse import issparse
-            if issparse(self.dist):
-                self.dist = self.dist.tocsr()
+            # If spatial radius is a number, we will be making searchlight
+            # patches based on distance computations. Alternatively, a list of
+            # predefined spatial patches may be provided, and we don't need
+            # `dist`.
+            if type(self.spatial_radius) in [float, int]:
+                if self.dist is None:
+                    raise ValueError('A spatial radius was requested, but no '
+                                     'distance information was specified '
+                                     '(=dist parameter).')
+                # Compressed Sparse Row format is optimal for our computations
+                from scipy.sparse import issparse
+                if issparse(self.dist):
+                    self.dist = self.dist.tocsr()
+                if self.sel_series is None:
+                    self.sel_series = np.arange(shape[self.series_dim])
+            else:
+                # Explicit spatial patches were provided
+                if self.sel_series is None:
+                    self.sel_series = np.arange(len(self.spatial_radius))
 
         # Will we be creating temporal searchlight patches?
         if temporal_radius is not None:
@@ -199,6 +222,7 @@ def __init__(self, shape, dist=None, spatial_radius=None,
             self._generator = iter([tuple(self.patch_template)])
 
     def __iter__(self):
+        """Get an iterator over the searchlight patches."""
         return self
 
     def __next__(self):
@@ -211,7 +235,11 @@ def _iter_spatio_temporal(self):
         patch = list(self.patch_template)  # Copy the template
         for series in self.sel_series:
             # Compute all spatial locations in the searchligh path.
-            spat_ind = _get_in_radius(self.dist, series, self.spatial_radius)
+            if type(self.spatial_radius) in [float, int]:
+                spat_ind = _get_in_radius(self.dist, series,
+                                          self.spatial_radius)
+            else:
+                spat_ind = self.spatial_radius[series]
             patch[self.series_dim] = spat_ind
             for sample in self.time_centers:
                 temp_ind = slice(sample - self.temporal_radius,
@@ -224,7 +252,11 @@ def _iter_spatial(self):
         logger.info('Creating spatial searchlight patches')
         patch = list(self.patch_template)  # Copy the template
         for series in self.sel_series:
-            spat_ind = _get_in_radius(self.dist, series, self.spatial_radius)
+            if type(self.spatial_radius) in [float, int]:
+                spat_ind = _get_in_radius(self.dist, series,
+                                          self.spatial_radius)
+            else:
+                spat_ind = self.spatial_radius[series]
             patch[self.series_dim] = spat_ind
             yield tuple(patch)
 
@@ -239,7 +271,7 @@ def _iter_temporal(self):
 
     @property
     def shape(self):
-        """Number of generated patches along multiple dimensions.
+        """Get the number of generated patches along multiple dimensions.
 
         This is useful for re-shaping the result obtained after consuming the
         this generator.