[ENH, MRG] Add interpolate bridged electrodes function

doc/changes/latest.inc

@@ -51,6 +51,8 @@ Enhancements
 
 - Add ``'voronoi'`` as an option for the ``image_interp`` argument in :func:`mne.viz.plot_topomap` to plot a topomap without interpolation using a Voronoi parcelation (:gh:`10571` by `Alex Rockhill`_)
 
+- Add :func:`mne.preprocessing.interpolate_bridged_electrodes` to use the spatially smeared signal to get a better interpolation rather than dropping those channels (:gh:`10587` by `Alex Rockhill`_)
+
 Bugs
 ~~~~
 - Make ``color`` parameter check in in :func:`mne.viz.plot_evoked_topo` consistent (:gh:`10217` by :newcontrib:`T. Wang` and `Stefan Appelhoff`_)

doc/preprocessing.rst

@@ -91,6 +91,7 @@ Projections:
    ica_find_ecg_events
    ica_find_eog_events
    infomax
+   interpolate_bridged_electrodes
    equalize_bads
    maxwell_filter
    maxwell_filter_prepare_emptyroom

examples/preprocessing/eeg_bridging.py

@@ -16,13 +16,14 @@
 bridging :footcite:`TenkeKayser2001`, which has been implemented in EEGLAB
 :footcite:`DelormeMakeig2004`. Unfortunately, there is not a lot to be
 done about electrode brigding once the data has been collected as far as
-preprocessing. Therefore, the recommendation is to check for electrode
-bridging early in data collection and address the problem. Or, if the data
-has already been collected, quantify the extent of the bridging so as not
-to introduce bias into the data from this effect and exclude subjects with
-bridging that might effect the outcome of a study. Preventing electrode
-bridging is ideal but awareness of the problem at least will mitigate its
-potential as a confound to a study. This tutorial follows
+preprocessing other than interpolating bridged channels. Therefore, our
+recommendation is to check for electrode bridging early in data collection
+and address the problem. Or, if the data has already been collected, quantify
+the extent of the bridging so as not to introduce bias into the data from this
+effect and exclude subjects with bridging that might effect the outcome of a
+study. Preventing electrode bridging is ideal but awareness of the problem at
+least will mitigate its potential as a confound to a study. This tutorial
+follows
 https://psychophysiology.cpmc.columbia.edu/software/eBridge/tutorial.html.
 """
 # Authors: Alex Rockhill <aprockhill@mailbox.org>
@@ -223,6 +224,23 @@
         raw_data[sub].info, bridged_idx, ed_matrix,
         title=f'Subject {sub} Bridged Electrodes', topomap_args=dict(vmax=5))
 
+# %%
+# For subjects with many bridged channels like Subject 6 shown in the example
+# above, it is advisable to exclude the subject. This because EEG recording
+# montage will not be comparable with the other subjects. And, if we tried to
+# interpole, the interpolation would depend on other channels which are also
+# bridged in that case. However, for subjects with only a few bridged channels,
+# those channels can be interpolated. Since the bridged data is still
+# biological (i.e. it is recording the subject's brain), it's just spatially
+# smeared, we can use :func:`mne.preprocessing.interpolate_bridged_electrodes`
+# to make a virtual channel midway between the two bridged channels
+# to aid in interpolation.
+
+# use subject 2, only one bridged electrode pair
+bridged_idx = ed_data[2][0]
+raw = mne.preprocessing.interpolate_bridged_electrodes(
+    raw_data[2], bridged_idx=bridged_idx)
+
 # %%
 # The Relationship Between Bridging and Impedances
 # ------------------------------------------------

mne/preprocessing/__init__.py

@@ -30,6 +30,6 @@
 from ._fine_cal import (compute_fine_calibration, read_fine_calibration,
                         write_fine_calibration)
 from .annotate_nan import annotate_nan
-from .interpolate import equalize_bads
+from .interpolate import equalize_bads, interpolate_bridged_electrodes
 from . import ieeg
 from ._css import cortical_signal_suppression

mne/preprocessing/interpolate.py

@@ -5,7 +5,12 @@
 from itertools import chain
 import numpy as np
 
-from ..epochs import BaseEpochs
+from ..utils import _validate_type
+from ..io import BaseRaw, RawArray
+from ..io.meas_info import create_info
+from ..epochs import BaseEpochs, EpochsArray
+from ..evoked import Evoked, EvokedArray
+from ..transforms import _sph_to_cart, _cart_to_sph
 
 
 def equalize_bads(insts, interp_thresh=1., copy=True):
@@ -65,3 +70,80 @@ def equalize_bads(insts, interp_thresh=1., copy=True):
         inst.info['bads'] = bads_keep
 
     return insts
+
+
+def interpolate_bridged_electrodes(inst, bridged_idx):
+    """Interpolate bridged electrode pairs.
+
+    Because bridged electrodes contain brain signal, it's just that the
+    signal is spatially smeared between the two electrodes, we can
+    make a virtual channel midway between the bridged pairs and use
+    that to aid in interpolation rather than completely discarding the
+    data from the two channels.
+
+    Parameters
+    ----------
+    inst : instance of Epochs, Evoked, or Raw
+        The data object with channels that are to be interpolated.
+    bridged_idx : list of tuple
+        The indices of channels marked as bridged with each bridged
+        pair stored as a tuple.
+
+    Returns
+    -------
+    inst : instance of Epochs, Evoked, or Raw
+        The modified data object.
+    """
+    _validate_type(inst, (BaseRaw, BaseEpochs, Evoked))
+    montage = inst.get_montage()
+    if montage is None:
+        raise RuntimeError('No channel positions found in ``inst``')
+    pos = montage.get_positions()
+    if pos['coord_frame'] != 'head':
+        raise RuntimeError('Montage channel positions must be in ``head``'
+                           'got {}'.format(pos['coord_frame']))
+    ch_pos = pos['ch_pos']
+    # store bads orig to put back at the end
+    bads_orig = inst.info['bads']
+    inst.info['bads'] = list()
+
+    # make virtual channels
+    virtual_chs = dict()
+    bads = set()
+    data = inst.get_data()
+    for idx0, idx1 in bridged_idx:
+        ch0 = inst.ch_names[idx0]
+        ch1 = inst.ch_names[idx1]
+        bads = bads.union([ch0, ch1])
+        # compute midway position in spherical coordinates in "head"
+        # (more accurate than cutting though the scalp by using cartesian)
+        pos0 = _cart_to_sph(ch_pos[ch0])
+        pos1 = _cart_to_sph(ch_pos[ch1])
+        pos_virtual = _sph_to_cart((pos0 + pos1) / 2)
+        virtual_info = create_info(
+            [f'{ch0}-{ch1} virtual'], inst.info['sfreq'], 'eeg')
+        virtual_info['chs'][0]['loc'][:3] = pos_virtual
+        if isinstance(inst, BaseRaw):
+            virtual_ch = RawArray(data[idx0:idx0 + 1], virtual_info,
+                                  first_samp=inst.first_samp)
+        elif isinstance(inst, BaseEpochs):
+            virtual_ch = EpochsArray(data[:, idx0:idx0 + 1], virtual_info,
+                                     tmin=inst.tmin)
+        else:  # evoked
+            virtual_ch = EvokedArray(data[idx0:idx0 + 1], virtual_info,
+                                     tmin=inst.tmin, nave=inst.nave,
+                                     kind=inst.kind)
+        virtual_chs[f'{ch0}-{ch1} virtual'] = virtual_ch
+
+    # add the virtual channels
+    inst.add_channels(list(virtual_chs.values()), force_update_info=True)
+
+    # use the virtual channels to interpolate
+    inst.info['bads'] = list(bads)
+    inst.interpolate_bads()
+
+    # drop virtual channels
+    inst.drop_channels(list(virtual_chs.keys()))
+
+    inst.info['bads'] = bads_orig
+    return inst

mne/preprocessing/tests/test_interpolate.py

@@ -1,8 +1,9 @@
 import os.path as op
+import numpy as np
 
 import pytest
 
-from mne.preprocessing.interpolate import equalize_bads
+from mne.preprocessing import equalize_bads, interpolate_bridged_electrodes
 from mne import io, pick_types, read_events, Epochs
 
 base_dir = op.join(op.dirname(__file__), '..', '..', 'io', 'tests', 'data')
@@ -58,3 +59,24 @@ def test_equalize_bads(interp_thresh, inst_type):
 
     for inst in insts_ok:
         assert set(inst.info['bads']) == set(bads_ok)
+
+
+def test_interpolate_bridged_electrodes():
+    """Test interpolate_bridged_electrodes function."""
+    raw, epochs, evoked = _load_data()
+    for inst in (raw, epochs, evoked):
+        idx0 = inst.ch_names.index('EEG 001')
+        idx1 = inst.ch_names.index('EEG 002')
+        ch_names_orig = inst.ch_names.copy()
+        bads_orig = inst.info['bads'].copy()
+        inst2 = inst.copy()
+        inst2.info['bads'] = ['EEG 001', 'EEG 002']
+        inst2.interpolate_bads()
+        data_interp_reg = inst2.get_data(picks=['EEG 001', 'EEG 002'])
+        inst = interpolate_bridged_electrodes(inst, [(idx0, idx1)])
+        data_interp = inst.get_data(picks=['EEG 001', 'EEG 002'])
+        assert not any(['virtual' in ch for ch in inst.ch_names])
+        assert inst.ch_names == ch_names_orig
+        assert inst.info['bads'] == bads_orig
+        # check closer to regular interpolation than original data
+        assert 1e-6 < np.mean(np.abs(data_interp - data_interp_reg)) < 5.4e-5

mne/transforms.py

@@ -721,8 +721,8 @@ def _cart_to_sph(cart):
     sph_pts : ndarray, shape (n_points, 3)
         Array containing points in spherical coordinates (rad, azimuth, polar)
     """
-    assert cart.ndim == 2 and cart.shape[1] == 3
     cart = np.atleast_2d(cart)
+    assert cart.ndim == 2 and cart.shape[1] == 3
     out = np.empty((len(cart), 3))
     out[:, 0] = np.sqrt(np.sum(cart * cart, axis=1))
     norm = np.where(out[:, 0] > 0, out[:, 0], 1)  # protect against / 0
@@ -746,8 +746,8 @@ def _sph_to_cart(sph_pts):
         Array containing points in Cartesian coordinates (x, y, z)
 
     """
-    assert sph_pts.ndim == 2 and sph_pts.shape[1] == 3
     sph_pts = np.atleast_2d(sph_pts)
+    assert sph_pts.ndim == 2 and sph_pts.shape[1] == 3
     cart_pts = np.empty((len(sph_pts), 3))
     cart_pts[:, 2] = sph_pts[:, 0] * np.cos(sph_pts[:, 2])
     xy = sph_pts[:, 0] * np.sin(sph_pts[:, 2])