17495122789

Committed 05 Sep 2025 01:50PM UTC coverage: 84.776% (+0.02%) from 84.761%

Build # 17495122789

Build Type

push

github

Committed by

JulietteFrancovich

Commit Message

Bump version: 1.8.1 → 1.8.2

Run Details

747 of 960 branches covered (77.81%)

Branch coverage included in aggregate %.

1 of 1 new or added line in 1 file covered. (100.0%)

4 existing lines in 1 file now uncovered.

2739 of 3152 relevant lines covered (86.9%)

0.87 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

89.87

/eitprocessing/features/pixel_breath.py

import itertools
import warnings
from collections.abc import Callable
from dataclasses import InitVar, dataclass, field
from typing import Final, Literal

import numpy as np
from scipy import signal

from eitprocessing.datahandling.breath import Breath
from eitprocessing.datahandling.continuousdata import ContinuousData
from eitprocessing.datahandling.eitdata import EITData
from eitprocessing.datahandling.intervaldata import IntervalData
from eitprocessing.datahandling.sequence import Sequence
from eitprocessing.features.breath_detection import BreathDetection

_SENTINEL_BREATH_DETECTION: Final = BreathDetection()
MAX_XCORR_LAG = 1
ALLOW_FRACTION_BREATHS_SKIPPED = 0.25


def _return_sentinel_breath_detection() -> BreathDetection:
    # Returns a sentinel of a BreathDetection, which only exists to signal that the default value for breath_detection
    # was used.
    return _SENTINEL_BREATH_DETECTION


@dataclass(kw_only=True)
class PixelBreath:
    """Algorithm for detecting timing of pixel breaths in pixel impedance data.

    This algorithm detects the position of start of inspiration, end of inspiration and
    end of expiration in pixel impedance data. It uses BreathDetection to find the global start and end
    of inspiration and expiration. These points are then used to find the start/end of pixel
    inspiration/expiration in pixel impedance data.

    Since this algorithm uses the previous and next global breath to determine the start and end of a pixel breath, the
    first and last global breaths can not be used to determine pixel breaths. They are always set to `None` in the
    return list. Breaths that could not properly be detected are set to `None` as well.

    Some pixel breaths may be phase shifted (inflation starts and ends later compared to others, e.g., due to pendelluft
    or late airway opening). Other pixel breaths may have a negative amplitude (impedance decreases during inspiration,
    e.g., due to pleural effusion or reconstruction artifacts). It is not always possible to determine whether a pixel
    is out of phase or has a negative amplitude. PixelBreath has three different phase correction modes. In 'negative
    amplitude' mode (default), pixels that have a decrease in amplitude between the start and end of globally defined
    inspiration, will have a negative amplitude and smaller phase shift. In 'phase shift' mode, all pixel breaths will
    have positive amplitudes, but can have large phase shifts. In 'none'/`None` mode, all pixels are assumed to be
    within rouglhy -90 to 90 degrees of phase. Note that the 'none' mode can lead to unexpected results, such as
    ultra-short (down to 2 frames) or very long breaths.

    Example:
    ```
    >>> pi = PixelBreath()
    >>> eit_data = sequence.eit_data['raw']
    >>> continuous_data = sequence.continuous_data['global_impedance_(raw)']
    >>> pixel_breaths = pi.find_pixel_breaths(eit_data, continuous_data, sequence)
    ```

    Args:
        breath_detection (BreathDetection): BreathDetection object to use for detecting breaths.
        phase_correction_mode: How to resolve pixels that are out-of-phase. Defaults to "negative amplitude".
    """

    breath_detection: BreathDetection = field(default_factory=_return_sentinel_breath_detection)
    breath_detection_kwargs: InitVar[dict | None] = None
    phase_correction_mode: Literal["negative amplitude", "phase shift", "none"] | None = "negative amplitude"

    def __post_init__(self, breath_detection_kwargs: dict | None):
        if breath_detection_kwargs is not None:
            if self.breath_detection is not _SENTINEL_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,
                stacklevel=2,
            )

    def find_pixel_breaths(  # noqa: C901, PLR0912, PLR0915
        self,
        eit_data: EITData,
        continuous_data: ContinuousData,
        sequence: Sequence | None = None,
        store: bool | None = None,
        result_label: str = "pixel_breaths",
    ) -> IntervalData:
        """Find pixel breaths in the data.

        This method finds the pixel start/end of inspiration/expiration based on the start/end of inspiration/expiration
        as detected in the continuous data.

        For most pixels, the start of a breath (start inspiration) is the valley between the middles (start of
        expiration) of the globally defined breaths on either side. The end of a pixel breath is the start of the next
        pixel breath. The middle of the pixel breath is the peak between the start and end of the pixel breath.

        If the pixel is out of phase or has negative amplitude, the definition of the breath depends on the phase
        correction mode. In 'negative amplitude' mode, the start of a breath is the peak between the middles of the
        globally defined breaths on either side, while the middle of the pixel breath is the valley of the start and end
        of the pixel breath. In 'phase shift' mode, first the phase shift between the pixel impedance and global
        impedance is determined as the highest crosscorrelation between the signals near a phase shift of 0. The start
        of breath is the valley between the phase shifted middles of the globally defined breaths on either side.

        Pixel breaths are constructed as a valley-peak-valley combination, representing the start of inspiration, the
        end of inspiration/start of expiration, and end of expiration.

        Args:
            eit_data: EITData to apply the algorithm to.
            continuous_data: ContinuousData to use for global breath detection.
            result_label: label of the returned IntervalData object, defaults to `'pixel_breaths'`.
            sequence: optional, Sequence that contains the object to detect pixel breaths in, and/or to store the result
            in.
            store: whether to store the result in the sequence, defaults to `True` if a Sequence if provided.

        Returns:
            An IntervalData object containing Breath objects.

        Raises:
            RuntimeError: If store is set to true but no sequence is provided.
            ValueError: If the provided sequence is not an instance of the Sequence dataclass.
        """
        if store is None and isinstance(sequence, Sequence):
            store = True

        if store and sequence is None:
            msg = "Can't store the result if no Sequence is provided."
            raise RuntimeError(msg)

        if store and not isinstance(sequence, Sequence):
            msg = "To store the result a Sequence dataclass must be provided."
            raise ValueError(msg)

        continuous_breaths = self.breath_detection.find_breaths(continuous_data)

        indices_breath_middles = np.searchsorted(
            eit_data.time,
            [breath.middle_time for breath in continuous_breaths.values],
        )

        _, n_rows, n_cols = eit_data.pixel_impedance.shape

        from eitprocessing.parameters.tidal_impedance_variation import TIV  # noqa: PLC0415

        pixel_tivs = TIV().compute_pixel_parameter(
            eit_data,
            continuous_data,
            sequence,
            tiv_method="inspiratory",
            tiv_timing="continuous",
            store=False,
        )  # Set store to false as to not save these pixel tivs as SparseData.

        pixel_tiv_with_continuous_data_timing = (
            np.empty((0, n_rows, n_cols)) if not len(pixel_tivs.values) else np.stack(pixel_tivs.values)
        )

        # Create a mask to detect slices that are entirely NaN
        all_nan_mask = np.isnan(pixel_tiv_with_continuous_data_timing).all(axis=0)

        # Initialize the mean_tiv_pixel array with NaNs
        mean_tiv_pixel = np.full((n_rows, n_cols), np.nan)

        # Only compute nanmean for slices that are not entirely NaN
        if not all_nan_mask.all():  # Check if there are any valid (non-all-NaN) slices
            mean_tiv_pixel[~all_nan_mask] = np.nanmean(pixel_tiv_with_continuous_data_timing[:, ~all_nan_mask], axis=0)

        time = eit_data.time
        pixel_impedance = eit_data.pixel_impedance

        pixel_breaths = np.full((len(continuous_breaths), n_rows, n_cols), None, dtype=object)

        lags = signal.correlation_lags(len(continuous_data), len(continuous_data), mode="same")

        match self.phase_correction_mode:
            case "negative amplitude":
                allow_negative_amplitude = True
                correct_for_phase_shift = None
            case "phase shift":
                allow_negative_amplitude = False
                correct_for_phase_shift = True
            case "none" | None:
                allow_negative_amplitude = False
                correct_for_phase_shift = False
            case _:
                msg = f"Unknown phase correction mode ({self.phase_correction_mode})."
                raise ValueError(msg)

        for row, col in itertools.product(range(n_rows), range(n_cols)):
            if all_nan_mask[row, col]:
                # pixel has no amplitude
                continue
            mean_tiv = mean_tiv_pixel[row, col]

            pi = np.copy(pixel_impedance[:, row, col])  # TODO: check whether copying is necessary
            if np.std(pi) == 0:
                # pixel has no amplitude
                continue

            if 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 = signal.detrend(cd, type="linear")

                if not np.all(np.isnan(pi)):
                    pi = signal.detrend(pi, type="linear")

                if correct_for_phase_shift:
                    # search for closest or 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)

                    # find the peaks in the cross correlation
                    peaks, _ = signal.find_peaks(xcorr[lag_range], height=0, prominence=xcorr.std())

                    # find the closest peak to zero lag
                    peak_lags = lags[lag_range][peaks]
                    peak_distances = np.abs(peak_lags)
                    min_peak_distance = np.min(peak_distances)
                    candidates = peak_lags[peak_distances == min_peak_distance]

                    # if there are multiple candidates, take the one with the highest cross correlation
                    if len(candidates) == 1:
                        lag = candidates[0]
                    elif len(candidates) == 2:  # noqa: PLR2004
                        # take the lag with the highest cross correlation
                        lag = candidates[np.argmax(xcorr[np.searchsorted(lags, candidates)])]
                    else:
                        msg = "Too many peaks found in cross correlation."
                        raise RuntimeError(msg)

                    # shift search area
                    lagged_indices_breath_middles = indices_breath_middles - lag
                else:
                    lagged_indices_breath_middles = indices_breath_middles

            skip = np.concatenate(
                (
                    np.flatnonzero(lagged_indices_breath_middles < 0),
                    np.flatnonzero(lagged_indices_breath_middles >= len(continuous_data)),
                )
            )
            lagged_indices_breath_middles = lagged_indices_breath_middles.clip(min=0, max=len(continuous_data) - 1)
            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)

            skip = np.concatenate((skip, np.flatnonzero(starts >= middles), np.flatnonzero(middles >= ends)))

            if len(starts) > len(skip) > len(starts) * ALLOW_FRACTION_BREATHS_SKIPPED:
                warnings.warn(
                    f"Skipping pixel ({row}, {col}) because more than half ({len(skip) / len(starts):.0%}) "
                    "of breaths skipped.",
                    RuntimeWarning,
                    stacklevel=2,
                )
                continue

            pixel_breaths[:, row, col] = self._construct_breaths(starts, middles, ends, time, skip=skip)

        intervals = [(breath.start_time, breath.end_time) for breath in continuous_breaths.values]

        pixel_breaths_container = IntervalData(
            label=result_label,
            name="Pixel in- and deflation timing as determined by Pixelbreath",
            unit=None,
            category="breath",
            intervals=intervals,
            values=list(
                pixel_breaths,
            ),  ## TODO: change back to pixel_breaths array when IntervalData works with 3D array
            derived_from=[eit_data],
        )
        if store:
            sequence.interval_data.add(pixel_breaths_container)

        return pixel_breaths_container

    def _construct_breaths(
        self, start: list[int], middle: list[int], end: list[int], time: np.ndarray, skip: np.ndarray | None = None
    ) -> list:
        skip_ = skip if skip is not None else np.array([], dtype=int)
        breaths = [
            Breath(time[s], time[m], time[e]) if i not in skip_ else None
            for i, (s, m, e) in enumerate(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]

    def _find_extreme_indices(
        self,
        pixel_impedance: np.ndarray,
        times: np.ndarray,
        row: int,
        col: int,
        function: Callable,
    ) -> np.ndarray:
        """Finds extreme indices in pixel impedance.

        This method divides the pixel impedance for a single pixel (selected using row and col) into smaller segments
        based on the `times` array. The times array consists of indices to divide these segments. The function iterates
        over each index in the times array to select consecutive segments of pixel impedance.
        For each segment, the method applies the `function` (either np.argmax or np.argmin) to extract an extreme value
        (local maximum or minimum).

        Args:
            pixel_impedance (np.ndarray): The pixel impedance array from which the function will extract values.
                Assumed to be 3-dimensional (e.g., time, rows, and columns).
            times (np.ndarray): 1D array of time indices. These times define the start
                and end of each segment in the pixel impedance.
            row (int): The row index in the pixel impedance
            col (int): The column index in the pixel impedance
            function (Callable): A function that is applied to each segment of data to find
                an extreme value (np.argmax or np.argmin)

        Returns:
            np.ndarray: An array of indices where the extreme values (based on the function)
            are located for each time segment.
        """
        return np.array(
            [function(pixel_impedance[t1:t2, row, col]) + t1 for t1, t2 in itertools.pairwise(times)],
        )

1	import itertools	1✔
2	import warnings	1✔
3	from collections.abc import Callable	1✔
4	from dataclasses import InitVar, dataclass, field	1✔
5	from typing import Final, Literal	1✔
6
7	import numpy as np	1✔
8	from scipy import signal	1✔
9
10	from eitprocessing.datahandling.breath import Breath	1✔
11	from eitprocessing.datahandling.continuousdata import ContinuousData	1✔
12	from eitprocessing.datahandling.eitdata import EITData	1✔
13	from eitprocessing.datahandling.intervaldata import IntervalData	1✔
14	from eitprocessing.datahandling.sequence import Sequence	1✔
15	from eitprocessing.features.breath_detection import BreathDetection	1✔
16
17	_SENTINEL_BREATH_DETECTION: Final = BreathDetection()	1✔
18	MAX_XCORR_LAG = 1	1✔
19	ALLOW_FRACTION_BREATHS_SKIPPED = 0.25	1✔
20
21
22	def _return_sentinel_breath_detection() -> BreathDetection:	1✔
23	# Returns a sentinel of a BreathDetection, which only exists to signal that the default value for breath_detection
24	# was used.
25	return _SENTINEL_BREATH_DETECTION	1✔
26
27
28	@dataclass(kw_only=True)	1✔
29	class PixelBreath:	1✔
30	"""Algorithm for detecting timing of pixel breaths in pixel impedance data.
31
32	This algorithm detects the position of start of inspiration, end of inspiration and
33	end of expiration in pixel impedance data. It uses BreathDetection to find the global start and end
34	of inspiration and expiration. These points are then used to find the start/end of pixel
35	inspiration/expiration in pixel impedance data.
36
37	Since this algorithm uses the previous and next global breath to determine the start and end of a pixel breath, the
38	first and last global breaths can not be used to determine pixel breaths. They are always set to `None` in the
39	return list. Breaths that could not properly be detected are set to `None` as well.
40
41	Some pixel breaths may be phase shifted (inflation starts and ends later compared to others, e.g., due to pendelluft
42	or late airway opening). Other pixel breaths may have a negative amplitude (impedance decreases during inspiration,
43	e.g., due to pleural effusion or reconstruction artifacts). It is not always possible to determine whether a pixel
44	is out of phase or has a negative amplitude. PixelBreath has three different phase correction modes. In 'negative
45	amplitude' mode (default), pixels that have a decrease in amplitude between the start and end of globally defined
46	inspiration, will have a negative amplitude and smaller phase shift. In 'phase shift' mode, all pixel breaths will
47	have positive amplitudes, but can have large phase shifts. In 'none'/`None` mode, all pixels are assumed to be
48	within rouglhy -90 to 90 degrees of phase. Note that the 'none' mode can lead to unexpected results, such as
49	ultra-short (down to 2 frames) or very long breaths.
50
51	Example:
52	```
53	>>> pi = PixelBreath()
54	>>> eit_data = sequence.eit_data['raw']
55	>>> continuous_data = sequence.continuous_data['global_impedance_(raw)']
56	>>> pixel_breaths = pi.find_pixel_breaths(eit_data, continuous_data, sequence)
57	```
58
59	Args:
60	breath_detection (BreathDetection): BreathDetection object to use for detecting breaths.
61	phase_correction_mode: How to resolve pixels that are out-of-phase. Defaults to "negative amplitude".
62	"""
63
64	breath_detection: BreathDetection = field(default_factory=_return_sentinel_breath_detection)	1✔
65	breath_detection_kwargs: InitVar[dict \| None] = None	1✔
66	phase_correction_mode: Literal["negative amplitude", "phase shift", "none"] \| None = "negative amplitude"	1✔
67
68	def __post_init__(self, breath_detection_kwargs: dict \| None):	1✔
69	if breath_detection_kwargs is not None:	1✔
70	if self.breath_detection is not _SENTINEL_BREATH_DETECTION:	1✔
71	msg = (	1✔
72	"`breath_detection_kwargs` is deprecated, and can't be used at the same time as `breath_detection`."
73	)
74	raise TypeError(msg)	1✔
75
76	self.breath_detection = BreathDetection(**breath_detection_kwargs)	1✔
77	warnings.warn(	1✔
78	"`breath_detection_kwargs` is deprecated and will be removed soon. "
79	"Replace with `breath_detection=BreathDetection(**breath_detection_kwargs)`.",
80	DeprecationWarning,
81	stacklevel=2,
82	)
83
84	def find_pixel_breaths( # noqa: C901, PLR0912, PLR0915	1✔
85	self,
86	eit_data: EITData,
87	continuous_data: ContinuousData,
88	sequence: Sequence \| None = None,
89	store: bool \| None = None,
90	result_label: str = "pixel_breaths",
91	) -> IntervalData:
92	"""Find pixel breaths in the data.
93
94	This method finds the pixel start/end of inspiration/expiration based on the start/end of inspiration/expiration
95	as detected in the continuous data.
96
97	For most pixels, the start of a breath (start inspiration) is the valley between the middles (start of
98	expiration) of the globally defined breaths on either side. The end of a pixel breath is the start of the next
99	pixel breath. The middle of the pixel breath is the peak between the start and end of the pixel breath.
100
101	If the pixel is out of phase or has negative amplitude, the definition of the breath depends on the phase
102	correction mode. In 'negative amplitude' mode, the start of a breath is the peak between the middles of the
103	globally defined breaths on either side, while the middle of the pixel breath is the valley of the start and end
104	of the pixel breath. In 'phase shift' mode, first the phase shift between the pixel impedance and global
105	impedance is determined as the highest crosscorrelation between the signals near a phase shift of 0. The start
106	of breath is the valley between the phase shifted middles of the globally defined breaths on either side.
107
108	Pixel breaths are constructed as a valley-peak-valley combination, representing the start of inspiration, the
109	end of inspiration/start of expiration, and end of expiration.
110
111	Args:
112	eit_data: EITData to apply the algorithm to.
113	continuous_data: ContinuousData to use for global breath detection.
114	result_label: label of the returned IntervalData object, defaults to `'pixel_breaths'`.
115	sequence: optional, Sequence that contains the object to detect pixel breaths in, and/or to store the result
116	in.
117	store: whether to store the result in the sequence, defaults to `True` if a Sequence if provided.
118
119	Returns:
120	An IntervalData object containing Breath objects.
121
122	Raises:
123	RuntimeError: If store is set to true but no sequence is provided.
124	ValueError: If the provided sequence is not an instance of the Sequence dataclass.
125	"""
126	if store is None and isinstance(sequence, Sequence):	1✔
127	store = True	1✔
128
129	if store and sequence is None:	1✔
130	msg = "Can't store the result if no Sequence is provided."	1✔
131	raise RuntimeError(msg)	1✔
132
133	if store and not isinstance(sequence, Sequence):	1✔
134	msg = "To store the result a Sequence dataclass must be provided."	1✔
135	raise ValueError(msg)	1✔
136
137	continuous_breaths = self.breath_detection.find_breaths(continuous_data)	1✔
138
139	indices_breath_middles = np.searchsorted(	1✔
140	eit_data.time,
141	[breath.middle_time for breath in continuous_breaths.values],
142	)
143
144	_, n_rows, n_cols = eit_data.pixel_impedance.shape	1✔
145
146	from eitprocessing.parameters.tidal_impedance_variation import TIV # noqa: PLC0415	1✔
147
148	pixel_tivs = TIV().compute_pixel_parameter(	1✔
149	eit_data,
150	continuous_data,
151	sequence,
152	tiv_method="inspiratory",
153	tiv_timing="continuous",
154	store=False,
155	) # Set store to false as to not save these pixel tivs as SparseData.
156
157	pixel_tiv_with_continuous_data_timing = (	1✔
158	np.empty((0, n_rows, n_cols)) if not len(pixel_tivs.values) else np.stack(pixel_tivs.values)
159	)
160
161	# Create a mask to detect slices that are entirely NaN
162	all_nan_mask = np.isnan(pixel_tiv_with_continuous_data_timing).all(axis=0)	1✔
163
164	# Initialize the mean_tiv_pixel array with NaNs
165	mean_tiv_pixel = np.full((n_rows, n_cols), np.nan)	1✔
166
167	# Only compute nanmean for slices that are not entirely NaN
168	if not all_nan_mask.all(): # Check if there are any valid (non-all-NaN) slices	1✔
169	mean_tiv_pixel[~all_nan_mask] = np.nanmean(pixel_tiv_with_continuous_data_timing[:, ~all_nan_mask], axis=0)	1✔
170
171	time = eit_data.time	1✔
172	pixel_impedance = eit_data.pixel_impedance	1✔
173
174	pixel_breaths = np.full((len(continuous_breaths), n_rows, n_cols), None, dtype=object)	1✔
175
176	lags = signal.correlation_lags(len(continuous_data), len(continuous_data), mode="same")	1✔
177
178	match self.phase_correction_mode:	1✔
179	case "negative amplitude":	1✔
180	allow_negative_amplitude = True	1✔
181	correct_for_phase_shift = None	1✔
182	case "phase shift":	1!
183	allow_negative_amplitude = False	1✔
184	correct_for_phase_shift = True	1✔
185	case "none" \| None:	×
186	allow_negative_amplitude = False	×
187	correct_for_phase_shift = False	×
188	case _:	×
189	msg = f"Unknown phase correction mode ({self.phase_correction_mode})."	×
190	raise ValueError(msg)	×
191
192	for row, col in itertools.product(range(n_rows), range(n_cols)):	1✔
193	if all_nan_mask[row, col]:	1✔
194	# pixel has no amplitude
195	continue	1✔
196	mean_tiv = mean_tiv_pixel[row, col]	1✔
197
198	pi = np.copy(pixel_impedance[:, row, col]) # TODO: check whether copying is necessary	1✔
199	if np.std(pi) == 0:	1✔
200	# pixel has no amplitude
201	continue	1✔
202
203	if allow_negative_amplitude and mean_tiv < 0:	1✔
204	start_func, middle_func = np.argmax, np.argmin	1✔
205	lagged_indices_breath_middles = indices_breath_middles	1✔
206	else:
207	start_func, middle_func = np.argmin, np.argmax	1✔
208
209	cd = np.copy(continuous_data.values)	1✔
210	cd = signal.detrend(cd, type="linear")	1✔
211
212	if not np.all(np.isnan(pi)):	1!
213	pi = signal.detrend(pi, type="linear")	1✔
214
215	if correct_for_phase_shift:	1✔
216	# search for closest or maximum cross correlation within MAX_XCORR_LAG times the average duration of
217	# a breath
218	xcorr = signal.correlate(cd, pi, mode="same")	1✔
219	max_lag = MAX_XCORR_LAG * np.mean(np.diff(indices_breath_middles))	1✔
220	lag_range = (lags >= -max_lag) & (lags <= max_lag)	1✔
221
222	# find the peaks in the cross correlation
223	peaks, _ = signal.find_peaks(xcorr[lag_range], height=0, prominence=xcorr.std())	1✔
224
225	# find the closest peak to zero lag
226	peak_lags = lags[lag_range][peaks]	1✔
227	peak_distances = np.abs(peak_lags)	1✔
228	min_peak_distance = np.min(peak_distances)	1✔
229	candidates = peak_lags[peak_distances == min_peak_distance]	1✔
230
231	# if there are multiple candidates, take the one with the highest cross correlation
232	if len(candidates) == 1:	1✔
233	lag = candidates[0]	1✔
234	elif len(candidates) == 2: # noqa: PLR2004	1!
235	# take the lag with the highest cross correlation
236	lag = candidates[np.argmax(xcorr[np.searchsorted(lags, candidates)])]	1✔
237	else:
UNCOV 238	msg = "Too many peaks found in cross correlation."	×
UNCOV 239	raise RuntimeError(msg)	×
240
241	# shift search area
242	lagged_indices_breath_middles = indices_breath_middles - lag	1✔
243	else:
244	lagged_indices_breath_middles = indices_breath_middles	1✔
245
246	skip = np.concatenate(	1✔
247	(
248	np.flatnonzero(lagged_indices_breath_middles < 0),
249	np.flatnonzero(lagged_indices_breath_middles >= len(continuous_data)),
250	)
251	)
252	lagged_indices_breath_middles = lagged_indices_breath_middles.clip(min=0, max=len(continuous_data) - 1)	1✔
253	outsides = self._find_extreme_indices(pixel_impedance, lagged_indices_breath_middles, row, col, start_func)	1✔
254	starts = outsides[:-1]	1✔
255	ends = outsides[1:]	1✔
256	middles = self._find_extreme_indices(pixel_impedance, outsides, row, col, middle_func)	1✔
257
258	skip = np.concatenate((skip, np.flatnonzero(starts >= middles), np.flatnonzero(middles >= ends)))	1✔
259
260	if len(starts) > len(skip) > len(starts) * ALLOW_FRACTION_BREATHS_SKIPPED:	1!
UNCOV 261	warnings.warn(	×
262	f"Skipping pixel ({row}, {col}) because more than half ({len(skip) / len(starts):.0%}) "
263	"of breaths skipped.",
264	RuntimeWarning,
265	stacklevel=2,
266	)
UNCOV 267	continue	×
268
269	pixel_breaths[:, row, col] = self._construct_breaths(starts, middles, ends, time, skip=skip)	1✔
270
271	intervals = [(breath.start_time, breath.end_time) for breath in continuous_breaths.values]	1✔
272
273	pixel_breaths_container = IntervalData(	1✔
274	label=result_label,
275	name="Pixel in- and deflation timing as determined by Pixelbreath",
276	unit=None,
277	category="breath",
278	intervals=intervals,
279	values=list(
280	pixel_breaths,
281	), ## TODO: change back to pixel_breaths array when IntervalData works with 3D array
282	derived_from=[eit_data],
283	)
284	if store:	1✔
285	sequence.interval_data.add(pixel_breaths_container)	1✔
286
287	return pixel_breaths_container	1✔
288
289	def _construct_breaths(	1✔
290	self, start: list[int], middle: list[int], end: list[int], time: np.ndarray, skip: np.ndarray \| None = None
291	) -> list:
292	skip_ = skip if skip is not None else np.array([], dtype=int)	1✔
293	breaths = [	1✔
294	Breath(time[s], time[m], time[e]) if i not in skip_ else None
295	for i, (s, m, e) in enumerate(zip(start, middle, end, strict=True))
296	]
297	# First and last breath are not detected by definition (need two breaths to find one breath)
298	return [None, *breaths, None]	1✔
299
300	def _find_extreme_indices(	1✔
301	self,
302	pixel_impedance: np.ndarray,
303	times: np.ndarray,
304	row: int,
305	col: int,
306	function: Callable,
307	) -> np.ndarray:
308	"""Finds extreme indices in pixel impedance.
309
310	This method divides the pixel impedance for a single pixel (selected using row and col) into smaller segments
311	based on the `times` array. The times array consists of indices to divide these segments. The function iterates
312	over each index in the times array to select consecutive segments of pixel impedance.
313	For each segment, the method applies the `function` (either np.argmax or np.argmin) to extract an extreme value
314	(local maximum or minimum).
315
316	Args:
317	pixel_impedance (np.ndarray): The pixel impedance array from which the function will extract values.
318	Assumed to be 3-dimensional (e.g., time, rows, and columns).
319	times (np.ndarray): 1D array of time indices. These times define the start
320	and end of each segment in the pixel impedance.
321	row (int): The row index in the pixel impedance
322	col (int): The column index in the pixel impedance
323	function (Callable): A function that is applied to each segment of data to find
324	an extreme value (np.argmax or np.argmin)
325
326	Returns:
327	np.ndarray: An array of indices where the extreme values (based on the function)
328	are located for each time segment.
329	"""
330	return np.array(	1✔
331	[function(pixel_impedance[t1:t2, row, col]) + t1 for t1, t2 in itertools.pairwise(times)],
332	)

EIT-ALIVE / eitprocessing / 17495122789

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous