CSD with working example and tests

added some tests WIP CSD MVP example working with example weird cherry-picking error again Update examples/preprocessing/plot_csd.py Co-Authored-By: Alexandre Gramfort <[email protected]> Update examples/preprocessing/plot_csd.py Co-Authored-By: Alexandre Gramfort <[email protected]> Update examples/preprocessing/plot_csd.py Co-Authored-By: Alexandre Gramfort <[email protected]> Update mne/preprocessing/_csd.py Co-Authored-By: Alexandre Gramfort <[email protected]> Update examples/preprocessing/plot_csd.py Co-Authored-By: Alexandre Gramfort <[email protected]> typos/small fixes incorporated Eric comments more Eric comments refactoring core comp per Eric realized I did that wrong Eric comments 3 plot fix Fix julian date conversions (mne-tools#6911) * Fix julian date conversions * remove tabs * use UTC and private functions. closes mne-tools#6901 [skip travis][skip azp] (mne-tools#6914) ENH: Add source space SNR estimate (mne-tools#6908) * source-snr-edit * FIX move estimate_snr to be a method * ENH: name change * Make it work * SNR works, now clean up code * fixed unnecessary stc copy, returns stc object * added plot source space SNR tutorial file * updated author names * fixed period for pydocsty and suggested documenation edit * edits suggested from @drammock * made edits suggested by reviewers * made edits suggested by reviewers * N -> n_channels * Fixed formatting of Reference in comment * fixed formatting in doc and in code * fixed snr_stc creation * further minor edits * fixed formatting issues and colormap choice in eg * removed colon after Reference * Corrected SNR calculation to take unwhitened sensor signal b instead of whitened, updated tutorial file * Minor changes to formatting * FIX: transparent * FIX: Rename * FIX? Added dSPM (mne-tools#6913) * Added dSPM * removed term signal covariance * Minor modification of dSPM definition * Corrected according to Mattis suggestions [MRG]: read *.dat electrode position files (mne-tools#6909) * ADD: channels.read_dig_dat() * DOC * DOC * TEST * DOC, __all__ * DOC: intro * DOC * pep & py3.5 compatibility * py3.5 MRG: allow users to annotate epochs with colors (mne-tools#6855) * ENH: copy over diff from autoreject * Make it work * More fixes * More fixes * Make horizontal patch * Cleanup * TST: add tests * Address Eric's comments * Address alex comments * DOC: update whats new * FIX travis + circle * FIX For ICA plot Fix infomax (mne-tools#6921) * fix extended infomax * update whats new * cleaner way, using setdefault DOC: Document authors better [ci skip] FIX: Comments and order umlauts MRG, MAINT: Simplify code for vol ravel (mne-tools#6920) * MAINT: Simplify code for vol ravel * STY: Flake ENH: Convert NIRS raw data to optical density (mne-tools#6827) * Initial framework for fnirs raw to optical density conversion * Add optical density to docs * Update mne/preprocessing/tests/test_optical_density.py Co-Authored-By: Eric Larson <[email protected]> * Update mne/preprocessing/tests/test_optical_density.py Co-Authored-By: Eric Larson <[email protected]> * Import testing * Put optical density docs in preprocessing section * Implement optical density conversion * Fix optical density code style issues * Add plotting of optical density data * Add defaults for optical density plotting * Update mne/preprocessing/optical_density.py Co-Authored-By: Eric Larson <[email protected]> * Handle negative intensities in optical density conversion * Add in -1* that was removed in previous commit * Add tests and fix imports for optical density * Add tests for optical density conversion * Use assert_allclose for testing optical density * Update test_optical_density.py * Update mne/preprocessing/_optical_density.py Co-Authored-By: Alexandre Gramfort <[email protected]> * Do not operate in place for optical density * Force data to load and modify test to include preload=False condition import fix tidying up weird file conflict errors more unexpected file changes fix weird off target changes doc docstring fixes docstyle reformatted plots changes formatting fix Update mne/preprocessing/_csd.py Co-Authored-By: Eric Larson <[email protected]> Update examples/preprocessing/plot_csd.py Co-Authored-By: Eric Larson <[email protected]> Update mne/preprocessing/tests/test_csd.py Co-Authored-By: Eric Larson <[email protected]> Update mne/preprocessing/_csd.py Co-Authored-By: Eric Larson <[email protected]> Update mne/preprocessing/_csd.py Co-Authored-By: Eric Larson <[email protected]> Update mne/preprocessing/_csd.py Co-Authored-By: Eric Larson <[email protected]> halfway through Eric suggested changes getting on plane WIP Eric comments almost done testing fixes updated fif constants csd formatting small change references
alexrockhill · Nov 4, 2019 · 1f0dbe1 · 1f0dbe1
1 parent daa162c
commit 1f0dbe1
Show file tree

Hide file tree

Showing 9 changed files with 457 additions and 11 deletions.
diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc
@@ -14,6 +14,8 @@ Current (0.20.dev0)
 
 Changelog
 ~~~~~~~~~
+- Add current source density :func:`mne.preprocessing.compute_current_source_density` to compute the surface Laplacian in order to reduce volume conduction in data by `Alex Rockhill`_
+
 - Add command :ref:`gen_mne_setup_source_space` to quickly set up bilateral hemisphere surface-based source space with subsampling by `Victor Ferat`_.
 
 - Add function :func:`mne.make_fixed_length_epochs` to segment raw into fixed length epochs by `Mohammad Daneshzand`_

diff --git a/doc/changes/names.inc b/doc/changes/names.inc
@@ -277,3 +277,5 @@
 .. _Victor Ferat: https://github.com/vferat
 
 .. _Demetres Kostas: https://github.com/kostasde
+
+.. _Alex Rockhill: https://github.com/alexrockhill/
diff --git a/doc/python_reference.rst b/doc/python_reference.rst
@@ -343,6 +343,7 @@ Projections:
 
    ICA
    Xdawn
+   compute_current_source_density
    compute_proj_ecg
    compute_proj_eog
    create_ecg_epochs

diff --git a/examples/preprocessing/plot_csd.py b/examples/preprocessing/plot_csd.py
@@ -0,0 +1,91 @@
+"""
+===============================================
+Analyze Data Using Current Source Density (CSD)
+===============================================
+
+This script shows an example of how to use CSD. CSD
+takes the spatial Laplacian of the sensor signal
+(derivative in both x and y). It does what a planar
+gradiometer does in MEG. Spatial derivative
+reduces point spread. CSD transformed data have a
+sharper or more distinct topography, reducing the
+negative impact of volume conduction.
+
+"""
+# Authors: Alex Rockhill <[email protected]>
+#
+# License: BSD (3-clause)
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import mne
+from mne.datasets import sample
+
+print(__doc__)
+
+data_path = sample.data_path()
+
+###############################################################################
+# Load sample subject data
+raw = mne.io.read_raw_fif(data_path + '/MEG/sample/sample_audvis_raw.fif',
+                          preload=True)
+raw = raw.pick_types(meg=False, eeg=True, eog=True, ecg=True, stim=True,
+                     exclude=raw.info['bads'])
+event_id = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3,
+            'visual/right': 4, 'smiley': 5, 'button': 32}
+events = mne.find_events(raw)
+epochs = mne.Epochs(raw, events, event_id=event_id, tmin=-0.2, tmax=.5,
+                    preload=True)
+evo = epochs['auditory'].average()
+
+###############################################################################
+# Let's look at the topography of CSD compared to average
+fig, axes = plt.subplots(1, 6)
+times = np.array([-0.1, 0., 0.05, 0.1, 0.15])
+evo.plot_topomap(times=times, cmap='Spectral_r', axes=axes[:5],
+                 outlines='skirt', contours=4, time_unit='s',
+                 colorbar=True, show=False, title='Average Reference')
+fig, axes = plt.subplots(1, 6)
+csd_evo = mne.preprocessing.compute_current_source_density(evo, copy=True)
+csd_evo.plot_topomap(times=times, axes=axes[:5], cmap='Spectral_r',
+                     outlines='skirt', contours=4, time_unit='s',
+                     colorbar=True, title='Current Source Density')
+
+###############################################################################
+# Plot evoked
+evo.plot_joint(title='Average Reference', show=False)
+csd_evo.plot_joint(title='Current Source Density')
+
+###############################################################################
+# Look at the effect of smoothing and spline flexibility
+fig, ax = plt.subplots(4, 4)
+fig.set_size_inches(10, 10)
+for i, lambda2 in enumerate([0, 1e-7, 1e-5, 1e-3]):
+    for j, m in enumerate([5, 4, 3, 2]):
+        this_csd_evo = \
+            mne.preprocessing.compute_current_source_density(evo, stiffness=m,
+                                                             lambda2=lambda2,
+                                                             copy=True)
+        this_csd_evo.plot_topomap(0.1, axes=ax[i, j],
+                                  outlines='skirt', contours=4, time_unit='s',
+                                  colorbar=False, show=False)
+        ax[i, j].set_title('m=%i, lambda=%s' % (m, lambda2))
+
+plt.show()
+
+# References
+# ----------
+#
+# [1] Perrin F, Bertrand O, Pernier J. "Scalp current density mapping:
+#     Value and estimation from potential data." IEEE Trans Biomed Eng.
+#     1987;34(4):283–288.
+#
+# [2] Perrin F, Pernier J, Bertrand O, Echallier JF. "Spherical splines
+#     for scalp potential and current density mapping."
+#     [Corrigenda EEG 02274, EEG Clin. Neurophysiol., 1990, 76, 565]
+#     Electroenceph Clin Neurophysiol. 1989;72(2):184–187.
+#
+# [3] Kayser J, Tenke CE. "On the benefits of using surface Laplacian
+#     (Current Source Density) methodology in electrophysiology."
+#     Int J Psychophysiol. 2015 Sep; 97(3): 171–173.
diff --git a/mne/channels/interpolation.py b/mne/channels/interpolation.py
@@ -13,6 +13,25 @@
 from ..forward import _map_meg_channels
 
 
+def _calc_h(cosang, stiffness=4, num_lterms=50):
+    """Calculate spherical spline h function between points on a sphere.
+
+    Parameters
+    ----------
+    cosang : array-like | float
+        cosine of angles between pairs of points on a spherical surface. This
+        is equivalent to the dot product of unit vectors.
+    stiffness : float
+        stiffnes of the spline. Also referred to as `m`.
+    num_lterms : int
+        number of Legendre terms to evaluate.
+    """
+    factors = [(2 * n + 1) /
+               (n ** (stiffness - 1) * (n + 1) ** (stiffness - 1) * 4 * np.pi)
+               for n in range(1, num_lterms + 1)]
+    return legval(cosang, [0] + factors)
+
+
 def _calc_g(cosang, stiffness=4, num_lterms=50):
     """Calculate spherical spline g function between points on a sphere.
 

diff --git a/mne/io/pick.py b/mne/io/pick.py
@@ -477,17 +477,18 @@ def pick_info(info, sel=(), copy=True, verbose=None):
     info['bads'] = [ch for ch in info['bads'] if ch in info['ch_names']]
 
     if 'comps' in info:
-        comps = deepcopy(info['comps'])
-        for c in comps:
-            row_idx = [k for k, n in enumerate(c['data']['row_names'])
-                       if n in info['ch_names']]
-            row_names = [c['data']['row_names'][i] for i in row_idx]
-            rowcals = c['rowcals'][row_idx]
-            c['rowcals'] = rowcals
-            c['data']['nrow'] = len(row_names)
-            c['data']['row_names'] = row_names
-            c['data']['data'] = c['data']['data'][row_idx]
-        info['comps'] = comps
+        if 'csd' not in info['comps']:  # added for csd 0.20
+            comps = deepcopy(info['comps'])
+            for c in comps:
+                row_idx = [k for k, n in enumerate(c['data']['row_names'])
+                           if n in info['ch_names']]
+                row_names = [c['data']['row_names'][i] for i in row_idx]
+                rowcals = c['rowcals'][row_idx]
+                c['rowcals'] = rowcals
+                c['data']['nrow'] = len(row_names)
+                c['data']['row_names'] = row_names
+                c['data']['data'] = c['data']['data'][row_idx]
+            info['comps'] = comps
     info._check_consistency()
     return info
 

diff --git a/mne/preprocessing/__init__.py b/mne/preprocessing/__init__.py
@@ -21,5 +21,6 @@
 from .stim import fix_stim_artifact
 from .maxwell import maxwell_filter
 from .xdawn import Xdawn
+from ._csd import compute_current_source_density
 from ._optical_density import optical_density
 from ._beer_lambert_law import beer_lambert_law
diff --git a/mne/preprocessing/_csd.py b/mne/preprocessing/_csd.py
@@ -0,0 +1,176 @@
+# Copyright 2003-2010 Jürgen Kayser <[email protected]>
+# Copyright 2017 Federico Raimondo <[email protected]> and
+#                Denis A. Engemann <[email protected]>
+#
+#
+# The original CSD Toolbox can be find at
+# http://psychophysiology.cpmc.columbia.edu/Software/CSDtoolbox/
+
+# Authors: Denis A. Engeman <[email protected]>
+#          Alex Rockhill <[email protected]>
+#
+# License: Relicensed under BSD (3-clause) and adapted with
+#          permission from authors of original GPL code
+
+import numpy as np
+
+from scipy.linalg import inv
+
+from mne import pick_types, pick_info
+from mne.utils import _validate_type, _check_ch_locs, _ensure_int
+from mne.io import BaseRaw
+from mne.io.constants import FIFF
+from mne.epochs import BaseEpochs
+from mne.evoked import Evoked
+from mne.bem import fit_sphere_to_headshape
+from mne.channels.interpolation import _calc_g, _calc_h
+
+
+def _prepare_G(G, lambda2):
+    G.flat[::len(G) + 1] += lambda2
+    # compute the CSD
+    Gi = inv(G)
+
+    TC = Gi.sum(0)
+    sgi = np.sum(TC)  # compute sum total
+
+    return Gi, TC, sgi
+
+
+def _compute_csd(G_precomputed, H, radius):
+    """Compute the CSD."""
+    n_channels = H.shape[0]
+    data = np.eye(n_channels)
+    mu = data.mean(0)
+    Z = data - mu
+
+    Gi, TC, sgi = G_precomputed
+
+    Cp2 = np.dot(Gi, Z)
+    c02 = np.sum(Cp2, axis=0) / sgi
+    C2 = Cp2 - np.dot(TC[:, None], c02[None, :])
+    X = np.dot(C2.T, H).T / radius ** 2
+    return X
+
+
+def compute_current_source_density(inst, lambda2=1e-5, stiffness=4,
+                                   n_legendre_terms=50, sphere='auto',
+                                   copy=True):
+    """Get the current source density (CSD) transformation.
+
+    Transformation based on spherical spline surface Laplacian.
+
+    .. note:: This function applies an average reference to the data.
+              Do not transform CSD data to source space.
+
+    Parameters
+    ----------
+    inst : instance of Raw, Epochs or Evoked
+        The data to be transformed.
+    lambda2 : float
+        Regularization parameter, produces smoothnes. Defaults to 1e-5.
+    stiffness : float
+        Stiffness of the spline. Also referred to as `m`.
+    n_legendre_terms : int
+        Number of Legendre terms to evaluate.
+    sphere : array-like, shape (4,)
+        The sphere, head-model of the form (x, y, z, r) where x, y, z
+        is the center of the sphere and r is the radius in meters.
+    copy : bool
+        Whether to overwrite instance data or create a copy.
+
+    Returns
+    -------
+    inst_csd : instance of Epochs or Evoked
+        The transformed data. Output type will match input type.
+    """
+    _validate_type(inst, (BaseEpochs, BaseRaw, Evoked), 'inst')
+
+    if inst.info['custom_ref_applied'] == -1:
+        raise ValueError('CSD already applied, should not be reapplied')
+
+    inst = inst.copy() if copy else inst
+
+    picks = pick_types(inst.info, meg=False, eeg=True, exclude='bads')
+
+    if len(picks) == 0:
+        raise ValueError('No EEG channels found.')
+
+    if lambda2 is None:
+        lambda2 = 1e-5
+
+    _validate_type(lambda2, 'numeric', 'lambda2')
+    if 0 > lambda2 or lambda2 > 1:
+        raise ValueError('lambda2 must be between 0 and 1, got %s' % lambda2)
+
+    _validate_type(stiffness, 'numeric', 'stiffness')
+    if stiffness < 0:
+        raise ValueError('stiffness must be non-negative got %s' % stiffness)
+
+    n_legendre_terms = _ensure_int(n_legendre_terms, 'n_legendre_terms')
+    if n_legendre_terms < 1:
+        raise ValueError('n_legendre_terms must be greater than 0, '
+                         'got %s' % n_legendre_terms)
+
+    if sphere == 'auto':
+        radius, origin_head, origin_device = fit_sphere_to_headshape(inst.info)
+        x, y, z = origin_head - origin_device
+    else:
+        _validate_type(sphere, tuple, 'sphere')
+        x, y, z, radius = sphere
+
+    _validate_type(x, 'numeric', 'x')
+    _validate_type(y, 'numeric', 'y')
+    _validate_type(z, 'numeric', 'z')
+    _validate_type(radius, 'numeric', 'radius')
+    if radius <= 0:
+        raise ValueError('sphere radius must be greater than 0, '
+                         'got %s' % radius)
+
+    _validate_type(copy, (bool), 'copy')
+
+    if not _check_ch_locs(inst.info['chs']):
+        raise ValueError('Zero or infinite position found in chs')
+
+    pos = np.array([inst.info['chs'][pick]['loc'][:3] for pick in picks])
+    pos -= (x, y, z)
+
+    G = _calc_g(np.dot(pos, pos.T), stiffness=stiffness,
+                num_lterms=n_legendre_terms)
+    H = _calc_h(np.dot(pos, pos.T), stiffness=stiffness,
+                num_lterms=n_legendre_terms)
+
+    G_precomputed = _prepare_G(G, lambda2)
+
+    trans_csd = _compute_csd(G_precomputed=G_precomputed,
+                             H=H, radius=radius)
+
+    if isinstance(inst, BaseEpochs):
+        for epo in inst._data:
+            epo[picks] = np.dot(trans_csd, epo[picks])
+    else:
+        inst._data = np.dot(trans_csd, inst._data[picks])
+
+    pick_info(inst.info, picks, copy=False)
+    inst.info['custom_ref_applied'] = -1
+    for ch in inst.info['chs']:
+        if ch['coil_type'] == FIFF.FIFFV_COIL_EEG:
+            ch.update(coil_type=FIFF.FIFFV_COIL_EEG_CSD,
+                      unit=FIFF.FIFF_UNIT_V_M2)
+    return inst
+
+# References
+# ----------
+#
+# [1] Perrin F, Bertrand O, Pernier J. "Scalp current density mapping:
+#     Value and estimation from potential data." IEEE Trans Biomed Eng.
+#     1987;34(4):283–288.
+#
+# [1] Perrin F, Pernier J, Bertrand O, Echallier JF. "Spherical splines
+#     for scalp potential and current density mapping."
+#     [Corrigenda EEG 02274, EEG Clin. Neurophysiol., 1990, 76, 565]
+#     Electroenceph Clin Neurophysiol. 1989;72(2):184–187.
+#
+# [2] Kayser J, Tenke CE. "On the benefits of using surface Laplacian
+#     (Current Source Density) methodology in electrophysiology."
+#     Int J Psychophysiol. 2015 Sep; 97(3): 171–173.