[WIP] Trying to allow plotting SEEG activity in volumetric source spa…

…ce (mne-tools#8402) * Trying to fix up source estimate for seeg. * Almost there? * Try swap out for volume source space. * Update mne/source_estimate.py Co-authored-by: Eric Larson <[email protected]> * Merging. * Update mne/source_estimate.py Co-authored-by: Eric Larson <[email protected]> * Update mne/source_estimate.py Co-authored-by: Eric Larson <[email protected]> * Update tutorials/misc/plot_seeg.py Co-authored-by: Eric Larson <[email protected]> * Update tutorials/misc/plot_ecog.py Co-authored-by: Eric Larson <[email protected]> * Update tutorials/misc/plot_seeg.py Co-authored-by: Eric Larson <[email protected]> * Update tutorials/misc/plot_seeg.py Co-authored-by: Eric Larson <[email protected]> * Update mne/source_estimate.py Co-authored-by: Eric Larson <[email protected]> * Update tutorials/misc/plot_seeg.py Co-authored-by: Eric Larson <[email protected]> * Removing nans and fixing units of coords. * Fixing example. * Fixing example. * Fixing flake. * Fixing flake. * Fixing flake. * FIX: Fix example * FIX: For real * FIX: Remove comments * FIX: Links * FIX: NaN * Adding subjects dir. * Fixing example. * Adding updated plot seeg. * Updating md5 and misc vs. * Fixing merge to add back larson's test. * ADding latest.inc entry. * Fixing up docstring. * Update tutorials/misc/plot_seeg.py Co-authored-by: Eric Larson <[email protected]> * Clean up test_source_estimate. * Update tutorials/misc/plot_seeg.py Co-authored-by: Eric Larson <[email protected]> * FIX: Order * FIX: No filename * FIX: Less memory * Updating docstring. * move up fs avg. * move up fs avg. * FIX: Space * move up fs avg. * Adding space back. * Fixing docstring. * Fixing flake. * ENH: Better thumb Co-authored-by: Eric Larson <[email protected]>
jordandc1 · Oct 26, 2020 · cdeb443 · cdeb443
1 parent 360e62d
commit cdeb443
Show file tree

Hide file tree

Showing 9 changed files with 316 additions and 36 deletions.
diff --git a/doc/Makefile b/doc/Makefile
@@ -5,7 +5,7 @@
 SPHINXOPTS    = -nWT --keep-going
 SPHINXBUILD   = sphinx-build
 PAPER         =
-MPROF         = SG_STAMP_STARTS=true mprof run --python sphinx
+MPROF         = SG_STAMP_STARTS=true mprof run -E --python sphinx
 
 # Internal variables.
 PAPEROPT_a4     = -D latex_paper_size=a4

diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc
@@ -27,6 +27,8 @@ Enhancements
 
 - Speed up :class:`mne.decoding.TimeDelayingRidge` with edge correction using Numba by `Eric Larson`_ (:gh:`8323`)
 
+- Add sEEG source visualization using :func:`mne.stc_near_sensors` and sEEG working tutorial by `Eric Larson`_ and `Adam Li`_ (:gh:`8402`)
+
 Bugs
 ~~~~
 - Fix bug with reading split files that have dashes in the filename **by new contributor** |Eduard Ort|_ (:gh:`8339`)

diff --git a/mne/datasets/utils.py b/mne/datasets/utils.py
@@ -240,7 +240,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.108', misc='0.6')
+    releases = dict(testing='0.108', misc='0.7')
     # And also update the "md5_hashes['testing']" variable below.
     # To update any other dataset, update the data archive itself (upload
     # an updated version) and update the md5 hash.
@@ -322,7 +322,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
             bst_raw='fa2efaaec3f3d462b319bc24898f440c',
             bst_resting='70fc7bf9c3b97c4f2eab6260ee4a0430'),
         fake='3194e9f7b46039bb050a74f3e1ae9908',
-        misc='e00808c3b05123059e2cf49ff276e919',
+        misc='2b2f2fec9d1197ed459117db1c6341ee',
         sample='12b75d1cb7df9dfb4ad73ed82f61094f',
         somato='32fd2f6c8c7eb0784a1de6435273c48b',
         spm='9f43f67150e3b694b523a21eb929ea75',

diff --git a/mne/source_estimate.py b/mne/source_estimate.py
@@ -17,6 +17,7 @@
 from .cov import Covariance
 from .evoked import _get_peak
 from .filter import resample
+from .io.constants import FIFF
 from .surface import (read_surface, _get_ico_surface, mesh_edges,
                       _project_onto_surface)
 from .source_space import (_ensure_src, _get_morph_src_reordering,
@@ -3176,13 +3177,13 @@ def extract_label_time_course(stcs, labels, src, mode='auto',
 
 @verbose
 def stc_near_sensors(evoked, trans, subject, distance=0.01, mode='sum',
-                     project=True, subjects_dir=None, verbose=None):
-    """Create a STC from ECoG sensor data.
+                     project=True, subjects_dir=None, src=None, verbose=None):
+    """Create a STC from ECoG and sEEG sensor data.
 
     Parameters
     ----------
     evoked : instance of Evoked
-        The evoked data. Must contain ECoG channels.
+        The evoked data. Must contain ECoG, or sEEG channels.
     %(trans)s
     subject : str
         The subject name.
@@ -3194,20 +3195,33 @@ def stc_near_sensors(evoked, trans, subject, distance=0.01, mode='sum',
         zero-order hold. See Notes.
     project : bool
         If True, project the electrodes to the nearest ``'pial`` surface
-        vertex before computing distances.
+        vertex before computing distances. Only used when doing a
+        surface projection.
     %(subjects_dir)s
+    src : instance of SourceSpaces
+        The source space.
+
+        .. warning:: If a surface source space is used, make sure that
+                     ``surf='pial'`` was used during construction.
     %(verbose)s
 
     Returns
     -------
     stc : instance of SourceEstimate
-        The surface source estimate.
+        The surface source estimate. If src is None, a surface source
+        estimate will be produced, and the number of vertices will equal
+        the number of pial-surface vertices that were close enough to
+        the sensors to take on a non-zero volue. If src is not None,
+        a surface, volume, or mixed source estimate will be produced
+        (depending on the kind of source space passed) and the
+        vertices will match those of src (i.e., there may be me
+        many all-zero values in stc.data).
 
     Notes
     -----
-    This function projects the ECoG sensors to the pial surface (if
-    ``project``), then the activation at each pial surface vertex is given
-    by the mode:
+    For surface projections, this function projects the ECoG sensors to
+    the pial surface (if ``project``), then the activation at each pial
+    surface vertex is given by the mode:
 
     - ``'sum'``
         Activation is the sum across each sensor weighted by the fractional
@@ -3224,37 +3238,84 @@ def stc_near_sensors(evoked, trans, subject, distance=0.01, mode='sum',
         The value is given by the value of the nearest sensor, up to a
         ``distance`` (beyond which it is zero).
 
-    .. versionadded:: 0.21
+    If creating a Volume STC, ``src`` must be passed in, and this
+    function will project sEEG sensors to nearby surrounding vertices.
+    Then the activation at each volume vertex is given by the mode
+    in the same way as ECoG surface projections.
+
+    .. versionadded:: 0.22
     """
     from scipy.spatial.distance import cdist, pdist
     from .evoked import Evoked
     _validate_type(evoked, Evoked, 'evoked')
     _validate_type(mode, str, 'mode')
+    _validate_type(src, (None, SourceSpaces), 'src')
     _check_option('mode', mode, ('sum', 'single', 'nearest'))
-    evoked = evoked.copy().pick_types(meg=False, ecog=True)
-    pos = np.array([ch['loc'][:3] for ch in evoked.info['chs']])
+
+    # create a copy of Evoked using ecog and seeg
+    evoked = evoked.copy().pick_types(ecog=True, seeg=True)
+
+    # get channel positions that will be used to pinpoint where
+    # in the Source space we will use the evoked data
+    pos = evoked._get_channel_positions()
+
+    # remove nan channels
+    nan_inds = np.where(np.isnan(pos).any(axis=1))[0]
+    nan_chs = [evoked.ch_names[idx] for idx in nan_inds]
+    evoked.drop_channels(nan_chs)
+    pos = [pos[idx] for idx in range(len(pos)) if idx not in nan_inds]
+
+    # coord_frame transformation from native mne "head" to MRI coord_frame
     trans, _ = _get_trans(trans, 'head', 'mri', allow_none=True)
+
+    # convert head positions -> coord_frame MRI
+    pos = apply_trans(trans, pos)
+
     subject = _check_subject(None, subject, False)
     subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
-    data, vertices = list(), list()
-    rrs = [read_surface(op.join(subjects_dir, subject,
-                                'surf', f'{hemi}.pial'))[0]
-           for hemi in ('lh', 'rh')]
-    offset = len(rrs[0])
-    rrs = np.concatenate(rrs)
-    rrs /= 1000.
-    pos = apply_trans(trans, pos)
+    if src is None:  # fake a full surface one
+        rrs = [read_surface(op.join(subjects_dir, subject,
+                                    'surf', f'{hemi}.pial'))[0]
+               for hemi in ('lh', 'rh')]
+        src = SourceSpaces([
+            dict(rr=rr / 1000., vertno=np.arange(len(rr)), type='surf',
+                 coord_frame=FIFF.FIFFV_COORD_MRI)
+            for rr in rrs])
+        del rrs
+        keep_all = False
+    else:
+        keep_all = True
+    # ensure it's a usable one
+    klass = dict(
+        surface=SourceEstimate,
+        volume=VolSourceEstimate,
+        mixed=MixedSourceEstimate,
+    )
+    _check_option('src.kind', src.kind, sorted(klass.keys()))
+    klass = klass[src.kind]
+    rrs = np.concatenate([s['rr'][s['vertno']] for s in src])
+    if src[0]['coord_frame'] == FIFF.FIFFV_COORD_HEAD:
+        rrs = apply_trans(trans, rrs)
+    # projection will only occur with surfaces
     logger.info(
-        f'Projecting {len(pos)} sensors onto {len(rrs)} vertices: {mode} mode')
-    if project:
+        f'Projecting data from {len(pos)} sensor{_pl(pos)} onto {len(rrs)} '
+        f'{src.kind} vertices: {mode} mode')
+    if project and src.kind == 'surface':
         logger.info('    Projecting electrodes onto surface')
         pos = _project_onto_surface(pos, dict(rr=rrs), project_rrs=True,
                                     method='nearest')[2]
+
     min_dist = pdist(pos).min() * 1000
     logger.info(
-        f'    Projected sensors have minimum distance {min_dist:0.1f} mm')
+        f'    Minimum {"projected " if project else ""}intra-sensor distance: '
+        f'{min_dist:0.1f} mm')
+
+    # compute pairwise distance between source space points and sensors
     dists = cdist(rrs, pos)
     assert dists.shape == (len(rrs), len(pos))
+
+    # only consider vertices within our "epsilon-ball"
+    # characterized by distance kwarg
     vertices = np.where((dists <= distance).any(-1))[0]
     logger.info(f'    {len(vertices)} / {len(rrs)} non-zero vertices')
     w = np.maximum(1. - dists[vertices] / distance, 0)
@@ -3269,9 +3330,20 @@ def stc_near_sensors(evoked, trans, subject, distance=0.01, mode='sum',
     if len(missing):
         warn(f'Channel{_pl(missing)} missing in STC: '
              f'{", ".join(evoked.ch_names[mi] for mi in missing)}')
-    data = w @ evoked.data
-    vertices = [vertices[vertices < offset],
-                vertices[vertices >= offset] - offset]
-    return SourceEstimate(
-        data, vertices, evoked.times[0], 1. / evoked.info['sfreq'],
-        subject=subject)
+
+    nz_data = w @ evoked.data
+    if not keep_all:
+        assert src.kind == 'surface'
+        data = nz_data
+        offset = len(src[0]['vertno'])
+        vertices = [vertices[vertices < offset],
+                    vertices[vertices >= offset] - offset]
+    else:
+        data = np.zeros(
+            (sum(len(s['vertno']) for s in src), len(evoked.times)),
+            dtype=nz_data.dtype)
+        data[vertices] = nz_data
+        vertices = [s['vertno'].copy() for s in src]
+
+    return klass(data, vertices, evoked.times[0], 1. / evoked.info['sfreq'],
+                 subject=subject, verbose=verbose)
diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py
@@ -26,15 +26,16 @@
                  spatial_src_adjacency, spatial_tris_adjacency, pick_info,
                  SourceSpaces, VolVectorSourceEstimate, read_trans, pick_types,
                  MixedVectorSourceEstimate, setup_volume_source_space,
-                 convert_forward_solution, pick_types_forward)
+                 convert_forward_solution, pick_types_forward,
+                 compute_source_morph)
 from mne.datasets import testing
 from mne.externals.h5io import write_hdf5
 from mne.fixes import fft, _get_img_fdata
 from mne.io import read_info
 from mne.io.constants import FIFF
 from mne.source_estimate import grade_to_tris, _get_vol_mask
 from mne.source_space import _get_src_nn
-from mne.transforms import apply_trans, invert_transform
+from mne.transforms import apply_trans, invert_transform, transform_surface_to
 from mne.minimum_norm import (read_inverse_operator, apply_inverse,
                               apply_inverse_epochs, make_inverse_operator)
 from mne.label import read_labels_from_annot, label_sign_flip
@@ -1474,6 +1475,7 @@ def test_stc_near_sensors(tmpdir):
     assert info['nchan'] == 33
     evoked = EvokedArray(np.eye(info['nchan']), info)
     trans = read_trans(fname_fwd)
+    assert trans['to'] == FIFF.FIFFV_COORD_HEAD
     this_dir = str(tmpdir)
     # testing does not have pial, so fake it
     os.makedirs(op.join(this_dir, 'sample', 'surf'))
@@ -1489,14 +1491,32 @@ def test_stc_near_sensors(tmpdir):
     with catch_logging() as log:
         stc = stc_near_sensors(evoked, **kwargs)
     log = log.getvalue()
-    assert 'minimum distance 7.' in log  # 7.4
+    assert 'Minimum projected intra-sensor distance: 7.' in log  # 7.4
     # this should be left-hemisphere dominant
     assert 5000 > len(stc.vertices[0]) > 4000
     assert 200 > len(stc.vertices[1]) > 100
     # and at least one vertex should have the channel values
     dists = cdist(stc.data, evoked.data)
     assert np.isclose(dists, 0., atol=1e-6).any(0).all()
 
+    src = read_source_spaces(fname_src)  # uses "white" but should be okay
+    for s in src:
+        transform_surface_to(s, 'head', trans, copy=False)
+    assert src[0]['coord_frame'] == FIFF.FIFFV_COORD_HEAD
+    stc_src = stc_near_sensors(evoked, src=src, **kwargs)
+    assert len(stc_src.data) == 7928
+    with pytest.warns(RuntimeWarning, match='not included'):  # some removed
+        stc_src_full = compute_source_morph(
+            stc_src, 'sample', 'sample', smooth=5, spacing=None,
+            subjects_dir=subjects_dir).apply(stc_src)
+    lh_idx = np.searchsorted(stc_src_full.vertices[0], stc.vertices[0])
+    rh_idx = np.searchsorted(stc_src_full.vertices[1], stc.vertices[1])
+    rh_idx += len(stc_src_full.vertices[0])
+    sub_data = stc_src_full.data[np.concatenate([lh_idx, rh_idx])]
+    assert sub_data.shape == stc.data.shape
+    corr = np.corrcoef(stc.data.ravel(), sub_data.ravel())[0, 1]
+    assert 0.6 < corr < 0.7
+
     # now single-weighting mode
     stc_w = stc_near_sensors(evoked, mode='single', **kwargs)
     assert_array_less(stc_w.data, stc.data + 1e-3)  # some tol
@@ -1518,3 +1538,13 @@ def test_stc_near_sensors(tmpdir):
         stc = stc_near_sensors(evoked, trans, 'sample', subjects_dir=this_dir,
                                project=False, distance=0.033)
     assert stc.data.any(0).sum() == len(evoked.ch_names) - 1
+
+    # and now with volumetric projection
+    src = read_source_spaces(fname_vsrc)
+    with catch_logging() as log:
+        stc_vol = stc_near_sensors(evoked, trans, 'sample', src=src,
+                                   subjects_dir=subjects_dir, verbose=True,
+                                   distance=0.033)
+    assert isinstance(stc_vol, VolSourceEstimate)
+    log = log.getvalue()
+    assert '4157 volume vertices' in log
diff --git a/mne/viz/_3d.py b/mne/viz/_3d.py
@@ -949,8 +949,11 @@ def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
         if other_bools[key] and len(picks):
             title = DEFAULTS["titles"][key] if key != 'fnirs' else 'fNIRS'
             if key != 'fnirs' or 'channels' in fnirs:
-                other_loc[key] = np.array([info['chs'][pick]['loc'][:3]
-                                           for pick in picks])
+                other_loc[key] = [
+                    info['chs'][pick]['loc'][:3] for pick in picks]
+                # deal with NaN
+                other_loc[key] = np.array([loc for loc in other_loc[key]
+                                           if np.isfinite(loc).all()], float)
                 logger.info(
                     f'Plotting {len(other_loc[key])} {title}'
                     f' location{_pl(other_loc[key])}')

diff --git a/mne/viz/_brain/_scraper.py b/mne/viz/_brain/_scraper.py
@@ -51,6 +51,7 @@ def __call__(self, block, block_vars, gallery_conf):
                                      ('time_viewer', False)]:
                     if key not in kwargs:
                         kwargs[key] = default
+                kwargs.pop('filename', None)  # always omit this one
                 if brain.time_viewer:
                     assert kwargs['time_viewer'], 'Must use time_viewer=True'
                 frames = brain._make_movie_frames(callback=None, **kwargs)

diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py
@@ -8,6 +8,15 @@
 MNE supports working with more than just MEG and EEG data. Here we show some
 of the functions that can be used to facilitate working with
 electrocorticography (ECoG) data.
+
+This example shows how to use:
+
+- ECoG data
+- channel locations in subject's MRI space
+- projection onto a surface
+
+For an example that involves sEEG data, channel locations in
+MNI space, or projection into a volume, see :ref:`tut_working_with_seeg`.
 """
 # Authors: Eric Larson <[email protected]>
 #          Chris Holdgraf <[email protected]>
-Original file line number
+Diff line change
@@ Expand Up / @@ -27,6 +27,8 @@ Enhancements @@
     - Speed up :class:`mne.decoding.TimeDelayingRidge` with edge correction using Numba by `Eric Larson`_ (:gh:`8323`)
+    - Add sEEG source visualization using :func:`mne.stc_near_sensors` and sEEG working tutorial by `Eric Larson`_ and `Adam Li`_ (:gh:`8402`)
     Bugs
     ~~~~
     - Fix bug with reading split files that have dashes in the filename **by new contributor** |Eduard Ort|_ (:gh:`8339`)
@@ Expand Down @@