diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index f779eb0a004..a0cc680a72a 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -41,6 +41,10 @@ Enhancements - Add keyboard shortcuts to toggle :meth:`mne.preprocessing.ICA.plot_properties` topomap channel types ('t') and power spectral density log-scale ('l') (:gh:`10557` by `Alex Rockhill`_) +- Add :func:`mne.preprocessing.compute_bridged_electrodes` to detect EEG electrodes with shared spatial sources due to a conductive medium connecting two or more electrodes, add :ref:`ex-eeg-bridging` for an example and :func:`mne.viz.plot_bridged_electrodes` to help visualize (:gh:`10571` by `Alex Rockhill`_) + +- 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`_) + Bugs ~~~~ - Fix bug in :func:`mne.io.read_raw_brainvision` when BrainVision data are acquired with the Brain Products "V-Amp" amplifier and disabled lowpass filter is marked with value ``0`` (:gh:`10517` by :newcontrib:`Alessandro Tonin`) diff --git a/doc/preprocessing.rst b/doc/preprocessing.rst index aec5eb5da10..25b9067b31a 100644 --- a/doc/preprocessing.rst +++ b/doc/preprocessing.rst @@ -76,6 +76,7 @@ Projections: annotate_nan compute_average_dev_head_t compute_current_source_density + compute_bridged_electrodes compute_fine_calibration compute_maxwell_basis compute_proj_ecg diff --git a/doc/references.bib b/doc/references.bib index 4e052285eec..66190aff1e7 100644 --- a/doc/references.bib +++ b/doc/references.bib @@ -362,6 +362,23 @@ @article{DarvasEtAl2006 year = {2006} } +@article{DelormeMakeig2004, + title = {{EEGLAB}: an open source toolbox for analysis of single-trial {EEG} dynamics including independent component analysis}, + volume = {134}, + issn = {0165-0270}, + shorttitle = {{EEGLAB}}, + doi = {10.1016/j.jneumeth.2003.10.009}, + language = {eng}, + number = {1}, + journal = {Journal of Neuroscience Methods}, + author = {Delorme, Arnaud and Makeig, Scott}, + month = mar, + year = {2004}, + pmid = {15102499}, + keywords = {Computer Simulation, Electroencephalography, Evoked Potentials, Software}, + pages = {9--21} +} + @article{DestrieuxEtAl2010, author = {Destrieux, Christophe and Fischl, Bruce and Dale, Anders and Halgren, Eric}, doi = {10.1016/j.neuroimage.2010.06.010}, @@ -579,6 +596,22 @@ @article{GramfortEtAl2013b year = {2013} } +@article{GreischarEtAl2004, + title = {Effects of electrode density and electrolyte spreading in dense array electroencephalographic recording}, + volume = {115}, + issn = {13882457}, + url = {https://linkinghub.elsevier.com/retrieve/pii/S1388245703003833}, + doi = {10.1016/j.clinph.2003.10.028}, + language = {en}, + number = {3}, + urldate = {2022-04-25}, + journal = {Clinical Neurophysiology}, + author = {Greischar, Lawrence L. and Burghy, Cory A. and van Reekum, Carien M. and Jackson, Daren C. and Pizzagalli, Diego A. and Mueller, Corrina and Davidson, Richard J.}, + month = mar, + year = {2004}, + pages = {710--720} +} + @article{GreveEtAl2013, author = {Greve, Douglas N. and {Van der Haegen}, Lise and Cai, Qing and Stufflebeam, Steven and Sabuncu, Mert R. and Fischl, Bruce and Brysbaert, Marc}, doi = {10.1162/jocn_a_00405}, @@ -1737,6 +1770,22 @@ @article{TauluSimola2006 year = {2006} } +@article{TenkeKayser2001, + title = {A convenient method for detecting electrolyte bridges in multichannel electroencephalogram and event-related potential recordings}, + volume = {112}, + issn = {1388-2457}, + doi = {10.1016/s1388-2457(00)00553-8}, + language = {eng}, + number = {3}, + journal = {Clinical Neurophysiology: Official Journal of the International Federation of Clinical Neurophysiology}, + author = {Tenke, C. E. and Kayser, J.}, + month = mar, + year = {2001}, + pmid = {11222978}, + keywords = {Algorithms, Artifacts, Brain, Electrodes, Electroencephalography, Electrolytes, Evoked Potentials, Humans, Scalp}, + pages = {545--550} +} + @article{TescheEtAl1995, author = {Tesche, Claudia D. and Uusitalo, Mikko A. and Ilmoniemi, Risto J. and Huotilainen, Minna and Kajola, Matti J. and Salonen, Oili L. M.}, doi = {10.1016/0013-4694(95)00064-6}, diff --git a/doc/visualization.rst b/doc/visualization.rst index 9f3b82753d1..080b8727ebf 100644 --- a/doc/visualization.rst +++ b/doc/visualization.rst @@ -24,6 +24,7 @@ Visualization mne_analyze_colormap plot_bem plot_brain_colorbar + plot_bridged_electrodes plot_chpi_snr plot_cov plot_channel_labels_circle diff --git a/examples/preprocessing/eeg_bridging.py b/examples/preprocessing/eeg_bridging.py new file mode 100644 index 00000000000..1c7f4ee36fa --- /dev/null +++ b/examples/preprocessing/eeg_bridging.py @@ -0,0 +1,262 @@ +# -*- coding: utf-8 -*- +""" +.. _ex-eeg-bridging: + +=============================================== +Identify EEG Electrodes Bridged by too much Gel +=============================================== + +Research-grade EEG often uses a gel based system, and when too much gel is +applied the gel conducting signal from the scalp to the electrode for one +electrode connects with the gel conducting signal from another electrode +"bridging" the two signals. This is undesirable because the signals from the +two (or more) electrodes are not as independent as they would otherwise be; +they are very similar to each other introducting additional +spatial smearing. An algorithm has been developed to detect electrode +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 +https://psychophysiology.cpmc.columbia.edu/software/eBridge/tutorial.html. +""" +# Authors: Alex Rockhill +# +# License: BSD-3-Clause + +# %% + +# sphinx_gallery_thumbnail_number = 2 + +import numpy as np +import matplotlib.pyplot as plt + +import mne + +print(__doc__) + +# %% +# Compute Electrical Distance Metric +# ---------------------------------- +# First, let's compute electrical distance metrics for a group of example +# subjects from the EEGBCI dataset in order to estimate electrode bridging. +# The electrical distance is just the variance of signals subtracted +# pairwise. Channels with activity that mirror another channel nearly +# exactly will have very low electrical distance. By inspecting the +# distribution of electrical distances, we can look for pairwise distances +# that are consistently near zero which are indicative of bridging. +# +# .. note:: It is likely to be sufficient to run this algorithm on a +# small portion (~3 minutes is probably plenty) of the data but +# that gel might settle over the course of a study causing more +# bridging so using the last segment of the data will +# give the most conservative estimate. + +montage = mne.channels.make_standard_montage('standard_1005') +ed_data = dict() # electrical distance/bridging data +raw_data = dict() # store infos for electrode positions +for sub in range(1, 11): + print(f'Computing electrode bridges for subject {sub}') + raw_fname = mne.datasets.eegbci.load_data(subject=sub, runs=(1,))[0] + raw = mne.io.read_raw(raw_fname, preload=True, verbose=False) + mne.datasets.eegbci.standardize(raw) # set channel names + raw.set_montage(montage, verbose=False) + raw_data[sub] = raw + ed_data[sub] = mne.preprocessing.compute_bridged_electrodes(raw) + + +# %% +# Examine an Electrical Distance Matrix +# ------------------------------------- +# Before we look at the electrical distance distributions across subjects, +# let's look at the distance matrix for one subject and try and understand +# how the algorithm works. We'll use subject 6 as it is a good example of +# bridging. In the zoomed out color scale version on the right, we can see +# that there is a distribution of electrical distances that are specific to +# that subject's head physiology/geometry and brain activity during the +# recording. On the right, when we clip the color range to zoom in, we can +# see several electrical distance outliers that are near zero; +# these indicate bridging. + +bridged_idx, ed_matrix = ed_data[6] + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) +fig.suptitle('Subject 6 Electrical Distance Matrix') + +# take median across epochs, only use upper triangular, lower is NaNs +ed_plot = np.zeros(ed_matrix.shape[1:]) * np.nan +triu_idx = np.triu_indices(ed_plot.shape[0], 1) +for idx0, idx1 in np.array(triu_idx).T: + ed_plot[idx0, idx1] = np.nanmedian(ed_matrix[:, idx0, idx1]) + +# plot full distribution color range +im1 = ax1.imshow(ed_plot, aspect='auto') +cax1 = fig.colorbar(im1, ax=ax1) +cax1.set_label(r'Electrical Distance ($\mu$$V^2$)') + +# plot zoomed in colors +im2 = ax2.imshow(ed_plot, aspect='auto', vmax=5) +cax2 = fig.colorbar(im2, ax=ax2) +cax2.set_label(r'Electrical Distance ($\mu$$V^2$)') +for ax in (ax1, ax2): + ax.set_xlabel('Channel Index') + ax.set_ylabel('Channel Index') + +fig.tight_layout() + +# %% +# Examine the Distribution of Electrical Distances +# ------------------------------------------------ +# Now let's plot a histogram of the electrical distance matrix. Note that the +# electrical distance matrix from the previous plot is upper triangular but +# does not include the diagonal. This means that the pairwise electrical +# distances are not computed between the same channel (which makes sense as +# the differences between a channel and itself would just be zero). The initial +# peak near zero therefore represents pairs of different channels that +# are nearly identical which is indicative of bridging. EEG recordings +# without bridged electrodes do not have a peak near zero. + +fig, ax = plt.subplots(figsize=(5, 5)) +fig.suptitle('Subject 6 Electrical Distance Matrix Distribution') +ax.hist(ed_matrix[~np.isnan(ed_matrix)], bins=np.linspace(0, 500, 51)) +ax.set_xlabel(r'Electrical Distance ($\mu$$V^2$)') +ax.set_ylabel('Count (channel pairs for all epochs)') + +# %% +# Plot Electrical Distances on a Topomap +# -------------------------------------- +# Now, let's look at the topography of the electrical distance matrix and +# see where our bridged channels are and check that their spatial +# arrangement makes sense. Here, we are looking at the minimum electrical +# distance for each channel and taking the median across all epochs +# (the raw data is epoched into 2 second non-overlapping intervals). +# This example is of the subject from the EEGBCI dataset with the most +# bridged channels so there are many light areas and red lines. They are +# generally grouped together and are biased toward horizontal connections +# (this may be because the EEG experimenter usually stands to the side and +# may have inserted the gel syringe tip in too far). + +mne.viz.plot_bridged_electrodes( + raw_data[6].info, bridged_idx, ed_matrix, + title='Subject 6 Bridged Electrodes', topomap_args=dict(vmax=5)) + +# %% +# Plot the Raw Voltage Time Series for Bridged Electrodes +# ------------------------------------------------------- +# Finally, let's do a sanity check and make sure that the bridged electrodes +# are indeed implausibly similar. We'll plot two bridged electrode pairs: +# F2-F4 and FC2-FC4, for subject 6 where they are bridged and subject 1 +# where they are not. As we can see, the pairs are nearly identical for +# subject 6 confirming that they are likely bridged. Interestingly, even +# though the two pairs are adjacent to each other, there are two distinctive +# pairs, meaning that it is unlikely that all four of these electrodes are +# bridged. + +raw = raw_data[6].copy().pick_channels(['F2', 'F4', 'FC2', 'FC4']) +raw.add_channels([mne.io.RawArray( + raw.get_data(ch1) - raw.get_data(ch2), + mne.create_info([f'{ch1}-{ch2}'], raw.info['sfreq'], 'eeg'), + raw.first_samp) for ch1, ch2 in [('F2', 'F4'), ('FC2', 'FC4')]]) +raw.plot(duration=20, scalings=dict(eeg=2e-4)) + +raw = raw_data[1].copy().pick_channels(['F2', 'F4', 'FC2', 'FC4']) +raw.add_channels([mne.io.RawArray( + raw.get_data(ch1) - raw.get_data(ch2), + mne.create_info([f'{ch1}-{ch2}'], raw.info['sfreq'], 'eeg'), + raw.first_samp) for ch1, ch2 in [('F2', 'F4'), ('FC2', 'FC4')]]) +raw.plot(duration=20, scalings=dict(eeg=2e-4)) + +# %% +# Compare Bridging Across Subjects in the EEGBCI Dataset +# ------------------------------------------------------- +# Now, let's look at the histograms of electrical distances for the whole +# EEGBCI dataset. As we can see in the zoomed in insert on the right, +# for subjects 6, 7 and 8 (and to a lesser extent 2 and 4), there is a +# different shape of the distribution of electrical distances around +# 0 :math:`{\mu}V^2` than for the other subjects. These subjects' +# distributions have a peak around 0 :math:`{\mu}V^2` distance +# and a trough around 5 :math:`{\mu}V^2` which is indicative of +# electrode bridging. The rest of the subjects' distributions increase +# monotonically, indicating normal spatial separation of sources. The +# large discrepancy in shapes of distributions is likely driven primarily by +# artifacts such as blinks which are an order of magnitude larger than +# neural data since this data has not been preprocessed but likely +# reflect neural or at least anatomical differences as well (i.e. the +# distance from the sensors to the brain). + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) +fig.suptitle('Electrical Distance Distribution for EEGBCI Subjects') +for ax in (ax1, ax2): + ax.set_ylabel('Count') + ax.set_xlabel(r'Electrical Distance ($\mu$$V^2$)') + +for sub, (bridged_idx, ed_matrix) in ed_data.items(): + # ed_matrix is upper triangular so exclude bottom half of NaNs + hist, edges = np.histogram(ed_matrix[~np.isnan(ed_matrix)].flatten(), + bins=np.linspace(0, 1000, 101)) + centers = (edges[1:] + edges[:-1]) / 2 + ax1.plot(centers, hist) + hist, edges = np.histogram(ed_matrix[~np.isnan(ed_matrix)].flatten(), + bins=np.linspace(0, 30, 21)) + centers = (edges[1:] + edges[:-1]) / 2 + ax2.plot(centers, hist, label=f'Sub {sub} #={len(bridged_idx)}') + +ax1.axvspan(0, 30, color='r', alpha=0.5) +ax2.legend(loc=(1.04, 0)) +fig.subplots_adjust(right=0.725, bottom=0.15, wspace=0.4) + +# %% +# For the group of subjects, let's look at their electrical distances +# and bridging. Especially since this is the same task, the lack of +# low electrical distances in many of the subjects is compelling +# evidence that the low electrical distance is caused by bridging +# and that it is avoidable given more judicious application of gel or +# other conductive electrolyte solution. + +for sub, (bridged_idx, ed_matrix) in ed_data.items(): + mne.viz.plot_bridged_electrodes( + raw_data[sub].info, bridged_idx, ed_matrix, + title=f'Subject {sub} Bridged Electrodes', topomap_args=dict(vmax=5)) + +# %% +# The Relationship Between Bridging and Impedances +# ------------------------------------------------ +# Electrode bridging is often brought about by inserting more gel in order +# to bring impendances down. Thus it can be helpful to compare bridging +# to impedances in the quest to be an ideal EEG technician! Low +# impedances lead to less noisy data and EEG without bridging is more +# spatially precise. Brain Imaging Data Structure (BIDS) recommends that +# impedances be stored in an EEG dataset in the `electrodes.tsv +# `_ file. +# Since the impedances are not stored for this dataset, we will fake +# them to demonstrate how they would be plotted. + +rng = np.random.default_rng(11) # seed for reproducibility +raw = raw_data[1] +impedances = rng.random((len(raw.ch_names,))) * 10 +fig, ax = plt.subplots(figsize=(5, 5)) +im, cn = mne.viz.plot_topomap(impedances, raw.info, axes=ax) +ax.set_title('Electrode Impendances') +cax = fig.colorbar(im, ax=ax) +cax.set_label(r'Impedance (k$\Omega$)') + +# %% +# Summary +# ------- +# In this example, we have shown a dataset where electrical bridging occurred +# during the EEG setup for several subjects. Hopefully this is convincing as to +# the importance of proper technique as well as checking your work to +# learn and improve as an EEG experimenter, and hopefully this tool will help +# us all collect better EEG data in the future. + +# %% +# References +# ---------- +# .. footbibliography:: diff --git a/examples/preprocessing/muscle_ica.py b/examples/preprocessing/muscle_ica.py index 9482e462b50..325c39c72b3 100644 --- a/examples/preprocessing/muscle_ica.py +++ b/examples/preprocessing/muscle_ica.py @@ -114,3 +114,8 @@ print(f'Manually found muscle artifact ICA components: {muscle_idx}\n' 'Automatically found muscle artifact ICA components: ' f'{muscle_idx_auto}') + +# %% +# References +# ---------- +# .. footbibliography:: diff --git a/mne/preprocessing/__init__.py b/mne/preprocessing/__init__.py index d12b7828354..efa3f15c7ef 100644 --- a/mne/preprocessing/__init__.py +++ b/mne/preprocessing/__init__.py @@ -23,12 +23,12 @@ compute_maxwell_basis, maxwell_filter_prepare_emptyroom) from .realign import realign_raw from .xdawn import Xdawn -from ._csd import compute_current_source_density +from ._csd import compute_current_source_density, compute_bridged_electrodes from . import nirs from .artifact_detection import (annotate_movement, compute_average_dev_head_t, annotate_muscle_zscore, annotate_break) from ._regress import regress_artifact -from ._fine_cal import (compute_fine_calibration, read_fine_calibration, +from ._fine_cal import (compute_fine_calibration, read_fine_calibration, write_fine_calibration) from .annotate_nan import annotate_nan from .interpolate import equalize_bads diff --git a/mne/preprocessing/_csd.py b/mne/preprocessing/_csd.py index a09db955206..988a3ffd718 100644 --- a/mne/preprocessing/_csd.py +++ b/mne/preprocessing/_csd.py @@ -12,13 +12,15 @@ # License: Relicensed under BSD-3-Clause and adapted with # permission from authors of original GPL code +import warnings import numpy as np from .. import pick_types -from ..utils import _validate_type, _ensure_int, _check_preload, verbose +from ..utils import (_validate_type, _ensure_int, _check_preload, verbose, + logger) from ..io import BaseRaw from ..io.constants import FIFF -from ..epochs import BaseEpochs +from ..epochs import BaseEpochs, make_fixed_length_epochs from ..evoked import Evoked from ..bem import fit_sphere_to_headshape from ..channels.interpolation import _calc_g, _calc_h @@ -175,3 +177,122 @@ def compute_current_source_density(inst, sphere='auto', lambda2=1e-5, inst.info['chs'][pick].update(coil_type=FIFF.FIFFV_COIL_EEG_CSD, unit=FIFF.FIFF_UNIT_V_M2) return inst + + +@verbose +def compute_bridged_electrodes(inst, lm_cutoff=16, epoch_threshold=0.5, + l_freq=0.5, h_freq=30, epoch_duration=2, + bw_method=None, verbose=None): + r"""Compute bridged EEG electrodes using the intrinsic Hjorth algorithm. + + First an electrical distance matrix is computed by taking the pairwise + variance. Then, a local maximum near 0 :math:`{\mu}V^2` and a + a local minimum below 5 :math:`{\mu}V^2` are found, the presence + of which is indicative of bridging. Finally, electrode distances below + the local minimum are marked as bridged as long as they happen on more + than the ``epoch_threshold`` proportion of epochs. + + Based on :footcite:`TenkeKayser2001,GreischarEtAl2004,DelormeMakeig2004` + and the `EEGLAB implementation + `_. + + Parameters + ---------- + inst : instance of Raw, Epochs or Evoked + The data to compute electrode bridging on. + lm_cutoff : float + The distance in :math:`{\mu}V^2` cutoff below which to + search for a local minimum (lm) indicative of bridging. + Defaults to 16 :math:`{\mu}V^2` to be conservative based + on the distributions in :footcite:`GreischarEtAl2004`. + epoch_threshold : float + The proportion of epochs with electrical distance less than + ``ed_threshold`` in order to consider the channel bridged. + The default is 0.5. + l_freq : float + The low cutoff frequency to use. Default is 0.5 Hz. + h_freq : float + The high cutoff frequency to use. Default is 30 Hz. + epoch_duration : float + The time in seconds to divide the raw into fixed-length epochs + to check for consistent bridging. Only used if ``inst`` is + :class:`mne.io.BaseRaw`. The default is 2 seconds. + bw_method : None + ``bw_method`` to pass to :class:`scipy.stats.gaussian_kde`. + %(verbose)s + + Returns + ------- + bridged_idx : list of tuple + The indices of channels marked as bridged with each bridged + pair stored as a tuple. + ed_matrix : ndarray of float, shape (n_epochs, n_channels, n_channels) + The electrical distance matrix for each pair of EEG electrodes. + + Notes + ----- + .. versionadded:: 1.1 + + References + ---------- + .. footbibliography:: + """ + from scipy.stats import gaussian_kde + from scipy.optimize import minimize_scalar + _check_preload(inst, 'Computing bridged electrodes') + inst = inst.copy() # don't modify original + picks = pick_types(inst.info, eeg=True) + if len(picks) == 0: + raise RuntimeError('No EEG channels found, cannot compute ' + 'electrode bridging') + # first, filter + inst.filter(l_freq=l_freq, h_freq=h_freq, picks=picks, verbose=False) + + if isinstance(inst, BaseRaw): + inst = make_fixed_length_epochs(inst, duration=epoch_duration, + preload=True, verbose=False) + + # standardize shape + data = inst.get_data(picks=picks) + if isinstance(inst, Evoked): + data = data[np.newaxis, ...] # expand evoked + + # next, compute electrical distance matrix, upper triangular + n_epochs = data.shape[0] + ed_matrix = np.zeros((n_epochs, picks.size, picks.size)) * np.nan + for i in range(picks.size): + for j in range(i + 1, picks.size): + ed_matrix[:, i, j] = np.var(data[:, i] - data[:, j], axis=1) + + # scale, fill in other half, diagonal + ed_matrix *= 1e12 # scale to muV**2 + + # initialize bridged indices + bridged_idx = list() + + # if not enough values below local minimum cutoff, return no bridges + ed_flat = ed_matrix[~np.isnan(ed_matrix)] + if ed_flat[ed_flat < lm_cutoff].size / n_epochs < epoch_threshold: + return bridged_idx, ed_matrix + + # kernel density estimation + kde = gaussian_kde(ed_flat[ed_flat < lm_cutoff]) + with warnings.catch_warnings(): + warnings.filterwarnings( + 'ignore', 'invalid value encountered in true_divide') + local_minimum = float(minimize_scalar( + lambda x: kde(x) if x < lm_cutoff and x > 0 else np.inf).x) + logger.info(f'Local minimum {local_minimum} found') + + # find electrodes that are below the cutoff local minimum on + # `epochs_threshold` proportion of epochs + for i in range(picks.size): + for j in range(i + 1, picks.size): + bridged_count = np.sum(ed_matrix[:, i, j] < local_minimum) + if bridged_count / n_epochs > epoch_threshold: + logger.info('Bridge detected between ' + f'{inst.ch_names[picks[i]]} and ' + f'{inst.ch_names[picks[j]]}') + bridged_idx.append((picks[i], picks[j])) + + return bridged_idx, ed_matrix diff --git a/mne/preprocessing/tests/test_csd.py b/mne/preprocessing/tests/test_csd.py index 9586a5ff479..4c339e132ca 100644 --- a/mne/preprocessing/tests/test_csd.py +++ b/mne/preprocessing/tests/test_csd.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -"""Test the compute_current_source_density function. +"""Test the current source density and related functions. For each supported file format, implement a test. """ @@ -23,7 +23,8 @@ from mne.utils import object_diff from mne.datasets import testing -from mne.preprocessing import compute_current_source_density +from mne.preprocessing import (compute_current_source_density, + compute_bridged_electrodes) data_path = op.join(testing.data_path(download=False), 'preprocessing') eeg_fname = op.join(data_path, 'test_eeg.mat') @@ -183,3 +184,29 @@ def test_csd_fif(): ch.update(coil_type=FIFF.FIFFV_COIL_EEG, unit=FIFF.FIFF_UNIT_V) raw_csd._data[pick] = raw._data[pick] assert object_diff(raw.info, raw_csd.info) == '' + + +def test_compute_bridged_electrodes(): + """Test computing bridged electrodes.""" + # test I/O + raw = read_raw_fif(raw_fname).load_data() + raw.pick_types(meg=True) + with pytest.raises(RuntimeError, match='No EEG channels found'): + bridged_idx, ed_matrix = compute_bridged_electrodes(raw) + + # test output + epoch_duration = 3 + raw = read_raw_fif(raw_fname).load_data() + idx0 = raw.ch_names.index('EEG 001') + idx1 = raw.ch_names.index('EEG 002') + raw._data[idx1] = raw._data[idx0] + bridged_idx, ed_matrix = compute_bridged_electrodes( + raw, epoch_duration=epoch_duration) + assert bridged_idx == [(idx0, idx1)] + picks = pick_types(raw.info, meg=False, eeg=True) + assert ed_matrix.shape == \ + (raw.times.size // (epoch_duration * raw.info['sfreq']), + picks.size, picks.size) + picks = list(picks) + assert np.all(ed_matrix[:, picks.index(idx0), picks.index(idx1)] == 0) + assert np.all(np.isnan(ed_matrix[0][np.tril_indices(len(picks), -1)])) diff --git a/mne/utils/docs.py b/mne/utils/docs.py index d97b543d9d2..4cf6557ebe5 100644 --- a/mne/utils/docs.py +++ b/mne/utils/docs.py @@ -1428,7 +1428,7 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): docdict['image_interp_topomap'] = """ image_interp : str The image interpolation to be used. All matplotlib options are - accepted. + accepted as well as ``'voronoi'``. """ docdict['include_tmax'] = """ diff --git a/mne/viz/__init__.py b/mne/viz/__init__.py index 1deffd19d70..2a7251ade2a 100644 --- a/mne/viz/__init__.py +++ b/mne/viz/__init__.py @@ -2,7 +2,8 @@ from .topomap import (plot_evoked_topomap, plot_projs_topomap, plot_arrowmap, plot_ica_components, plot_tfr_topomap, plot_topomap, - plot_epochs_psd_topomap, plot_layout) + plot_epochs_psd_topomap, plot_layout, + plot_bridged_electrodes) from .topo import plot_topo_image_epochs, iter_topography from .utils import (tight_layout, mne_analyze_colormap, compare_fiff, ClickableImage, add_background_image, plot_sensors, diff --git a/mne/viz/tests/test_topomap.py b/mne/viz/tests/test_topomap.py index b9b79137bbb..67392cc8733 100644 --- a/mne/viz/tests/test_topomap.py +++ b/mne/viz/tests/test_topomap.py @@ -30,7 +30,8 @@ from mne.viz import plot_evoked_topomap, plot_projs_topomap, topomap from mne.viz.topomap import (_get_pos_outlines, _onselect, plot_topomap, - plot_arrowmap, plot_psds_topomap) + plot_arrowmap, plot_psds_topomap, + plot_bridged_electrodes) from mne.viz.utils import _find_peaks, _fake_click from mne.utils import requires_sklearn, check_version @@ -647,8 +648,8 @@ def test_plot_topomap_cnorm(): else: from matplotlib.colors import DivergingNorm as TwoSlopeNorm - np.random.seed(42) - v = np.random.uniform(low=-1, high=2.5, size=64) + rng = np.random.default_rng(42) + v = rng.uniform(low=-1, high=2.5, size=64) v[:3] = [-1, 0, 2.5] montage = make_standard_montage("biosemi64") @@ -667,3 +668,25 @@ def test_plot_topomap_cnorm(): msg = "vmax=2.5 is implicitly defined by cnorm, ignoring vmax=10.*" with pytest.warns(RuntimeWarning, match=msg): plot_topomap(v, info, vmax=10, cnorm=cnorm) + + +def test_plot_bridged_electrodes(): + """Test plotting of bridged electrodes.""" + rng = np.random.default_rng(42) + montage = make_standard_montage("biosemi64") + info = create_info(montage.ch_names, 256, "eeg").set_montage("biosemi64") + bridged_idx = [(0, 1), (2, 3)] + n_epochs = 10 + ed_matrix = np.zeros((n_epochs, len(info.ch_names), + len(info.ch_names))) * np.nan + triu_idx = np.triu_indices(len(info.ch_names), 1) + for i in range(n_epochs): + ed_matrix[i][triu_idx] = rng.random() + rng.random(triu_idx[0].size) + fig = plot_bridged_electrodes(info, bridged_idx, ed_matrix, + topomap_args=dict(names=info.ch_names, + vmax=1, show_names=True)) + # two bridged lines plus head outlines + assert len(fig.axes[0].lines) == 6 + + with pytest.raises(RuntimeError, match='Expected'): + plot_bridged_electrodes(info, bridged_idx, np.zeros((5, 6, 7))) diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py index 42fd393fad4..ad345c45a82 100644 --- a/mne/viz/topomap.py +++ b/mne/viz/topomap.py @@ -737,9 +737,7 @@ def plot_topomap(data, pos, vmin=None, vmax=None, cmap=None, sensors=True, If an array, the values represent the levels for the contours. The values are in µV for EEG, fT for magnetometers and fT/m for gradiometers. Defaults to 6. - image_interp : str - The image interpolation to be used. All matplotlib options are - accepted. + %(image_interp_topomap)s show : bool Show figure if True. onselect : callable | None @@ -834,6 +832,41 @@ def _setup_interp(pos, res, extrapolate, sphere, outlines, border): return extent, Xi, Yi, interp +_VORONOI_CIRCLE_RES = 100 + + +def _voronoi_topomap(data, pos, outlines, ax, cmap, norm, extent, res): + """Make a Voronoi diagram on a topomap.""" + from scipy.spatial import Voronoi + # we need an image axis object so first empty image to plot over + im = ax.imshow(np.zeros((res, res)) * np.nan, cmap=cmap, + origin='lower', aspect='equal', extent=extent, + norm=norm) + rx, ry = outlines['clip_radius'] + cx, cy = outlines.get('clip_origin', (0., 0.)) + # add points on the circle to make boundaries, expand out to + # ensure regions extend to the edge of the topomap + vor = Voronoi(np.concatenate( + [pos, [(rx * 1.5 * np.cos(2 * np.pi / _VORONOI_CIRCLE_RES * t), + ry * 1.5 * np.sin(2 * np.pi / _VORONOI_CIRCLE_RES * t)) + for t in range(_VORONOI_CIRCLE_RES)]])) + for point_idx, region_idx in enumerate( + vor.point_region[:-_VORONOI_CIRCLE_RES]): + if -1 in vor.regions[region_idx]: + continue + polygon = list() + for i in vor.regions[region_idx]: + x, y = vor.vertices[i] + if (x - cx)**2 / rx**2 + (y - cy)**2 / ry**2 < 1: + polygon.append((x, y)) + else: + x *= rx / np.linalg.norm(vor.vertices[i]) + y *= ry / np.linalg.norm(vor.vertices[i]) + polygon.append((x, y)) + ax.fill(*zip(*polygon), color=cmap(norm(data[point_idx]))) + return im + + def _get_patch(outlines, extrapolate, interp, ax): from matplotlib import patches clip_radius = outlines['clip_radius'] @@ -931,7 +964,7 @@ def _plot_topomap(data, pos, vmin=None, vmax=None, cmap=None, sensors=True, norm = min(data) >= 0 vmin, vmax = _setup_vmin_vmax(data, vmin, vmax, norm) if cmap is None: - cmap = 'Reds' if norm else 'RdBu_r' + cmap = plt.get_cmap('Reds' if norm else 'RdBu_r') outlines = _make_head_outlines(sphere, pos, outlines, (0., 0.)) assert isinstance(outlines, dict) @@ -950,11 +983,17 @@ def _plot_topomap(data, pos, vmin=None, vmax=None, cmap=None, sensors=True, # plot outline patch_ = _get_patch(outlines, extrapolate, interp, ax) - # plot interpolated map + # get colormap normalization if cnorm is None: cnorm = Normalize(vmin=vmin, vmax=vmax) - im = ax.imshow(Zi, cmap=cmap, origin='lower', aspect='equal', - extent=extent, interpolation=image_interp, norm=cnorm) + + # plot interpolated map + if image_interp == 'voronoi': + im = _voronoi_topomap(data, pos=pos, outlines=outlines, ax=ax, + cmap=cmap, norm=cnorm, extent=extent, res=res) + else: + im = ax.imshow(Zi, cmap=cmap, origin='lower', aspect='equal', + extent=extent, interpolation=image_interp, norm=cnorm) # gh-1432 had a workaround for no contours here, but we'll remove it # because mpl has probably fixed it @@ -1130,9 +1169,7 @@ def plot_ica_components(ica, picks=None, ch_type=None, res=64, values for the contour thresholds (may sometimes be inaccurate, use array for accuracy). If an array, the values represent the levels for the contours. Defaults to 6. - image_interp : str - The image interpolation to be used. All matplotlib options are - accepted. + %(image_interp_topomap)s inst : Raw | Epochs | None To be able to see component properties after clicking on component topomap you need to pass relevant data - instances of Raw or Epochs @@ -2619,3 +2656,71 @@ def plot_arrowmap(data, info_from, info_to=None, scale=3e-10, vmin=None, plt_show(show) return fig + + +@fill_doc +def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, + topomap_args=None): + """Topoplot electrode distance matrix with bridged electrodes connected. + + Parameters + ---------- + %(info_not_none)s + bridged_idx : list of tuple + The indices of channels marked as bridged with each bridged + pair stored as a tuple. + ed_matrix : ndarray of float, shape (n_channels, n_channels) + The electrical distance matrix for each pair of EEG electrodes. + title : str + A title to add to the plot. + topomap_args : dict | None + Arguments to pass to :func:`mne.viz.plot_topomap`. + + Returns + ------- + fig : instance of matplotlib.figure.Figure + The topoplot figure handle. + """ + import matplotlib.pyplot as plt + if topomap_args is None: + topomap_args = dict() + else: + topomap_args = topomap_args.copy() # don't change original + topomap_args.setdefault('image_interp', 'voronoi') + topomap_args.setdefault('names', info.ch_names) + topomap_args.setdefault('show_names', True) + topomap_args.setdefault('contours', False) + if 'axes' not in topomap_args: + fig, ax = plt.subplots() + topomap_args['axes'] = ax + else: + fig = None + # handle colorbar here instead of in plot_topomap + colorbar = topomap_args.pop('colorbar') if \ + 'colorbar' in topomap_args else True + # use sphere to find positions + sphere = topomap_args['sphere'] if 'sphere' in topomap_args else None + picks = pick_types(info, eeg=True) + if ed_matrix.shape[1:] != (picks.size, picks.size): + raise RuntimeError( + f'Expected {(ed_matrix.shape[0], picks.size, picks.size)} ' + f'shaped `ed_matrix`, got {ed_matrix.shape}') + # fill in lower triangular + ed_matrix = ed_matrix.copy() + tril_idx = np.tril_indices(picks.size) + for epo_idx in range(ed_matrix.shape[0]): + ed_matrix[epo_idx][tril_idx] = ed_matrix[epo_idx].T[tril_idx] + elec_dists = np.median(np.nanmin(ed_matrix, axis=1), axis=0) + pos = _find_topomap_coords(info, picks, sphere=sphere) + im, cn = plot_topomap(elec_dists, info, **topomap_args) + fig = im.figure if fig is None else fig + # add bridged connections + for idx0, idx1 in bridged_idx: + im.axes.plot([pos[idx0][0], pos[idx1][0]], + [pos[idx0][1], pos[idx1][1]], color='r') + if title is not None: + im.axes.set_title(title) + if colorbar: + cax = fig.colorbar(im) + cax.set_label(r'Electrical Distance ($\mu$$V^2$)') + return fig