From 8cd5020f6e0226e77c6b9ada9d2afb3193b310ca Mon Sep 17 00:00:00 2001 From: Alex Date: Tue, 26 Apr 2022 09:10:17 -0700 Subject: [PATCH 01/12] wip, save commit with instrinsic hjorth calculation --- doc/references.bib | 49 ++++++++++++ examples/preprocessing/eeg_bridging.py | 82 ++++++++++++++++++++ mne/preprocessing/__init__.py | 4 +- mne/preprocessing/_csd.py | 101 ++++++++++++++++++++++++- mne/preprocessing/tests/test_csd.py | 6 +- 5 files changed, 238 insertions(+), 4 deletions(-) create mode 100644 examples/preprocessing/eeg_bridging.py 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/examples/preprocessing/eeg_bridging.py b/examples/preprocessing/eeg_bridging.py new file mode 100644 index 00000000000..15f34c2741f --- /dev/null +++ b/examples/preprocessing/eeg_bridging.py @@ -0,0 +1,82 @@ +# -*- coding: utf-8 -*- +""" +.. _ex-eeg-briding: + +=============================================== +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 as the signals from the two +(or more) electrodes are not as independent as they would otherwise be; +spatial smearing is caused wherein the signals are more similar to each other +been developed to detect electrode bridging :footcite:`TenkeKayser2001`, which +was previously 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. +""" +# Authors: Alex Rockhill +# +# License: BSD-3-Clause + +# %% + +# sphinx_gallery_thumbnail_number = 2 + +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 +meg_path = data_path / 'MEG' / 'sample' +raw = mne.io.read_raw_fif(meg_path / 'sample_audvis_raw.fif') +raw = raw.pick_types(eeg=True, exclude=raw.info['bads']).load_data() + +# %% +# Compute the bridged channels and plot them on a topomap. + +bridged_idx, scores = mne.preprocessing.compute_bridged_electrodes(raw) +mne.viz.plot_bridged_electrodes(bridged_idx, scores, raw.info, + title='Bridged Electrodes') + +# %% +# 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. Impedances can be stored with an EEG dataset +# like in the :ref:`electrodes-tsv` Brain Imaging Data Structure (BIDS) +# file. Since the impedances are not stored for this dataset, we will fake +# them to demonstrate how they would be plotted. + +np.random.seed(11) # seed for reproducibility +impedances = np.random.random((len(raw.ch_names,))) +fig, ax = plt.subplots(figsize=(5, 5)) +im, cn = mne.viz.plot_topomap(impedances, raw.info, axes=ax) +ax.set_title('Electrode Impendances Audio/Visual Task') +cax = fig.colorbar(im, ax=ax) +cax.set_label(r'Impedance (k$\Omega$)') + + +# %% +# References +# ---------- +# .. footbibliography:: +# +# +# .. _electrodes-tsv: https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/03-electroencephalography.html#electrodes-description-_electrodestsv # noqa E501 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..4163425718c 100644 --- a/mne/preprocessing/_csd.py +++ b/mne/preprocessing/_csd.py @@ -18,7 +18,7 @@ from ..utils import _validate_type, _ensure_int, _check_preload, verbose 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 +175,102 @@ 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, ed_threshold=3, epoch_threshold=0.5, + l_freq=0.5, h_freq=30, epoch_duration=2, + verbose=None): + """Compute bridged EEG electrodes using the intrinsic Hjorth algorithm. + + Based on :footcite:`TenkeKayser2001,GreischarEtAl2004,DelormeMakeig2004`. + Original implementation + https://psychophysiology.cpmc.columbia.edu/software/eBridge/index.html. + + Parameters + ---------- + inst : instance of Raw, Epochs or Evoked + The data to compute electrode bridging on. + ed_threshold : float + The bridged electrode threshold in μV:sup:`2`. The default is + 3 μV:sup:`2` as used 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. + %(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_channels, n_channels) + The electrical distance matrix for each pair of EEG electrodes. + + Notes + ----- + .. versionadded:: 1.1 + + References + ---------- + .. footbibliography:: + """ + _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) + + if isinstance(inst, BaseRaw): + inst = make_fixed_length_epochs(inst, duration=epoch_duration, + preload=True) + + # standardize shape + data = inst.get_data() + if isinstance(inst, Evoked): + data = data.reshape((1,) + data.shape) # expand evoked + + # next, compute electrical distance matrix + n_epochs = data.shape[0] + ed_matrix = np.zeros((n_epochs, picks.size, picks.size)) + for i, idx0 in enumerate(picks): + for j, idx1 in enumerate(picks[:i + 1]): + ed_matrix[:, i, j] = np.var(data[:, idx0] - data[:, idx1], axis=1) + + # scale, fill in other half, diagonal + ed_matrix *= 1e12 # scale to muV**2 + ed_matrix += np.swapaxes(ed_matrix, 2, 1) + for i in range(n_epochs): + np.fill_diagonal(ed_matrix[i], np.inf) + + ''' + ed_matrix = np.median(ed_matrix, axis=0) + np.fill_diagonal(ed_matrix, 1) + # now, compute weighting factor + weights = (1 / ed_matrix) / np.max(1 / ed_matrix, axis=0) + + # finally compute intrinsic hjorth by subtracting scale nearest + # electical distance neighbor + hjorth = np.empty(data.shape) + for i, idx0 in enumerate(picks): + idx1 = np.argmax(weights[idx0]) + hjorth[i] = data[idx0] - data[idx1] * weights[idx0, idx1] + hjorth *= 1e12 # scale to μV**2 + ''' + + bridged_idx = [idx for i, idx in enumerate(picks) if + np.sum(np.min(ed_matrix[:, i], axis=0) < ed_threshold) > + n_epochs * epoch_threshold] + return bridged_idx, ed_matrix diff --git a/mne/preprocessing/tests/test_csd.py b/mne/preprocessing/tests/test_csd.py index 9586a5ff479..d29694cadb3 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. """ @@ -183,3 +183,7 @@ 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.""" From b86d695097fdc02cdf9eb24579f7c43caee8d78a Mon Sep 17 00:00:00 2001 From: Alex Date: Tue, 26 Apr 2022 16:24:49 -0700 Subject: [PATCH 02/12] working version but a bit off --- doc/changes/latest.inc | 2 + doc/preprocessing.rst | 1 + doc/visualization.rst | 1 + examples/preprocessing/eeg_bridging.py | 65 +++++++++++++----- examples/preprocessing/muscle_ica.py | 5 ++ mne/preprocessing/_csd.py | 95 +++++++++++++++++--------- mne/viz/__init__.py | 3 +- mne/viz/topomap.py | 49 +++++++++++++ 8 files changed, 171 insertions(+), 50 deletions(-) diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index f779eb0a004..e4938de166b 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -41,6 +41,8 @@ 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:`XX` 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/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 index 15f34c2741f..2229903ac4c 100644 --- a/examples/preprocessing/eeg_bridging.py +++ b/examples/preprocessing/eeg_bridging.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- """ -.. _ex-eeg-briding: +.. _ex-eeg-bridging: =============================================== Identify EEG Electrodes Bridged by too much Gel @@ -22,6 +22,8 @@ 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 # @@ -35,37 +37,68 @@ import matplotlib.pyplot as plt import mne -from mne.datasets import sample print(__doc__) -data_path = sample.data_path() - # %% -# Load sample subject data -meg_path = data_path / 'MEG' / 'sample' -raw = mne.io.read_raw_fif(meg_path / 'sample_audvis_raw.fif') -raw = raw.pick_types(eeg=True, exclude=raw.info['bads']).load_data() +# Let's look at the histograms of electrical distances for the EEGBCI dataset. +# As we can see, for subjects 6, 7 and 8 (and to a lesser extent 4), there is +# a different shape of the distribution of electrical distances than for the +# other subjects. These subjects' distributions have a peak around 0 +# :math:`{\\mu}`V:sup:`2` distance and a trough around 5 +# :math:`{\\mu}`V:sup:`2` which is indicative of electrode bridging. +# The rest of the subjects' distributions increase monotonically, +# indicating normal spatial separation of sources. + +fig, ax = plt.subplots(figsize=(6, 4)) +ax.set_title('Electrical Distance Distribution for EEGBCI Subjects') +ax.set_ylabel('Count') +ax.set_xlabel(r'Electrical Distance ($\mu$$V^2$)') +for sub in range(1, 11): + print(f'Computing electrode bridges for subject {sub}') + raw = mne.io.read_raw(mne.datasets.eegbci.load_data( + subject=sub, runs=(1,))[0], preload=True, verbose=False) + mne.datasets.eegbci.standardize(raw) # set channel names + montage = mne.channels.make_standard_montage('standard_1005') + raw.set_montage(montage, verbose=False) + + bridged_idx, ed_matrix = mne.preprocessing.compute_bridged_electrodes(raw) + # 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, 100, 51)) + ax.plot((edges[1:] + edges[:-1]) / 2, hist, + label=f'Sub {sub} #={len(bridged_idx)}') + +ax.legend(loc=(1.04, 0)) +fig.subplots_adjust(right=0.725, bottom=0.15) # %% -# Compute the bridged channels and plot them on a topomap. +# Let's look at one subject with bridged electrodes and plot what's +# going on to try and understand electrode bridging better. + +raw = mne.io.read_raw(mne.datasets.eegbci.load_data( + subject=6, runs=(1,))[0], preload=True) +mne.datasets.eegbci.standardize(raw) # set channel names +montage = mne.channels.make_standard_montage('standard_1005') +raw.set_montage(montage) -bridged_idx, scores = mne.preprocessing.compute_bridged_electrodes(raw) -mne.viz.plot_bridged_electrodes(bridged_idx, scores, raw.info, - title='Bridged Electrodes') +bridged_idx, ed_matrix = mne.preprocessing.compute_bridged_electrodes(raw) +mne.viz.plot_bridged_electrodes( + raw.info, bridged_idx, ed_matrix, title='Bridged Electrodes', + topomap_args=dict(names=raw.ch_names, show_names=True)) # %% # 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. Impedances can be stored with an EEG dataset -# like in the :ref:`electrodes-tsv` Brain Imaging Data Structure (BIDS) -# file. Since the impedances are not stored for this dataset, we will fake +# spatially precise. Impedances are recommended to be stored in an EEG dataset +# in the :ref:`electrodes-tsv` file within th eBrain Imaging Data Structure +# (BIDS). Since the impedances are not stored for this dataset, we will fake # them to demonstrate how they would be plotted. np.random.seed(11) # seed for reproducibility -impedances = np.random.random((len(raw.ch_names,))) +impedances = np.random.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 Audio/Visual Task') 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/_csd.py b/mne/preprocessing/_csd.py index 4163425718c..1e865e6fdf3 100644 --- a/mne/preprocessing/_csd.py +++ b/mne/preprocessing/_csd.py @@ -15,7 +15,8 @@ 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, make_fixed_length_epochs @@ -178,11 +179,19 @@ def compute_current_source_density(inst, sphere='auto', lambda2=1e-5, @verbose -def compute_bridged_electrodes(inst, ed_threshold=3, epoch_threshold=0.5, - l_freq=0.5, h_freq=30, epoch_duration=2, - verbose=None): +def compute_bridged_electrodes(inst, ed_threshold=3, lm_cutoff=5, + n_bins=20, ed_max=12, epoch_threshold=0.5, + l_freq=0.5, h_freq=30, + epoch_duration=2, verbose=None): """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:sup:`2` and a + a local minimum below 5 :math:`{\\mu}`V:sup:`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`. Original implementation https://psychophysiology.cpmc.columbia.edu/software/eBridge/index.html. @@ -192,8 +201,19 @@ def compute_bridged_electrodes(inst, ed_threshold=3, epoch_threshold=0.5, inst : instance of Raw, Epochs or Evoked The data to compute electrode bridging on. ed_threshold : float - The bridged electrode threshold in μV:sup:`2`. The default is - 3 μV:sup:`2` as used in :footcite:`GreischarEtAl2004`. + The bridged electrode threshold in :math:`{\\mu}`V:sup:`2`. The + default is 3 :math:`{\\mu}`V:sup:`2` as used in + :footcite:`GreischarEtAl2004`. + lm_cutoff : float + The distance in :math:`{\\mu}`V:sup:`2` cutoff below which to + search for a local minimum (lm) indicative of bridging. + Defaults to 5 :math:`{\\mu}`V:sup:`2`. + n_bins : int + The number of bins to use to create a histogram of the electrical + distance matrix. + ed_max : float + The maximum electrical distance to consider when making the + histogram of electrical distances. epoch_threshold : float The proportion of epochs with electrical distance less than ``ed_threshold`` in order to consider the channel bridged. @@ -231,46 +251,55 @@ def compute_bridged_electrodes(inst, ed_threshold=3, epoch_threshold=0.5, raise RuntimeError('No EEG channels found, cannot compute ' 'electrode bridging') # first, filter - inst.filter(l_freq=l_freq, h_freq=h_freq, picks=picks) + 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) + preload=True, verbose=False) # standardize shape data = inst.get_data() if isinstance(inst, Evoked): data = data.reshape((1,) + data.shape) # expand evoked - # next, compute electrical distance matrix + # next, compute electrical distance matrix, upper triangular n_epochs = data.shape[0] - ed_matrix = np.zeros((n_epochs, picks.size, picks.size)) + ed_matrix = np.zeros((n_epochs, picks.size, picks.size)) * np.nan for i, idx0 in enumerate(picks): - for j, idx1 in enumerate(picks[:i + 1]): - ed_matrix[:, i, j] = np.var(data[:, idx0] - data[:, idx1], axis=1) + for j, idx1 in enumerate(picks[i + 1:]): + ed_matrix[:, i, i + 1 + j] = \ + np.var(data[:, idx0] - data[:, idx1], axis=1) # scale, fill in other half, diagonal ed_matrix *= 1e12 # scale to muV**2 - ed_matrix += np.swapaxes(ed_matrix, 2, 1) - for i in range(n_epochs): - np.fill_diagonal(ed_matrix[i], np.inf) - - ''' - ed_matrix = np.median(ed_matrix, axis=0) - np.fill_diagonal(ed_matrix, 1) - # now, compute weighting factor - weights = (1 / ed_matrix) / np.max(1 / ed_matrix, axis=0) - - # finally compute intrinsic hjorth by subtracting scale nearest - # electical distance neighbor - hjorth = np.empty(data.shape) + + # initialize bridged indices + bridged_idx = list() + + # compute histogram + bins, edges = np.histogram(ed_matrix[~np.isnan(ed_matrix)], + bins=np.linspace(0, ed_max, n_bins)) + centers = (edges[1:] + edges[:-1]) / 2 + + # check if has peak near 0, if not, no bridges + if np.diff(bins[centers <= ed_threshold]).sum() >= 0: + return bridged_idx, ed_matrix + + # find local minimum, also indicative of bridging + centers_interp = np.linspace(0, lm_cutoff, n_bins) + bins_interp = np.interp(centers_interp, centers, bins) + local_minimum = centers_interp[np.argmin(bins_interp)] + + # find electrodes that are below the cutoff local minumum on + # `epochs_threshold` proportion of epochs + check_bridged = ed_matrix < local_minimum for i, idx0 in enumerate(picks): - idx1 = np.argmax(weights[idx0]) - hjorth[i] = data[idx0] - data[idx1] * weights[idx0, idx1] - hjorth *= 1e12 # scale to μV**2 - ''' - - bridged_idx = [idx for i, idx in enumerate(picks) if - np.sum(np.min(ed_matrix[:, i], axis=0) < ed_threshold) > - n_epochs * epoch_threshold] + for j, idx1 in enumerate(picks[i + 1:]): + if np.sum(check_bridged[:, i, i + 1 + j]) / \ + n_epochs > epoch_threshold: + logger.info('Bridge detected between ' + f'{inst.ch_names[picks[idx0]]} and ' + f'{inst.ch_names[picks[idx1]]}') + bridged_idx.append((idx0, idx1)) + return bridged_idx, ed_matrix 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/topomap.py b/mne/viz/topomap.py index 42fd393fad4..4da70e3d502 100644 --- a/mne/viz/topomap.py +++ b/mne/viz/topomap.py @@ -2619,3 +2619,52 @@ def plot_arrowmap(data, info_from, info_to=None, scale=3e-10, vmin=None, plt_show(show) return fig + + +def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, + colorbar=True, sphere=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. + colorbar : bool + Whether to add a colorbar. + %(sphere_topomap)s + topomap_args : dict | None + Arguments to pass to :func:`mne.viz.plot_topomap`. + + Returns + ------- + fig : instance of matplotlib.figure.Figure + The topoplot figure handle. + """ + picks = pick_types(info, eeg=True) + # fill in lower triangular + 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) + if topomap_args is None: + topomap_args = dict() + if 'vmax' not in topomap_args: + topomap_args['vmax'] = max([elec_dists[idx] for idxs in bridged_idx + for idx in idxs]) * 100 + im, cn = plot_topomap(elec_dists, info, **topomap_args) + # add bridged connections + pos = _find_topomap_coords(info, picks, sphere=sphere) + 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 = im.figure.colorbar(im) + cax.set_label(r'Electrical Distance ($\mu$$V^2$)') From 2f0821efa1d560efd862fbe4d3ba03e973ab43b8 Mon Sep 17 00:00:00 2001 From: Alex Date: Wed, 27 Apr 2022 13:50:50 -0700 Subject: [PATCH 03/12] wip, save Voroni not working --- doc/changes/latest.inc | 2 +- examples/preprocessing/eeg_bridging.py | 170 ++++++++++++++++++------- mne/preprocessing/_csd.py | 30 ++--- mne/viz/topomap.py | 59 +++++++-- 4 files changed, 185 insertions(+), 76 deletions(-) diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index e4938de166b..55d8c2085c9 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -41,7 +41,7 @@ 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:`XX` 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`_) Bugs ~~~~ diff --git a/examples/preprocessing/eeg_bridging.py b/examples/preprocessing/eeg_bridging.py index 2229903ac4c..6782ca93421 100644 --- a/examples/preprocessing/eeg_bridging.py +++ b/examples/preprocessing/eeg_bridging.py @@ -9,20 +9,20 @@ 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 as the signals from the two -(or more) electrodes are not as independent as they would otherwise be; -spatial smearing is caused wherein the signals are more similar to each other -been developed to detect electrode bridging :footcite:`TenkeKayser2001`, which -was previously 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 +"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 more similar to each other than they would otherwise be causing +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 @@ -41,19 +41,22 @@ print(__doc__) # %% -# Let's look at the histograms of electrical distances for the EEGBCI dataset. -# As we can see, for subjects 6, 7 and 8 (and to a lesser extent 4), there is -# a different shape of the distribution of electrical distances than for the -# other subjects. These subjects' distributions have a peak around 0 -# :math:`{\\mu}`V:sup:`2` distance and a trough around 5 -# :math:`{\\mu}`V:sup:`2` which is indicative of electrode bridging. -# The rest of the subjects' distributions increase monotonically, -# indicating normal spatial separation of sources. - -fig, ax = plt.subplots(figsize=(6, 4)) -ax.set_title('Electrical Distance Distribution for EEGBCI Subjects') -ax.set_ylabel('Count') -ax.set_xlabel(r'Electrical Distance ($\mu$$V^2$)') +# 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. + +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 = mne.io.read_raw(mne.datasets.eegbci.load_data( @@ -61,31 +64,108 @@ mne.datasets.eegbci.standardize(raw) # set channel names montage = mne.channels.make_standard_montage('standard_1005') raw.set_montage(montage, verbose=False) + raw_data[sub] = raw + ed_data[sub] = mne.preprocessing.compute_bridged_electrodes(raw) - bridged_idx, ed_matrix = mne.preprocessing.compute_bridged_electrodes(raw) - # 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, 100, 51)) - ax.plot((edges[1:] + edges[:-1]) / 2, hist, - label=f'Sub {sub} #={len(bridged_idx)}') -ax.legend(loc=(1.04, 0)) -fig.subplots_adjust(right=0.725, bottom=0.15) +# %% +# 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 zoom in, we can see several electical +# distance outliers that are near zero; these indicate bridging. + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) +fig.suptitle('Subject 6 Electrical Distance Matrix') +bridged_idx, ed_matrix = ed_data[6] +im1 = ax1.imshow(np.nanmedian(ed_matrix, axis=0)) # take median across epochs +cax1 = fig.colorbar(im1, ax=ax1) +cax1.set_label(r'Electrical Distance ($\mu$$V^2$)') +# zoomed in colors +im2 = ax2.imshow(np.nanmedian(ed_matrix, axis=0), 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() # %% -# Let's look at one subject with bridged electrodes and plot what's -# going on to try and understand electrode bridging better. +# Now let's plot a histogram of the electrical distance matrix. Note that the +# electrical distance matrix is upper triangular but does not include the +# diagonal from the previous plot. This means that the pairwise electrical +# distances are not computed between the same channel (which makes sense +# the differences between a channel and itself would just be zero). The initial +# peak near zero is therefore represents pairs of different channels with +# that are nearly identical which is indicative of bridging. EEG recordings +# without ridged electrodes do not have a peak near zero. -raw = mne.io.read_raw(mne.datasets.eegbci.load_data( - subject=6, runs=(1,))[0], preload=True) -mne.datasets.eegbci.standardize(raw) # set channel names -montage = mne.channels.make_standard_montage('standard_1005') -raw.set_montage(montage) +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') + +# %% +# 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. -bridged_idx, ed_matrix = mne.preprocessing.compute_bridged_electrodes(raw) mne.viz.plot_bridged_electrodes( - raw.info, bridged_idx, ed_matrix, title='Bridged Electrodes', - topomap_args=dict(names=raw.ch_names, show_names=True)) + raw_data[6], bridged_idx, ed_matrix, + title=f'Subject 6 Bridged Electrodes', + topomap_args=dict(names=raw_data[6].ch_names, axes=ax, + image_interp='voroni', vmax=0.05, show_names=True)) + +# %% +# Let's look + +# %% +# Let's look at the histograms of electrical distances for the 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:sup:`2` +# than for the other subjects. These subjects' distributions have a peak around +# 0 :math:`{\\mu}`V:sup:`2` distance and a trough around +# 5 :math:`{\\mu}`V:sup:`2` which is indicative of electrode bridging. +# The rest of the subjects' distributions increase monotonically, +# indicating normal spatial separation of sources. + +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) + +# %% +# Let's look at topoplots with bridged electrodes to try and +# understand electrode bridging better. + +for sub, (bridged_idx, ed_matrix) in ed_data.items(): + fig, ax = plt.subplots() + mne.viz.plot_bridged_electrodes( + infos[sub], bridged_idx, ed_matrix, + title=f'Subject {sub} Bridged Electrodes', + topomap_args=dict(names=infos[sub].ch_names, axes=ax, + image_interp='voroni', vmax=0.05, show_names=True)) # %% # Electrode bridging is often brought about by inserting more gel in order diff --git a/mne/preprocessing/_csd.py b/mne/preprocessing/_csd.py index 1e865e6fdf3..c999dcc3839 100644 --- a/mne/preprocessing/_csd.py +++ b/mne/preprocessing/_csd.py @@ -179,8 +179,8 @@ def compute_current_source_density(inst, sphere='auto', lambda2=1e-5, @verbose -def compute_bridged_electrodes(inst, ed_threshold=3, lm_cutoff=5, - n_bins=20, ed_max=12, epoch_threshold=0.5, +def compute_bridged_electrodes(inst, lm_cutoff=12, kde_args=None, + epoch_threshold=0.5, l_freq=0.5, h_freq=30, epoch_duration=2, verbose=None): """Compute bridged EEG electrodes using the intrinsic Hjorth algorithm. @@ -200,20 +200,15 @@ def compute_bridged_electrodes(inst, ed_threshold=3, lm_cutoff=5, ---------- inst : instance of Raw, Epochs or Evoked The data to compute electrode bridging on. - ed_threshold : float - The bridged electrode threshold in :math:`{\\mu}`V:sup:`2`. The - default is 3 :math:`{\\mu}`V:sup:`2` as used in - :footcite:`GreischarEtAl2004`. lm_cutoff : float The distance in :math:`{\\mu}`V:sup:`2` cutoff below which to search for a local minimum (lm) indicative of bridging. - Defaults to 5 :math:`{\\mu}`V:sup:`2`. - n_bins : int - The number of bins to use to create a histogram of the electrical - distance matrix. - ed_max : float - The maximum electrical distance to consider when making the - histogram of electrical distances. + Defaults to 8 :math:`{\\mu}`V:sup:`2` based on + :footcite:`GreischarEtAl2004`. + kde_args : dict | None + Arguments to pass to :class:`sklearn.neighbors.KernelDensity` to + modify the estimation of distribution of electrical distances + from which the local minimum is estimated. epoch_threshold : float The proportion of epochs with electrical distance less than ``ed_threshold`` in order to consider the channel bridged. @@ -244,6 +239,7 @@ def compute_bridged_electrodes(inst, ed_threshold=3, lm_cutoff=5, ---------- .. footbibliography:: """ + from sklearn.neighbors import KernelDensity _check_preload(inst, 'Computing bridged electrodes') inst = inst.copy() # don't modify original picks = pick_types(inst.info, eeg=True) @@ -278,14 +274,10 @@ def compute_bridged_electrodes(inst, ed_threshold=3, lm_cutoff=5, # compute histogram bins, edges = np.histogram(ed_matrix[~np.isnan(ed_matrix)], - bins=np.linspace(0, ed_max, n_bins)) + bins=np.linspace(0, lm_cutoff, n_bins)) centers = (edges[1:] + edges[:-1]) / 2 - # check if has peak near 0, if not, no bridges - if np.diff(bins[centers <= ed_threshold]).sum() >= 0: - return bridged_idx, ed_matrix - - # find local minimum, also indicative of bridging + # find local minimum centers_interp = np.linspace(0, lm_cutoff, n_bins) bins_interp = np.interp(centers_interp, centers, bins) local_minimum = centers_interp[np.argmin(bins_interp)] diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py index 4da70e3d502..a308f1a3696 100644 --- a/mne/viz/topomap.py +++ b/mne/viz/topomap.py @@ -2621,8 +2621,41 @@ def plot_arrowmap(data, info_from, info_to=None, scale=3e-10, vmin=None, return fig +def _voroni_topomap(data, pos, info, sphere, ch_type, outlines, ax, cmap): + """Make a Voroni diagram on a topomap.""" + from scipy.spatial import Voronoi + clip_origin = _adjust_meg_sphere(sphere, info, ch_type)[1] + outlines = _make_head_outlines( + sphere, pos, outlines, clip_origin) + rx, ry = outlines['clip_radius'] + cx, cy = clip_origin + vor = Voronoi(np.concatenate([pos, [[1, 1], [1, -1], [-1, 1], [1, -1]]])) + for point_idx, region_idx in enumerate(vor.point_region[:-4]): + polygon = list() + if -1 in vor.regions[region_idx]: + continue + for i in vor.regions[region_idx]: + if i == -1: + # outer bound, take each other point and project to edge + for j in vor.regions[region_idx]: + if j != -1: + x, y = vor.vertices[j] + x *= rx / np.linalg.norm(vor.vertices[j]) + y *= ry / np.linalg.norm(vor.vertices[j]) + polygon.append((x, y)) + else: + 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(data[point_idx])) + + def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, - colorbar=True, sphere=None, topomap_args=None): + topomap_args=None): """Topoplot electrode distance matrix with bridged electrodes connected. Parameters @@ -2635,9 +2668,6 @@ def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, The electrical distance matrix for each pair of EEG electrodes. title : str A title to add to the plot. - colorbar : bool - Whether to add a colorbar. - %(sphere_topomap)s topomap_args : dict | None Arguments to pass to :func:`mne.viz.plot_topomap`. @@ -2646,20 +2676,27 @@ def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, fig : instance of matplotlib.figure.Figure The topoplot figure handle. """ + # 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) # fill in lower triangular 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) - if topomap_args is None: - topomap_args = dict() - if 'vmax' not in topomap_args: - topomap_args['vmax'] = max([elec_dists[idx] for idxs in bridged_idx - for idx in idxs]) * 100 - im, cn = plot_topomap(elec_dists, info, **topomap_args) - # add bridged connections pos = _find_topomap_coords(info, picks, sphere=sphere) + if 'image_interp' not in topomap_args: # default is Vononi + im, cn = plot_topomap(np.zeros((picks.size,)), info, **topomap_args) + # plot Vonorni Diagram over topomap + _voroni_topomap(elec_dists, pos, info, sphere, 'eeg', + topomap_args.get('outlines', 'head'), im.axes, + im.cmap) + else: + im, cn = plot_topomap(elec_dists, info, **topomap_args) + # 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') From a81afe467728ba5495b67241efb5f9f588cb02f7 Mon Sep 17 00:00:00 2001 From: Alex Date: Wed, 27 Apr 2022 18:07:25 -0700 Subject: [PATCH 04/12] working version --- examples/preprocessing/eeg_bridging.py | 96 +++++++++++++++++++------- mne/preprocessing/_csd.py | 60 ++++++++-------- mne/preprocessing/tests/test_csd.py | 22 +++++- mne/viz/tests/test_topomap.py | 26 ++++++- mne/viz/topomap.py | 54 ++++++++------- 5 files changed, 175 insertions(+), 83 deletions(-) diff --git a/examples/preprocessing/eeg_bridging.py b/examples/preprocessing/eeg_bridging.py index 6782ca93421..7e14ff053cb 100644 --- a/examples/preprocessing/eeg_bridging.py +++ b/examples/preprocessing/eeg_bridging.py @@ -75,12 +75,12 @@ # 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 zoom in, we can see several electical +# recording. On the right, when we 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') -bridged_idx, ed_matrix = ed_data[6] im1 = ax1.imshow(np.nanmedian(ed_matrix, axis=0)) # take median across epochs cax1 = fig.colorbar(im1, ax=ax1) cax1.set_label(r'Electrical Distance ($\mu$$V^2$)') @@ -112,27 +112,61 @@ # %% # 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. - +# 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). + +fig, ax = plt.subplots() mne.viz.plot_bridged_electrodes( - raw_data[6], bridged_idx, ed_matrix, + raw_data[6].info, bridged_idx, ed_matrix, title=f'Subject 6 Bridged Electrodes', topomap_args=dict(names=raw_data[6].ch_names, axes=ax, - image_interp='voroni', vmax=0.05, show_names=True)) + vmax=5, show_names=True)) # %% -# Let's look +# 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)) # %% -# Let's look at the histograms of electrical distances for the 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:sup:`2` -# than for the other subjects. These subjects' distributions have a peak around -# 0 :math:`{\\mu}`V:sup:`2` distance and a trough around -# 5 :math:`{\\mu}`V:sup:`2` which is indicative of electrode bridging. -# The rest of the subjects' distributions increase monotonically, -# indicating normal spatial separation of sources. +# 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:sup:`2` than for the other subjects. These subjects' +# distributions have a peak around 0 :math:`{\\mu}`V:sup:`2` distance +# and a trough around 5 :math:`{\\mu}`V:sup:`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') @@ -156,35 +190,47 @@ fig.subplots_adjust(right=0.725, bottom=0.15, wspace=0.4) # %% -# Let's look at topoplots with bridged electrodes to try and -# understand electrode bridging better. +# 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(): fig, ax = plt.subplots() mne.viz.plot_bridged_electrodes( - infos[sub], bridged_idx, ed_matrix, + raw_data[sub].info, bridged_idx, ed_matrix, title=f'Subject {sub} Bridged Electrodes', - topomap_args=dict(names=infos[sub].ch_names, axes=ax, - image_interp='voroni', vmax=0.05, show_names=True)) + topomap_args=dict(names=raw_data[sub].ch_names, axes=ax, + vmax=5, show_names=True)) # %% # 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. Impedances are recommended to be stored in an EEG dataset -# in the :ref:`electrodes-tsv` file within th eBrain Imaging Data Structure -# (BIDS). Since the impedances are not stored for this dataset, we will fake +# spatially precise. Brain Imaging Data Structure (BIDS) recommendeds that +# impedances be stored in an EEG dataset in the :ref:`electrodes-tsv` file. +# Since the impedances are not stored for this dataset, we will fake # them to demonstrate how they would be plotted. np.random.seed(11) # seed for reproducibility impedances = np.random.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 Audio/Visual Task') +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 diff --git a/mne/preprocessing/_csd.py b/mne/preprocessing/_csd.py index c999dcc3839..54524fdbdc5 100644 --- a/mne/preprocessing/_csd.py +++ b/mne/preprocessing/_csd.py @@ -179,10 +179,9 @@ def compute_current_source_density(inst, sphere='auto', lambda2=1e-5, @verbose -def compute_bridged_electrodes(inst, lm_cutoff=12, kde_args=None, - epoch_threshold=0.5, - l_freq=0.5, h_freq=30, - epoch_duration=2, verbose=None): +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): """Compute bridged EEG electrodes using the intrinsic Hjorth algorithm. First an electrical distance matrix is computed by taking the pairwise @@ -203,12 +202,8 @@ def compute_bridged_electrodes(inst, lm_cutoff=12, kde_args=None, lm_cutoff : float The distance in :math:`{\\mu}`V:sup:`2` cutoff below which to search for a local minimum (lm) indicative of bridging. - Defaults to 8 :math:`{\\mu}`V:sup:`2` based on - :footcite:`GreischarEtAl2004`. - kde_args : dict | None - Arguments to pass to :class:`sklearn.neighbors.KernelDensity` to - modify the estimation of distribution of electrical distances - from which the local minimum is estimated. + Defaults to 16 :math:`{\\mu}`V:sup:`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. @@ -221,6 +216,8 @@ def compute_bridged_electrodes(inst, lm_cutoff=12, kde_args=None, 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 :ref:`scipy.stats.gaussian_kde`. %(verbose)s Returns @@ -239,7 +236,8 @@ def compute_bridged_electrodes(inst, lm_cutoff=12, kde_args=None, ---------- .. footbibliography:: """ - from sklearn.neighbors import KernelDensity + 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) @@ -254,17 +252,16 @@ def compute_bridged_electrodes(inst, lm_cutoff=12, kde_args=None, preload=True, verbose=False) # standardize shape - data = inst.get_data() + data = inst.get_data(picks=picks) if isinstance(inst, Evoked): data = data.reshape((1,) + data.shape) # 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, idx0 in enumerate(picks): - for j, idx1 in enumerate(picks[i + 1:]): - ed_matrix[:, i, i + 1 + j] = \ - np.var(data[:, idx0] - data[:, idx1], axis=1) + 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 @@ -272,26 +269,25 @@ def compute_bridged_electrodes(inst, lm_cutoff=12, kde_args=None, # initialize bridged indices bridged_idx = list() - # compute histogram - bins, edges = np.histogram(ed_matrix[~np.isnan(ed_matrix)], - bins=np.linspace(0, lm_cutoff, n_bins)) - centers = (edges[1:] + edges[:-1]) / 2 + # if not enough values below local minimum cutoff, return no bridges + if ed_matrix[ed_matrix < lm_cutoff].size / n_epochs < epoch_threshold: + return bridged_idx, ed_matrix - # find local minimum - centers_interp = np.linspace(0, lm_cutoff, n_bins) - bins_interp = np.interp(centers_interp, centers, bins) - local_minimum = centers_interp[np.argmin(bins_interp)] + # kernel density estimation + kde = gaussian_kde(ed_matrix[ed_matrix < lm_cutoff]) + 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 minumum on + # find electrodes that are below the cutoff local minimum on # `epochs_threshold` proportion of epochs check_bridged = ed_matrix < local_minimum - for i, idx0 in enumerate(picks): - for j, idx1 in enumerate(picks[i + 1:]): - if np.sum(check_bridged[:, i, i + 1 + j]) / \ - n_epochs > epoch_threshold: + for i in range(picks.size): + for j in range(i + 1, picks.size): + if np.sum(check_bridged[:, i, j]) / n_epochs > epoch_threshold: logger.info('Bridge detected between ' - f'{inst.ch_names[picks[idx0]]} and ' - f'{inst.ch_names[picks[idx1]]}') - bridged_idx.append((idx0, idx1)) + 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 d29694cadb3..9c75b1a6352 100644 --- a/mne/preprocessing/tests/test_csd.py +++ b/mne/preprocessing/tests/test_csd.py @@ -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') @@ -187,3 +188,22 @@ def test_csd_fif(): 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 + 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) + assert bridged_idx == [(idx0, idx1)] + picks = pick_types(raw.info, meg=False, eeg=True) + assert ed_matrix.shape == (raw.times.size // (2 * 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/viz/tests/test_topomap.py b/mne/viz/tests/test_topomap.py index b9b79137bbb..e174f9ed2a9 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 @@ -667,3 +668,26 @@ 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.""" + np.random.seed(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] = \ + np.random.random() + np.random.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 a308f1a3696..f340f551cbd 100644 --- a/mne/viz/topomap.py +++ b/mne/viz/topomap.py @@ -2621,39 +2621,36 @@ def plot_arrowmap(data, info_from, info_to=None, scale=3e-10, vmin=None, return fig -def _voroni_topomap(data, pos, info, sphere, ch_type, outlines, ax, cmap): +def _voroni_topomap(data, pos, info, sphere, ch_type, outlines, ax, cmap, + norm): """Make a Voroni diagram on a topomap.""" from scipy.spatial import Voronoi + sphere = _check_sphere(sphere) clip_origin = _adjust_meg_sphere(sphere, info, ch_type)[1] outlines = _make_head_outlines( sphere, pos, outlines, clip_origin) rx, ry = outlines['clip_radius'] cx, cy = clip_origin - vor = Voronoi(np.concatenate([pos, [[1, 1], [1, -1], [-1, 1], [1, -1]]])) - for point_idx, region_idx in enumerate(vor.point_region[:-4]): - polygon = list() + # add faroff points in a circle + vor = Voronoi(np.concatenate([pos, [(np.cos(2 * np.pi / 100 * t), + np.sin(2 * np.pi / 100 * t)) + for t in range(101)]])) + for point_idx, region_idx in enumerate(vor.point_region[:-101]): if -1 in vor.regions[region_idx]: continue + polygon = list() for i in vor.regions[region_idx]: - if i == -1: - # outer bound, take each other point and project to edge - for j in vor.regions[region_idx]: - if j != -1: - x, y = vor.vertices[j] - x *= rx / np.linalg.norm(vor.vertices[j]) - y *= ry / np.linalg.norm(vor.vertices[j]) - polygon.append((x, y)) + x, y = vor.vertices[i] + if (x - cx)**2 / rx**2 + (y - cy)**2 / ry**2 < 1: + polygon.append((x, y)) else: - 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(data[point_idx])) + 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]))) +@fill_doc def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, topomap_args=None): """Topoplot electrode distance matrix with bridged electrodes connected. @@ -2676,26 +2673,34 @@ def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, fig : instance of matplotlib.figure.Figure The topoplot figure handle. """ + if topomap_args is None: + topomap_args = dict() # 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) - if 'image_interp' not in topomap_args: # default is Vononi + if 'image_interp' in topomap_args: + im, cn = plot_topomap(elec_dists, info, **topomap_args) + else: # default is Vononi im, cn = plot_topomap(np.zeros((picks.size,)), info, **topomap_args) # plot Vonorni Diagram over topomap _voroni_topomap(elec_dists, pos, info, sphere, 'eeg', topomap_args.get('outlines', 'head'), im.axes, - im.cmap) - else: - im, cn = plot_topomap(elec_dists, info, **topomap_args) + im.cmap, im.norm) + # add bridged connections for idx0, idx1 in bridged_idx: im.axes.plot([pos[idx0][0], pos[idx1][0]], @@ -2705,3 +2710,4 @@ def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, if colorbar: cax = im.figure.colorbar(im) cax.set_label(r'Electrical Distance ($\mu$$V^2$)') + return im.figure From 4c0fd897427b55a14ef87b8f19401daaaa017d90 Mon Sep 17 00:00:00 2001 From: Alex Date: Wed, 27 Apr 2022 18:13:22 -0700 Subject: [PATCH 05/12] fix flake fix runtime warning divide by zero and 3.7 nan handling fix warning filter fix pydocstyle fix links, refs fix forgot to delete one more error fixed try to fix microvolts squared that's not going to work, trying something else --- examples/preprocessing/eeg_bridging.py | 19 ++++++++++------ mne/preprocessing/_csd.py | 30 +++++++++++++++----------- 2 files changed, 31 insertions(+), 18 deletions(-) diff --git a/examples/preprocessing/eeg_bridging.py b/examples/preprocessing/eeg_bridging.py index 7e14ff053cb..06052b53d8c 100644 --- a/examples/preprocessing/eeg_bridging.py +++ b/examples/preprocessing/eeg_bridging.py @@ -33,6 +33,7 @@ # sphinx_gallery_thumbnail_number = 2 +import warnings import numpy as np import matplotlib.pyplot as plt @@ -81,7 +82,12 @@ bridged_idx, ed_matrix = ed_data[6] fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) fig.suptitle('Subject 6 Electrical Distance Matrix') -im1 = ax1.imshow(np.nanmedian(ed_matrix, axis=0)) # take median across epochs +with warnings.catch_warnings(): + # lower triangular matrix is NaNs so we should to suppress the warning + warnings.filterwarnings( + action='ignore', message='All-NaN slice encountered') + # take median across epochs + im1 = ax1.imshow(np.nanmedian(ed_matrix, axis=0)) cax1 = fig.colorbar(im1, ax=ax1) cax1.set_label(r'Electrical Distance ($\mu$$V^2$)') # zoomed in colors @@ -124,7 +130,7 @@ fig, ax = plt.subplots() mne.viz.plot_bridged_electrodes( raw_data[6].info, bridged_idx, ed_matrix, - title=f'Subject 6 Bridged Electrodes', + title='Subject 6 Bridged Electrodes', topomap_args=dict(names=raw_data[6].ch_names, axes=ax, vmax=5, show_names=True)) @@ -211,11 +217,15 @@ # 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) recommendeds that -# impedances be stored in an EEG dataset in the :ref:`electrodes-tsv` file. +# 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. np.random.seed(11) # seed for reproducibility +raw = raw_data[1] impedances = np.random.random((len(raw.ch_names,))) * 10 fig, ax = plt.subplots(figsize=(5, 5)) im, cn = mne.viz.plot_topomap(impedances, raw.info, axes=ax) @@ -236,6 +246,3 @@ # References # ---------- # .. footbibliography:: -# -# -# .. _electrodes-tsv: https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/03-electroencephalography.html#electrodes-description-_electrodestsv # noqa E501 diff --git a/mne/preprocessing/_csd.py b/mne/preprocessing/_csd.py index 54524fdbdc5..b33487bc2a1 100644 --- a/mne/preprocessing/_csd.py +++ b/mne/preprocessing/_csd.py @@ -12,6 +12,7 @@ # 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 @@ -182,18 +183,18 @@ def compute_current_source_density(inst, sphere='auto', lambda2=1e-5, 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): - """Compute bridged EEG electrodes using the intrinsic Hjorth algorithm. + 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:sup:`2` and a - a local minimum below 5 :math:`{\\mu}`V:sup:`2` are found, the presence + variance. Then, a local maximum near 0 :math:`\mu`V:sup:`2` and a + a local minimum below 5 :math:`\mu`V:sup:`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`. - Original implementation - https://psychophysiology.cpmc.columbia.edu/software/eBridge/index.html. + Based on :footcite:`TenkeKayser2001,GreischarEtAl2004,DelormeMakeig2004` + and the `EEGLAB implementation + `_. Parameters ---------- @@ -217,7 +218,7 @@ def compute_bridged_electrodes(inst, lm_cutoff=16, epoch_threshold=0.5, 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 :ref:`scipy.stats.gaussian_kde`. + ``bw_method`` to pass to :class:`scipy.stats.gaussian_kde`. %(verbose)s Returns @@ -270,18 +271,23 @@ def compute_bridged_electrodes(inst, lm_cutoff=16, epoch_threshold=0.5, bridged_idx = list() # if not enough values below local minimum cutoff, return no bridges - if ed_matrix[ed_matrix < lm_cutoff].size / n_epochs < epoch_threshold: + 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_matrix[ed_matrix < lm_cutoff]) - local_minimum = float(minimize_scalar( - lambda x: kde(x) if x < lm_cutoff and x > 0 else np.inf).x) + 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 - check_bridged = ed_matrix < local_minimum + check_bridged = np.logical_and(~np.isnan(ed_matrix), + ed_matrix < local_minimum) for i in range(picks.size): for j in range(i + 1, picks.size): if np.sum(check_bridged[:, i, j]) / n_epochs > epoch_threshold: From 8f85d1b4c83d344bf3fd576971ee0baa929dfb9a Mon Sep 17 00:00:00 2001 From: Alex Date: Thu, 28 Apr 2022 10:09:42 -0700 Subject: [PATCH 06/12] alex review --- examples/preprocessing/eeg_bridging.py | 29 +++++++++++++------------- mne/preprocessing/_csd.py | 9 ++++---- mne/preprocessing/tests/test_csd.py | 9 +++++--- mne/viz/tests/test_topomap.py | 9 ++++---- mne/viz/topomap.py | 12 +++++------ 5 files changed, 34 insertions(+), 34 deletions(-) diff --git a/examples/preprocessing/eeg_bridging.py b/examples/preprocessing/eeg_bridging.py index 06052b53d8c..354c4baa46d 100644 --- a/examples/preprocessing/eeg_bridging.py +++ b/examples/preprocessing/eeg_bridging.py @@ -33,7 +33,6 @@ # sphinx_gallery_thumbnail_number = 2 -import warnings import numpy as np import matplotlib.pyplot as plt @@ -82,16 +81,16 @@ bridged_idx, ed_matrix = ed_data[6] fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) fig.suptitle('Subject 6 Electrical Distance Matrix') -with warnings.catch_warnings(): - # lower triangular matrix is NaNs so we should to suppress the warning - warnings.filterwarnings( - action='ignore', message='All-NaN slice encountered') - # take median across epochs - im1 = ax1.imshow(np.nanmedian(ed_matrix, axis=0)) +# 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]) +im1 = ax1.imshow(ed_plot, aspect='auto') cax1 = fig.colorbar(im1, ax=ax1) cax1.set_label(r'Electrical Distance ($\mu$$V^2$)') # zoomed in colors -im2 = ax2.imshow(np.nanmedian(ed_matrix, axis=0), vmax=5) +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): @@ -103,11 +102,11 @@ # Now let's plot a histogram of the electrical distance matrix. Note that the # electrical distance matrix is upper triangular but does not include the # diagonal from the previous plot. This means that the pairwise electrical -# distances are not computed between the same channel (which makes sense +# 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 is therefore represents pairs of different channels with -# that are nearly identical which is indicative of bridging. EEG recordings -# without ridged electrodes do not have a peak near zero. +# 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') @@ -216,7 +215,7 @@ # 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) recommendeds that +# spatially precise. Brain Imaging Data Structure (BIDS) recommends that # impedances be stored in an EEG dataset in the `electrodes.tsv # epoch_threshold: + 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]]}') diff --git a/mne/preprocessing/tests/test_csd.py b/mne/preprocessing/tests/test_csd.py index 9c75b1a6352..4c339e132ca 100644 --- a/mne/preprocessing/tests/test_csd.py +++ b/mne/preprocessing/tests/test_csd.py @@ -195,15 +195,18 @@ def test_compute_bridged_electrodes(): 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) + 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 // (2 * raw.info['sfreq']), - picks.size, picks.size) + 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/viz/tests/test_topomap.py b/mne/viz/tests/test_topomap.py index e174f9ed2a9..67392cc8733 100644 --- a/mne/viz/tests/test_topomap.py +++ b/mne/viz/tests/test_topomap.py @@ -648,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") @@ -672,7 +672,7 @@ def test_plot_topomap_cnorm(): def test_plot_bridged_electrodes(): """Test plotting of bridged electrodes.""" - np.random.seed(42) + 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)] @@ -681,8 +681,7 @@ def test_plot_bridged_electrodes(): 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] = \ - np.random.random() + np.random.random(triu_idx[0].size) + 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)) diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py index f340f551cbd..1bd60fd9863 100644 --- a/mne/viz/topomap.py +++ b/mne/viz/topomap.py @@ -2621,9 +2621,9 @@ def plot_arrowmap(data, info_from, info_to=None, scale=3e-10, vmin=None, return fig -def _voroni_topomap(data, pos, info, sphere, ch_type, outlines, ax, cmap, - norm): - """Make a Voroni diagram on a topomap.""" +def _voronoi_topomap(data, pos, info, sphere, ch_type, outlines, ax, cmap, + norm): + """Make a Voronoi diagram on a topomap.""" from scipy.spatial import Voronoi sphere = _check_sphere(sphere) clip_origin = _adjust_meg_sphere(sphere, info, ch_type)[1] @@ -2697,9 +2697,9 @@ def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, else: # default is Vononi im, cn = plot_topomap(np.zeros((picks.size,)), info, **topomap_args) # plot Vonorni Diagram over topomap - _voroni_topomap(elec_dists, pos, info, sphere, 'eeg', - topomap_args.get('outlines', 'head'), im.axes, - im.cmap, im.norm) + _voronoi_topomap(elec_dists, pos, info, sphere, 'eeg', + topomap_args.get('outlines', 'head'), im.axes, + im.cmap, im.norm) # add bridged connections for idx0, idx1 in bridged_idx: From 88998145bb055df6e9ef2af301512c2da0210042 Mon Sep 17 00:00:00 2001 From: Alex Date: Thu, 28 Apr 2022 14:12:00 -0700 Subject: [PATCH 07/12] a few last aethetic fixes (muV^2) and add voronoi as option for image interpolation --- examples/preprocessing/eeg_bridging.py | 30 ++++--- mne/preprocessing/_csd.py | 8 +- mne/utils/docs.py | 2 +- mne/viz/topomap.py | 114 ++++++++++++++----------- 4 files changed, 88 insertions(+), 66 deletions(-) diff --git a/examples/preprocessing/eeg_bridging.py b/examples/preprocessing/eeg_bridging.py index 354c4baa46d..caf1f02b326 100644 --- a/examples/preprocessing/eeg_bridging.py +++ b/examples/preprocessing/eeg_bridging.py @@ -41,6 +41,8 @@ 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 @@ -69,6 +71,8 @@ # %% +# 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 @@ -99,6 +103,8 @@ fig.tight_layout() # %% +# Examine the Distibution of Electrical Distances +# ----------------------------------------------- # Now let's plot a histogram of the electrical distance matrix. Note that the # electrical distance matrix is upper triangular but does not include the # diagonal from the previous plot. This means that the pairwise electrical @@ -115,6 +121,8 @@ ax.set_ylabel('Count') # %% +# 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 @@ -126,14 +134,13 @@ # (this may be because the EEG experimenter usually stands to the side and # may have inserted the gel syringe tip in too far). -fig, ax = plt.subplots() mne.viz.plot_bridged_electrodes( raw_data[6].info, bridged_idx, ed_matrix, - title='Subject 6 Bridged Electrodes', - topomap_args=dict(names=raw_data[6].ch_names, axes=ax, - vmax=5, show_names=True)) + 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 @@ -158,13 +165,15 @@ raw.plot(duration=20, scalings=dict(eeg=2e-4)) # %% +# Compare Bridging Accross 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:sup:`2` than for the other subjects. These subjects' -# distributions have a peak around 0 :math:`{\\mu}`V:sup:`2` distance -# and a trough around 5 :math:`{\\mu}`V:sup:`2` which is indicative of +# 0 :math:`\muV^2` than for the other subjects. These subjects' +# distributions have a peak around 0 :math:`\muV^2` distance +# and a trough around 5 :math:`\muV^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 @@ -203,14 +212,13 @@ # other conductive electrolyte solution. for sub, (bridged_idx, ed_matrix) in ed_data.items(): - fig, ax = plt.subplots() mne.viz.plot_bridged_electrodes( raw_data[sub].info, bridged_idx, ed_matrix, - title=f'Subject {sub} Bridged Electrodes', - topomap_args=dict(names=raw_data[sub].ch_names, axes=ax, - vmax=5, show_names=True)) + 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 diff --git a/mne/preprocessing/_csd.py b/mne/preprocessing/_csd.py index d4058014b47..67564dce126 100644 --- a/mne/preprocessing/_csd.py +++ b/mne/preprocessing/_csd.py @@ -186,8 +186,8 @@ def compute_bridged_electrodes(inst, lm_cutoff=16, epoch_threshold=0.5, 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:sup:`2` and a - a local minimum below 5 :math:`\mu`V:sup:`2` are found, the presence + variance. Then, a local maximum near 0 :math:`\muV^2` and a + a local minimum below 5 :math:`\muV^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. @@ -201,9 +201,9 @@ def compute_bridged_electrodes(inst, lm_cutoff=16, epoch_threshold=0.5, inst : instance of Raw, Epochs or Evoked The data to compute electrode bridging on. lm_cutoff : float - The distance in :math:`{\\mu}`V:sup:`2` cutoff below which to + The distance in :math:`\muV^2` cutoff below which to search for a local minimum (lm) indicative of bridging. - Defaults to 16 :math:`{\\mu}`V:sup:`2` to be conservative based + Defaults to 16 :math:`\muV^2` to be conservative based on the distributions in :footcite:`GreischarEtAl2004`. epoch_threshold : float The proportion of epochs with electrical distance less than 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/topomap.py b/mne/viz/topomap.py index 1bd60fd9863..a0a30795a2b 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 emtpy 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 @@ -2621,35 +2658,6 @@ def plot_arrowmap(data, info_from, info_to=None, scale=3e-10, vmin=None, return fig -def _voronoi_topomap(data, pos, info, sphere, ch_type, outlines, ax, cmap, - norm): - """Make a Voronoi diagram on a topomap.""" - from scipy.spatial import Voronoi - sphere = _check_sphere(sphere) - clip_origin = _adjust_meg_sphere(sphere, info, ch_type)[1] - outlines = _make_head_outlines( - sphere, pos, outlines, clip_origin) - rx, ry = outlines['clip_radius'] - cx, cy = clip_origin - # add faroff points in a circle - vor = Voronoi(np.concatenate([pos, [(np.cos(2 * np.pi / 100 * t), - np.sin(2 * np.pi / 100 * t)) - for t in range(101)]])) - for point_idx, region_idx in enumerate(vor.point_region[:-101]): - 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]))) - - @fill_doc def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, topomap_args=None): @@ -2673,8 +2681,21 @@ def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, fig : instance of matplotlib.figure.Figure The topoplot figure handle. """ + import matplotlib.pyplot as plt if topomap_args is None: topomap_args = dict() + if 'image_interp' not in topomap_args: + topomap_args['image_interp'] = 'voronoi' + if 'names' not in topomap_args: + topomap_args['names'] = info.ch_names + if 'show_names' not in topomap_args: + topomap_args['show_names'] = True + if 'contours' not in topomap_args: + topomap_args['contours'] = False + if 'axes' not in topomap_args: + fig, ax = plt.subplots() + else: + fig = None # handle colorbar here instead of in plot_topomap colorbar = topomap_args.pop('colorbar') if \ 'colorbar' in topomap_args else True @@ -2692,15 +2713,8 @@ def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, 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) - if 'image_interp' in topomap_args: - im, cn = plot_topomap(elec_dists, info, **topomap_args) - else: # default is Vononi - im, cn = plot_topomap(np.zeros((picks.size,)), info, **topomap_args) - # plot Vonorni Diagram over topomap - _voronoi_topomap(elec_dists, pos, info, sphere, 'eeg', - topomap_args.get('outlines', 'head'), im.axes, - im.cmap, im.norm) - + 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]], @@ -2708,6 +2722,6 @@ def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, if title is not None: im.axes.set_title(title) if colorbar: - cax = im.figure.colorbar(im) + cax = fig.colorbar(im) cax.set_label(r'Electrical Distance ($\mu$$V^2$)') - return im.figure + return fig From 8a8ed9cfc24b9c5e6c350a661a28ba0d2fb002af Mon Sep 17 00:00:00 2001 From: Alex Date: Thu, 28 Apr 2022 14:13:55 -0700 Subject: [PATCH 08/12] update latest --- doc/changes/latest.inc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 55d8c2085c9..a0cc680a72a 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -43,6 +43,8 @@ Enhancements - 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`) From c1f62694c209edf360fa43f0cc9323d8a33ebe74 Mon Sep 17 00:00:00 2001 From: Alex Date: Thu, 28 Apr 2022 14:16:43 -0700 Subject: [PATCH 09/12] fix flake --- examples/preprocessing/eeg_bridging.py | 6 +++--- mne/viz/topomap.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/preprocessing/eeg_bridging.py b/examples/preprocessing/eeg_bridging.py index caf1f02b326..0569c3226fb 100644 --- a/examples/preprocessing/eeg_bridging.py +++ b/examples/preprocessing/eeg_bridging.py @@ -103,8 +103,8 @@ fig.tight_layout() # %% -# Examine the Distibution of Electrical Distances -# ----------------------------------------------- +# Examine the Distribution of Electrical Distances +# ------------------------------------------------ # Now let's plot a histogram of the electrical distance matrix. Note that the # electrical distance matrix is upper triangular but does not include the # diagonal from the previous plot. This means that the pairwise electrical @@ -165,7 +165,7 @@ raw.plot(duration=20, scalings=dict(eeg=2e-4)) # %% -# Compare Bridging Accross Subjects in the EEGBCI Dataset +# 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, diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py index a0a30795a2b..d3927ae049f 100644 --- a/mne/viz/topomap.py +++ b/mne/viz/topomap.py @@ -838,7 +838,7 @@ def _setup_interp(pos, res, extrapolate, sphere, outlines, border): 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 emtpy image to plot over + # 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) From 254db4e9231715e4630ec535dba7400de8aa3cd1 Mon Sep 17 00:00:00 2001 From: Alex Date: Thu, 28 Apr 2022 19:45:13 -0700 Subject: [PATCH 10/12] fix muV2 --- examples/preprocessing/eeg_bridging.py | 6 +++--- mne/preprocessing/_csd.py | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/examples/preprocessing/eeg_bridging.py b/examples/preprocessing/eeg_bridging.py index 0569c3226fb..956dc6c5027 100644 --- a/examples/preprocessing/eeg_bridging.py +++ b/examples/preprocessing/eeg_bridging.py @@ -171,9 +171,9 @@ # 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:`\muV^2` than for the other subjects. These subjects' -# distributions have a peak around 0 :math:`\muV^2` distance -# and a trough around 5 :math:`\muV^2` which is indicative of +# 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 diff --git a/mne/preprocessing/_csd.py b/mne/preprocessing/_csd.py index 67564dce126..988a3ffd718 100644 --- a/mne/preprocessing/_csd.py +++ b/mne/preprocessing/_csd.py @@ -186,8 +186,8 @@ def compute_bridged_electrodes(inst, lm_cutoff=16, epoch_threshold=0.5, 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:`\muV^2` and a - a local minimum below 5 :math:`\muV^2` are found, the presence + 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. @@ -201,9 +201,9 @@ def compute_bridged_electrodes(inst, lm_cutoff=16, epoch_threshold=0.5, inst : instance of Raw, Epochs or Evoked The data to compute electrode bridging on. lm_cutoff : float - The distance in :math:`\muV^2` cutoff below which to + The distance in :math:`{\mu}V^2` cutoff below which to search for a local minimum (lm) indicative of bridging. - Defaults to 16 :math:`\muV^2` to be conservative based + 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 From 3f510cf562ebdda5a31069d7c6e6710d9ab65c6b Mon Sep 17 00:00:00 2001 From: Alex Date: Fri, 29 Apr 2022 08:02:29 -0700 Subject: [PATCH 11/12] mikolaj review --- examples/preprocessing/eeg_bridging.py | 23 +++++++++++++++-------- mne/viz/topomap.py | 15 +++++++-------- 2 files changed, 22 insertions(+), 16 deletions(-) diff --git a/examples/preprocessing/eeg_bridging.py b/examples/preprocessing/eeg_bridging.py index 956dc6c5027..91037f70cda 100644 --- a/examples/preprocessing/eeg_bridging.py +++ b/examples/preprocessing/eeg_bridging.py @@ -11,7 +11,7 @@ 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 more similar to each other than they would otherwise be causing +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 @@ -57,14 +57,14 @@ # 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 = mne.io.read_raw(mne.datasets.eegbci.load_data( - subject=sub, runs=(1,))[0], preload=True, verbose=False) + 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 - montage = mne.channels.make_standard_montage('standard_1005') raw.set_montage(montage, verbose=False) raw_data[sub] = raw ed_data[sub] = mne.preprocessing.compute_bridged_electrodes(raw) @@ -79,27 +79,34 @@ # 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 zoom in, we can see several electrical -# distance outliers that are near zero; these indicate bridging. +# 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$)') -# zoomed in colors + +# 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() # %% @@ -118,7 +125,7 @@ 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') +ax.set_ylabel('Count (channel pairs for all epochs)') # %% # Plot Electrical Distances on a Topomap diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py index d3927ae049f..ad345c45a82 100644 --- a/mne/viz/topomap.py +++ b/mne/viz/topomap.py @@ -2684,16 +2684,15 @@ def plot_bridged_electrodes(info, bridged_idx, ed_matrix, title=None, import matplotlib.pyplot as plt if topomap_args is None: topomap_args = dict() - if 'image_interp' not in topomap_args: - topomap_args['image_interp'] = 'voronoi' - if 'names' not in topomap_args: - topomap_args['names'] = info.ch_names - if 'show_names' not in topomap_args: - topomap_args['show_names'] = True - if 'contours' not in topomap_args: - topomap_args['contours'] = False + 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 From 1ba55d513e06f7adc5afc349bf55de72a60a6d88 Mon Sep 17 00:00:00 2001 From: Alex Date: Fri, 29 Apr 2022 12:31:38 -0700 Subject: [PATCH 12/12] clarify wording --- examples/preprocessing/eeg_bridging.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/preprocessing/eeg_bridging.py b/examples/preprocessing/eeg_bridging.py index 91037f70cda..1c7f4ee36fa 100644 --- a/examples/preprocessing/eeg_bridging.py +++ b/examples/preprocessing/eeg_bridging.py @@ -113,8 +113,8 @@ # Examine the Distribution of Electrical Distances # ------------------------------------------------ # Now let's plot a histogram of the electrical distance matrix. Note that the -# electrical distance matrix is upper triangular but does not include the -# diagonal from the previous plot. This means that the pairwise electrical +# 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