add pixel analysis tool to close #27, and add Python 3.12 support to c…

…lose #26

.github/workflows/format.yaml

@@ -7,7 +7,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.10", "3.11"]
+        python-version: ["3.10", "3.11", "3.12"]
         os: [ubuntu-22.04]
     runs-on: ${{ matrix.os }}
 

.github/workflows/lint.yaml

@@ -7,7 +7,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.10", "3.11"]
+        python-version: ["3.10", "3.11", "3.12"]
         os: [ubuntu-22.04]
     runs-on: ${{ matrix.os }}
 

.github/workflows/tests.yaml

@@ -7,7 +7,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.10", "3.11"]
+        python-version: ["3.10", "3.11", "3.12"]
         os: [ubuntu-latest, ubuntu-20.04, macos-latest, windows-latest]
     runs-on: ${{ matrix.os }}
 

CHANGELOG.md

@@ -4,6 +4,16 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [v0.5.0] - 2024-11-21
+
+### Added
+* `ImagingChannel` wrapper type
+* Cross-chanel signal extraction/analysis tool for images; see [Issue 27](https://github.com/gerlichlab/gertils/issues/27).
+
+### Changed
+* Depend on v2.2.1 of `numpydoc_decorator` directly from PyPI, rather than our custom release from earlier.
+* Support Python 3.12
+
 ## [v0.4.4] - 2024-04-19
 
 ### Fixed

README.md

@@ -9,6 +9,7 @@ These tools are organised by use case at the module level; that is, tools that a
 - [geometry](./gertils/geometry.py) -- tools for working with entities in space
 - [gpu](./gertils/gpu.py) -- tools for running computations on GPUs, especially with TensorFlow
 - [pathtools](./gertils/pathtools.py) -- tools for working with filesystem paths generally
+- [pixel_value_statistics](./gertils/pixel_value_statistics.py) -- tools for computing statistics of pixel values
 - [types](./gertils/pathtools.py) -- data types for working with genome biology, especially imaging
 - [zarr_tools](./gertils/zarr_tools.py) -- functions and types for working with ZARR-stored data
 

gertils/__init__.py

@@ -12,3 +12,7 @@
     find_single_path_by_fov,
     get_fov_sort_key,
 )
+from .pixel_value_statistics import (
+    RegionalPixelStatistics,
+    compute_pixel_statistics,
+)

gertils/pixel_value_statistics.py

@@ -0,0 +1,144 @@
+"""Tools for taking statistics over pixel values"""
+
+import dataclasses
+import logging
+from collections.abc import Iterable
+from typing import TypeAlias
+
+import numpy as np
+import numpy.typing as npt
+from numpydoc_decorator import doc  # type: ignore[import]
+
+from .geometry import ImagePoint3D, ZCoordinate
+from .types import ImagingChannel
+
+__all__ = ["RegionalPixelStatistics", "compute_pixel_statistics"]
+
+Numeric: TypeAlias = (
+    float | int | np.float16 | np.float32 | np.float64 | np.int8 | np.int16 | np.int32 | np.int64
+)
+PixelValue: TypeAlias = np.uint8 | np.uint16
+
+
+@doc(
+    summary="Store a handful of pixel value stats for a particular ROI.",
+    parameters=dict(
+        mean_value="The mean pixel intensity in a ROI",
+        sigma_value="The standard deviation of pixel intensity in a ROI",
+        min_value="The minimum pixel value in a ROI",
+        med_value="The median pixel value in a ROI",
+        max_value="The maximum pixel value in a ROI",
+        proj_mean="Mean of values in the max-z-projection",
+        proj_sigma="Standard deviation of values in the max-z-projection",
+        proj_min="Minimum of values in the max-z-projection",
+        proj_med="Median of values in the max-z-projection",
+        proj_max="Maximum of values in the max-z-projection",
+        center_mean="Mean of values in the central z-slice",
+        center_sigma="Standard deviation of values in the central z-slice",
+        center_min="Minimum of values in the central z-slice",
+        center_med="Median of values in the central z-slice",
+        center_max="Maximum of values in the central z-slice",
+    ),
+)
+@dataclasses.dataclass(kw_only=True, frozen=True)
+class RegionalPixelStatistics:  # noqa: D101
+    mean_value: float
+    sigma_value: float
+    min_value: float
+    med_value: float
+    max_value: float
+    proj_mean: float
+    proj_sigma: float
+    proj_min: float
+    proj_med: float
+    proj_max: float
+    center_mean: float
+    center_sigma: float
+    center_min: float
+    center_med: float
+    center_max: float
+
+    @property
+    def to_dict(self) -> dict[str, float]:  # noqa: D102
+        return dataclasses.asdict(self)
+
+    @classmethod
+    def from_image(
+        cls, img: npt.NDArray[PixelValue], central_z: ZCoordinate
+    ) -> "RegionalPixelStatistics":
+        """Compute stats for the given region (defined by whole given image)."""
+        if len(img.shape) != 3:  # noqa: PLR2004
+            raise ValueError(f"To build {cls.__name__}, image must be 3D, not {len(img.shape)}D")
+        round_z: int = int(round(central_z))
+        if round_z < 0:
+            raise ValueError(
+                f"Cannot extract pixel values from negative z-slice. ({round_z}, from {central_z})"
+            )
+        if round_z == img.shape[0] and central_z < img.shape[0]:
+            logging.warning(
+                f"Rounding central_z down from {central_z} to comply with z-depth of {img.shape[0]}"  # noqa: G004
+            )
+            round_z = int(central_z)
+        elif round_z >= img.shape[0]:
+            raise ValueError(
+                f"Cannot extract pixel values from z-slice ({round_z}, from {central_z}) for image with {img.shape[0]} z-slice(s)."
+            )
+
+        central_plane_img = img[round_z]
+        maxproj = np.max(img, axis=0)
+
+        return cls(
+            mean_value=img.mean(),
+            sigma_value=img.std(),
+            min_value=img.min(),
+            med_value=np.median(img),  # type: ignore[arg-type]
+            max_value=img.max(),
+            proj_mean=maxproj.mean(),
+            proj_sigma=maxproj.std(),
+            proj_min=maxproj.min(),
+            proj_med=np.median(maxproj),
+            proj_max=maxproj.max(),
+            center_mean=central_plane_img.mean(),
+            center_sigma=central_plane_img.std(),
+            center_min=central_plane_img.min(),
+            center_med=np.median(central_plane_img),
+            center_max=central_plane_img.max(),
+        )
+
+
+@doc(
+    summary="Compute statistics over pixels from multiple channels, in a region centered on a point.",
+    parameters=dict(
+        img="Image in which to measure pixels",
+        pt="Center of region to measure",
+        channels="Channels of image in which to measure pixels",
+        diameter="Size (width and height) of region around point in which to measure pixels",
+        signal_column="Name for the field/column in which to store channel from which pixels were taken",
+    ),
+    returns="List of records, each mapping key/field to value",
+)
+def compute_pixel_statistics(  # noqa: D103
+    img: npt.NDArray[PixelValue],
+    pt: ImagePoint3D,
+    *,
+    channels: Iterable[ImagingChannel],
+    diameter: int,
+    signal_column: str,
+) -> list[dict[str, Numeric]]:
+    left: int = round(pt.x - diameter / 2)
+    right: int = left + diameter
+    top: int = round(pt.y - diameter / 2)
+    bottom: int = top + diameter
+    bounds: dict[str, int] = {
+        "y_min_px": top,
+        "y_max_px": bottom,
+        "x_min_px": left,
+        "x_max_px": right,
+    }
+    # Build up records, e.g. rows of data table/frame
+    result: list[dict[str, Numeric]] = []
+    for ch in channels:
+        subimg = img[ch.get, :, max(0, top) : bottom, max(0, left) : right]
+        stats = RegionalPixelStatistics.from_image(subimg, central_z=pt.z)
+        result.append({signal_column: ch.get, **bounds, **stats.to_dict})
+    return result

gertils/types.py

@@ -19,7 +19,7 @@
 
 @doc(
     summary="Wrap an int as a 1-based field of view (FOV).",
-    parameters=dict(get="The value to wrape"),
+    parameters=dict(get="The value to wrap"),
     raises=dict(
         TypeError="If the given value to wrap isn't an integer",
         ValueError="If the given value to wrap isn't positive",
@@ -38,6 +38,18 @@ def __post_init__(self) -> None:
             raise ValueError(f"1-based FOV view must be positive int; got {self.get}")
 
 
+@doc(summary="An imaging channel", parameters=dict(get="Wrapped value"))
+@dataclass(frozen=True, order=True)
+class ImagingChannel:  # noqa: D101
+    get: int
+
+    def __post_init__(self) -> None:
+        if not isinstance(self.get, int):
+            raise TypeError(f"Imaging channel must be integer, not {type(self.get).__name__}")
+        if self.get < 0:
+            raise ValueError(f"Imaging channel must be nonnegative; got {self.get}")
+
+
 @doc(
     summary="Wrap an int as a 1-based nucleus number.",
     parameters=dict(get="The value to wrap"),

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "gertils"
-version = "0.4.4"
+version = "0.5.0"
 description = "General computing utilities for data processing and analysis in the Gerlich group at IMBA"
 authors = [
     "Vince Reuter <vince.reuter@gmail.com>",
@@ -13,6 +13,7 @@ classifiers = [
     "Programming Language :: Python",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
+    "Programming Language :: Python :: 3.12",
     "Operating System :: POSIX :: Linux",
     "Operating System :: MacOS",
     "Operating System :: Microsoft :: Windows",
@@ -28,10 +29,10 @@ build-backend = "poetry.core.masonry.api"
 
 [tool.poetry.dependencies]
 # These are the main runtime dependencies.
-python = ">=3.10.0,<3.12"
+python = ">=3.10.0,<3.13"
 dask = "^2023.5.1"
 numpy = "^1.24.2"
-numpydoc_decorator = { git = "https://github.com/vreuter/numpydoc_decorator", tag = "v2.2.1" }
+numpydoc_decorator = "^2.2.1"
 zarr = "^2.4.12"
 
 [tool.poetry.group.dev]

tests/test_gertils_top_level.py

@@ -17,6 +17,7 @@
 @pytest.mark.parametrize(
     ("member_name", "expected_presence"),
     [(name, True) for name in PATHTOOLS_PUBLIC_MEMBERS]
+    + [(name, True) for name in ["RegionalPixelStatistics", "compute_pixel_statistics"]]
     + [(name, False) for name in envs_module.__all__ + COLLECTIONS_PUBLIC_MEMBERS],
 )
 def test_import_visibility(member_name, expected_presence):