Extend artifact detection: DetectAndRemoveArtifacts preprocessor and signed saturation (#4539)

Co-authored-by: Olivier Winter <olivier.winter@hotmail.fr>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: 3 people

doc/modules/preprocessing.rst

@@ -149,6 +149,37 @@ can also be obtained from the pipeline object directly:
     dict_used_to_make_pipeline = preprocessing_pipeline.preprocessor_dict
 
 
+Some preprocessing steps, such as :code:`detect_and_remove_artifacts`, allow you to specify an input recording
+and optionally another recording to perform some computation (e.g., detect artifacts on the output of a previous
+preprocessor, but remove them on the the parent preprocessor). In this case, the string "pipeline[preprocessor_name]"
+can be used in the dictionary to specify that the recording argument for this step should be the output of a previous
+preprocessor in the same pipeline. For example, if we want to use the output of the "bandpass_filter" step as the
+recording to detect artifacts, we can specify it as follows:
+
+.. code-block:: python
+
+    preprocessing_dict = {
+        'bandpass_filter': {'freq_min': 250},
+        'common_reference': {'operator': 'median', 'reference': 'global'},
+        'detect_and_remove_artifacts': {'recording_to_detect': 'pipeline[bandpass_filter]'},
+    }
+
+This will detect artifacts on the output of the "bandpass_filter" step, but the artifacts will be removed on the output
+of the "common_reference" step (since the parent recording for "detect_and_remove_artifacts" is by default the output of
+the previous step in the pipeline, which is "common_reference" in this case).
+To specify the "raw" recording, i.e., the input to the pipeline, we can use "pipeline[raw]".
+For example, if we want to detect artifacts on the raw recording, we can specify it as follows:
+
+
+.. code-block:: python
+
+    preprocessing_dict = {
+        'bandpass_filter': {'freq_min': 250},
+        'common_reference': {'operator': 'median', 'reference': 'global'},
+        'detect_and_remove_artifacts': {'recording_to_detect': 'pipeline[raw]'},
+    }
+
+
 Impact on recording dtype
 -------------------------
 

src/spikeinterface/extractors/neoextractors/openephys.py

@@ -15,7 +15,10 @@
 
 import probeinterface
 
-from spikeinterface.extractors.neuropixels_utils import get_neuropixels_sample_shifts_from_probe
+from spikeinterface.extractors.neuropixels_utils import (
+    get_neuropixels_sample_shifts_from_probe,
+    compute_saturation_threshold_from_probe,
+)
 from spikeinterface.extractors.neoextractors.neobaseextractor import NeoBaseRecordingExtractor, NeoBaseEventExtractor
 
 
@@ -338,6 +341,11 @@ def __init__(
                     if sample_shifts is not None:
                         self.set_property("inter_sample_shift", sample_shifts)
 
+                    # add saturation levels if available
+                    saturation_threshold_uV = compute_saturation_threshold_from_probe(probe, oe_stream_name)
+                    if saturation_threshold_uV is not None:
+                        self.annotate(saturation_threshold_uV=saturation_threshold_uV)
+
             # folder_path can point to different levels of the OE folder structure
             # (root, record node, experiment, or recording). We need to find the root folder
             # in order to load the sync timestamps and set them as times to the recording.

src/spikeinterface/extractors/neoextractors/spikegadgets.py

@@ -56,9 +56,8 @@ def __init__(
         )
         self._kwargs.update(dict(file_path=str(Path(file_path).absolute()), stream_id=stream_id))
 
-        probegroup = None  # TODO remove once probeinterface is updated to 0.2.22 in the pyproject.toml
-        if packaging.version.parse(probeinterface.__version__) > packaging.version.parse("0.2.21"):
-            probegroup = probeinterface.read_spikegadgets(file_path, raise_error=False)
+        probegroup = probeinterface.read_spikegadgets(file_path, raise_error=False)
+        # TODO: add adc sample shifts and saturation levels if available in the probe metadata
 
         if probegroup is not None:
             self.set_probegroup(probegroup, in_place=True)

src/spikeinterface/extractors/neoextractors/spikeglx.py

@@ -4,7 +4,10 @@
 import probeinterface
 
 from spikeinterface.core.core_tools import define_function_from_class
-from spikeinterface.extractors.neuropixels_utils import get_neuropixels_sample_shifts_from_probe
+from spikeinterface.extractors.neuropixels_utils import (
+    get_neuropixels_sample_shifts_from_probe,
+    compute_saturation_threshold_from_probe,
+)
 from spikeinterface.extractors.neoextractors.neobaseextractor import NeoBaseRecordingExtractor, NeoBaseEventExtractor
 
 
@@ -91,6 +94,11 @@ def __init__(
             sample_shifts = get_neuropixels_sample_shifts_from_probe(probe)
             if sample_shifts is not None:
                 self.set_property("inter_sample_shift", sample_shifts)
+
+            # add saturation levels if available
+            saturation_threshold_uV = compute_saturation_threshold_from_probe(probe, self.stream_id)
+            if saturation_threshold_uV is not None:
+                self.annotate(saturation_threshold_uV=saturation_threshold_uV)
         else:
             warning_message = (
                 "Unable to find a corresponding metadata file for the recording. "

src/spikeinterface/extractors/neuropixels_utils.py

@@ -53,6 +53,43 @@ def get_neuropixels_sample_shifts_from_probe(probe: Probe) -> np.ndarray:
     return sample_shifts
 
 
+def compute_saturation_threshold_from_probe(probe: Probe, stream_name: str) -> float:
+    """
+    Compute the saturation threshold in microvolts for a Neuropixels probe based on the probe information.
+
+    Parameters
+    ----------
+    probe : Probe
+        The probe object containing channel and ADC information.
+    stream_name : str
+        The name of the stream. If it contains "lf" (or "LF"), the function will look for the
+        lf gain and saturation in the probe annotations.
+
+    Returns
+    -------
+    saturation_threshold_uV : float
+        The saturation threshold in microvolts, or None if it cannot be computed.
+    """
+    adc_range_vpp = probe.annotations.get("adc_range_vpp", None)
+    gain = None
+    if "lf" in stream_name.lower():
+        gain = probe.annotations.get("lf_gain", None)
+    else:
+        gain = probe.annotations.get("ap_gain", None)
+
+    if adc_range_vpp is not None and gain is not None:
+        saturation_threshold_uV = (adc_range_vpp / 2) / gain * 1e6
+        return saturation_threshold_uV
+
+    warnings.warn(
+        "Unable to compute saturation threshold from the Neuropixels probe metadata. "
+        "The saturation threshold will not be loaded. ",
+        UserWarning,
+        stacklevel=2,
+    )
+    return None
+
+
 def synchronize_neuropixel_streams(recording_ref, recording_other):
     """
     Use the last "sync" channel from spikeglx or openephys neuropixels to synchronize