Add position scale

src/funtracks/annotators/_regionprops_annotator.py

@@ -28,6 +28,7 @@
 
     from funtracks.actions import BasicAction
     from funtracks.data_model import Tracks
+    from funtracks.data_model.tracks import PositionUnits
 
     IntensityImage: TypeAlias = np.ndarray | da.Array
 
@@ -40,26 +41,27 @@
 
 
 def _centroid(mask: Mask, spacing: tuple[float, ...] | None) -> list[float]:
-    """Centroid in world units, read directly from the mask array.
+    """Centroid read directly from the mask array, in pixel or world coordinates.
 
-    Equivalent to ``ExtendedRegionProperties.centroid`` (``(local_centroid + bbox_min)
-    * spacing``) but skips the skimage regionprops machinery (find_objects, region
-    caching), which is wasted overhead when only the centroid is needed.
+    Equivalent to ``ExtendedRegionProperties.centroid_pixel``/``.centroid`` but
+    skips the skimage regionprops machinery (find_objects, region caching), which
+    is wasted overhead when only the centroid is needed.
 
     Args:
         mask: A Mask object representing one detection.
-        spacing: Voxel spacing per spatial dimension, or None for unit spacing.
+        spacing: Voxel spacing per spatial dimension to scale the centroid by, or
+            None to return it in pixel coordinates.
 
     Returns:
         The centroid coordinates, one float per spatial dimension.
     """
     arr = mask.mask
     bbox_min = mask.bbox[: arr.ndim]
     local = np.array([idx.mean() for idx in np.nonzero(arr)])
-    world = local + bbox_min
+    pixel = local + bbox_min
     if spacing is not None:
-        world = world * np.asarray(spacing)
-    return [float(v) for v in world]
+        pixel = pixel * np.asarray(spacing)
+    return [float(v) for v in pixel]
 
 
 def _bbox_slicing(mask: Mask) -> tuple[slice, ...]:
@@ -146,6 +148,12 @@ class RegionpropsAnnotator(GraphAnnotator):
     - perimeter/surface area
     - mean intensity (requires an intensity image, see ``set_intensity_images``)
 
+    The centroid is stored in pixel or world coordinates according to
+    ``tracks.position_units`` (pixel by default, matching how the segmentation
+    itself is indexed). All size/shape measurements are always computed with
+    ``tracks.scale`` as the voxel spacing and are therefore in world units,
+    regardless of ``position_units``.
+
     Defaults to computing all features, but individual ones can be turned off by changing
     the self.include value at the corresponding index to the feature in self.features.
     """
@@ -211,6 +219,7 @@ def __init__(
         specs = RegionpropsAnnotator._define_features(
             tracks.ndim,
             self.channel_names,
+            position_units=tracks.position_units,
         )
         # Replace the default position spec with the requested key(s): one spec
         # holding the whole coordinate, or one scalar spec per axis. Keyed off the
@@ -355,6 +364,7 @@ def _define_features(
         cls,
         ndim: int,
         channel_names: Sequence[str] | None = None,
+        position_units: PositionUnits = "pixel",
     ) -> list[FeatureSpec]:
         """Define all supported regionprops features along with keys and function names.
 
@@ -365,6 +375,10 @@ def _define_features(
             ndim: Total number of dimensions including time (3 or 4)
             channel_names: Display names for the intensity channels, one per channel.
                 Controls how many values the intensity feature holds.
+            position_units: Whether the position feature should be computed (and
+                recomputed) in pixel or world coordinates. World positions are
+                ``scale_dependent``: they need re-deriving whenever ``tracks.scale``
+                changes, the same as area/perimeter/etc. Pixel positions do not.
 
         Returns:
             list[FeatureSpec]: List of feature specifications with key, feature,
@@ -374,8 +388,16 @@ def _define_features(
         # Default to 3D when ndim is None to enable matching all position columns
         axis_names = ["z", "y", "x"] if ndim is None or ndim == 4 else ["y", "x"]
 
+        position_feature = Position(axes=axis_names)
+        if position_units == "world":
+            position_feature["scale_dependent"] = True
+
         return [
-            FeatureSpec(DEFAULT_POS_KEY, Position(axes=axis_names), "centroid"),
+            FeatureSpec(
+                DEFAULT_POS_KEY,
+                position_feature,
+                "centroid" if position_units == "world" else "centroid_pixel",
+            ),
             FeatureSpec(DEFAULT_AREA_KEY, Area(ndim=ndim), "area"),
             FeatureSpec(
                 DEFAULT_INTENSITY_KEY, Intensity(channel_names), "intensity_mean"
@@ -468,6 +490,9 @@ def compute(self, feature_keys: list[str] | None = None) -> None:
         # regionprops pass anyway and its centroid comes for free, so pos goes through
         # the normal path with everything else.
         fast_pos = set(keys_to_compute) == set(self.pos_keys)
+        # Only scale the fast-path centroid when positions are stored in world units;
+        # centroid_pixel is the pixel-coordinate path and must stay unscaled.
+        pos_spacing = spacing if self.tracks.position_units == "world" else None
 
         node_ids = [
             node_id for node_id in self.graph.node_ids() if self.graph.has_node(node_id)
@@ -491,7 +516,7 @@ def compute(self, feature_keys: list[str] | None = None) -> None:
             mask = self.graph.nodes[node_id]["mask"]
             all_node_ids.append(node_id)
             if fast_pos:
-                centroid = _centroid(mask, spacing)
+                centroid = _centroid(mask, pos_spacing)
                 for key in keys_to_compute:
                     all_values[key].append(self._pos_value(key, centroid))
                 continue

src/funtracks/annotators/_regionprops_extended.py

@@ -39,6 +39,16 @@ def centroid(self):
             )
         return raw
 
+    @property
+    def centroid_pixel(self):
+        """Global centroid in pixel coordinates, ignoring the spacing."""
+        local = np.array(
+            [idx.mean() for idx in np.nonzero(self._label_image == self.label)]
+        )
+        if self._bbox_min is not None:
+            local = local + self._bbox_min
+        return tuple(float(v) for v in local)
+
     @property
     def axes(self):
         """
@@ -233,7 +243,7 @@ def regionprops_extended(
 
     region = mask.regionprops(spacing=spacing, intensity_image=intensity_image)
 
-    bbox_min = np.array(mask.bbox[: mask._mask.ndim]) if spacing is not None else None
+    bbox_min = np.array(mask.bbox[: mask._mask.ndim])
 
     extended_region = ExtendedRegionProperties(
         slice=region.slice,

src/funtracks/candidate_graph/compute_graph.py

@@ -18,19 +18,20 @@ def compute_graph_from_seg(
     backend: str = "memory",
 ) -> td.graph.BaseGraph:
     """Construct a candidate graph from a segmentation array. Nodes are placed at the
-    centroid of each segmentation and edges are added for all nodes in adjacent frames
-    within max_edge_distance.
+    centroid of each segmentation (in pixel coordinates) and edges are added for all
+    nodes in adjacent frames within max_edge_distance.
 
     Args:
         segmentation (np.ndarray): A numpy array with integer labels and dimensions
             (t, [z], y, x).
-        max_edge_distance (float): Maximum distance that objects can travel between
-            frames. All nodes with centroids within this distance in adjacent frames
-            will by connected with a candidate edge.
+        max_edge_distance (float): Maximum distance in world units that objects can
+            travel between frames. All nodes with centroids within this distance in
+            adjacent frames will by connected with a candidate edge.
         iou (bool, optional): Whether to include IOU on the candidate graph.
             Defaults to False.
         scale (list[float] | None, optional): The scale of the segmentation data.
-            Will be used to rescale the point locations and attribute computations.
+            Node positions stay in pixel coordinates; the scale is used for
+            measurements in world units (area) and for the distance computations.
             Defaults to None, which implies the data is isotropic.
         t_start (int, optional): The time value to assign to the first frame of the
             segmentation. Frame i will get t = t_start + i. Useful when the
@@ -61,6 +62,7 @@ def compute_graph_from_seg(
         max_edge_distance=max_edge_distance,
         node_frame_dict=node_frame_dict,
         iou_dict=iou_dict,
+        scale=scale,
     )
 
     logger.info("Candidate edges: %d", cand_graph.num_edges())
@@ -83,28 +85,28 @@ def compute_graph_from_points_list(
 
     Args:
         points_list (np.ndarray): An NxD numpy array with N points and D
-            (3 or 4) dimensions. Dimensions should be in order  (t, [z], y, x).
-        max_edge_distance (float): Maximum distance that objects can travel between
-            frames. All nodes with centroids within this distance in adjacent frames
-            will by connected with a candidate edge.
-        scale (list[float] | None, optional): Amount to scale the points in each
-            dimension. Only needed if the provided points are in "voxel" coordinates
-            instead of world coordinates. Defaults to None, which implies the data is
-            isotropic.
+            (3 or 4) dimensions, in pixel coordinates. Dimensions should be in
+            order (t, [z], y, x).
+        max_edge_distance (float): Maximum distance in world units that objects can
+            travel between frames. All nodes with centroids within this distance in
+            adjacent frames will by connected with a candidate edge.
+        scale (list[float] | None, optional): The scale of the data in each dimension
+            (including time). The points are stored unchanged (in pixel coordinates);
+            the scale is only used to measure distances in world units. Defaults to
+            None, which implies the data is isotropic.
         backend (str, optional): Graph backend, "memory" or "sql". Defaults to "memory".
 
     Returns:
         td.graph.BaseGraph: A candidate graph that can be passed to the motile solver.
     """
     # add nodes
-    cand_graph, node_frame_dict = nodes_from_points_list(
-        points_list, scale=scale, backend=backend
-    )
+    cand_graph, node_frame_dict = nodes_from_points_list(points_list, backend=backend)
     logger.info("Candidate nodes: %d", cand_graph.num_nodes())
     # add edges
     add_cand_edges(
         cand_graph,
         max_edge_distance=max_edge_distance,
         node_frame_dict=node_frame_dict,
+        scale=scale,
     )
     return cand_graph

src/funtracks/candidate_graph/utils.py

@@ -27,8 +27,8 @@ def nodes_from_segmentation(
 
     Each node will have the following attributes:
         - t
-        - pos
-        - area
+        - pos (in pixel coordinates)
+        - area (in world units, i.e. scaled by the voxel size)
         - mask (if mask=True): cropped boolean mask (tracksdata Mask object)
         - bbox (if mask=True): bounding box as [min_0, ..., max_0, ...] int array
 
@@ -39,7 +39,8 @@ def nodes_from_segmentation(
             with funtracks.utils.ensure_unique_labels before calling this function.
         scale (list[float] | None, optional): The scale of the segmentation data in all
             dimensions (including time, which should have a dummy 1 value).
-            Will be used to rescale the point locations and attribute computations.
+            Positions and masks stay in pixel coordinates; the scale is only applied
+            to measurements that are expressed in world units (here: area).
             Defaults to None, which implies the data is isotropic.
         mask (bool, optional): Whether to include mask and bbox attributes for each
             node. Uses regionprop.image (already computed by regionprops at no extra
@@ -67,6 +68,9 @@ def nodes_from_segmentation(
             f"Scale {scale} should have {segmentation.ndim} dims"
         )
 
+    # Voxel volume, used to express the pixel-count area in world units.
+    voxel_volume = float(np.prod(scale[1:]))
+
     node_attributes = ["pos", "area", "mask", "bbox"] if mask else ["pos", "area"]
     cand_graph = create_empty_graph(
         node_attributes=node_attributes,
@@ -83,7 +87,9 @@ def nodes_from_segmentation(
     for t in tqdm(range(len(segmentation)), desc="Extracting nodes from segmentation"):
         segs = segmentation[t]
         nodes_in_frame: list[int] = []
-        props = regionprops(segs, spacing=tuple(scale[1:]))
+        # No spacing: the centroid must stay in pixel coordinates. The area is converted
+        # to world coordinates
+        props = regionprops(segs)
         for regionprop in props:
             node_id = regionprop.label
             if node_id in seen_ids:
@@ -92,7 +98,7 @@ def nodes_from_segmentation(
             attrs: dict[str, Any] = {
                 "t": t + t_start,
                 "pos": list(regionprop.centroid),
-                "area": float(regionprop.area),
+                "area": float(regionprop.area) * voxel_volume,
             }
             if mask:
                 # regionprop.image is the cropped boolean mask within the bbox,
@@ -117,21 +123,20 @@ def nodes_from_segmentation(
 
 def nodes_from_points_list(
     points_list: np.ndarray,
-    scale: list[float] | None = None,
     backend: str = "memory",
 ) -> tuple[td.graph.BaseGraph, dict[int, list[Any]]]:
     """Extract candidate nodes from a list of points. Uses the index of the
     point in the list as its unique id.
     Returns a tracksdata graph with only nodes, and also a dictionary from frames to
     node_ids for efficient edge adding.
 
+    The points are stored as provided: funtracks keeps positions in pixel
+    coordinates, so the caller should pass pixel coordinates here.
+
     Args:
         points_list (np.ndarray): An NxD numpy array with N points and D
-            (3 or 4) dimensions. Dimensions should be in order (t, [z], y, x).
-        scale (list[float] | None, optional): Amount to scale the points in each
-            dimension (including time). Only needed if the provided points are in
-            "voxel" coordinates instead of world coordinates. Defaults to None, which
-            implies the data is isotropic.
+            (3 or 4) dimensions, in pixel coordinates. Dimensions should be in
+            order (t, [z], y, x).
         backend (str, optional): Graph backend, "memory" or "sql". Defaults to "memory".
 
     Returns:
@@ -142,12 +147,6 @@ def nodes_from_points_list(
 
     ndim = points_list.shape[1]
 
-    if scale is not None:
-        assert len(scale) == ndim, (
-            f"Cannot scale points with {ndim} dims by factor {scale}"
-        )
-        points_list = points_list * np.array(scale)
-
     cand_graph = create_empty_graph(
         node_attributes=["pos"],
         position_attrs=["pos"],
@@ -193,7 +192,11 @@ def _compute_node_frame_dict(cand_graph: td.graph.BaseGraph) -> dict[int, list[A
     return node_frame_dict
 
 
-def create_kdtree(cand_graph: td.graph.BaseGraph, node_ids: list[Any]) -> KDTree:
+def create_kdtree(
+    cand_graph: td.graph.BaseGraph,
+    node_ids: list[Any],
+    spatial_scale: np.ndarray | None = None,
+) -> KDTree:
     """Create a kdtree with the given nodes from the candidate graph.
     Will fail if provided node ids are not in the candidate graph.
 
@@ -203,6 +206,10 @@ def create_kdtree(cand_graph: td.graph.BaseGraph, node_ids: list[Any]) -> KDTree
             include in the KDTree. Useful for limiting to one time frame.
             Must be a list (not a generic iterable) to preserve order for
             correct mapping of query_ball_tree results back to node ids.
+        spatial_scale (np.ndarray | None, optional): Per-spatial-dimension voxel
+            size. Positions are stored in pixel coordinates; multiplying by this
+            puts the tree in world units so distance queries are in world units
+            too. Defaults to None (positions used as-is).
 
     Returns:
         KDTree: A KDTree containing the positions of the given nodes,
@@ -213,7 +220,9 @@ def create_kdtree(cand_graph: td.graph.BaseGraph, node_ids: list[Any]) -> KDTree
     # graph attribute schema on every call, which dominates runtime at scale.
     df = cand_graph.filter(node_ids=node_ids).node_attrs(attr_keys=["node_id", "pos"])
     pos_by_id = dict(zip(df["node_id"].to_list(), df["pos"].to_list(), strict=True))
-    positions = [pos_by_id[nid] for nid in node_ids]
+    positions = np.asarray([pos_by_id[nid] for nid in node_ids], dtype=float)
+    if spatial_scale is not None:
+        positions = positions * spatial_scale
     return KDTree(positions)
 
 
@@ -222,6 +231,7 @@ def add_cand_edges(
     max_edge_distance: float,
     node_frame_dict: None | dict[int, list[Any]] = None,
     iou_dict: dict[int, dict[int, float]] | None = None,
+    scale: list[float] | None = None,
 ) -> None:
     """Add candidate edges to a candidate graph by connecting all nodes in adjacent
     frames that are closer than max_edge_distance.
@@ -230,8 +240,8 @@ def add_cand_edges(
         cand_graph (td.graph.BaseGraph): Candidate graph with only nodes populated.
             Will be modified in-place to add edges.
         max_edge_distance (float): Maximum distance that objects can travel between
-            frames. All nodes within this distance in adjacent frames will by connected
-            with a candidate edge.
+            frames, in world units. All nodes within this distance in adjacent frames
+            will by connected with a candidate edge.
         node_frame_dict (dict[int, list[Any]] | None, optional): A mapping from frames
             to node ids. If not provided, it will be computed from cand_graph. Defaults
             to None.
@@ -240,6 +250,10 @@ def add_cand_edges(
             provided, the iou value is included directly in each edge dict passed to
             bulk_add_edges, avoiding a separate per-edge update pass after insertion.
             Defaults to None (no iou attribute added).
+        scale (list[float] | None, optional): The scale of the data in all dimensions
+            (including time). Node positions are in pixel coordinates, so they are
+            multiplied by the spatial part of the scale before measuring distances,
+            keeping max_edge_distance in world units. Defaults to None (isotropic).
     """
     logger.info("Extracting candidate edges")
     if not node_frame_dict:
@@ -248,14 +262,16 @@ def add_cand_edges(
     if iou_dict is not None:
         cand_graph.add_edge_attr_key("iou", default_value=0.0, dtype=pl.Float64)
 
+    spatial_scale = None if scale is None else np.asarray(scale[1:], dtype=float)
+
     frames = sorted(node_frame_dict.keys())
     prev_node_ids = node_frame_dict[frames[0]]
-    prev_kdtree = create_kdtree(cand_graph, prev_node_ids)
+    prev_kdtree = create_kdtree(cand_graph, prev_node_ids, spatial_scale)
     for frame in tqdm(frames, desc="Adding candidate edges"):
         if frame + 1 not in node_frame_dict:
             continue
         next_node_ids = node_frame_dict[frame + 1]
-        next_kdtree = create_kdtree(cand_graph, next_node_ids)
+        next_kdtree = create_kdtree(cand_graph, next_node_ids, spatial_scale)
 
         matched_indices = prev_kdtree.query_ball_tree(next_kdtree, max_edge_distance)