From 88eb400d8ded6e1c88b88c190d7ceb797b35d27a Mon Sep 17 00:00:00 2001 From: AnniekStok Date: Fri, 2 Oct 2026 14:55:04 +0200 Subject: [PATCH] add point_intensity_annotator for measuring intensity in disk/ball regions when only points are available --- src/funtracks/annotators/__init__.py | 2 + src/funtracks/annotators/_intensity_images.py | 257 +++++++++++++++++ .../annotators/_point_intensity_annotator.py | 268 ++++++++++++++++++ .../annotators/_regionprops_annotator.py | 207 +------------- src/funtracks/data_model/tracks.py | 76 ++++- .../test_point_intensity_annotator.py | 165 +++++++++++ .../annotators/test_regionprops_annotator.py | 6 +- 7 files changed, 781 insertions(+), 200 deletions(-) create mode 100644 src/funtracks/annotators/_intensity_images.py create mode 100644 src/funtracks/annotators/_point_intensity_annotator.py create mode 100644 tests/annotators/test_point_intensity_annotator.py diff --git a/src/funtracks/annotators/__init__.py b/src/funtracks/annotators/__init__.py index f605dbde..229f9a3c 100644 --- a/src/funtracks/annotators/__init__.py +++ b/src/funtracks/annotators/__init__.py @@ -1,6 +1,7 @@ from ._annotator_registry import AnnotatorRegistry from ._edge_annotator import EdgeAnnotator from ._graph_annotator import GraphAnnotator +from ._point_intensity_annotator import PointIntensityAnnotator from ._regionprops_annotator import RegionpropsAnnotator from ._track_annotator import TrackAnnotator @@ -8,6 +9,7 @@ "AnnotatorRegistry", "EdgeAnnotator", "GraphAnnotator", + "PointIntensityAnnotator", "RegionpropsAnnotator", "TrackAnnotator", ] diff --git a/src/funtracks/annotators/_intensity_images.py b/src/funtracks/annotators/_intensity_images.py new file mode 100644 index 00000000..2de709ae --- /dev/null +++ b/src/funtracks/annotators/_intensity_images.py @@ -0,0 +1,257 @@ +"""Intensity image handling shared by the annotators that measure mean intensity. + +Mean intensity is measured in a region per node: the segmentation mask when there is +one (``RegionpropsAnnotator``), or a disk/sphere around the point otherwise +(``PointIntensityAnnotator``). In both cases the region is converted to a Mask. All +functions related to validating, reading, and stacking the intensity images lives here. +""" + +from __future__ import annotations + +from collections.abc import Sequence +from typing import TYPE_CHECKING, TypeAlias + +import numpy as np +from tracksdata.nodes import Mask + +from funtracks.features import Intensity + +if TYPE_CHECKING: + import dask.array as da + + from funtracks.data_model import Tracks + from funtracks.features import Feature + + IntensityImage: TypeAlias = np.ndarray | da.Array + + +def _bbox_slicing(mask: Mask) -> tuple[slice, ...]: + """The spatial slices covering a mask's bounding box, one per spatial axis.""" + ndim = mask.mask.ndim + bbox = mask.bbox + return tuple(slice(bbox[i], bbox[i + ndim]) for i in range(ndim)) + + +def _as_intensity_image(crops: list[np.ndarray]) -> np.ndarray: + """Combine one crop per channel into the intensity image skimage expects. + + Several channels are stacked on a trailing axis, which skimage reads as a + multichannel intensity image and answers with one mean per channel. + """ + return crops[0] if len(crops) == 1 else np.stack(crops, axis=-1) + + +class _FrameCache: + """Serves bounding box crops for many nodes out of one materialized time point. + + ``compute`` walks nodes in time order, so holding the current frame lets every node + in it be cropped from memory. + + Only the current time point is held (one frame per channel), and the cache is local + to a single ``compute`` call, so nothing is retained afterwards. + """ + + def __init__(self, intensity_images: list[IntensityImage] | None): + self._images = intensity_images + self._time: int | None = None + self._frames: list[np.ndarray] = [] + + def crop(self, mask: Mask, time: int | None) -> np.ndarray | None: + """The intensity image for one mask, read from the cached time point.""" + if self._images is None or time is None: + return None + if time != self._time: + self._frames = [np.asarray(image[time]) for image in self._images] + self._time = time + slicing = _bbox_slicing(mask) + return _as_intensity_image([frame[slicing] for frame in self._frames]) + + +def _intensity_crop( + intensity_images: list[IntensityImage] | None, mask: Mask, time: int | None +) -> np.ndarray | None: + """Crop the intensity image(s) of one time point to a mask's bounding box. + + For updating a single node, where caching a whole frame would be wasted. + + Args: + intensity_images: The images to crop, one per channel. + mask: The mask defining the bounding box to crop to. + time: The time point to take the intensity frame from. + + Returns: + The cropped intensity image, shaped like the mask with a trailing channel + axis when there is more than one channel, or None if no image is set. + """ + if intensity_images is None or time is None: + return None + # Slice the time point and the bounding box in a single indexing operation, so + # a store that supports it fetches only the box instead of the whole frame. + slicing = (time, *_bbox_slicing(mask)) + return _as_intensity_image([np.asarray(image[slicing]) for image in intensity_images]) + + +def _times_in_order(tracks: Tracks, node_ids: list[int]) -> tuple[list[int], list[int]]: + """The nodes sorted by time, with their times, so each frame is read only once. + + The times are fetched in one bulk query rather than one graph lookup per node. + """ + in_time_order = sorted( + zip(node_ids, tracks.get_times(node_ids), strict=True), + key=lambda pair: pair[1], + ) + return ( + [node_id for node_id, _ in in_time_order], + [time for _, time in in_time_order], + ) + + +class IntensityImagesMixin: + """Holds the raw images an annotator measures mean intensity in. + + Expects the host annotator to provide ``tracks``, ``all_features``, + ``intensity_key`` and ``compute``, as ``GraphAnnotator`` subclasses do. + """ + + tracks: Tracks + all_features: dict[str, tuple[Feature, bool]] + intensity_key: str + intensity_images: list[IntensityImage] | None + channel_names: list[str] | None + + if TYPE_CHECKING: + # Provided by the annotator; declared only for type checking, so that it + # cannot shadow the annotator's own compute in the method resolution order. + def compute(self, feature_keys: list[str] | None = None) -> None: ... + + def _validate_intensity_images( + self, + tracks: Tracks, + intensity_images: Sequence[IntensityImage] | None, + channel_names: Sequence[str] | None, + ) -> None: + """Validate and store the intensity images and channel names. + + With a segmentation, every image must match its shape. Without one, every image + must have the dimensions of the tracks (time first), and all must share a shape. + + Args: + tracks: The tracks, used to validate the image shapes. + intensity_images: raw images to measure intensity on. + channel_names: Optional display names. + + Raises: + TypeError: If a single array is passed instead of a sequence of them. + ValueError: If the image shapes do not fit the tracks or each other, or if + the number of channel names does not match the number of channels. + """ + if hasattr(intensity_images, "shape"): + raise TypeError( + "intensity_images takes one image per channel: pass [image], not a " + "bare array" + ) + if not intensity_images: + self.intensity_images = None + self.channel_names = None + return + + images = list(intensity_images) + if tracks.segmentation is not None: + seg_shape = tuple(tracks.segmentation.shape) + for image in images: + if tuple(image.shape) != seg_shape: + raise ValueError( + f"Intensity image shape {tuple(image.shape)} does not match " + f"the segmentation shape {seg_shape}" + ) + else: + first_shape = tuple(images[0].shape) + for image in images: + if len(image.shape) != tracks.ndim: + raise ValueError( + f"Intensity image shape {tuple(image.shape)} does not have " + f"the {tracks.ndim} dimensions of the tracks (t, [z], y, x)" + ) + if tuple(image.shape) != first_shape: + raise ValueError( + f"Intensity images have different shapes: " + f"{tuple(image.shape)} and {first_shape}" + ) + + num_channels = len(images) + if channel_names is None: + names = ( + None + if num_channels == 1 + else [f"channel_{i}" for i in range(num_channels)] + ) + else: + if len(channel_names) != num_channels: + raise ValueError( + f"Got {len(channel_names)} channel names for {num_channels} " + "intensity channels" + ) + names = list(channel_names) + + self.intensity_images = images + self.channel_names = names + + def _intensity_feature(self) -> Feature: + """The intensity Feature for the current channels.""" + return Intensity(self.channel_names) + + def set_intensity_images( + self, + intensity_images: Sequence[IntensityImage] | None, + channel_names: Sequence[str] | None = None, + ) -> None: + """Attach (or clear) the intensity images used to compute the intensity feature. + + ``Tracks`` builds this annotator before any raw image is known, so this is the + normal way to supply them. If the intensity feature is already enabled, it is + brought up to date here: re-registered when the number of channels changed + (the column holds one value per channel), recomputed when the images changed, + and left alone when only the channel names differ. + + Args: + intensity_images: Raw images, one per channel, each shaped like the + segmentation, or ``(t, [z], y, x)`` without one. Pass None (or an empty + list) to clear. + channel_names: Display names, one per intensity image. + """ + previous_images = self.intensity_images + previous_feature, included = self.all_features[self.intensity_key] + + self._validate_intensity_images(self.tracks, intensity_images, channel_names) + + # Rebuild the intensity Feature: its num_values follows the channel count. + feature = self._intensity_feature() + self.all_features[self.intensity_key] = (feature, included) + + if not included or self.intensity_key not in self.tracks.features: + return + + if feature["num_values"] != previous_feature["num_values"]: + # The column shape changed, so it has to be dropped and rebuilt + self.tracks.disable_features([self.intensity_key]) + self.tracks.enable_features([self.intensity_key]) + return + + self.tracks.features[self.intensity_key] = feature + if not self._is_same_image(previous_images, self.intensity_images): + self.compute([self.intensity_key]) + + @staticmethod + def _is_same_image( + previous: list[IntensityImage] | None, current: list[IntensityImage] | None + ) -> bool: + """Whether two intensity inputs are the very same images, channel for channel. + + Compared by identity: renaming a channel should not trigger a recompute, but + swapping in a different image (even an equal-looking one) should. + """ + if previous is None or current is None: + return previous is current + return len(previous) == len(current) and all( + before is after for before, after in zip(previous, current, strict=True) + ) diff --git a/src/funtracks/annotators/_point_intensity_annotator.py b/src/funtracks/annotators/_point_intensity_annotator.py new file mode 100644 index 00000000..2e66bccc --- /dev/null +++ b/src/funtracks/annotators/_point_intensity_annotator.py @@ -0,0 +1,268 @@ +from __future__ import annotations + +import warnings +from collections.abc import Sequence +from typing import TYPE_CHECKING, Any + +import numpy as np +from tracksdata.nodes import Mask + +from funtracks.actions.add_delete_node import AddNode +from funtracks.actions.update_node_attrs import UpdateNodeAttrs +from funtracks.features import Feature, Intensity + +from ._graph_annotator import GraphAnnotator +from ._intensity_images import ( + IntensityImagesMixin, + _FrameCache, + _intensity_crop, + _times_in_order, +) +from ._regionprops_annotator import DEFAULT_INTENSITY_KEY + +if TYPE_CHECKING: + from funtracks.actions import BasicAction + from funtracks.data_model import Tracks + + from ._intensity_images import IntensityImage + +DEFAULT_POINT_DIAMETER = 5.0 + + +def _ball_mask( + position: Sequence[float], + diameter: float, + spacing: np.ndarray, + frame_shape: tuple[int, ...], +) -> Mask | None: + """The disk (2D) or sphere (3D) of pixels around a point, clipped to the frame. + + A pixel belongs to the ball if its center lies within ``diameter / 2`` of the point, + measured in world units, so anisotropic spacing gives an ellipse/ellipsoid in pixel + space. When the ball is too small to contain a pixel center, the pixel nearest to + the point is used, so that every point inside the frame gets a measurement. + + Args: + position: The point in world units, one value per spatial axis. + diameter: The ball diameter in world units. + spacing: Pixel spacing per spatial axis. + frame_shape: The spatial shape of the image the ball is clipped to. + + Returns: + The ball as a Mask, or None if it does not overlap the frame at all. + """ + center = np.asarray(position, dtype=float) / spacing + shape = np.asarray(frame_shape) + radius = diameter / 2 + radii = radius / spacing + + lo = np.maximum(np.ceil(center - radii).astype(int), 0) + hi = np.minimum(np.floor(center + radii).astype(int) + 1, shape) + if np.all(hi > lo): + grids = np.ogrid[ + tuple(slice(start, stop) for start, stop in zip(lo, hi, strict=True)) + ] + dist2 = sum( + ((grid - c) * s) ** 2 + for grid, c, s in zip(grids, center, spacing, strict=True) + ) + inside = dist2 <= radius**2 + if inside.any(): + return Mask(inside, bbox=np.concatenate([lo, hi])) + + nearest = np.round(center).astype(int) + if np.any(nearest < 0) or np.any(nearest >= shape): + return None + return Mask( + np.ones((1,) * len(nearest), dtype=bool), + bbox=np.concatenate([nearest, nearest + 1]), + ) + + +def _mean_intensity(mask: Mask, intensity_image: np.ndarray) -> Any: + """The mean intensity inside a mask, one value per channel. + + Matches what ``RegionpropsAnnotator`` stores for ``intensity_mean``: a float for a + single channel, a list of floats (from the trailing channel axis) for several. + """ + mean = intensity_image[mask.mask].mean(axis=0) + if np.ndim(mean) == 0: + return float(mean) + return [float(v) for v in mean] + + +class PointIntensityAnnotator(IntensityImagesMixin, GraphAnnotator): + """Measures mean intensity around each point, for tracks without a segmentation. + + Intensity is measured in a disk (2D) or sphere (3D) of a user-chosen diameter, in + world units, centered on the node position. + """ + + @classmethod + def can_annotate(cls, tracks) -> bool: + """Point-based intensity is only for tracks without a segmentation. + + With a segmentation, ``RegionpropsAnnotator`` measures inside the mask instead. + """ + return tracks.segmentation is None + + def __init__( + self, + tracks: Tracks, + intensity_images: Sequence[IntensityImage] | None = None, + channel_names: Sequence[str] | None = None, + diameter: float = DEFAULT_POINT_DIAMETER, + ): + """ + Args: + tracks: The tracks to compute intensity for. + intensity_images: Optional raw images to measure intensity in, one per + channel, each shaped ``(t, [z], y, x)``. + channel_names: Optional display names, one per intensity image. + diameter: Diameter of the disk/sphere to measure in, in world units + """ + self.intensity_key = DEFAULT_INTENSITY_KEY + self.diameter = self._validate_diameter(diameter) + + self.intensity_images: list[IntensityImage] | None = None + self.channel_names: list[str] | None = None + self._validate_intensity_images(tracks, intensity_images, channel_names) + + super().__init__(tracks, {self.intensity_key: Intensity(self.channel_names)}) + + @classmethod + def get_available_features( + cls, ndim: int = 3, channel_names: Sequence[str] | None = None + ) -> dict[str, Feature]: + """Get all features that can be computed by this annotator. + + Args: + ndim: Total number of dimensions including time (3 or 4). Unused, for + consistency with the other annotators. + channel_names: Display names for the intensity channels. + + Returns: + Dictionary mapping feature keys to Feature definitions. + """ + return {DEFAULT_INTENSITY_KEY: Intensity(channel_names)} + + @staticmethod + def _validate_diameter(diameter: float) -> float: + diameter = float(diameter) + if not np.isfinite(diameter) or diameter <= 0: + raise ValueError(f"Diameter must be a positive number, got {diameter}") + return diameter + + def set_diameter(self, diameter: float) -> None: + """Change the disk/sphere diameter, recomputing intensity if it is enabled. + + Args: + diameter: The new diameter, in world units. + """ + diameter = self._validate_diameter(diameter) + if diameter == self.diameter: + return + self.diameter = diameter + if self.intensity_key in self.features: + self.compute([self.intensity_key]) + + @property + def _spacing(self) -> np.ndarray: + """Pixel spacing per spatial axis, to turn world positions into pixels.""" + if self.tracks.scale is None: + return np.ones(self.tracks.ndim - 1) + return np.asarray(self.tracks.scale[1:], dtype=float) + + def _images_to_measure(self) -> list[IntensityImage] | None: + """The images to measure in, or None (with a warning) if there are none.""" + if self.intensity_images is not None: + return self.intensity_images + warnings.warn( + f"Cannot compute {self.intensity_key!r}: no intensity image is set on " + "the PointIntensityAnnotator. Call set_intensity_images() first.", + stacklevel=3, + ) + return None + + def compute(self, feature_keys: list[str] | None = None) -> None: + """Compute mean intensity around every node and add it to the tracks. + + Args: + feature_keys: Optional list of specific feature keys to compute. + If None, computes all currently active features. Keys not in + self.features (not enabled) are ignored. + """ + if not self._filter_feature_keys(feature_keys): + return + images = self._images_to_measure() + if images is None: + return + + node_ids = [ + node_id for node_id in self.graph.node_ids() if self.graph.has_node(node_id) + ] + if not node_ids: + return + # Walk the nodes in time order so that each frame is read once and serves + # every node in it. + node_ids, times = _times_in_order(self.tracks, node_ids) + positions = self.tracks.get_positions(node_ids) + spacing = self._spacing + frame_shape = tuple(images[0].shape[1:]) + frames = _FrameCache(images) + + values: list[Any] = [] + for position, time in zip(positions, times, strict=True): + mask = _ball_mask(position, self.diameter, spacing, frame_shape) + if mask is None: + values.append(None) + continue + values.append(_mean_intensity(mask, frames.crop(mask, time))) + + self.tracks._set_nodes_attrs(node_ids, {self.intensity_key: values}) + + def _compute_node(self, node: int, images: list[IntensityImage]) -> None: + """Measure the intensity around a single node.""" + mask = _ball_mask( + self.tracks.get_position(node), + self.diameter, + self._spacing, + tuple(images[0].shape[1:]), + ) + value = None + if mask is not None: + time = self.tracks.get_time(node) + value = _mean_intensity(mask, _intensity_crop(images, mask, time)) + self.tracks._set_node_attr(node, self.intensity_key, value) + + def update(self, action: BasicAction) -> None: + """Remeasure a node when it is added or moved. + + Args: + action (BasicAction): The action that triggered this update + """ + moved = isinstance(action, UpdateNodeAttrs) and self._moves_node(action) + if not (isinstance(action, AddNode) or moved): + return + if not self._filter_feature_keys(None): + return + images = self._images_to_measure() + if images is not None: + self._compute_node(action.node, images) + + def _moves_node(self, action: UpdateNodeAttrs) -> bool: + """Whether an attribute update changes the node position.""" + position_key = self.tracks.features.position_key + pos_keys = [position_key] if isinstance(position_key, str) else position_key + return any(key in action.new_attrs for key in pos_keys or []) + + def change_key(self, old_key: str, new_key: str) -> None: + """Rename a feature key, keeping the intensity key in sync. + + Args: + old_key: Existing key to rename. + new_key: New key to replace it with. + """ + super().change_key(old_key, new_key) + if old_key == self.intensity_key: + self.intensity_key = new_key diff --git a/src/funtracks/annotators/_regionprops_annotator.py b/src/funtracks/annotators/_regionprops_annotator.py index bf1ca1ac..fbe0142f 100644 --- a/src/funtracks/annotators/_regionprops_annotator.py +++ b/src/funtracks/annotators/_regionprops_annotator.py @@ -2,7 +2,7 @@ import warnings from collections.abc import Sequence -from typing import TYPE_CHECKING, Any, NamedTuple, TypeAlias +from typing import TYPE_CHECKING, Any, NamedTuple import numpy as np from tracksdata.nodes import Mask @@ -21,15 +21,19 @@ ) from ._graph_annotator import GraphAnnotator +from ._intensity_images import ( + IntensityImagesMixin, + _FrameCache, + _intensity_crop, + _times_in_order, +) from ._regionprops_extended import regionprops_extended if TYPE_CHECKING: - import dask.array as da - from funtracks.actions import BasicAction from funtracks.data_model import Tracks - IntensityImage: TypeAlias = np.ndarray | da.Array + from ._intensity_images import IntensityImage DEFAULT_POS_KEY = "pos" DEFAULT_AREA_KEY = "area" @@ -62,48 +66,6 @@ def _centroid(mask: Mask, spacing: tuple[float, ...] | None) -> list[float]: return [float(v) for v in world] -def _bbox_slicing(mask: Mask) -> tuple[slice, ...]: - """The spatial slices covering a mask's bounding box, one per spatial axis.""" - ndim = mask.mask.ndim - bbox = mask.bbox - return tuple(slice(bbox[i], bbox[i + ndim]) for i in range(ndim)) - - -def _as_intensity_image(crops: list[np.ndarray]) -> np.ndarray: - """Combine one crop per channel into the intensity image skimage expects. - - Several channels are stacked on a trailing axis, which skimage reads as a - multichannel intensity image and answers with one mean per channel. - """ - return crops[0] if len(crops) == 1 else np.stack(crops, axis=-1) - - -class _FrameCache: - """Serves bounding box crops for many nodes out of one materialized time point. - - ``compute`` walks nodes in time order, so holding the current frame lets every node - in it be cropped from memory. - - Only the current time point is held (one frame per channel), and the cache is local - to a single ``compute`` call, so nothing is retained afterwards. - """ - - def __init__(self, intensity_images: list[IntensityImage] | None): - self._images = intensity_images - self._time: int | None = None - self._frames: list[np.ndarray] = [] - - def crop(self, mask: Mask, time: int | None) -> np.ndarray | None: - """The intensity image for one mask, read from the cached time point.""" - if self._images is None or time is None: - return None - if time != self._time: - self._frames = [np.asarray(image[time]) for image in self._images] - self._time = time - slicing = _bbox_slicing(mask) - return _as_intensity_image([frame[slicing] for frame in self._frames]) - - def _to_attr_value(value: Any) -> Any: """Convert a regionprops value into something the graph backend can store. @@ -135,7 +97,7 @@ class FeatureSpec(NamedTuple): regionprops_attr: str -class RegionpropsAnnotator(GraphAnnotator): +class RegionpropsAnnotator(IntensityImagesMixin, GraphAnnotator): """A graph annotator using regionprops to extract node features from segmentations. The possible features include: @@ -239,117 +201,6 @@ def __init__( # Build regionprops name mapping from specs self.regionprops_names = {spec.key: spec.regionprops_attr for spec in specs} - def _validate_intensity_images( - self, - tracks: Tracks, - intensity_images: Sequence[IntensityImage] | None, - channel_names: Sequence[str] | None, - ) -> None: - """Validate and store the intensity images and channel names. - - Args: - tracks: The tracks, used to validate the image shapes against the - segmentation. - intensity_images: raw images to measure intensity on. - channel_names: Optional display names. - - Raises: - TypeError: If a single array is passed instead of a sequence of them. - ValueError: If an image does not match the segmentation shape, or if the - number of channel names does not match the number of channels. - """ - if hasattr(intensity_images, "shape"): - raise TypeError( - "intensity_images takes one image per channel: pass [image], not a " - "bare array" - ) - if not intensity_images: - self.intensity_images = None - self.channel_names = None - return - - images = list(intensity_images) - seg_shape = tracks.segmentation.shape if tracks.segmentation is not None else None - if seg_shape is not None: - for image in images: - if tuple(image.shape) != tuple(seg_shape): - raise ValueError( - f"Intensity image shape {tuple(image.shape)} does not match " - f"the segmentation shape {tuple(seg_shape)}" - ) - - num_channels = len(images) - if channel_names is None: - names = ( - None - if num_channels == 1 - else [f"channel_{i}" for i in range(num_channels)] - ) - else: - if len(channel_names) != num_channels: - raise ValueError( - f"Got {len(channel_names)} channel names for {num_channels} " - "intensity channels" - ) - names = list(channel_names) - - self.intensity_images = images - self.channel_names = names - - def set_intensity_images( - self, - intensity_images: Sequence[IntensityImage] | None, - channel_names: Sequence[str] | None = None, - ) -> None: - """Attach (or clear) the intensity images used to compute the intensity feature. - - ``Tracks`` builds this annotator before any raw image is known, so this is the - normal way to supply them. If the intensity feature is already enabled, it is - brought up to date here: re-registered when the number of channels changed - (the column holds one value per channel), recomputed when the images changed, - and left alone when only the channel names differ. - - Args: - intensity_images: See ``__init__``. Pass None (or an empty list) to clear. - channel_names: See ``__init__``. - """ - previous_images = self.intensity_images - previous_feature, included = self.all_features[self.intensity_key] - - self._validate_intensity_images(self.tracks, intensity_images, channel_names) - - # Rebuild the intensity Feature: its num_values follows the channel count. - feature = Intensity(self.channel_names) - self.all_features[self.intensity_key] = (feature, included) - - if not included or self.intensity_key not in self.tracks.features: - return - - if feature["num_values"] != previous_feature["num_values"]: - # The column shape changed, so it has to be dropped and rebuilt - self.tracks.disable_features([self.intensity_key]) - self.tracks.enable_features([self.intensity_key]) - return - - self.tracks.features[self.intensity_key] = feature - if not self._is_same_image(previous_images, self.intensity_images): - self.compute([self.intensity_key]) - - @staticmethod - def _is_same_image( - previous: list[IntensityImage] | None, current: list[IntensityImage] | None - ) -> bool: - """Whether two intensity inputs are the very same images, channel for channel. - - Compared by identity: renaming a channel should not trigger a recompute, but - swapping in a different image (even an equal-looking one) should. - """ - if previous is None or current is None: - return previous is current - return len(previous) == len(current) and all( - before is after for before, after in zip(previous, current, strict=True) - ) - @classmethod def _define_features( cls, @@ -405,29 +256,6 @@ def get_available_features( specs = RegionpropsAnnotator._define_features(ndim, channel_names) return {spec.key: spec.feature for spec in specs} - def _intensity_crop(self, mask: Mask, time: int | None) -> np.ndarray | None: - """Crop the intensity image(s) of one time point to a mask's bounding box. - - skimage requires the intensity image to match the shape of the label image, - which in this case is the bbox-sized mask array rather than the full frame. - - Args: - mask: The mask defining the bounding box to crop to. - time: The time point to take the intensity frame from. - - Returns: - The cropped intensity image, shaped like the mask with a trailing channel - axis when there is more than one channel, or None if no image is set. - """ - if self.intensity_images is None or time is None: - return None - # Slice the time point and the bounding box in a single indexing operation, so - # a store that supports it fetches only the box instead of the whole frame. - slicing = (time, *_bbox_slicing(mask)) - return _as_intensity_image( - [np.asarray(image[slicing]) for image in self.intensity_images] - ) - def compute(self, feature_keys: list[str] | None = None) -> None: """Compute the currently included features and add them to the tracks. @@ -474,17 +302,11 @@ def compute(self, feature_keys: list[str] | None = None) -> None: ] # Times are only needed for intensity, and stay None otherwise so that the # frame cache has nothing to read. - times: list[int | None] = [None] * len(node_ids) + times: Sequence[int | None] = [None] * len(node_ids) if self.intensity_key in keys_to_compute: - # Fetch the times in one bulk query rather than one graph lookup per node, - # and walk the nodes in time order so that each frame is read once and - # serves every node in it. - in_time_order = sorted( - zip(node_ids, self.tracks.get_times(node_ids), strict=True), - key=lambda pair: pair[1], - ) - node_ids = [node_id for node_id, _ in in_time_order] - times = [time for _, time in in_time_order] + # Walk the nodes in time order so that each frame is read once and serves + # every node in it. + node_ids, times = _times_in_order(self.tracks, node_ids) frames = _FrameCache(self.intensity_images) for node_id, time in zip(node_ids, times, strict=True): @@ -568,8 +390,9 @@ def _regionprops_update( spacing = None if self.tracks.scale is None else tuple(self.tracks.scale[1:]) if time is None and self.intensity_key in feature_keys: time = self.tracks.get_time(node_id) + intensity_image = _intensity_crop(self.intensity_images, mask, time) for region in regionprops_extended( - mask, spacing=spacing, intensity_image=self._intensity_crop(mask, time) + mask, spacing=spacing, intensity_image=intensity_image ): # Skip labels that aren't nodes in the graph (e.g., unselected detections) if not self.graph.has_node(node_id): diff --git a/src/funtracks/data_model/tracks.py b/src/funtracks/data_model/tracks.py index 8c1666fc..60367018 100644 --- a/src/funtracks/data_model/tracks.py +++ b/src/funtracks/data_model/tracks.py @@ -41,6 +41,7 @@ from funtracks.annotators import ( AnnotatorRegistry, GraphAnnotator, + PointIntensityAnnotator, RegionpropsAnnotator, ) @@ -379,6 +380,7 @@ def _get_annotators(self) -> AnnotatorRegistry: Creates annotators conditionally: - RegionpropsAnnotator: Only if segmentation is provided + - PointIntensityAnnotator: Only if segmentation is NOT provided - EdgeAnnotator: Only if segmentation is provided - TrackAnnotator: Always (every Tracks has track ids) @@ -391,6 +393,7 @@ def _get_annotators(self) -> AnnotatorRegistry: from funtracks.annotators import ( AnnotatorRegistry, EdgeAnnotator, + PointIntensityAnnotator, RegionpropsAnnotator, ) @@ -411,6 +414,17 @@ def _get_annotators(self) -> AnnotatorRegistry: ) ) + # PointIntensityAnnotator: measures around the points when there is no + # segmentation to measure inside + if PointIntensityAnnotator.can_annotate(self): + annotator_list.append( + PointIntensityAnnotator( + self, + intensity_images=self._intensity_images, + channel_names=self._channel_names, + ) + ) + # EdgeAnnotator: requires segmentation if EdgeAnnotator.can_annotate(self): annotator_list.append(EdgeAnnotator(self)) @@ -938,6 +952,25 @@ def regionprops_annotator(self) -> RegionpropsAnnotator | None: return annotator return None + @property + def point_intensity_annotator(self) -> PointIntensityAnnotator | None: + """The registered PointIntensityAnnotator, or None when there is a + segmentation.""" + from funtracks.annotators import PointIntensityAnnotator + + for annotator in self.annotators: + if isinstance(annotator, PointIntensityAnnotator): + return annotator + return None + + @property + def intensity_annotator( + self, + ) -> RegionpropsAnnotator | PointIntensityAnnotator | None: + """The annotator that computes the "intensity" feature: inside the masks when + there is a segmentation, around the points otherwise.""" + return self.regionprops_annotator or self.point_intensity_annotator + def set_intensity_images( self, intensity_images: Sequence[IntensityImage] | None, @@ -949,25 +982,56 @@ def set_intensity_images( carry. Call this before enable_features(["intensity"]); if intensity is already enabled, its values are recomputed here. + With a segmentation, intensity is the mean inside each mask. Without one, it + is the mean in a disk/sphere around each point (see set_intensity_diameter). + Args: intensity_images: Raw images, one per channel, each shaped like the - segmentation (t, [z], y, x). Pass None or an empty list to clear. + segmentation (t, [z], y, x), or with the dimensions of the tracks if + there is no segmentation. Pass None or an empty list to clear. channel_names: Display names, one per intensity image. Raises: - ValueError: If there is no segmentation (and hence no RegionpropsAnnotator), - or the images do not match the segmentation shape. + ValueError: If there is no annotator to compute intensity with, or the + images do not fit the tracks. """ - annotator = self.regionprops_annotator + annotator = self.intensity_annotator if annotator is None: raise ValueError( - "Cannot set intensity images: this Tracks has no segmentation, so " - "there is no RegionpropsAnnotator to compute intensity with." + "Cannot set intensity images: there is no annotator to compute " + "intensity with." ) annotator.set_intensity_images(intensity_images, channel_names) self._intensity_images = intensity_images self._channel_names = channel_names + def set_intensity_diameter(self, diameter: float) -> None: + """Set the diameter of the disk/sphere that point intensity is measured in. + + Only applies to tracks without a segmentation, where intensity is measured + around each point. Recomputes the intensity feature if it is enabled. + + Args: + diameter: The diameter in world units (the units of the node positions). + + Raises: + ValueError: If the tracks have a segmentation (intensity is then measured + inside the masks), or the diameter is not positive. + """ + annotator = self.point_intensity_annotator + if annotator is None: + raise ValueError( + "Cannot set an intensity diameter: this Tracks has a segmentation, so " + "intensity is measured inside the masks." + ) + annotator.set_diameter(diameter) + + @property + def intensity_diameter(self) -> float | None: + """The point intensity diameter, or None when there is a segmentation.""" + annotator = self.point_intensity_annotator + return None if annotator is None else annotator.diameter + def enable_features(self, feature_keys: list[str], recompute: bool = True) -> None: """Enable multiple features for computation efficiently. diff --git a/tests/annotators/test_point_intensity_annotator.py b/tests/annotators/test_point_intensity_annotator.py new file mode 100644 index 00000000..2975a86e --- /dev/null +++ b/tests/annotators/test_point_intensity_annotator.py @@ -0,0 +1,165 @@ +import numpy as np +import pytest + +from funtracks.annotators import PointIntensityAnnotator +from funtracks.annotators._point_intensity_annotator import _ball_mask +from funtracks.data_model import Tracks +from funtracks.user_actions import UserAddNode, UserUpdateNodeAttrs + +track_attrs = {"time_attr": "t", "tracklet_attr": "track_id"} + + +def _image_shape(ndim: int) -> tuple[int, ...]: + """Image shape matching the extent of the get_graph fixture positions.""" + return (5, 100, 100) if ndim == 3 else (5, 100, 100, 100) + + +def _frame_index_image(ndim: int) -> np.ndarray: + """Intensity image whose value is the time index everywhere in the frame.""" + shape = _image_shape(ndim) + frames = np.arange(shape[0], dtype=np.float32).reshape( + (-1,) + (1,) * (len(shape) - 1) + ) + return np.broadcast_to(frames, shape) + + +def _x_gradient_image(ndim: int) -> np.ndarray: + """Intensity image whose value is the x index, so a centered ball averages to x.""" + shape = _image_shape(ndim) + return np.broadcast_to(np.arange(shape[-1], dtype=np.float32), shape) + + +class TestBallMask: + def test_disk_pixels(self): + # Pixel centers within 2.5 of (50, 50): 1 + 4 + 4 + 4 + 8 + mask = _ball_mask([50, 50], 5, np.ones(2), (100, 100)) + assert mask.mask.sum() == 21 + assert list(mask.bbox) == [48, 48, 53, 53] + + def test_sphere_pixels(self): + mask = _ball_mask([50, 50, 50], 3, np.ones(3), (100, 100, 100)) + # Radius 1.5: center, 6 face neighbours and 12 edge neighbours + assert mask.mask.sum() == 19 + + def test_anisotropic_spacing(self): + """The diameter is in world units, so coarser axes span fewer pixels.""" + mask = _ball_mask([50, 50], 9, np.array([4.0, 1.0]), (100, 100)) + rows, cols = np.nonzero(mask.mask) + assert np.ptp(rows) + 1 == 2 # rows 12 and 13 in 12.5 +- 1.125 + assert np.ptp(cols) + 1 == 9 # columns 46 to 54 in 50 +- 4.5 + + def test_clipped_to_frame(self): + mask = _ball_mask([0, 0], 5, np.ones(2), (100, 100)) + assert list(mask.bbox[:2]) == [0, 0] + assert mask.mask.sum() == 8 # the quarter disk inside the frame + + def test_tiny_diameter_uses_nearest_pixel(self): + mask = _ball_mask([10.4, 20.6], 0.1, np.ones(2), (100, 100)) + assert mask.mask.sum() == 1 + assert list(mask.bbox) == [10, 21, 11, 22] + + def test_outside_frame(self): + assert _ball_mask([-10, 50], 5, np.ones(2), (100, 100)) is None + assert _ball_mask([50, 200], 5, np.ones(2), (100, 100)) is None + + +@pytest.mark.parametrize("ndim", [3, 4]) +class TestPointIntensityAnnotator: + def test_registered_without_seg(self, get_graph, ndim): + tracks = Tracks(get_graph(ndim), ndim=ndim, **track_attrs) + assert tracks.regionprops_annotator is None + assert isinstance(tracks.point_intensity_annotator, PointIntensityAnnotator) + assert tracks.intensity_annotator is tracks.point_intensity_annotator + assert "intensity" in tracks.get_available_features() + + def test_not_registered_with_seg(self, get_graph, ndim): + tracks = Tracks(get_graph(ndim, with_seg=True), ndim=ndim, **track_attrs) + assert tracks.point_intensity_annotator is None + assert tracks.intensity_diameter is None + with pytest.raises(ValueError, match="has a segmentation"): + tracks.set_intensity_diameter(3) + + def test_compute(self, get_graph, ndim): + """Each node measures in its own frame.""" + tracks = Tracks(get_graph(ndim), ndim=ndim, **track_attrs) + tracks.set_intensity_images([_frame_index_image(ndim)]) + tracks.enable_features(["intensity"]) + + for node in tracks.graph_full.node_ids(): + assert tracks.get_node_attr(node, "intensity") == pytest.approx( + tracks.get_time(node) + ) + + def test_centered_on_position(self, get_graph, ndim): + tracks = Tracks(get_graph(ndim), ndim=ndim, **track_attrs) + tracks.set_intensity_images([_x_gradient_image(ndim)]) + tracks.set_intensity_diameter(7) + tracks.enable_features(["intensity"]) + # Node 1 sits at x=50, far from the frame edge, so the ball is symmetric + assert tracks.get_node_attr(1, "intensity") == pytest.approx(50) + + def test_scale(self, get_graph, ndim): + """Positions are in world units, so they are divided by the scale to find the + pixels to measure.""" + scale = [1.0] + [2.0] * (ndim - 1) + tracks = Tracks(get_graph(ndim), ndim=ndim, scale=scale, **track_attrs) + tracks.set_intensity_images([_x_gradient_image(ndim)]) + tracks.enable_features(["intensity"]) + # Node 1 at x=50 (world) is pixel 25 + assert tracks.get_node_attr(1, "intensity") == pytest.approx(25) + + def test_multichannel(self, get_graph, ndim): + tracks = Tracks(get_graph(ndim), ndim=ndim, **track_attrs) + tracks.set_intensity_images( + [_frame_index_image(ndim), _x_gradient_image(ndim)], + channel_names=["time", "x"], + ) + tracks.enable_features(["intensity"]) + assert tracks.features["intensity"]["num_values"] == 2 + assert list(tracks.get_node_attr(1, "intensity")) == pytest.approx([0, 50]) + + def test_set_diameter_recomputes(self, get_graph, ndim): + tracks = Tracks(get_graph(ndim), ndim=ndim, **track_attrs) + tracks.set_intensity_images([_x_gradient_image(ndim)]) + tracks.set_intensity_diameter(1) + tracks.enable_features(["intensity"]) + # Node 4 sits at x=1.5, so a larger ball is clipped by the frame edge and + # leans towards higher x + small = tracks.get_node_attr(4, "intensity") + tracks.set_intensity_diameter(9) + assert tracks.intensity_diameter == 9 + assert tracks.get_node_attr(4, "intensity") > small + + def test_invalid_diameter(self, get_graph, ndim): + tracks = Tracks(get_graph(ndim), ndim=ndim, **track_attrs) + for diameter in [0, -1, np.nan]: + with pytest.raises(ValueError, match="positive"): + tracks.set_intensity_diameter(diameter) + + def test_update_on_add_and_move(self, get_graph, ndim): + tracks = Tracks(get_graph(ndim), ndim=ndim, **track_attrs) + tracks.set_intensity_images([_x_gradient_image(ndim)]) + tracks.enable_features(["intensity"]) + + pos = [50.0] * (ndim - 3) + [50.0, 30.0] + UserAddNode(tracks, node=10, attributes={"pos": pos, "t": 3, "track_id": 10}) + assert tracks.get_node_attr(10, "intensity") == pytest.approx(30) + + new_pos = [50.0] * (ndim - 3) + [50.0, 70.0] + UserUpdateNodeAttrs(tracks, node=10, attrs={"pos": new_pos}) + assert tracks.get_node_attr(10, "intensity") == pytest.approx(70) + + def test_wrong_image_ndim(self, get_graph, ndim): + tracks = Tracks(get_graph(ndim), ndim=ndim, **track_attrs) + with pytest.raises(ValueError, match="dimensions of the tracks"): + tracks.set_intensity_images([np.zeros((5,) * (ndim + 1))]) + + def test_mismatched_channel_shapes(self, get_graph, ndim): + tracks = Tracks(get_graph(ndim), ndim=ndim, **track_attrs) + with pytest.raises(ValueError, match="different shapes"): + tracks.set_intensity_images([np.zeros((5,) * ndim), np.zeros((6,) * ndim)]) + + def test_no_image_warns(self, get_graph, ndim): + tracks = Tracks(get_graph(ndim), ndim=ndim, **track_attrs) + with pytest.warns(UserWarning, match="no intensity image"): + tracks.enable_features(["intensity"]) diff --git a/tests/annotators/test_regionprops_annotator.py b/tests/annotators/test_regionprops_annotator.py index 3a82cfe0..c8f0c04b 100644 --- a/tests/annotators/test_regionprops_annotator.py +++ b/tests/annotators/test_regionprops_annotator.py @@ -404,8 +404,10 @@ def test_set_intensity_images_without_seg(self, get_graph, ndim): graph = get_graph(ndim, with_seg=False) tracks = Tracks(graph, ndim=ndim, **track_attrs) assert tracks.regionprops_annotator is None - with pytest.raises(ValueError, match="no segmentation"): - tracks.set_intensity_images([_frame_index_image(ndim)]) + # Without a segmentation, intensity is measured around the points instead + raw = _frame_index_image(ndim) + tracks.set_intensity_images([raw]) + assert tracks.point_intensity_annotator.intensity_images[0] is raw def test_missing_seg(self, get_graph, ndim): """Test that RegionpropsAnnotator gracefully handles missing segmentation."""