diff --git a/.github/workflows/release_pypi.yml b/.github/workflows/release_pypi.yml index e378911d8..07bc812ff 100644 --- a/.github/workflows/release_pypi.yml +++ b/.github/workflows/release_pypi.yml @@ -27,27 +27,26 @@ jobs: needs: [build] runs-on: ubuntu-latest if: github.event_name == 'workflow_dispatch' + permissions: + id-token: write steps: - uses: actions/download-artifact@v4 with: name: artifact path: dist - - uses: pypa/gh-action-pypi-publish@v1.4.2 + - uses: pypa/gh-action-pypi-publish@release/v1 with: - user: __token__ - password: ${{ secrets.TEST_PYPI_TOKEN_ALIVE }} - repository_url: https://test.pypi.org/legacy/ + repository-url: https://test.pypi.org/legacy/ upload_pypi: needs: [build] runs-on: ubuntu-latest if: github.event_name == 'release' && github.event.action == 'published' + permissions: + id-token: write steps: - uses: actions/download-artifact@v4 with: name: artifact path: dist - - uses: pypa/gh-action-pypi-publish@v1.4.2 - with: - user: __token__ - password: ${{ secrets.PYPI_TOKEN_ALIVE }} + - uses: pypa/gh-action-pypi-publish@release/v1 diff --git a/eitprocessing/filters/butterworth_filters.py b/eitprocessing/filters/butterworth_filters.py index 459e2b5a9..f7a74609a 100644 --- a/eitprocessing/filters/butterworth_filters.py +++ b/eitprocessing/filters/butterworth_filters.py @@ -52,7 +52,7 @@ class ButterworthFilter(TimeDomainFilter): """ filter_type: Literal["lowpass", "highpass", "bandpass", "bandstop"] - cutoff_frequency: float | tuple[float] + cutoff_frequency: float | tuple[float, float] order: int sample_frequency: float ignore_max_order: InitVar[bool] = False diff --git a/eitprocessing/parameters/odcl.py b/eitprocessing/parameters/odcl.py new file mode 100644 index 000000000..6e8516699 --- /dev/null +++ b/eitprocessing/parameters/odcl.py @@ -0,0 +1,142 @@ +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING, Literal + +import numpy as np + +from eitprocessing.datahandling.pixelmap import ODCLMap +from eitprocessing.parameters.tidal_impedance_variation import TIV +from eitprocessing.roi import PixelMask + +if TYPE_CHECKING: + from eitprocessing.datahandling.eitdata import EITData + + +@dataclass(frozen=True, kw_only=True) +class OdclPeepTrial: + """Calculate relative overdistention and collapse from a series of PEEP steps. + + Provided a list of EITData objects from a decremental PEEP trial, calculates the relative overdistention and + collapse according to Costa et al. (2009). + """ + + compliance_method: Literal["constant driving pressure"] = "constant driving pressure" + + def apply( + self, + eit_datasets: list[EITData], + ) -> OdclPeepTrialResult: + """Apply the algorithm to the provided datasets. + + The provided datasets should be in decremental PEEP order, i.e. in order of measurement. + + Returns an object containing an ODCLMap, overdistention an collapse value per provided PEEP step. + + Arguments: + eit_datasets (list or EITData): EITData objects + + Returns: + OdclPeepTrialResult: object containing results and plotting functionality + odcl_maps: list of ODCLMap + overdistention: array of overdistention values + collapse: array of collapse values + plotting: plotting functionality + + """ + n_steps = len(eit_datasets) + continuous_datasets = [eit_data.get_summed_impedance() for eit_data in eit_datasets] + + mask = PixelMask( + ~np.all( + [np.all(np.isnan(eit_data.pixel_impedance), axis=0) for eit_data in eit_datasets], + axis=0, + ) + ) # infer mask from all-NaN pixels + + tiv = TIV() + tiv_per_breath = [ + tiv.compute_pixel_parameter( + eit_data, + continuous_data=continuous_data, + sequence=None, + store=False, + ) + for eit_data, continuous_data in zip(eit_datasets, continuous_datasets, strict=True) + ] + + mean_tiv_per_step = np.asarray([np.nanmean(maps.values, axis=0) for maps in tiv_per_breath]) + masked_tiv_per_step = mask.apply(mean_tiv_per_step) + + if self.compliance_method == "constant driving pressure": + compliance = masked_tiv_per_step + else: + msg = f"Unknown compliance method '{self.compliance_method}'" + raise ValueError(msg) + + compliance = np.clip(compliance, 0, np.inf) + + per_pixel_max_compliance_idx = np.argmax(compliance, axis=0)[np.newaxis] + + max_compliance_per_pixel = np.take_along_axis( + compliance, + per_pixel_max_compliance_idx, + axis=0, + ) + + normalized_max_compliance = max_compliance_per_pixel / np.nanmax(max_compliance_per_pixel) + relative_compliance = compliance / max_compliance_per_pixel + relative_compliance_change = (1 - relative_compliance) * normalized_max_compliance + + steps = np.arange(n_steps)[:, np.newaxis, np.newaxis] + + overdistention = np.where( + steps < per_pixel_max_compliance_idx, + relative_compliance_change, + 0, + ) + collapse = np.where( + steps > per_pixel_max_compliance_idx, + relative_compliance_change, + 0, + ) + + total_overdistention = np.nansum(overdistention, axis=(1, 2)) / np.nansum(normalized_max_compliance) + total_collapse = np.nansum(collapse, axis=(1, 2)) / np.nansum(normalized_max_compliance) + + odcl_maps = [ODCLMap(map_) for map_ in mask.apply(overdistention - collapse)] + + return OdclPeepTrialResult( + odcl_maps=odcl_maps, + overdistention=total_overdistention, + collapse=total_collapse, + ) + + +@dataclass(kw_only=True, frozen=True) +class OdclPeepTrialResult: + """Result object for the OdclPeepTrial algorithm. + + Contains an ODCLMap, overdistention value and collapse value per PEEP step. + + Attributes: + odcl_maps: list of ODCLMap pixel maps + overdistention: array of floats + collapse: array of floats + plotting: plotting functionality + plot_overdistention_collapse: + plot the overdistention and collapse values in a line plot + plot_odcl_maps: + show the ODCL maps in a grid + """ + + odcl_maps: list[ODCLMap] + overdistention: np.ndarray + collapse: np.ndarray + + @property + def plotting(self): + """Plotting functionality for ODCL PEEP Trial results.""" + from eitprocessing.plotting.odcl import OdclPeepTrialPlotting + + return OdclPeepTrialPlotting(self) diff --git a/eitprocessing/plotting/odcl.py b/eitprocessing/plotting/odcl.py new file mode 100644 index 000000000..74716f08a --- /dev/null +++ b/eitprocessing/plotting/odcl.py @@ -0,0 +1,124 @@ +from dataclasses import dataclass + +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.axes import Axes +from matplotlib.figure import Figure +from matplotlib.ticker import PercentFormatter + +from eitprocessing.parameters.odcl import OdclPeepTrialResults + + +@dataclass(frozen=True) +class OdclPeepTrialPlotting: + """Plotting functionality for OdclPeepTrialResults.""" + + results: OdclPeepTrialResults + + def plot_overdistention_collapse( + self, + *, + ax: Axes | None = None, + labels: list[str] | None = None, + xlabel: str | None = None, + ) -> tuple[Figure, Axes]: + """Plot overdistention and collapse values as line plots. + + If an axis is provided, the plots are created there. Otherwise, a new figure and axis is created. + + Labels can be provided for the steps, which will be shown as tick labels. + The label for the x-axis can also be provided. + + Returns the figure and axis. + """ + if ax is None: + fig, ax = plt.subplots(layout="constrained") + else: + fig = None + + x = np.arange(1, len(self.overdistention) + 1) + ax.plot( + x, + self.overdistention, + label="Overdistention", + marker="o", + lw=3, + color="darkorange", + ) + ax.plot(x, self.collapse, label="Collapse", marker="s", lw=3, color="#707070") + ax.set(xticks=x) + ax.yaxis.set_major_formatter(PercentFormatter(xmax=1)) + ax.grid() + if fig: + fig.legend(loc="outside upper center", ncol=2, frameon=False) + else: + ax.legend(loc="upper center", ncol=2) + + ndigits = int(np.log10(len(x))) + 1 + if labels is None: + labels = [f"#{x_:0{ndigits}d}" for x_ in x] + + xlabel = xlabel or "PEEP steps (decremental)" + + ax.set(xticklabels=labels, xlabel=xlabel) + return fig, ax + + def plot_odcl_maps( + self, + *, + labels: list[str] | None = None, + ncols: int = 4, + figsize: tuple[float, float] | None = None, + ) -> tuple[Figure, np.ndarray]: + """Create plots of the ODCL maps. + + Generates a n_rows x n_cols grid of axes fitting all PEEP steps. Provided labels are set as axis titles. + + Returns the figure and axes. + """ + n_steps = len(self.odcl_maps) + n_rows = int(np.ceil(n_steps / ncols)) + fig, axes = plt.subplots( + n_rows, + ncols, + layout="compressed", + sharex=True, + sharey=True, + figsize=figsize or (2 * ncols + 1, 2.4 * n_rows), + ) + + for axis, map_, od, cl in zip( + axes.flatten(), + self.odcl_maps, + self.overdistention, + self.collapse, + strict=False, + ): + map_.plotting.imshow(ax=axis) + axis.set_xlabel(f"OD: {od:.0%}, CL: {cl:.0%}") + + n_digits = int(np.log10(n_steps)) + 1 + if labels is None: + labels = [f"#{x + 1:0{n_digits}d}" for x in range(n_steps)] + for axis, label in zip(axes.flatten(), labels, strict=False): + axis.set(title=label) + + images = [axis.images[0] for axis in axes.flat if axis.images] + old_colorbars = [image.colorbar for image in images if image.colorbar is not None] + old = old_colorbars[0] + + fig.colorbar( + images[0], + ax=axes, + format=old.formatter, + extend=old.extend, + label=old.ax.get_ylabel(), + ) + + for colorbar in old_colorbars: + colorbar.remove() + + for axis in axes.flatten()[len(self.odcl_maps) :]: + axis.remove() + + return fig, axes diff --git a/eitprocessing/roi/__init__.py b/eitprocessing/roi/__init__.py index 592dafda6..3abb0cfcf 100644 --- a/eitprocessing/roi/__init__.py +++ b/eitprocessing/roi/__init__.py @@ -4,21 +4,8 @@ module is `PixelMask`. Any type of region of interest selection results in a `PixelMask` object. A mask can be applied to any pixel dataset (EITData, PixelMap) with the same shape. -Several default masks have been predefined. NB: the right side of the patient is to the left side of the EIT image and -vice versa. - -- `VENTRAL_MASK` includes only the first 16 rows; -- `DORSAL_MASK` includes only the last 16 rows; -- `ANATOMICAL_RIGHT_MASK` includes only the first 16 columns; -- `ANATOMICAL_LEFT_MASK` includes only the last 16 columns; -- `QUADRANT_1_MASK` includes the top right quadrant; -- `QUADRANT_2_MASK` includes the top left quadrant; -- `QUADRANT_3_MASK` includes the bottom right quadrant; -- `QUADRANT_4_MASK` includes the bottom left quadrant; -- `LAYER_1_MASK` includes only the first 8 rows; -- `LAYER_2_MASK` includes only the second set of 8 rows; -- `LAYER_3_MASK` includes only the third set of 8 rows; -- `LAYER_4_MASK` includes only the last 8 rows. +Predefined default masks can be created using `get_geometric_mask()`. These masks are based on common anatomical regions +of interest. """ from __future__ import annotations @@ -296,7 +283,7 @@ def __eq__(self, other: object) -> bool: return self.mask.shape == other.mask.shape and np.array_equal(self.mask, other.mask, equal_nan=True) -def get_geometric_mask(mask: str, shape: tuple[int, int] = (32, 32)) -> PixelMask: +def get_geometric_mask(mask: str, shape: tuple[int, int] = (32, 32)) -> PixelMask: # noqa: PLR0911 """Get a geometric mask by name. Masks can be generated for appropriates shapes, provided by the shape` argument, a tuple of two integers @@ -306,6 +293,7 @@ def get_geometric_mask(mask: str, shape: tuple[int, int] = (32, 32)) -> PixelMas The function accepts both full names (e.g., "layer 1") and abbreviations (e.g., "L1"). The following masks are available: + - "global" or "G": all pixels in the EIT image. - "ventral" or "V": the first half rows of the EIT image. - "dorsal" or "D": the last half rows of the EIT image. - "anatomical right" or "R": the first half columns of the EIT image. @@ -332,23 +320,6 @@ def get_geometric_mask(mask: str, shape: tuple[int, int] = (32, 32)) -> PixelMas ValueError: If an unknown mask name is provided. ValueError: If the shape is not compatible with the requested mask. """ - - def _check_dimensions( - name: str, shape: tuple[int, int], *, height_divisor: int | None = None, width_divisor: int | None = None - ) -> None: - total_height, total_width = shape - if isinstance(height_divisor, int) and total_height % height_divisor != 0: - msg = ( - f"Shape {shape} is not compatible with a {name} mask. " - "The height must be a multiple of {height_divisor}." - ) - raise ValueError(msg) - if isinstance(width_divisor, int) and total_width % width_divisor != 0: - msg = ( - f"Shape {shape} is not compatible with a {name} mask. The width must be a multiple of {width_divisor}." - ) - raise ValueError(msg) - total_height, total_width = shape n_layers_quadrants = 4 @@ -357,6 +328,8 @@ def _check_dimensions( mask = re.sub(r"^L([1-4]{1})$", r"layer \1", mask) match mask.split(" "): + case ["global"] | ["G"]: + return PixelMask(np.ones(shape), label="global") case ["ventral"] | ["V"]: _check_dimensions("ventral", shape, height_divisor=2) height = total_height // 2 @@ -387,7 +360,7 @@ def _check_dimensions( case ["layer", num_str] if num_str.isdigit() and 1 <= int(num_str) <= n_layers_quadrants: num = int(num_str) - _check_dimensions("layer", shape, height_divisor=4) + _check_dimensions("layer", shape, height_divisor=n_layers_quadrants) height = total_height // n_layers_quadrants return PixelMask( np.concatenate( @@ -411,3 +384,17 @@ def _check_dimensions( case _: msg = f"Unknown mask name: {mask}." raise ValueError(msg) + + +def _check_dimensions( + name: str, shape: tuple[int, int], *, height_divisor: int | None = None, width_divisor: int | None = None +) -> None: + total_height, total_width = shape + if isinstance(height_divisor, int) and total_height % height_divisor != 0: + msg = ( + f"Shape {shape} is not compatible with a {name} mask. The height must be a multiple of {{height_divisor}}." + ) + raise ValueError(msg) + if isinstance(width_divisor, int) and total_width % width_divisor != 0: + msg = f"Shape {shape} is not compatible with a {name} mask. The width must be a multiple of {width_divisor}." + raise ValueError(msg) diff --git a/pyproject.toml b/pyproject.toml index b8c989308..73aaeec30 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -110,6 +110,7 @@ eitprocessing = ["config/*.yaml", "py.typed"] output-format = "concise" line-length = 120 extend-include = ["*.ipynb"] +target-version = "py310" [tool.ruff.lint] select = ["ALL"] diff --git a/tests/test_pixelmask.py b/tests/test_pixelmask.py index a69e16698..eb2d5e5a2 100644 --- a/tests/test_pixelmask.py +++ b/tests/test_pixelmask.py @@ -220,6 +220,14 @@ def test_pixelmask_subtract(): assert np.array_equal(pm3.mask, np.array([[np.nan, np.nan, 1, np.nan], [0.1, np.nan, np.nan, 0.5]]), equal_nan=True) +@pytest.mark.parametrize("shape", [(32, 1), (64, 1), (16, 1), (100, 1), (4, 1)]) +def test_predefined_global_mask(shape: tuple[int, int]): + global_mask = get_geometric_mask("global", shape) + assert global_mask.label == "global" + assert global_mask.shape == shape + assert np.all(global_mask.mask == 1.0) + + @pytest.mark.parametrize("shape", [(32, 1), (64, 1), (16, 1), (100, 1), (4, 1)]) def test_predefined_layer_masks(shape: tuple[int, int]): quarter_height = shape[0] // 4