Introduce lower bounds for ExG peak detection thresholds

doc/changes/latest.inc

@@ -297,6 +297,8 @@ Bugs
 
 - Fix bug where ``backend='notebook'`` could not be used in :meth:`mne.SourceEstimate.plot` (:gh:`9305` by `Jean-Remi King`_)
 
+- `mne.preprocessing.find_eog_events`, `mne.preprocessing.find_ecg_events`, `mne.preprocessing.create_eog_epochs`, and `mne.preprocessing.create_ecg_epochs` now have better handling for situations where the signal is flat or the provided threshold value is too small due to floating point imprecisions (:gh:`xxx` by `Richard Höchenberger`_)
+
 - `mne.preprocessing.compute_proj_eog` and `mne.preprocessing.compute_proj_ecg` now return empty lists if no EOG or ECG events, respectively, could be found. Previously, we'd return ``None`` in these situations, which does not match the documented behavior of returning a list of projectors (:gh:`9277` by `Richard Höchenberger`_)
 
 API changes

doc/conf.py

@@ -793,6 +793,7 @@ def reset_warnings(gallery_conf, fname):
         'ignore', message='.*mne-realtime.*', category=DeprecationWarning)
 
 
+
 reset_warnings(None, None)
 
 

mne/preprocessing/_peak_finder.py

@@ -50,11 +50,26 @@ def peak_finder(x0, thresh=None, extrema=1, verbose=None):
     s = x0.size
 
     if x0.ndim >= 2 or s == 0:
-        raise ValueError('The input data must be a non empty 1D vector')
+        raise ValueError('The input data must be a non-empty 1D vector')
+
+    if x0.dtype.kind in 'uib':  # integer dtype
+        x0 = x0.astype(float)
+
+    eps = 1e2 * np.finfo(x0.dtype).eps
+
+    if thresh is not None:
+        assert thresh >= 0.
 
     if thresh is None:
         thresh = (np.max(x0) - np.min(x0)) / 4
+        if thresh < eps:
+            thresh = eps
         logger.debug('Peak finder automatic threshold: %0.2g' % (thresh,))
+    elif thresh < eps:
+        raise ValueError(
+            f'The selected threshold value is too small for peak detection '
+            f'due to inevitable floating-point imprecisions.\nPlease pick a '
+            f'value >= {eps:g} (You passed thresh={thresh})')
 
     assert extrema in [-1, 1]
 
@@ -160,19 +175,18 @@ def peak_finder(x0, thresh=None, extrema=1, verbose=None):
             peak_mags = []
             peak_inds = []
 
-    # Change sign of data if was finding minima
-    if extrema < 0:
-        peak_mags *= -1.0
-
     # ensure output type array
     if not isinstance(peak_inds, np.ndarray):
         peak_inds = np.atleast_1d(peak_inds).astype('int64')
 
     if not isinstance(peak_mags, np.ndarray):
         peak_mags = np.atleast_1d(peak_mags).astype('float64')
 
-    # Plot if no output desired
-    if len(peak_inds) == 0:
+    # Change sign of data if was finding minima
+    if extrema < 0 and peak_mags.size > 0:
+        peak_mags *= -1.0
+
+    if peak_inds.size == 0:
         logger.info('No significant peaks found')
     else:
         logger.info('Found %d significant peak%s'

mne/preprocessing/ecg.py

@@ -31,7 +31,7 @@ def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3,
         Sampling rate
     ecg : array
         ECG signal
-    thresh_value : float | str
+    thresh_value : float | 'auto'
         qrs detection threshold. Can also be "auto" for automatic
         selection of threshold.
     levels : float
@@ -65,18 +65,32 @@ def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3,
 
     n_points = len(ecg_abs)
 
+    # get peak amplitudes in the 1st, 2nd, and 3rd second of the recording
     maxpt = np.empty(3)
     maxpt[0] = np.max(ecg_abs[:init])
     maxpt[1] = np.max(ecg_abs[init:init * 2])
     maxpt[2] = np.max(ecg_abs[init * 2:init * 3])
 
     init_max = np.mean(maxpt)
+    eps = 1e2 * np.finfo(ecg_abs.dtype).eps
+
+    if init_max < eps:
+        warn('The first 3 seconds of your ECG data appear to be flat. Please '
+             'inspect your data to ensure this is intentional.')
+        init_max = eps
 
     if thresh_value == 'auto':
         thresh_runs = np.arange(0.3, 1.1, 0.05)
     elif isinstance(thresh_value, str):
         raise ValueError('threshold value must be "auto" or a float')
     else:
+        if init_max * thresh_value < eps:
+            raise ValueError(
+                f'The selected threshold value is too small for peak '
+                f'detection with the provided data, due to inevitable '
+                f'floating-point imprecisions.\n'
+                f'Please pick a value >= {eps/init_max:g} '
+                f'(You passed thresh={thresh_value})')
         thresh_runs = [thresh_value]
 
     # Try a few thresholds (or just one)
@@ -100,7 +114,7 @@ def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3,
             else:
                 ii += 1
 
-        if len(rms) == 0:
+        if not rms:
             rms.append(0.0)
             time.append(0.0)
         time = np.array(time)
@@ -151,7 +165,7 @@ def find_ecg_events(raw, event_id=999, ch_name=None, tstart=0.0,
     %(ecg_ch_name)s
     %(ecg_tstart)s
     %(ecg_filter_freqs)s
-    qrs_threshold : float | str
+    qrs_threshold : float | 'auto'
         Between 0 and 1. qrs detection threshold. Can also be "auto" to
         automatically choose the threshold that generates a reasonable
         number of heartbeats (40-160 beats / min).

mne/preprocessing/ica.py

@@ -1246,7 +1246,7 @@ def find_bads_ecg(self, inst, ch_name=None, threshold='auto', start=None,
             The name of the channel to use for ECG peak detection.
             The argument is mandatory if the dataset contains no ECG
             channels.
-        threshold : float | str
+        threshold : float | 'auto'
             The value above which a feature is classified as outlier. If 'auto'
             and method is 'ctps', automatically compute the threshold. If
             'auto' and method is 'correlation', defaults to 3.0. The default

mne/preprocessing/tests/test_ecg.py

@@ -86,12 +86,20 @@ def test_find_ecg():
     assert 'ECG-SYN' not in raw.ch_names
     assert 'ECG-SYN' not in ecg_epochs.ch_names
 
-    # Test behavior if no peaks can be found -> achieve this by providing
-    # all-zero'd data
-    raw._data[ecg_idx] = 0.
+    # test behavior if no peaks can be found -> achieve this by setting a high
+    # threshold
     ecg_events, _, average_pulse, ecg = find_ecg_events(
-        raw, ch_name=raw.ch_names[ecg_idx], return_ecg=True
+        raw, ch_name=raw.ch_names[ecg_idx], return_ecg=True, qrs_threshold=1e5
     )
     assert ecg_events.size == 0
     assert average_pulse == 0
+
+    # test behavior for an entirely flat channel
+    raw._data[ecg_idx] = 0.
+    with pytest.warns(RuntimeWarning, match='first 3 seconds.*flat'):
+        ecg_events, _, average_pulse, ecg = find_ecg_events(
+            raw, ch_name=raw.ch_names[ecg_idx], return_ecg=True
+        )
+    assert ecg_events.size == 0
+    assert average_pulse == 0
     assert np.allclose(ecg, np.zeros_like(ecg))

mne/preprocessing/tests/test_eog.py

@@ -1,5 +1,6 @@
 import os.path as op
 
+import numpy as np
 import pytest
 
 from mne import Annotations
@@ -39,3 +40,14 @@ def test_find_eog():
 
     events = find_eog_events(raw, ch_name=['EEG 060', 'EOG 061'])
     assert len(events) == 4
+
+    # test too low threshold
+    eps = np.finfo(raw.get_data().dtype).eps
+    with pytest.raises(ValueError, match='threshold value is too small'):
+        find_eog_events(raw, ch_name=None, thresh=eps)
+
+    # test with an entirely flat signal
+    eog = raw.copy().pick_types(meg=False, eog=True).load_data()
+    eog._data[:] = 0.
+    events = find_eog_events(eog, ch_name=None)
+    assert events.size == 0

mne/preprocessing/tests/test_ssp.py

@@ -63,11 +63,7 @@ def test_compute_proj_ecg(short_raw, average):
         projs, events, drop_log = compute_proj_ecg(
             raw, n_mag=2, n_grad=2, n_eeg=2, ch_name='MEG 1531', bads=[],
             average=average, avg_ref=True, no_proj=True, l_freq=None,
-            h_freq=None, tmax=dur_use, return_drop_log=True,
-            # XXX can be removed once
-            # XXX https://github.com/mne-tools/mne-python/issues/9273
-            # XXX has been resolved:
-            qrs_threshold=1e-15)
+            h_freq=None, tmax=dur_use, return_drop_log=True)
     assert projs == []
     assert len(events) == len(drop_log)