[ENH] Add tutorial on time-frequency source estimation with STC viewer GUI

doc/changes/latest.inc

@@ -30,6 +30,7 @@ Enhancements
 - Add :func:`mne.beamformer.apply_dics_tfr_epochs` to apply a DICS beamformer to time-frequency resolved epochs (:gh:`11096` by `Alex Rockhill`_)
 - Check whether head radius (estimated from channel positions) is correct when reading EEGLAB data with :func:`~mne.io.read_raw_eeglab` and :func:`~mne.read_epochs_eeglab`. If head radius is not within likely values, warn informing about possible units mismatch and the new ``montage_units`` argument (:gh:`11283` by `Mikołaj Magnuski`_).
 - Add support for a callable passed in ``combine`` for `mne.time_frequency.AverageTFR.plot` and `mne.time_frequency.AverageTFR.plot_joint` (:gh:`11329` by `Mathieu Scheltienne`_)
+- Add to :ref:`ex-source-loc-methods` how to apply inverse methods to time-frequency resolved epochs and use :func:`mne.gui.view_vol_stc` to view the output (:gh:`10920` by `Alex Rockhill`_)
 
 Bugs
 ~~~~

doc/conf.py

@@ -272,6 +272,7 @@
     # unlinkable
     'CoregistrationUI',
     'IntracranialElectrodeLocator',
+    'VolSourceEstimateViewer',
     'mne_qt_browser.figure.MNEQtBrowser',
 }
 numpydoc_validate = True

doc/mri.rst

@@ -34,5 +34,6 @@ Step by step instructions for using :func:`gui.coregistration`:
    transforms.apply_volume_registration
    transforms.compute_volume_registration
    vertex_to_mni
+   gui.view_vol_stc
    warp_montage_volume
    coreg.Coregistration

examples/inverse/evoked_ers_source_power.py

@@ -13,6 +13,7 @@
 """
 # Authors: Luke Bloy <luke.bloy@gmail.com>
 #          Eric Larson <larson.eric.d@gmail.com>
+#          Alex Rockhill <aprockhill@mailbox.org>
 #
 # License: BSD-3-Clause
 
@@ -22,7 +23,7 @@
 import mne
 from mne.cov import compute_covariance
 from mne.datasets import somato
-from mne.time_frequency import csd_morlet
+from mne.time_frequency import csd_tfr
 from mne.beamformer import (make_dics, apply_dics_csd, make_lcmv,
                             apply_lcmv_cov)
 from mne.minimum_norm import (make_inverse_operator, apply_inverse_cov)
@@ -31,8 +32,10 @@
 
 # %%
 # Reading the raw data and creating epochs:
+
 data_path = somato.data_path()
 subject = '01'
+subjects_dir = data_path / 'derivatives' / 'freesurfer' / 'subjects'
 task = 'somato'
 raw_fname = (data_path / 'sub-{}'.format(subject) / 'meg' /
              'sub-{}_task-{}_meg.fif'.format(subject, task))
@@ -53,15 +56,13 @@
                     preload=True, decim=3)
 
 # Read forward operator and point to freesurfer subject directory
-fname_fwd = (data_path / 'derivatives' / 'sub-{}'.format(subject) /
+fwd_fname = (data_path / 'derivatives' / 'sub-{}'.format(subject) /
              'sub-{}_task-{}-fwd.fif'.format(subject, task))
-subjects_dir = data_path / 'derivatives' / 'freesurfer' / 'subjects'
-
-fwd = mne.read_forward_solution(fname_fwd)
+fwd = mne.read_forward_solution(fwd_fname)
 
 # %%
-# Compute covariances
-# -------------------
+# Compute covariances and cross-spectral density
+# ----------------------------------------------
 # ERS activity starts at 0.5 seconds after stimulus onset. Because these
 # data have been processed by MaxFilter directly (rather than MNE-Python's
 # version), we have to be careful to compute the rank with a more conservative
@@ -70,17 +71,32 @@
 # will be correctly preserved.
 
 rank = mne.compute_rank(epochs, tol=1e-6, tol_kind='relative')
-active_win = (0.5, 1.5)
-baseline_win = (-1, 0)
-baseline_cov = compute_covariance(epochs, tmin=baseline_win[0],
-                                  tmax=baseline_win[1], method='shrunk',
+win_active = (0.5, 1.5)
+win_baseline = (-1, 0)
+cov_baseline = compute_covariance(epochs, tmin=win_baseline[0],
+                                  tmax=win_baseline[1], method='shrunk',
                                   rank=rank, verbose=True)
-active_cov = compute_covariance(epochs, tmin=active_win[0], tmax=active_win[1],
+cov_active = compute_covariance(epochs, tmin=win_active[0], tmax=win_active[1],
                                 method='shrunk', rank=rank, verbose=True)
 
-# Weighted averaging is already in the addition of covariance objects.
-common_cov = baseline_cov + active_cov
-mne.viz.plot_cov(baseline_cov, epochs.info)
+# Weighted averaging is already in the addition of covariance objects
+cov_common = cov_baseline + cov_active
+cov_baseline.plot(epochs.info)
+
+# compute cross-spectral density matrices
+freqs = np.logspace(np.log10(12), np.log10(30), 9)
+
+# time-frequency decomposition
+epochs_tfr = mne.time_frequency.tfr_morlet(
+    epochs, freqs=freqs, n_cycles=freqs / 2, return_itc=False,
+    average=False, output='complex')
+epochs_tfr.decimate(20)  # decimate for speed
+
+csd = csd_tfr(epochs_tfr, tmin=-1, tmax=1.5)
+csd_baseline = csd_tfr(epochs_tfr, tmin=win_baseline[0], tmax=win_baseline[1])
+csd_ers = csd_tfr(epochs_tfr, tmin=win_active[0], tmax=win_active[1])
+
+csd_baseline.plot()
 
 # %%
 # Compute some source estimates
@@ -90,13 +106,7 @@
 # See :ref:`ex-inverse-source-power` for more information about DICS.
 
 
-def _gen_dics(active_win, baseline_win, epochs):
-    freqs = np.logspace(np.log10(12), np.log10(30), 9)
-    csd = csd_morlet(epochs, freqs, tmin=-1, tmax=1.5, decim=20)
-    csd_baseline = csd_morlet(epochs, freqs, tmin=baseline_win[0],
-                              tmax=baseline_win[1], decim=20)
-    csd_ers = csd_morlet(epochs, freqs, tmin=active_win[0], tmax=active_win[1],
-                         decim=20)
+def _gen_dics(csd, ers_csd, csd_baseline, fwd):
     filters = make_dics(epochs.info, fwd, csd.mean(), pick_ori='max-power',
                         reduce_rank=True, real_filter=True, rank=rank)
     stc_base, freqs = apply_dics_csd(csd_baseline.mean(), filters)
@@ -106,30 +116,30 @@ def _gen_dics(active_win, baseline_win, epochs):
 
 
 # generate lcmv source estimate
-def _gen_lcmv(active_cov, baseline_cov, common_cov):
+def _gen_lcmv(active_cov, cov_baseline, common_cov, fwd):
     filters = make_lcmv(epochs.info, fwd, common_cov, reg=0.05,
                         noise_cov=None, pick_ori='max-power')
-    stc_base = apply_lcmv_cov(baseline_cov, filters)
-    stc_act = apply_lcmv_cov(active_cov, filters)
+    stc_base = apply_lcmv_cov(cov_baseline, filters)
+    stc_act = apply_lcmv_cov(cov_active, filters)
     stc_act /= stc_base
     return stc_act
 
 
 # generate mne/dSPM source estimate
-def _gen_mne(active_cov, baseline_cov, common_cov, fwd, info, method='dSPM'):
-    inverse_operator = make_inverse_operator(info, fwd, common_cov)
-    stc_act = apply_inverse_cov(active_cov, info, inverse_operator,
+def _gen_mne(cov_active, cov_baseline, cov_common, fwd, info, method='dSPM'):
+    inverse_operator = make_inverse_operator(info, fwd, cov_common)
+    stc_act = apply_inverse_cov(cov_active, info, inverse_operator,
                                 method=method, verbose=True)
-    stc_base = apply_inverse_cov(baseline_cov, info, inverse_operator,
+    stc_base = apply_inverse_cov(cov_baseline, info, inverse_operator,
                                  method=method, verbose=True)
     stc_act /= stc_base
     return stc_act
 
 
 # Compute source estimates
-stc_dics = _gen_dics(active_win, baseline_win, epochs)
-stc_lcmv = _gen_lcmv(active_cov, baseline_cov, common_cov)
-stc_dspm = _gen_mne(active_cov, baseline_cov, common_cov, fwd, epochs.info)
+stc_dics = _gen_dics(csd, csd_ers, csd_baseline, fwd)
+stc_lcmv = _gen_lcmv(cov_active, cov_baseline, cov_common, fwd)
+stc_dspm = _gen_mne(cov_active, cov_baseline, cov_common, fwd, epochs.info)
 
 # %%
 # Plot source estimates
@@ -153,3 +163,56 @@ def _gen_mne(active_cov, baseline_cov, common_cov, fwd, info, method='dSPM'):
 brain_dspm = stc_dspm.plot(
     hemi='rh', subjects_dir=subjects_dir, subject=subject,
     time_label='dSPM source power in the 12-30 Hz frequency band')
+
+# %%
+# Use volume source estimate with time-frequency resolution
+# ---------------------------------------------------------
+
+# make a volume source space
+surface = subjects_dir / subject / 'bem' / 'inner_skull.surf'
+vol_src = mne.setup_volume_source_space(
+    subject=subject, subjects_dir=subjects_dir, surface=surface,
+    pos=10, add_interpolator=False)  # just for speed!
+
+conductivity = (0.3,)  # one layer for MEG
+model = mne.make_bem_model(subject=subject, ico=3,  # just for speed
+                           conductivity=conductivity,
+                           subjects_dir=subjects_dir)
+bem = mne.make_bem_solution(model)
+
+trans = fwd['info']['mri_head_t']
+vol_fwd = mne.make_forward_solution(
+    raw.info, trans=trans, src=vol_src, bem=bem, meg=True, eeg=True,
+    mindist=5.0, n_jobs=1, verbose=True)
+
+# Compute source estimate using MNE solver
+snr = 3.0
+lambda2 = 1.0 / snr ** 2
+method = 'MNE'  # use MNE method (could also be dSPM or sLORETA)
+
+# make a different inverse operator for each frequency so as to properly
+# whiten the sensor data
+inverse_operator = list()
+for freq_idx in range(epochs_tfr.freqs.size):
+    # for each frequency, compute a separate covariance matrix
+    cov_baseline = csd_baseline.get_data(index=freq_idx, as_cov=True)
+    cov_baseline['data'] = cov_baseline['data'].real  # only normalize by real
+    # then use that covariance matrix as normalization for the inverse
+    # operator
+    inverse_operator.append(mne.minimum_norm.make_inverse_operator(
+        epochs.info, vol_fwd, cov_baseline))
+
+# finally, compute the stcs for each epoch and frequency
+stcs = mne.minimum_norm.apply_inverse_tfr_epochs(
+    epochs_tfr, inverse_operator, lambda2, method=method,
+    pick_ori='vector')
+
+# %%
+# Plot volume source estimates
+# ----------------------------
+
+viewer = mne.gui.view_vol_stc(stcs, subject=subject, subjects_dir=subjects_dir,
+                              src=vol_src, inst=epochs_tfr)
+viewer.go_to_max()  # show the maximum intensity source vertex
+viewer.set_cmap(vmin=0.25, vmid=0.8)
+viewer.set_3d_view(azimuth=40, elevation=35, distance=350)

mne/gui/__init__.py

@@ -1,9 +1,11 @@
 """Convenience functions for opening GUIs."""
 
 # Authors: Christian Brodbeck <christianbrodbeck@nyu.edu>
+#          Alex Rockhill <aprockhill@mailbox.org>
 #
 # License: BSD-3-Clause
 
+import numpy as np
 from ..utils import verbose, get_config, warn
 
 
@@ -247,6 +249,97 @@ def locate_ieeg(info, trans, aligned_ct, subject=None, subjects_dir=None,
     return gui
 
 
+@verbose
+def view_vol_stc(stcs, freq_first=True, subject=None, subjects_dir=None,
+                 src=None, inst=None, show=True, block=False, verbose=None):
+    """View a volume time and/or frequency source time course estimate.
+
+    Parameters
+    ----------
+    stcs : list of list | generator
+        The source estimates. List of lists or generators for epochs
+        and frequencies (i.e. using
+        :func:`mne.minimum_norm.apply_inverse_tfr_epochs` or
+        :func:`mne.beamformer.apply_dics_tfr_epochs`-- in this case
+        use ``freq_first=False``). Lists of source estimates across
+        frequencies (e.g. :func::func:`mne.beamformer.apply_dics_csd`)
+        and lists of source estimates across epochs
+        (e.g. :func:`mne.minimum_norm.apply_inverse_epochs` and
+        :func:`mne.beamformer.apply_dics_epochs`--in these
+        case use ``freq_first=False``) are also allowed. Single
+        source estimates (e.g. :func:`mne.minimum_norm.apply_inverse`
+        and :func:`mne.beamformer.apply_dics`) are also allowed
+        (``freq_first`` will not be used in this case).
+    freq_first : bool
+        If frequencies are the outer list of ``stcs`` use ``True``.
+    %(subject)s
+    %(subjects_dir)s
+    src : instance of SourceSpaces
+        The volume source space for the ``stc``.
+    inst : EpochsTFR | AverageTFR | None
+        The time-frequency or data object to use to plot topography.
+    show : bool
+        Show the GUI if True.
+    block : bool
+        Whether to halt program execution until the figure is closed.
+    %(verbose)s
+
+    Returns
+    -------
+    gui : instance of VolSourceEstimateViewer
+        The graphical user interface (GUI) window.
+    """
+    from ..viz.backends._utils import _qt_app_exec
+    from ._vol_stc import VolSourceEstimateViewer, COMPLEX_DTYPE, RANGE_VALUE
+    from qtpy.QtWidgets import QApplication
+    # get application
+    app = QApplication.instance()
+    if app is None:
+        app = QApplication(['Source Estimate Viewer'])
+
+    # cast to integers to lower memory usage, use custom complex data
+    # type if necessary
+    data = list()
+    # can be generator, compute using first stc object, just a general
+    # rescaling of data, does not need to be precise
+    scalar = None
+    for inner_stcs in (stcs if np.iterable(stcs) else [stcs]):
+        inner_data = list()
+        for stc in (inner_stcs if np.iterable(inner_stcs) else [inner_stcs]):
+            if np.iscomplexobj(stc.data):
+                if scalar is None:
+                    # this is an order of magnitude approximation,
+                    # larger stcs will have some clipping
+                    scalar = (RANGE_VALUE - 1) / stc.data.real.max()
+                stc_data = np.zeros(stc.data.shape, COMPLEX_DTYPE)
+                stc_data['re'] = np.clip(stc.data.real * scalar,
+                                         -RANGE_VALUE, RANGE_VALUE - 1)
+                stc_data['im'] = np.clip(stc.data.imag * scalar,
+                                         -RANGE_VALUE, RANGE_VALUE - 1)
+                inner_data.append(stc_data)
+            else:
+                if scalar is None:
+                    scalar = (RANGE_VALUE - 1) / stc.data.max() / 10
+                inner_data.append(np.clip(stc.data * scalar,
+                                          -RANGE_VALUE, RANGE_VALUE - 1
+                                          ).astype(np.int16))
+        data.append(inner_data)
+    data = np.array(data)
+    if data.ndim == 4:  # scalar solution, add dimension at the end
+        data = data[:, :, :, None]
+
+    # move frequencies to penultimate
+    data = data.transpose((1, 2, 3, 0, 4) if freq_first else (0, 2, 3, 1, 4))
+
+    gui = VolSourceEstimateViewer(
+        data, subject=subject, subjects_dir=subjects_dir,
+        src=src, inst=inst, show=show,
+        verbose=verbose)
+    if block:
+        _qt_app_exec(app)
+    return gui
+
+
 class _GUIScraper(object):
     """Scrape GUI outputs."""
 
@@ -256,11 +349,13 @@ def __repr__(self):
     def __call__(self, block, block_vars, gallery_conf):
         from ._ieeg_locate import IntracranialElectrodeLocator
         from ._coreg import CoregistrationUI
+        from ._vol_stc import VolSourceEstimateViewer
         from sphinx_gallery.scrapers import figure_rst
         from qtpy import QtGui
         for gui in block_vars['example_globals'].values():
             if (isinstance(gui, (IntracranialElectrodeLocator,
-                                 CoregistrationUI)) and
+                                 CoregistrationUI,
+                                 VolSourceEstimateViewer)) and
                     not getattr(gui, '_scraped', False) and
                     gallery_conf['builder_name'] == 'html'):
                 gui._scraped = True  # monkey-patch but it's easy enough