WIP: Plot glass brain

.travis.yml

@@ -14,41 +14,41 @@ matrix:
         # Full (Linux, 2.7) first half
         - os: linux
           env: CONDA_DEPENDENCIES="numpy scipy matplotlib pandas scikit-learn h5py pillow statsmodels mayavi pytest"
-               PIP_DEPENDENCIES="pysurfer nitime faulthandler joblib nibabel codecov nose-timer nose-faulthandler pytest-sugar pytest-cov pytest-attrib"
+               PIP_DEPENDENCIES="pysurfer nitime nilearn faulthandler joblib nibabel codecov nose-timer nose-faulthandler pytest-sugar pytest-cov pytest-attrib"
                MNE_DIRS=". beamformer channels commands connectivity datasets decoding forward gui inverse_sparse io"
                OPTION=""
         # Full (Linux, 2.7) second half
         - os: linux
           env: CONDA_DEPENDENCIES="numpy scipy matplotlib pandas scikit-learn h5py pillow statsmodels mayavi pytest"
-               PIP_DEPENDENCIES="pysurfer nitime faulthandler joblib nibabel codecov nose-timer nose-faulthandler pytest-sugar pytest-cov pytest-attrib"
+               PIP_DEPENDENCIES="pysurfer nitime nilearn faulthandler joblib nibabel codecov nose-timer nose-faulthandler pytest-sugar pytest-cov pytest-attrib"
                MNE_DIRS="minimum_norm preprocessing realtime simulation stats time_frequency viz"
                OPTION=""
 
         # OSX first half
         - os: osx
           env: CONDA_DEPENDENCIES="numpy scipy matplotlib pandas scikit-learn h5py pillow statsmodels mayavi pytest"
-               PIP_DEPENDENCIES="pysurfer nitime faulthandler joblib nibabel codecov nose-timer nose-faulthandler pytest-sugar pytest-cov pytest-attrib"
+               PIP_DEPENDENCIES="pysurfer nitime nilearn faulthandler joblib nibabel codecov nose-timer nose-faulthandler pytest-sugar pytest-cov pytest-attrib"
                MNE_DIRS=". beamformer channels commands connectivity datasets decoding forward gui inverse_sparse io"
                OPTION=""
         # OSX second half
         - os: osx
           env: CONDA_DEPENDENCIES="numpy scipy matplotlib pandas scikit-learn h5py pillow statsmodels mayavi pytest"
-               PIP_DEPENDENCIES="pysurfer nitime faulthandler joblib nibabel codecov nose-timer nose-faulthandler pytest-sugar pytest-cov pytest-attrib"
+               PIP_DEPENDENCIES="pysurfer nitime nilearn faulthandler joblib nibabel codecov nose-timer nose-faulthandler pytest-sugar pytest-cov pytest-attrib"
                MNE_DIRS="minimum_norm preprocessing realtime simulation stats time_frequency viz"
                OPTION=""
 
         # Py3k + non-default stim channel first half
         - os: linux
           env: PYTHON_VERSION=3.6 TEST_LOCATION=install MNE_STIM_CHANNEL=STI101
                CONDA_DEPENDENCIES="numpy scipy matplotlib pandas scikit-learn h5py pillow statsmodels nose pytest"
-               PIP_DEPENDENCIES="nitime joblib nibabel codecov nose-timer pytest-sugar pytest-cov pytest-attrib"
+               PIP_DEPENDENCIES="nitime nilearn joblib nibabel codecov nose-timer pytest-sugar pytest-cov pytest-attrib"
                MNE_DIRS=". beamformer channels commands connectivity datasets decoding forward gui inverse_sparse io"
                OPTION=""
         # Py3k + non-default stim channel second half
         - os: linux
           env: PYTHON_VERSION=3.6 TEST_LOCATION=install MNE_STIM_CHANNEL=STI101
                CONDA_DEPENDENCIES="numpy scipy matplotlib pandas scikit-learn h5py pillow statsmodels nose pytest"
-               PIP_DEPENDENCIES="nitime joblib nibabel codecov nose-timer pytest-sugar pytest-cov pytest-attrib"
+               PIP_DEPENDENCIES="nitime nilearn joblib nibabel codecov nose-timer pytest-sugar pytest-cov pytest-attrib"
                MNE_DIRS="minimum_norm preprocessing realtime simulation stats time_frequency viz"
                OPTION=""
 

doc/python_reference.rst

@@ -244,6 +244,7 @@ Visualization
    plot_evoked_field
    plot_evoked_white
    plot_filter
+   plot_stc_glass_brain
    plot_head_positions
    plot_ideal_filter
    plot_compare_evokeds

doc/whats_new.rst

@@ -87,6 +87,8 @@ Changelog
 
     - :class:`mne.decoding.SlidingEstimator` and :class:`mne.decoding.GeneralizingEstimator` now accept **fit_params at fitting by `Jean-Remi King`_
 
+    - Add :func:`mne.viz.plot_stc_glass_brain` for plotting source estimates on glass brain by `Jaakko Leppakangas`_
+
 BUG
 ~~~
 

examples/inverse/plot_compute_mne_inverse_volume.py

@@ -25,6 +25,8 @@
 data_path = sample.data_path()
 fname_inv = data_path + '/MEG/sample/sample_audvis-meg-vol-7-meg-inv.fif'
 fname_evoked = data_path + '/MEG/sample/sample_audvis-ave.fif'
+subjects_dir = data_path + '/subjects'
+fname_t1 = subjects_dir + '/sample/mri/T1.mgz'
 
 snr = 3.0
 lambda2 = 1.0 / snr ** 2
@@ -39,16 +41,22 @@
 stc = apply_inverse(evoked, inverse_operator, lambda2, method)
 stc.crop(0.0, 0.2)
 
-# Export result as a 4D nifti object
+# Export result as a 4D NIfTI object
 img = stc.as_volume(src,
                     mri_resolution=False)  # set True for full MRI resolution
 
-# Save it as a nifti file
+# Save it as a NIfTI file:
 # nib.save(img, 'mne_%s_inverse.nii.gz' % method)
 
-t1_fname = data_path + '/subjects/sample/mri/T1.mgz'
+###############################################################################
+# Plot with our wrapper to :func:`nilearn.plotting.plot_glass_brain`:
 
-# Plotting with nilearn ######################################################
-plot_stat_map(index_img(img, 61), t1_fname, threshold=8.,
+stc.plot_glass_brain(subject='sample', subjects_dir=subjects_dir,
+                     src=src, vmin=8, vmax=26)
+
+###############################################################################
+# Plot directly with :func:`nilearn.plotting.plot_stat_map`:
+
+plot_stat_map(index_img(img, 61), fname_t1, threshold=8.,
               title='%s (t=%.1f s.)' % (method, stc.times[61]))
 plt.show()

mne/datasets/utils.py

@@ -228,7 +228,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
     path = _get_path(path, key, name)
     # To update the testing or misc dataset, push commits, then make a new
     # release on GitHub. Then update the "releases" variable:
-    releases = dict(testing='0.36', misc='0.3')
+    releases = dict(testing='0.37', misc='0.3')
     # And also update the "hashes['testing']" variable below.
 
     # To update any other dataset, update the data archive itself (upload
@@ -285,7 +285,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
         sample='1d5da3a809fded1ef5734444ab5bf857',
         somato='f3e3a8441477bb5bacae1d0c6e0964fb',
         spm='ecce87351d88def59d3d4cdc561e2a60',
-        testing='03e84fe9a50dfeb04ab5f138ceed441d',
+        testing='19b828839f51fb33d2d8f6f46c87978e',
         multimodal='26ec847ae9ab80f58f204d09e2c08367',
         visual_92_categories='46c7e590f4a48596441ce001595d5e58',
         mtrf='273a390ebbc48da2c3184b01a82e4636',

mne/source_estimate.py

@@ -23,7 +23,7 @@
                            _ensure_src_subject, SourceSpaces)
 from .utils import (get_subjects_dir, _check_subject, logger, verbose,
                     _time_mask, warn as warn_, copy_function_doc_to_method_doc)
-from .viz import plot_source_estimates
+from .viz import plot_source_estimates, plot_stc_glass_brain
 from .io.base import ToDataFrameMixin, TimeMixin
 
 from .externals.six import string_types
@@ -963,6 +963,22 @@ def transform(self, func, idx=None, tmin=None, tmax=None, copy=False):
 
         return stcs
 
+    @copy_function_doc_to_method_doc(plot_stc_glass_brain)
+    def plot_glass_brain(self, subject=None, subjects_dir=None, src=None,
+                         initial_time=0., display_mode='ortho', colorbar=False,
+                         axes=None, title=None, threshold='auto',
+                         annotate=True, black_bg=False, cmap=None, alpha=0.7,
+                         vmin=None, vmax=None, plot_abs=True,
+                         symmetric_cbar="auto", show=True):
+        return plot_stc_glass_brain(self, subject, subjects_dir=subjects_dir,
+                                    src=src, initial_time=initial_time,
+                                    display_mode=display_mode,
+                                    colorbar=colorbar, axes=axes, title=title,
+                                    threshold=threshold, annotate=annotate,
+                                    black_bg=black_bg, cmap=cmap, alpha=alpha,
+                                    vmin=vmin, vmax=vmax, plot_abs=plot_abs,
+                                    symmetric_cbar=symmetric_cbar, show=show)
+
 
 def _center_of_mass(vertices, values, hemi, surf, subject, subjects_dir,
                     restrict_vertices):
@@ -2689,6 +2705,9 @@ def save_stc_as_volume(fname, stc, src, dest='mri', mri_resolution=False):
     if not isinstance(stc, VolSourceEstimate):
         raise Exception('Only volume source estimates can be saved as '
                         'volumes')
+    if not isinstance(src, SourceSpaces) or src.kind != 'volume':
+        raise TypeError('src must be a volume source space, got %s'
+                        % (repr(src),))
 
     n_times = stc.data.shape[1]
     shape = src[0]['shape']

mne/source_space.py

@@ -1158,14 +1158,8 @@ def vertex_to_mni(vertices, hemis, subject, subjects_dir=None, mode=None,
 
 
 @verbose
-def _read_talxfm(subject, subjects_dir, mode=None, verbose=None):
-    """Read MNI transform from FreeSurfer talairach.xfm file.
-
-    Adapted from freesurfer m-files. Altered to deal with Norig
-    and Torig correctly.
-    """
-    if mode is not None and mode not in ['nibabel', 'freesurfer']:
-        raise ValueError('mode must be "nibabel" or "freesurfer"')
+def _read_talairach(subject, subjects_dir, verbose=None):
+    """Read MNI transform from FreeSurfer talairach.xfm file."""
     fname = op.join(subjects_dir, subject, 'mri', 'transforms',
                     'talairach.xfm')
     # read the RAS to MNI transform from talairach.xfm
@@ -1195,6 +1189,19 @@ def _read_talxfm(subject, subjects_dir, mode=None, verbose=None):
         else:
             raise ValueError('failed to find \'Linear_Transform\' string in '
                              'xfm file:\n%s' % fname)
+    return xfm
+
+
+@verbose
+def _read_talxfm(subject, subjects_dir, mode=None, verbose=None):
+    """Read MRI-MNI transform from FreeSurfer talairach.xfm file.
+
+    Adapted from freesurfer m-files. Altered to deal with Norig
+    and Torig correctly.
+    """
+    if mode is not None and mode not in ['nibabel', 'freesurfer']:
+        raise ValueError('mode must be "nibabel" or "freesurfer"')
+    xfm = _read_talairach(subject, subjects_dir, verbose=verbose)
 
     # Setup the RAS to MNI transform
     ras_mni_t = {'from': FIFF.FIFFV_MNE_COORD_RAS,

mne/utils.py

@@ -1097,6 +1097,7 @@ def wrapper(func):
 requires_nitime = partial(requires_module, name='nitime')
 requires_h5py = partial(requires_module, name='h5py')
 requires_numpydoc = partial(requires_module, name='numpydoc')
+requires_nilearn = partial(requires_module, name='nilearn')
 
 
 def check_version(library, min_version):

mne/viz/_3d.py

@@ -27,7 +27,7 @@
 from ..io.pick import pick_types
 from ..io.constants import FIFF
 from ..io.meas_info import read_fiducials
-from ..source_space import SourceSpaces, _create_surf_spacing, _check_spacing
+from ..source_space import (SourceSpaces, _create_surf_spacing, _check_spacing)
 
 from ..surface import (_get_head_surface, get_meg_helmet_surf, read_surface,
                        transform_surface_to, _project_onto_surface,

mne/viz/__init__.py

@@ -23,3 +23,4 @@
                   _plot_sources_raw, _plot_sources_epochs, plot_ica_properties)
 from .montage import plot_montage
 from .decoding import plot_gat_matrix, plot_gat_times
+from .source_estimate import plot_stc_glass_brain

mne/viz/source_estimate.py

@@ -0,0 +1,168 @@
+# -*- coding: utf-8 -*-
+# Authors: Jaakko Leppakangas <jaeilepp@gmail.com>
+#
+# License: Simplified BSD
+
+"""Functions to plot source estimate data."""
+from __future__ import print_function
+
+import os.path as op
+
+import numpy as np
+
+from ..source_space import _read_talairach
+from ..surface import _compute_nearest, read_surface
+from ..utils import has_nibabel, get_subjects_dir, _check_subject
+from ..transforms import apply_trans
+from .utils import plt_show
+
+
+def plot_stc_glass_brain(stc, subject, subjects_dir=None, src=None,
+                         initial_time=0., display_mode='ortho', colorbar=False,
+                         axes=None, title=None, threshold='auto',
+                         annotate=True, black_bg=False, cmap=None, alpha=0.7,
+                         vmin=None, vmax=None, plot_abs=True,
+                         symmetric_cbar="auto", show=True):
+    """Plot stc on glass brain.
+
+    Parameters
+    ----------
+    stc : instance of SourceEstimate | instance of VolSourceEstimate
+        Data to plot.
+    subject : str | None
+        The subject name corresponding to FreeSurfer environment
+        variable SUBJECT.  If None stc.subject will be used. If that
+        is None, the environment will be used.
+    subjects_dir : None | str
+        The path to the freesurfer subjects reconstructions.
+        It corresponds to Freesurfer environment variable SUBJECTS_DIR.
+        The default is None.
+    src : instance of SourceSpace | None
+        The source space. Only needed for VolumeSourceEstimate. Has no effect
+        when using instance of SourceEstimate.
+    initial_time : float
+        The time instance to plot in seconds. Defaults to 0.
+    display_mode : string
+        Choose the direction of the cuts: 'x' - sagittal, 'y' - coronal,
+        'z' - axial, 'l' - sagittal left hemisphere only,
+        'r' - sagittal right hemisphere only, 'ortho' - three cuts are
+        performed in orthogonal directions. Possible values are: 'ortho',
+        'x', 'y', 'z', 'xz', 'yx', 'yz', 'l', 'r', 'lr', 'lzr', 'lyr',
+        'lzry', 'lyrz'. Defaults to 'ortho'.
+    colorbar : bool
+        If True, display a colorbar on the right of the plots. Defaults to
+        False.
+    axes : matplotlib axes or 4 tuple of float: (xmin, ymin, width, height)
+        The axes, or the coordinates, in matplotlib figure space,
+        of the axes used to display the plot. If None, the complete
+        figure is used.
+    title : string
+        The title displayed on the figure.
+    threshold : float | None | 'auto'
+        The value used to threshold the image: values below the threshold (in
+        absolute value) are plotted as transparent. If None, the image is not
+        hresholded. If float, If auto, the threshold is determined magically by
+        analysis of the image.
+    annotate : bool
+        If True, positions and left/right annotation are added to the plot.
+    black_bg : bool
+        If True, the background of the image is set to be black. If you wish to
+        save figures with a black background, you will need to pass
+        "facecolor='k', edgecolor='k'" to matplotlib.pyplot.savefig.
+    cmap : matplotlib colormap
+        The colormap for specified image
+    alpha : float
+        Alpha transparency for the brain schematics.
+    vmax : float
+        Upper bound for plotting.
+    vmin : float
+        Lower bound for plotting.
+    plot_abs : bool
+        If True (default), maximum intensity projection of the absolute value
+        will be used (rendering positive and negative values in the same
+        manner). If False, the sign of the maximum intensity will be
+        represented with different colors. See `nilearn documentation`_
+        for examples.
+    symmetric_cbar : bool | 'auto'
+        Specifies whether the colorbar should range from -vmax to vmax
+        or from vmin to vmax. Setting to 'auto' will select the latter if
+        the range of the whole image is either positive or negative.
+        Note: The colormap will always be set to range from -vmax to vmax.
+    show : bool
+        Whether to show the figure.
+
+    Returns
+    -------
+    fig : instance of matplotlib.Figure
+        The matplotlib figure.
+
+    Notes
+    -----
+    This function requires :mod:`nilearn` and :mod:`nibabel`.
+
+    .. _nilearn documentation: http://nilearn.github.io/auto_examples/01_plotting/plot_demo_glass_brain_extensive.html
+    """  # noqa: E501
+    from scipy import sparse
+    if has_nibabel(vox2ras_tkr=True):
+        import nibabel as nib
+    else:
+        raise ImportError('This function requires nibabel.')
+    try:
+        from nilearn import plotting
+    except ImportError:
+        raise ImportError('This function requires nilearn.')
+    from ..source_estimate import VolSourceEstimate, SourceEstimate
+
+    subject = _check_subject(stc.subject, subject, True)
+    subjects_dir = get_subjects_dir(subjects_dir=subjects_dir,
+                                    raise_error=True)
+    time_idx = np.argmin(np.abs(stc.times - initial_time))
+    xfm = _read_talairach(subject, subjects_dir)
+    if isinstance(stc, VolSourceEstimate):
+        if src is None:
+            raise ValueError('src cannot be None when plotting volume source '
+                             'estimate.')
+        img = stc.as_volume(src)
+        affine = np.dot(xfm, img.affine)
+        data = img.get_data()[:, :, :, time_idx]
+        img = nib.Nifti1Image(data, affine)
+    elif isinstance(stc, SourceEstimate):
+        n_vertices = sum(len(v) for v in stc.vertices)
+        offset = 0
+        surf_to_mri = 0.
+        for hi, hemi in enumerate(['lh', 'rh']):
+            ribbon = nib.load(op.join(subjects_dir, subject, 'mri',
+                                      '%s.ribbon.mgz' % hemi))
+            xfm = ribbon.header.get_vox2ras_tkr()
+            mri_data = ribbon.get_data()
+            ijk = np.array(np.where(mri_data)).T
+            xyz = apply_trans(xfm, ijk) / 1000.
+            row_ind = np.where(mri_data.ravel())[0]
+            data = np.ones(row_ind.size)
+            rr = read_surface(op.join(subjects_dir, subject, 'surf',
+                                      '%s.white' % hemi))[0]
+            rr /= 1000.
+            rr = rr[stc.vertices[hi]]
+            col_ind = _compute_nearest(rr, xyz) + offset
+            surf_to_mri = surf_to_mri + sparse.csr_matrix(
+                (data, (row_ind, col_ind)), shape=(mri_data.size, n_vertices))
+            offset += len(stc.vertices[hi])
+
+        data = surf_to_mri.dot(stc.data[:, time_idx])
+        data.shape = mri_data.shape
+        xfm = _read_talairach(subject, subjects_dir)
+
+        affine = np.dot(xfm, ribbon.header.get_vox2ras())
+        img = nib.Nifti1Image(data, affine)
+
+    else:
+        ValueError('Glass brain plotting is only supported for SourceEstimate '
+                   'and VolSourceEstimate. Got %s.' % type(stc))
+    fig = plotting.plot_glass_brain(
+        img, display_mode=display_mode, colorbar=colorbar, axes=axes,
+        title=title, threshold=threshold, annotate=annotate,
+        black_bg=black_bg, cmap=cmap, alpha=alpha, vmax=vmax, vmin=vmin,
+        plot_abs=plot_abs, symmetric_cbar=symmetric_cbar)
+    # resampling_interpolation='nearest')
+    plt_show(show)
+    return fig