diff --git a/src/tracksdata/nodes/_regionprops.py b/src/tracksdata/nodes/_regionprops.py index 11ad7261..201f764c 100644 --- a/src/tracksdata/nodes/_regionprops.py +++ b/src/tracksdata/nodes/_regionprops.py @@ -1,4 +1,4 @@ -from collections.abc import Callable +from collections.abc import Callable, Sequence from functools import partial from typing import Any @@ -6,7 +6,12 @@ import polars as pl from numpy.typing import NDArray from polars.datatypes import numpy_char_code_to_dtype -from skimage.measure._regionprops import RegionProperties, regionprops +from skimage.measure._regionprops import ( + PROPS, + RegionProperties, + _require_intensity_image, + regionprops, +) from typing_extensions import override from tracksdata.constants import DEFAULT_ATTR_KEYS @@ -37,6 +42,21 @@ class RegionPropsNodes(BaseNodesOperator): Physical spacing between pixels. If provided, affects distance-based measurements. Should be (row_spacing, col_spacing) for 2D or (depth_spacing, row_spacing, col_spacing) for 3D. + separate_arrays : bool, optional + If True, array-like properties (e.g. ``inertia_tensor`` or multi-channel + ``intensity_mean``) are flattened into multiple scalar attributes instead + of being stored as a single array attribute. The new attributes are named + ``_`` (e.g. ``intensity_mean_0``, ``intensity_mean_1``) using + the same convention as ``node_attrs(unpack=True)``. This makes individual + components filterable. Defaults to False. + channel_names : Sequence[str] | None, optional + Names of the channels of a multi-channel `intensity_image`, whose channels + are expected on its last axis. When provided, intensity-dependent properties + (e.g. ``intensity_mean``, ``centroid_weighted``) are split along the channel + axis into one attribute per channel, named ``_`` + (e.g. ``intensity_mean_dapi``) instead of a single array attribute. + Custom callables in `extra_properties` are never split, since their output + layout is unknown. Defaults to None (no splitting). Attributes ---------- @@ -44,6 +64,10 @@ class RegionPropsNodes(BaseNodesOperator): List of additional properties to compute. _spacing : tuple[float, float] | None Physical spacing between pixels. + _separate_arrays : bool + Whether array-like properties are flattened into scalar attributes. + _channel_names : list[str] | None + Names used to split intensity properties along the channel axis. Examples -------- @@ -86,12 +110,26 @@ def custom_property(region): labels_series = np.random.randint(0, 10, (10, 100, 100)) node_op.add_nodes(graph, labels=labels_series) ``` + + Name the channels of a multi-channel intensity image: + + ```python + node_op = RegionPropsNodes( + extra_properties=["intensity_mean"], + channel_names=["dapi", "gfp"], + ) + # intensity_image has shape (t, y, x, 2) + node_op.add_nodes(graph, labels=labels, intensity_image=intensity) + # nodes have `intensity_mean_dapi` and `intensity_mean_gfp` attributes + ``` """ def __init__( self, extra_properties: list[str | Callable[[RegionProperties], Any]] | None = None, spacing: tuple[float, float] | None = None, + separate_arrays: bool = False, + channel_names: Sequence[str] | None = None, ): super().__init__() self._extra_properties = extra_properties or [] @@ -102,6 +140,54 @@ def __init__( if "bbox" in self._extra_properties: raise ValueError("`bbox` is not supported as an extra property. It's already included by default.") self._spacing = spacing + self._separate_arrays = separate_arrays + self._channel_names = self._validate_channel_names(channel_names) + + @staticmethod + def _validate_channel_names(channel_names: Sequence[str] | None) -> list[str] | None: + """ + Validate and normalize the `channel_names` argument. + + Parameters + ---------- + channel_names : Sequence[str] | None + The channel names provided by the user. + + Returns + ------- + list[str] | None + The channel names as a list, or None when not provided. + + Raises + ------ + ValueError + If the names are empty, not strings, or not unique. + """ + if channel_names is None: + return None + + channel_names = list(channel_names) + if len(channel_names) == 0: + raise ValueError("`channel_names` must not be empty, use `None` for single-channel intensity images.") + + non_str = [name for name in channel_names if not isinstance(name, str)] + if non_str: + raise ValueError(f"`channel_names` must be strings, got {non_str}.") + + if len(set(channel_names)) != len(channel_names): + raise ValueError(f"`channel_names` must be unique, got {channel_names}.") + + return channel_names + + @staticmethod + def _is_intensity_prop(prop: str) -> bool: + """ + Whether a regionprops property name is computed from the intensity image. + + Such properties gain a trailing channel axis for multi-channel intensity + images, hence they are the ones split by `channel_names`. + """ + return PROPS.get(prop, prop) in _require_intensity_image def _axis_names(self, labels: NDArray[np.integer]) -> list[str]: """ @@ -124,6 +210,61 @@ def _axis_names(self, labels: NDArray[np.integer]) -> list[str]: else: raise ValueError(f"`labels` must be 't + 2D' or 't + 3D', got '{labels.ndim}' dimensions.") + def _attr_items(self, key: str, value: Any, split_channels: bool = False) -> list[tuple[str, Any]]: + """ + Normalize a single property value into one or more node attribute items. + + Tuple/list/array-like numeric values are converted to numpy arrays so they + are stored consistently as fixed-shape array attributes. When + ``separate_arrays`` is enabled, such values are instead flattened into + scalar attributes named ``_`` (row-major), matching the + ``node_attrs(unpack=True)`` naming convention. + + When ``split_channels`` is set, the last axis is first split into one item + per channel, named ``_``; the remainder of each channel's + value is then normalized as above. + + Parameters + ---------- + key : str + The base attribute name. + value : Any + The property value returned by regionprops or a custom callable. + split_channels : bool + Whether the last axis of `value` is a channel axis to be split using + `channel_names`. + + Returns + ------- + list[tuple[str, Any]] + The (name, value) pairs to add to the node attributes. + + Raises + ------ + ValueError + If the channel axis length does not match the number of channel names. + """ + if isinstance(value, np.ndarray | tuple | list): + arr = np.asarray(value) + if arr.dtype.kind in "biufc" and arr.ndim >= 1: + if split_channels: + if arr.shape[-1] != len(self._channel_names): + raise ValueError( + f"Property '{key}' has {arr.shape[-1]} channels, " + f"but {len(self._channel_names)} `channel_names` were provided: {self._channel_names}." + ) + items = [] + for i, name in enumerate(self._channel_names): + channel_value = arr[..., i] + # 0-dim arrays are stored as scalars, not as shape-() array attributes + items.extend(self._attr_items(f"{key}_{name}", channel_value[()])) + return items + if self._separate_arrays: + return [("_".join([key, *map(str, idx)]), arr[idx]) for idx in np.ndindex(arr.shape)] + return [(key, arr)] + + return [(key, value)] + def _init_node_attrs(self, graph: BaseGraph, node_attrs: dict[str, Any]) -> None: """ Initialize the node attributes for the graph. @@ -150,6 +291,9 @@ def attr_keys(self) -> list[str]: Returns only the keys for extra_properties. The centroid coordinates (x, y, z) and mask are always included but not listed here. + When `channel_names` is set, intensity-dependent properties are listed + once per channel. Note that `separate_arrays` further splits array-valued + properties into per-index keys that are not listed here. Returns ------- @@ -164,7 +308,15 @@ def attr_keys(self) -> list[str]: print(keys) # ['area', 'perimeter'] ``` """ - return [prop.__name__ if callable(prop) else prop for prop in self._extra_properties] + keys = [] + for prop in self._extra_properties: + if callable(prop): + keys.append(prop.__name__) + elif self._channel_names is not None and self._is_intensity_prop(prop): + keys.extend(f"{prop}_{name}" for name in self._channel_names) + else: + keys.append(prop) + return keys @override def add_nodes( @@ -200,7 +352,14 @@ def add_nodes( intensity_image : NDArray | None, optional Intensity image(s) corresponding to the labels. Used for computing intensity-based properties. Must have the same shape as labels - (excluding the label values). + (excluding the label values), plus a trailing channel axis when + `channel_names` is used. + + Raises + ------ + ValueError + If `channel_names` was provided but `intensity_image` is missing or + does not have a matching trailing channel axis. Examples -------- @@ -229,6 +388,16 @@ def add_nodes( node_op.add_nodes(graph, labels=labels, t=0, intensity_image=fluorescence_image) ``` """ + if self._channel_names is not None: + if intensity_image is None: + raise ValueError("`channel_names` was provided but `intensity_image` is None.") + if intensity_image.ndim != labels.ndim + 1 or intensity_image.shape[-1] != len(self._channel_names): + raise ValueError( + f"`intensity_image` must have shape '{(*labels.shape, len(self._channel_names))}' to match " + f"`labels` plus the {len(self._channel_names)} channels of " + f"`channel_names` {self._channel_names}, got '{intensity_image.shape}'." + ) + if "shape" not in graph.metadata: graph.metadata.update(shape=labels.shape) @@ -302,9 +471,11 @@ def _nodes_per_time( for prop in self._extra_properties: if callable(prop): - attrs[prop.__name__] = prop(obj) + key, value, split_channels = prop.__name__, prop(obj), False else: - attrs[prop] = getattr(obj, prop) + key, value = prop, getattr(obj, prop) + split_channels = self._channel_names is not None and self._is_intensity_prop(prop) + attrs.update(self._attr_items(key, value, split_channels=split_channels)) attrs[DEFAULT_ATTR_KEYS.MASK] = Mask(obj.image, obj.bbox) attrs[DEFAULT_ATTR_KEYS.BBOX] = np.asarray(obj.bbox, dtype=int) diff --git a/src/tracksdata/nodes/_test/test_regionprops.py b/src/tracksdata/nodes/_test/test_regionprops.py index 567c62e0..31073c6b 100644 --- a/src/tracksdata/nodes/_test/test_regionprops.py +++ b/src/tracksdata/nodes/_test/test_regionprops.py @@ -1,7 +1,9 @@ import numpy as np +import polars as pl import pytest from skimage.measure._regionprops import RegionProperties +from tracksdata.attrs import NodeAttr from tracksdata.constants import DEFAULT_ATTR_KEYS from tracksdata.graph import RustWorkXGraph from tracksdata.nodes import Mask, RegionPropsNodes @@ -334,3 +336,216 @@ def test_regionprops_multiprocessing_isolation() -> None: """Test that multiprocessing options don't affect subsequent tests.""" # Verify default n_workers is 1 assert get_options().n_workers == 1 + + +def test_regionprops_multichannel_intensity_array() -> None: + """Multi-channel intensity props are stored as fixed-shape array attributes (#195).""" + graph = RustWorkXGraph() + + labels = np.array([[[1, 1, 0], [1, 0, 2], [0, 2, 2]]], dtype=np.int32) + intensity = np.zeros((1, 3, 3, 2), dtype=np.float32) + intensity[..., 0] = [[10, 20, 0], [30, 0, 40], [0, 50, 60]] + intensity[..., 1] = [[1, 2, 0], [3, 0, 4], [0, 5, 6]] + + operator = RegionPropsNodes(extra_properties=["intensity_max"]) + operator.add_nodes(graph, labels=labels, intensity_image=intensity) + + nodes_df = graph.node_attrs() + assert isinstance(nodes_df.schema["intensity_max"], pl.Array) + assert nodes_df["intensity_max"].dtype.shape == (2,) + + +def test_regionprops_tuple_property_stored_as_array() -> None: + """Tuple-returning props (e.g. centroid_weighted) are normalized to array attributes (#191).""" + graph = RustWorkXGraph() + + labels = np.array([[[1, 1, 0], [1, 0, 2], [0, 2, 2]]], dtype=np.int32) + intensity = np.array([[[10, 20, 0], [30, 0, 40], [0, 50, 60]]], dtype=np.float32) + + operator = RegionPropsNodes(extra_properties=["centroid_weighted"]) + operator.add_nodes(graph, labels=labels, intensity_image=intensity) + + nodes_df = graph.node_attrs() + # tuple props must become fixed-shape arrays (not pl.List) so they are unpackable + assert isinstance(nodes_df.schema["centroid_weighted"], pl.Array) + unpacked = graph.node_attrs(unpack=True) + assert "centroid_weighted_0" in unpacked.columns + assert "centroid_weighted_1" in unpacked.columns + + +def test_regionprops_separate_arrays() -> None: + """`separate_arrays=True` flattens array props into filterable scalar columns (#269).""" + graph = RustWorkXGraph() + + labels = np.array([[[1, 1, 0], [1, 0, 2], [0, 2, 2]]], dtype=np.int32) + intensity = np.zeros((1, 3, 3, 2), dtype=np.float32) + intensity[..., 0] = [[10, 20, 0], [30, 0, 40], [0, 50, 60]] + intensity[..., 1] = [[1, 2, 0], [3, 0, 4], [0, 5, 6]] + + operator = RegionPropsNodes(extra_properties=["intensity_max", "inertia_tensor"], separate_arrays=True) + operator.add_nodes(graph, labels=labels, intensity_image=intensity) + + nodes_df = graph.node_attrs() + # 1D property -> single index suffix; 2D property -> row-major double index suffix + for col in ["intensity_max_0", "intensity_max_1", "inertia_tensor_0_0", "inertia_tensor_1_1"]: + assert col in nodes_df.columns + assert nodes_df[col].dtype == pl.Float64 + + # the array column itself must not exist when separated + assert "intensity_max" not in nodes_df.columns + + # separated columns are now filterable + subgraph = graph.filter(NodeAttr("intensity_max_0") > 30) + filtered = subgraph.node_attrs() + assert len(filtered) == 1 + assert filtered["intensity_max_0"][0] == 60.0 + + +def _multichannel_data() -> tuple[np.ndarray, np.ndarray]: + labels = np.array([[[1, 1, 0], [1, 0, 2], [0, 2, 2]]], dtype=np.int32) + intensity = np.zeros((1, 3, 3, 2), dtype=np.float32) + intensity[..., 0] = [[10, 20, 0], [30, 0, 40], [0, 50, 60]] + intensity[..., 1] = [[1, 2, 0], [3, 0, 4], [0, 5, 6]] + return labels, intensity + + +def test_regionprops_channel_names() -> None: + """`channel_names` splits intensity props into one attribute per named channel.""" + graph = RustWorkXGraph() + labels, intensity = _multichannel_data() + + operator = RegionPropsNodes(extra_properties=["intensity_max", "area"], channel_names=["dapi", "gfp"]) + operator.add_nodes(graph, labels=labels, intensity_image=intensity) + + nodes_df = graph.node_attrs() + # the packed array column is replaced by one scalar column per channel + assert "intensity_max" not in nodes_df.columns + assert nodes_df["intensity_max_dapi"].dtype == pl.Float64 + assert nodes_df["intensity_max_gfp"].dtype == pl.Float64 + # non-intensity properties are untouched + assert "area" in nodes_df.columns + + assert sorted(nodes_df["intensity_max_dapi"]) == [30.0, 60.0] + assert sorted(nodes_df["intensity_max_gfp"]) == [3.0, 6.0] + + # per-channel columns are filterable + filtered = graph.filter(NodeAttr("intensity_max_gfp") > 4).node_attrs() + assert len(filtered) == 1 + assert filtered["intensity_max_dapi"][0] == 60.0 + + +def test_regionprops_channel_names_multi_axis_property() -> None: + """Channel splitting keeps the remaining axes of multi-axis intensity props as arrays.""" + graph = RustWorkXGraph() + labels, intensity = _multichannel_data() + + operator = RegionPropsNodes(extra_properties=["centroid_weighted"], channel_names=["dapi", "gfp"]) + operator.add_nodes(graph, labels=labels, intensity_image=intensity) + + nodes_df = graph.node_attrs() + for name in ["dapi", "gfp"]: + col = f"centroid_weighted_{name}" + assert isinstance(nodes_df.schema[col], pl.Array) + # (y, x) is kept as a fixed-shape array + assert nodes_df.schema[col].shape == (2,) + + +def test_regionprops_channel_names_with_separate_arrays() -> None: + """`separate_arrays` further splits the per-channel values into scalars.""" + graph = RustWorkXGraph() + labels, intensity = _multichannel_data() + + operator = RegionPropsNodes( + extra_properties=["intensity_max", "centroid_weighted"], + channel_names=["dapi", "gfp"], + separate_arrays=True, + ) + operator.add_nodes(graph, labels=labels, intensity_image=intensity) + + nodes_df = graph.node_attrs() + for col in [ + "intensity_max_dapi", + "intensity_max_gfp", + "centroid_weighted_dapi_0", + "centroid_weighted_dapi_1", + "centroid_weighted_gfp_0", + "centroid_weighted_gfp_1", + ]: + assert col in nodes_df.columns + assert nodes_df[col].dtype == pl.Float64 + + +def test_regionprops_channel_names_custom_property_not_split() -> None: + """Custom callables keep their own layout, they are not split per channel.""" + graph = RustWorkXGraph() + labels, intensity = _multichannel_data() + + def max_per_channel(region: RegionProperties) -> np.ndarray: + return np.asarray(region.intensity_max) + + operator = RegionPropsNodes(extra_properties=[max_per_channel], channel_names=["dapi", "gfp"]) + operator.add_nodes(graph, labels=labels, intensity_image=intensity) + + nodes_df = graph.node_attrs() + assert isinstance(nodes_df.schema["max_per_channel"], pl.Array) + assert nodes_df.schema["max_per_channel"].shape == (2,) + + +def test_regionprops_channel_names_attr_keys() -> None: + """attr_keys lists intensity properties once per channel.""" + operator = RegionPropsNodes(extra_properties=["area", "intensity_mean"], channel_names=["dapi", "gfp"]) + assert operator.attr_keys() == ["area", "intensity_mean_dapi", "intensity_mean_gfp"] + + +@pytest.mark.parametrize( + "channel_names,match", + [ + ([], "must not be empty"), + (["dapi", "dapi"], "must be unique"), + (["dapi", 0], "must be strings"), + ], +) +def test_regionprops_channel_names_invalid(channel_names: list, match: str) -> None: + """Invalid `channel_names` are rejected at construction time.""" + with pytest.raises(ValueError, match=match): + RegionPropsNodes(channel_names=channel_names) + + +def test_regionprops_channel_names_shape_mismatch() -> None: + """The intensity image must have a trailing axis matching `channel_names`.""" + graph = RustWorkXGraph() + labels, intensity = _multichannel_data() + + operator = RegionPropsNodes(extra_properties=["intensity_max"], channel_names=["dapi", "gfp", "rfp"]) + with pytest.raises(ValueError, match="must have shape"): + operator.add_nodes(graph, labels=labels, intensity_image=intensity) + + # single-channel intensity image, no trailing channel axis + with pytest.raises(ValueError, match="must have shape"): + operator.add_nodes(graph, labels=labels, intensity_image=intensity[..., 0]) + + with pytest.raises(ValueError, match="`intensity_image` is None"): + operator.add_nodes(graph, labels=labels) + + +def test_regionprops_separate_arrays_matches_unpack() -> None: + """`separate_arrays=True` column names match `node_attrs(unpack=True)`.""" + labels = np.array([[[1, 1, 0], [1, 0, 2], [0, 2, 2]]], dtype=np.int32) + intensity = np.zeros((1, 3, 3, 2), dtype=np.float32) + intensity[..., 0] = [[10, 20, 0], [30, 0, 40], [0, 50, 60]] + intensity[..., 1] = [[1, 2, 0], [3, 0, 4], [0, 5, 6]] + + extra = ["intensity_max", "inertia_tensor"] + + sep_graph = RustWorkXGraph() + RegionPropsNodes(extra_properties=extra, separate_arrays=True).add_nodes( + sep_graph, labels=labels, intensity_image=intensity + ) + + packed_graph = RustWorkXGraph() + RegionPropsNodes(extra_properties=extra).add_nodes(packed_graph, labels=labels, intensity_image=intensity) + + def _prop_cols(df: pl.DataFrame) -> set[str]: + return {c for c in df.columns if c.startswith(("intensity_max", "inertia_tensor"))} + + assert _prop_cols(sep_graph.node_attrs()) == _prop_cols(packed_graph.node_attrs(unpack=True))