Keep coordinates in pixel units

src/funtracks/annotators/_regionprops_annotator.py

@@ -31,27 +31,23 @@
 DEFAULT_PERIMETER_KEY = "perimeter"
 
 
-def _centroid(mask: Mask, spacing: tuple[float, ...] | None) -> list[float]:
-    """Centroid in world units, read directly from the mask array.
+def _centroid(mask: Mask) -> list[float]:
+    """Centroid in pixel coordinates, read directly from the mask array.
 
-    Equivalent to ``ExtendedRegionProperties.centroid`` (``(local_centroid + bbox_min)
-    * spacing``) but skips the skimage regionprops machinery (find_objects, region
+    Equivalent to ``ExtendedRegionProperties.centroid_pixel`` (``local_centroid +
+    bbox_min``) 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.
 
     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
-    if spacing is not None:
-        world = world * np.asarray(spacing)
-    return [float(v) for v in world]
+    return [float(v) for v in local + bbox_min]
 
 
 class FeatureSpec(NamedTuple):
@@ -78,6 +74,10 @@ class RegionpropsAnnotator(GraphAnnotator):
     - circularity/sphericity
     - perimeter/surface area
 
+    The centroid is stored in pixel coordinates, matching how the segmentation
+    itself is indexed. All size/shape measurements are computed with
+    ``tracks.scale`` as the voxel spacing and are therefore in world 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.
     """
@@ -148,7 +148,7 @@ def _define_features(
         axis_names = ["z", "y", "x"] if ndim is None or ndim == 4 else ["y", "x"]
 
         return [
-            FeatureSpec(DEFAULT_POS_KEY, Position(axes=axis_names), "centroid"),
+            FeatureSpec(DEFAULT_POS_KEY, Position(axes=axis_names), "centroid_pixel"),
             FeatureSpec(DEFAULT_AREA_KEY, Area(ndim=ndim), "area"),
             # TODO: Add in intensity when image is passed
             # FeatureSpec("intensity", Intensity(ndim=ndim), "intensity"),
@@ -207,7 +207,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:
-                all_values[self.pos_key].append(_centroid(mask, spacing))
+                all_values[self.pos_key].append(_centroid(mask))
                 continue
             (region,) = regionprops_extended(mask, spacing=spacing)
             for key in keys_to_compute:

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):
         """
@@ -228,7 +238,7 @@ def regionprops_extended(
 
     region = mask.regionprops(spacing=spacing)
 
-    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

@@ -17,19 +17,20 @@ def compute_graph_from_seg(
     t_start: int = 0,
 ) -> 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
@@ -59,6 +60,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())
@@ -80,25 +82,27 @@ 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.
 
     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)
+    cand_graph, node_frame_dict = nodes_from_points_list(points_list)
     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

@@ -26,8 +26,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
 
@@ -38,7 +38,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
@@ -65,6 +66,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,
@@ -80,7 +84,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:
@@ -89,7 +95,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,
@@ -114,20 +120,19 @@ def nodes_from_segmentation(
 
 def nodes_from_points_list(
     points_list: np.ndarray,
-    scale: list[float] | None = None,
 ) -> 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).
 
     Returns:
         tuple[td.graph.BaseGraph, dict[int, list[Any]]]: A candidate graph with only
@@ -137,12 +142,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"],
@@ -187,7 +186,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.
 
@@ -197,6 +200,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,
@@ -207,7 +214,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)
 
 
@@ -216,6 +225,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.
@@ -224,8 +234,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.
@@ -234,6 +244,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:
@@ -242,14 +256,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)