Fix pixel breath negative TIV

eitprocessing/features/pixel_breath.py

@@ -1,8 +1,11 @@
 import itertools
+import warnings
 from collections.abc import Callable
-from dataclasses import dataclass, field
+from dataclasses import InitVar, dataclass, field
+from typing import Final
 
 import numpy as np
+from scipy import signal
 
 from eitprocessing.datahandling.breath import Breath
 from eitprocessing.datahandling.continuousdata import ContinuousData
@@ -11,6 +14,15 @@
 from eitprocessing.datahandling.sequence import Sequence
 from eitprocessing.features.breath_detection import BreathDetection
 
+_SENTINAL_BREATH_DETECTION: Final = BreathDetection()
+MAX_XCORR_LAG = 0.75
+
+
+def _return_sentinal_breath_detection() -> BreathDetection:
+    # Returns a sential of a BreathDetection, which only exists to signal that the default value for breath_detection
+    # was used.
+    return _SENTINAL_BREATH_DETECTION
+
 
 @dataclass
 class PixelBreath:
@@ -30,20 +42,31 @@ class PixelBreath:
     ```
 
     Args:
-    breath_detection_kwargs (dict): A dictionary of keyword arguments for breath detection.
-        The available keyword arguments are:
-        minimum_duration: minimum expected duration of breaths, defaults to 2/3 of a second
-        averaging_window_duration: duration of window used for averaging the data, defaults to 15 seconds
-        averaging_window_function: function used to create a window for averaging the data, defaults to np.blackman
-        amplitude_cutoff_fraction: fraction of the median amplitude below which breaths are removed, defaults to 0.25
-        invalid_data_removal_window_length: window around invalid data in which breaths are removed, defaults to 0.5
-        invalid_data_removal_percentile: the nth percentile of values used to remove outliers, defaults to 5
-        invalid_data_removal_multiplier: the multiplier used to remove outliers, defaults to 4
+        breath_detection (BreathDetection): BreathDetection object to use for detecing breaths.
+        allow_negative_amplitude (bool): whether to asume out-of-phase pixels have negative amplitude instead.
     """
 
-    breath_detection_kwargs: dict = field(default_factory=dict)
-
-    def find_pixel_breaths(
+    breath_detection: BreathDetection = field(default_factory=_return_sentinal_breath_detection)
+    breath_detection_kwargs: InitVar[dict | None] = None
+    allow_negative_amplitude: bool = True
+    correct_for_phase_shift: bool = True
+
+    def __post_init__(self, breath_detection_kwargs: dict | None):
+        if breath_detection_kwargs is not None:
+            if self.breath_detection is not _SENTINAL_BREATH_DETECTION:
+                msg = (
+                    "`breath_detection_kwargs` is deprecated, and can't be used at the same time as `breath_detection`."
+                )
+                raise TypeError(msg)
+
+            self.breath_detection = BreathDetection(**breath_detection_kwargs)
+            warnings.warn(
+                "`breath_detection_kwargs` is deprecated and will be removed soon. "
+                "Replace with `breath_detection=BreathDetection(**breath_detection_kwargs)`.",
+                DeprecationWarning,
+            )
+
+    def find_pixel_breaths(  # noqa: C901, PLR0912, PLR0915
         self,
         eit_data: EITData,
         continuous_data: ContinuousData,
@@ -103,9 +126,7 @@ def find_pixel_breaths(
             msg = "To store the result a Sequence dataclass must be provided."
             raise ValueError(msg)
 
-        bd_kwargs = self.breath_detection_kwargs.copy()
-        breath_detection = BreathDetection(**bd_kwargs)
-        continuous_breaths = breath_detection.find_breaths(continuous_data)
+        continuous_breaths = self.breath_detection.find_breaths(continuous_data)
 
         indices_breath_middles = np.searchsorted(
             eit_data.time,
@@ -144,28 +165,57 @@ def find_pixel_breaths(
 
         pixel_breaths = np.full((len(continuous_breaths), n_rows, n_cols), None)
 
+        lags = signal.correlation_lags(len(continuous_data), len(continuous_data), mode="same")
+
         for row, col in itertools.product(range(n_rows), range(n_cols)):
             mean_tiv = mean_tiv_pixel[row, col]
 
-            if np.any(pixel_impedance[:, row, col] > 0):
-                if mean_tiv < 0:
-                    start_func, middle_func = np.argmax, np.argmin
-                else:
-                    start_func, middle_func = np.argmin, np.argmax
-
-                outsides = self._find_extreme_indices(pixel_impedance, indices_breath_middles, row, col, start_func)
-                starts = outsides[:-1]
-                ends = outsides[1:]
-                middles = self._find_extreme_indices(pixel_impedance, outsides, row, col, middle_func)
-                # TODO discuss; this block of code is implemented to prevent noisy pixels from breaking the code.
-                # Quick solve is to make entire breath object None if any breath in a pixel does not have
-                # consecutive start, middle and end.
-                # However, this might cause problems elsewhere.
-                if (starts >= middles).any() or (middles >= ends).any():
-                    pixel_breath = None
+            if np.std(pixel_impedance[:, row, col]) == 0:
+                # pixel has no amplitude
+                continue
+
+            if self.allow_negative_amplitude and mean_tiv < 0:
+                start_func, middle_func = np.argmax, np.argmin
+                lagged_indices_breath_middles = indices_breath_middles
+            else:
+                start_func, middle_func = np.argmin, np.argmax
+
+                cd = np.copy(continuous_data.values)
+                cd -= np.nanmean(cd)
+                pi = np.copy(pixel_impedance[:, row, col])
+                pi -= np.nanmean(pixel_impedance[:, row, col])
+
+                if self.correct_for_phase_shift:
+                    # search for maximum cross correlation within MAX_XCORR_LAG times the average
+                    # duration of a breath
+                    xcorr = signal.correlate(cd, pi, mode="same")
+                    max_lag = MAX_XCORR_LAG * np.mean(np.diff(indices_breath_middles))
+                    lag_range = (lags > -max_lag) & (lags < max_lag)
+                    # TODO: if this does not work, implement robust peak detection
+                    lag = lags[lag_range][np.argmax(xcorr[lag_range])]
+                    # positive lag: pixel inflates later than summed
+
+                    # shift search area
+                    lagged_indices_breath_middles = indices_breath_middles - lag
+                    lagged_indices_breath_middles = lagged_indices_breath_middles[
+                        (lagged_indices_breath_middles >= 0) & (lagged_indices_breath_middles < len(cd))
+                    ]
                 else:
-                    pixel_breath = self._construct_breaths(starts, middles, ends, time)
-                pixel_breaths[:, row, col] = pixel_breath
+                    lagged_indices_breath_middles = indices_breath_middles
+
+            outsides = self._find_extreme_indices(pixel_impedance, lagged_indices_breath_middles, row, col, start_func)
+            starts = outsides[:-1]
+            ends = outsides[1:]
+            middles = self._find_extreme_indices(pixel_impedance, outsides, row, col, middle_func)
+            # TODO discuss; this block of code is implemented to prevent noisy pixels from breaking the code.
+            # Quick solve is to make entire breath object None if any breath in a pixel does not have
+            # consecutive start, middle and end.
+            # However, this might cause problems elsewhere.
+            if (starts >= middles).any() or (middles >= ends).any():
+                pixel_breath = None
+            else:
+                pixel_breath = self._construct_breaths(starts, middles, ends, time)
+            pixel_breaths[:, row, col] = pixel_breath
 
         intervals = [(breath.start_time, breath.end_time) for breath in continuous_breaths.values]
 
@@ -178,15 +228,14 @@ def find_pixel_breaths(
             values=list(
                 pixel_breaths,
             ),  ## TODO: change back to pixel_breaths array when IntervalData works with 3D array
-            parameters=self.breath_detection_kwargs,
             derived_from=[eit_data],
         )
         if store:
             sequence.interval_data.add(pixel_breaths_container)
 
         return pixel_breaths_container
 
-    def _construct_breaths(self, start: list, middle: list, end: list, time: np.ndarray) -> list:
+    def _construct_breaths(self, start: list[int], middle: list[int], end: list[int], time: np.ndarray) -> list:
         breaths = [Breath(time[s], time[m], time[e]) for s, m, e in zip(start, middle, end, strict=True)]
         # First and last breath are not detected by definition (need two breaths to find one breath)
         return [None, *breaths, None]
@@ -222,5 +271,5 @@ def _find_extreme_indices(
             are located for each time segment.
         """
         return np.array(
-            [function(pixel_impedance[times[i] : times[i + 1], row, col]) + times[i] for i in range(len(times) - 1)],
+            [function(pixel_impedance[t1:t2, row, col]) + t1 for t1, t2 in itertools.pairwise(times)],
         )

eitprocessing/parameters/eeli.py

@@ -1,5 +1,6 @@
-from dataclasses import dataclass, field
-from typing import Literal, get_args
+import warnings
+from dataclasses import InitVar, dataclass, field
+from typing import Final, Literal, get_args
 
 import numpy as np
 
@@ -10,15 +11,38 @@
 from eitprocessing.features.breath_detection import BreathDetection
 from eitprocessing.parameters import ParameterCalculation
 
+_SENTINAL_BREATH_DETECTION: Final = BreathDetection()
+
+
+def _return_sentinal_breath_detection() -> BreathDetection:
+    # Returns a sential of a BreathDetection, which only exists to signal that the default value for breath_detection
+    # was used.
+    return _SENTINAL_BREATH_DETECTION
+
 
 @dataclass
 class EELI(ParameterCalculation):
     """Compute the end-expiratory lung impedance (EELI) per breath."""
 
+    breath_detection: BreathDetection = field(default_factory=_return_sentinal_breath_detection)
     method: Literal["breath_detection"] = "breath_detection"
-    breath_detection_kwargs: dict = field(default_factory=dict)
+    breath_detection_kwargs: InitVar[dict | None] = None
+
+    def __post_init__(self, breath_detection_kwargs: dict | None):
+        if breath_detection_kwargs is not None:
+            if self.breath_detection is not _SENTINAL_BREATH_DETECTION:
+                msg = (
+                    "`breath_detection_kwargs` is deprecated, and can't be used at the same time as `breath_detection`."
+                )
+                raise TypeError(msg)
+
+            self.breath_detection = BreathDetection(**breath_detection_kwargs)
+            warnings.warn(
+                "`breath_detection_kwargs` is deprecated and will be removed soon. "
+                "Replace with `breath_detection=BreathDetection(**breath_detection_kwargs)`.",
+                DeprecationWarning,
+            )
 
-    def __post_init__(self):
         _methods = get_args(EELI.__dataclass_fields__["method"].type)
         if self.method not in _methods:
             msg = f"Method {self.method} is not valid. Use any of {', '.join(_methods)}"
@@ -66,9 +90,7 @@ def compute_parameter(
 
         check_category(continuous_data, "impedance", raise_=True)
 
-        bd_kwargs = self.breath_detection_kwargs.copy()
-        breath_detection = BreathDetection(**bd_kwargs)
-        breaths = breath_detection.find_breaths(continuous_data)
+        breaths = self.breath_detection.find_breaths(continuous_data)
 
         if not len(breaths):
             time = np.array([], dtype=float)

eitprocessing/parameters/tidal_impedance_variation.py

@@ -1,7 +1,9 @@
 import itertools
-from dataclasses import dataclass, field
+import sys
+import warnings
+from dataclasses import InitVar, dataclass, field
 from functools import singledispatchmethod
-from typing import Literal, NoReturn
+from typing import Final, Literal, NoReturn
 
 import numpy as np
 
@@ -15,19 +17,57 @@
 from eitprocessing.features.pixel_breath import PixelBreath
 from eitprocessing.parameters import ParameterCalculation
 
+_SENTINAL_PIXEL_BREATH: Final = PixelBreath()
+_SENTINAL_BREATH_DETECTION: Final = BreathDetection()
+
+
+def _return_sentinal_pixel_breath() -> PixelBreath:
+    # Returns a sential of a PixelBreath, which only exists to signal that the default value for pixel_breath was used.
+    return _SENTINAL_PIXEL_BREATH
+
+
+def _return_sentinal_breath_detection() -> BreathDetection:
+    # Returns a sential of a BreathDetection, which only exists to signal that the default value for breath_detection
+    # was used.
+    return _SENTINAL_BREATH_DETECTION
+
 
 @dataclass
 class TIV(ParameterCalculation):
     """Compute the tidal impedance variation (TIV) per breath."""
 
     method: Literal["extremes"] = "extremes"
-    breath_detection_kwargs: dict = field(default_factory=dict)
+    breath_detection: BreathDetection = field(default_factory=_return_sentinal_breath_detection)
+    breath_detection_kwargs: InitVar[dict | None] = None
 
-    def __post_init__(self) -> None:
+    # The default is a sentinal that will be replaced in __post_init__
+    pixel_breath: PixelBreath = field(default_factory=_return_sentinal_pixel_breath)
+
+    def __post_init__(self, breath_detection_kwargs: dict | None) -> None:
         if self.method != "extremes":
             msg = f"Method {self.method} is not implemented. The method must be 'extremes'."
             raise NotImplementedError(msg)
 
+        if breath_detection_kwargs is not None:
+            if self.breath_detection is not _SENTINAL_BREATH_DETECTION:
+                msg = (
+                    "`breath_detection_kwargs` is deprecated, and can't be used at the same time as `breath_detection`."
+                )
+                raise TypeError(msg)
+
+            self.breath_detection = BreathDetection(**breath_detection_kwargs)
+            warnings.warn(
+                "`breath_detection_kwargs` is deprecated and will be removed soon. "
+                "Replace with `breath_detection=BreathDetection(**breath_detection_kwargs)`.",
+                DeprecationWarning,
+            )
+
+        if self.pixel_breath is _SENTINAL_PIXEL_BREATH:
+            # If no value was provided at initialization, PixelBreath should use the same BreathDetection object as TIV.
+            # However, a default factory cannot be used because it can't access self.breath_detection. The sentinal
+            # object is replaced here (only if pixel_breath was not provided) with the correct BreathDetection object.
+            self.pixel_breath = PixelBreath(breath_detection=self.breath_detection)
+
     @singledispatchmethod
     def compute_parameter(
         self,
@@ -101,7 +141,6 @@ def compute_continuous_parameter(
             category="impedance difference",
             time=[breath.middle_time for breath in breaths.values if breath is not None],
             description="Tidal impedance variation determined on continuous data",
-            parameters=self.breath_detection_kwargs,
             derived_from=[continuous_data],
             values=tiv_values,
         )
@@ -213,7 +252,6 @@ def compute_pixel_parameter(
             category="impedance difference",
             time=list(all_pixels_breath_timings),
             description="Tidal impedance variation determined on pixel impedance",
-            parameters=self.breath_detection_kwargs,
             derived_from=[eit_data],
             values=list(all_pixels_tiv_values.astype(float)),
         )
@@ -224,9 +262,7 @@ def compute_pixel_parameter(
         return tiv_container
 
     def _detect_breaths(self, data: ContinuousData) -> IntervalData:
-        bd_kwargs = self.breath_detection_kwargs.copy()
-        breath_detection = BreathDetection(**bd_kwargs)
-        return breath_detection.find_breaths(data)
+        return self.breath_detection.find_breaths(data)
 
     def _detect_pixel_breaths(
         self,
@@ -235,9 +271,7 @@ def _detect_pixel_breaths(
         sequence: Sequence,
         store: bool,
     ) -> IntervalData:
-        bd_kwargs = self.breath_detection_kwargs.copy()
-        pi = PixelBreath(breath_detection_kwargs=bd_kwargs)
-        return pi.find_pixel_breaths(
+        return self.pixel_breath.find_pixel_breaths(
             eit_data,
             continuous_data,
             result_label="pixel breaths",
@@ -274,6 +308,13 @@ def _calculate_tiv_values(
             mean_outer_values = data[[start_indices, end_indices]].mean(axis=0)
             end_inspiratory_values = data[middle_indices]
             tiv_values = end_inspiratory_values - mean_outer_values
+        else:
+            msg = f"`tiv_method` ({tiv_method}) not valid."
+            exc = ValueError(msg)
+            if sys.version_info >= (3, 11):
+                exc.add_note("Valid value for `tiv_method` are 'inspiratory', 'expiratory' and 'mean'.")
+            raise exc
+
         if tiv_timing == "pixel":
             tiv_values = [None, *tiv_values, None]
 

tests/test_eeli.py

@@ -93,6 +93,16 @@ def test_eeli_values(repeat_n: int):  # noqa: ARG001
     assert np.array_equal(eeli_values, valley_values[1:])
 
 
+def test_bd_init():
+    assert EELI(breath_detection_kwargs={"minimum_duration": 0}) == EELI(
+        breath_detection=BreathDetection(minimum_duration=0)
+    )
+    with pytest.warns(DeprecationWarning):
+        EELI(breath_detection_kwargs={"minimum_duration": 0})
+    with pytest.raises(TypeError):
+        EELI(breath_detection_kwargs={"minimum_duration": 0}, breath_detection=BreathDetection(minimum_duration=0))
+
+
 def test_with_data(draeger1: Sequence, pytestconfig: pytest.Config):
     if pytestconfig.getoption("--cov"):
         pytest.skip("Skip with option '--cov' so other tests can cover 100%.")