BUG: use fixed normalized CTPS cutoff for automatic ECG detection

doc/changes/dev/14279.bugfix.rst

@@ -0,0 +1 @@
+Use a fixed normalized Kuiper-index cutoff of 0.3 for automatic CTPS ECG detection in :meth:`mne.preprocessing.ICA.find_bads_ecg`, independent of sampling frequency and trial count, by :newcontrib:`Amirali Moradniaei`.

doc/changes/names.inc

@@ -20,6 +20,7 @@
 .. _Alexandra Corneyllie: https://github.com/Lx37
 .. _Alexandre Barachant: https://alexandre.barachant.org
 .. _Aman Srivastava: https://github.com/aman-coder03
+.. _Amirali Moradniaei: https://github.com/AtomicGlance
 .. _Ana Radanovic: https://github.com/anaradanovic
 .. _Andrea Brovelli: https://brovelli.github.io/
 .. _Andreas Højlund: https://github.com/ahoejlund

mne/channels/_dig_montage_utils.py

@@ -19,7 +19,9 @@ def _read_dig_montage_egi(
             "hsp, hpi, elp, point_names, fif must all be None if egi is not None"
         )
     _check_fname(fname, overwrite="read", must_exist=True)
-    defusedxml = _soft_import("defusedxml", "reading EGI montages")
+    _soft_import("defusedxml", "reading EGI montages")
+    import defusedxml.ElementTree
+
     root = defusedxml.ElementTree.parse(fname).getroot()
     ns = root.tag[root.tag.index("{") : root.tag.index("}") + 1]
     sensors = root.find(f"{ns}sensorLayout/{ns}sensors")

mne/channels/_standard_montage_utils.py

@@ -423,7 +423,9 @@ def _read_brainvision(fname, head_size):
     # standard electrode positions: X-axis from T7 to T8, Y-axis from Oz to
     # Fpz, Z-axis orthogonal from XY-plane through Cz, fit to a sphere if
     # idealized (when radius=1), specified in millimeters
-    defusedxml = _soft_import("defusedxml", "reading BrainVision montages")
+    _soft_import("defusedxml", "reading BrainVision montages")
+    import defusedxml.ElementTree
+
     root = defusedxml.ElementTree.parse(fname).getroot()
     ch_names = [s.text for s in root.findall("./Electrode/Name")]
     theta = [float(s.text) for s in root.findall("./Electrode/Theta")]

mne/io/nedf/nedf.py

@@ -55,7 +55,9 @@ def _parse_nedf_header(header):
     n_samples : int
         The number of data samples.
     """
-    defusedxml = _soft_import("defusedxml", "reading NEDF data")
+    _soft_import("defusedxml", "reading NEDF data")
+    import defusedxml.ElementTree  # ty: ignore[unresolved-import]
+
     info = {}
     # nedf files have three accelerometer channels sampled at 100Hz followed
     # by five EEG samples + TTL trigger sampled at 500Hz

mne/preprocessing/ica.py

@@ -4,7 +4,6 @@
 # Copyright the MNE-Python contributors.
 
 import json
-import math
 import warnings
 from collections import namedtuple
 from collections.abc import Sequence
@@ -1617,30 +1616,6 @@ def _find_bads_ch(
 
         return labels, scores
 
-    def _get_ctps_threshold(self, pk_threshold=20):
-        """Automatically decide the threshold of Kuiper index for CTPS method.
-
-        This function finds the threshold of Kuiper index based on the
-        threshold of pk. Kuiper statistic that minimizes the difference between
-        pk and the pk threshold (defaults to 20 :footcite:`DammersEtAl2008`)
-        is returned. It is assumed that the data are appropriately filtered and
-        bad data are rejected at least based on peak-to-peak amplitude
-        when/before running the ICA decomposition on data.
-
-        References
-        ----------
-        .. footbibliography::
-        """
-        N = self.info["sfreq"]
-        Vs = np.arange(1, 100) / 100
-        C = math.sqrt(N) + 0.155 + 0.24 / math.sqrt(N)
-        # in formula (13), when k gets large, only k=1 matters for the
-        # summation. k*V*C thus becomes V*C
-        Pks = 2 * (4 * (Vs * C) ** 2 - 1) * (np.exp(-2 * (Vs * C) ** 2))
-        # NOTE: the threshold of pk is transformed to Pk for comparison
-        # pk = -log10(Pk)
-        return Vs[np.argmin(np.abs(Pks - 10 ** (-pk_threshold)))]
-
     @verbose
     def find_bads_ecg(
         self,
@@ -1714,9 +1689,9 @@ def find_bads_ecg(
         The ``threshold``, ``method``, and ``measure`` parameters interact in
         the following ways:
 
-        - If ``method='ctps'``, ``threshold`` refers to the significance value
-          of a Kuiper statistic, and ``threshold='auto'`` will compute the
-          threshold automatically based on the sampling frequency.
+        - If ``method='ctps'``, ``threshold`` refers to the maximum normalized
+          Kuiper index across time, and ``threshold='auto'`` sets the threshold
+          to 0.3, independent of the sampling frequency and number of trials.
         - If ``method='correlation'`` and ``measure='correlation'``,
           ``threshold`` refers to the Pearson correlation value, and
           ``threshold='auto'`` sets the threshold to 0.9.
@@ -1747,9 +1722,6 @@ def find_bads_ecg(
             ecg = inst.ch_names[idx_ecg]
 
         if method == "ctps":
-            if threshold == "auto":
-                threshold = self._get_ctps_threshold()
-                logger.info(f"Using threshold: {threshold:.2f} for CTPS ECG detection")
             if isinstance(inst, BaseRaw):
                 sources = self.get_sources(
                     create_ecg_epochs(
@@ -1771,6 +1743,9 @@ def find_bads_ecg(
                 sources = self.get_sources(inst).get_data(copy=False)
             else:
                 raise ValueError("With `ctps` only Raw and Epochs input is supported")
+            if threshold == "auto":
+                threshold = 0.3
+                logger.info(f"Using threshold: {threshold:.2f} for CTPS ECG detection")
             _, p_vals, _ = ctps(sources)
             scores = p_vals.max(-1)
             ecg_idx = np.where(scores >= threshold)[0]

mne/preprocessing/tests/test_ica.py

@@ -762,6 +762,24 @@ def short_raw_epochs():
     return raw, epochs, epochs_eog
 
 
+@pytest.mark.parametrize("n_epochs", [2, 3])
+@pytest.mark.parametrize("sfreq", [50, 100])
+def test_ctps_auto_threshold(short_raw_epochs, n_epochs, sfreq):
+    """Check auto uses the same normalized cutoff for different inputs."""
+    _, epochs, _ = short_raw_epochs
+    epochs = epochs[:n_epochs].copy().resample(sfreq)
+    ica = _ICA(n_components=2, rng=0)
+    with _baseline_corrected:
+        ica.fit(epochs)
+    with catch_logging(True) as log:
+        auto_idx, auto_scores = ica.find_bads_ecg(epochs, threshold="auto")
+    manual_idx, manual_scores = ica.find_bads_ecg(epochs, threshold=0.3)
+    assert "Using threshold: 0.30" in log.getvalue()
+    assert_array_equal(auto_scores, manual_scores)
+    assert_array_equal(auto_idx, manual_idx)
+    assert_array_equal(sorted(auto_idx), np.flatnonzero(auto_scores >= 0.3))
+
+
 @pytest.mark.slowtest
 @pytest.mark.parametrize("method", ["picard", "fastica"])
 def test_ica_additional(method, tmp_path, short_raw_epochs):
@@ -793,8 +811,6 @@ def test_ica_additional(method, tmp_path, short_raw_epochs):
         ica.fit(raw, np.arange(1, 6))
     _assert_ica_attributes(ica, raw.get_data(np.arange(1, 6)))
 
-    # check Kuiper index threshold
-    assert_allclose(ica._get_ctps_threshold(), 0.5)
     with pytest.raises(TypeError, match="str or numeric"):
         ica.find_bads_ecg(raw, threshold=None)
     with pytest.warns(RuntimeWarning, match="is longer than the signal"):
@@ -806,7 +822,11 @@ def test_ica_additional(method, tmp_path, short_raw_epochs):
         )
     # check passing a ch_name to find_bads_ecg
     with pytest.warns(RuntimeWarning, match="longer"):
-        _, scores_1 = ica.find_bads_ecg(raw, threshold="auto")
+        auto_idx, scores_1 = ica.find_bads_ecg(raw, threshold="auto")
+    with pytest.warns(RuntimeWarning, match="longer"):
+        manual_idx, manual_scores = ica.find_bads_ecg(raw, threshold=0.3)
+    assert_array_equal(scores_1, manual_scores)
+    assert_array_equal(auto_idx, manual_idx)
     with pytest.warns(RuntimeWarning, match="longer"):
         _, scores_2 = ica.find_bads_ecg(raw, raw.ch_names[1], threshold="auto")
     assert scores_1[0] != scores_2[0]

tutorials/preprocessing/40_artifact_correction_ica.py

@@ -444,8 +444,8 @@
 # necessary to pass a specific channel name.
 # `~mne.preprocessing.ICA.find_bads_ecg` also has two options for its
 # ``method`` parameter: ``'ctps'`` (cross-trial phase statistics
-# :footcite:`DammersEtAl2008`) and
-# ``'correlation'`` (Pearson correlation between data and ECG channel).
+# :footcite:`DammersEtAl2008`) and ``'correlation'`` (Pearson correlation
+# between data and ECG channel).
 
 ica.exclude = []
 # find which ICs match the ECG pattern