From c265f85ee0330b197bb7ef756f9d53610e67f60a Mon Sep 17 00:00:00 2001
From: AnniekStok <anniek.stokkermans@gmail.com>
Date: Mon, 27 Jul 2026 17:16:07 +0200
Subject: [PATCH 1/5] scale world coordinates to pixel coordinates when
 exporting to geff, and scale back at import

---
 .../import_export/_tracks_builder.py          | 29 +++++++++
 src/funtracks/import_export/geff/_export.py   | 25 +++++++-
 src/funtracks/import_export/geff/_import.py   | 63 ++++++++++++++++++-
 tests/import_export/test_import_from_geff.py  |  9 ++-
 4 files changed, 119 insertions(+), 7 deletions(-)

diff --git a/src/funtracks/import_export/_tracks_builder.py b/src/funtracks/import_export/_tracks_builder.py
index 7ed3ce75..29d9494b 100644
--- a/src/funtracks/import_export/_tracks_builder.py
+++ b/src/funtracks/import_export/_tracks_builder.py
@@ -678,6 +678,24 @@ def enable_features(
         if static_features:
             tracks.features.update(static_features)
 
+    def _resolve_import_scale(self, scale: list[float] | None) -> list[float] | None:
+        """Return the scale to use for this import.
+
+        By default the caller-provided scale is used as-is (e.g. the scale a user
+        entered in a scale widget for a CSV import). Format-specific builders can
+        override this to source the scale elsewhere.
+        """
+        return scale
+
+    def _scale_to_world_coords(self, scale: list[float] | None) -> None:
+        """Rescale the loaded position columns in ``self.in_memory_geff``.
+
+        No-op by default. Formats that store positions in a different coordinate
+        space than funtracks' world coordinates (e.g. GEFF stores pixels) override
+        this to convert using the resolved import ``scale``.
+        """
+        return
+
     def build(
         self,
         source: Path | pd.DataFrame,
@@ -740,11 +758,22 @@ def build(
         # Validate node_name_map is complete and valid
         self.validate_name_map(has_segmentation=segmentation is not None)
 
+        # Resolve the scale to use for this import. By default the caller-provided
+        # scale is used (e.g. from a scale widget for CSV imports); format-specific
+        # builders may override this (the GEFF builder reads the authoritative
+        # scale from the geff metadata).
+        scale = self._resolve_import_scale(scale)
+
         # 1. Load source data to InMemoryGeff
         self.load_source(source, self.node_name_map)
         if self.in_memory_geff is None:
             raise ValueError("load_source() must populate self.in_memory_geff")
 
+        # Rescale loaded positions into world coordinates if the format requires
+        # it (GEFF stores pixel coordinates). Uses the same resolved scale that is
+        # attached to the resulting Tracks, so the two stay consistent.
+        self._scale_to_world_coords(scale)
+
         # 2. Combine multi-value feature columns
         self._combine_multi_value_props(
             self.in_memory_geff["node_props"], self.node_name_map
diff --git a/src/funtracks/import_export/geff/_export.py b/src/funtracks/import_export/geff/_export.py
index ad496bcd..9121db14 100644
--- a/src/funtracks/import_export/geff/_export.py
+++ b/src/funtracks/import_export/geff/_export.py
@@ -6,6 +6,7 @@
 )
 
 import geff_spec
+import numpy as np
 import polars as pl
 import tracksdata as td
 from geff_spec import GeffMetadata
@@ -170,6 +171,22 @@ def _build_geff_metadata(
             }
         )
 
+    # Convert world coordinates to pixel coordinates by dividing each spatial
+    # position attribute by its axis scale, to align with the geff format, which
+    # stores pixel coordinates together with a per-axis scale in the metadata.
+    # `graph` is a detached copy (see split_position_attr), so this does not
+    # mutate the caller's live graph.
+    pos_pixel_coords: dict[str, np.ndarray] = {}
+    node_df = graph.node_attrs()
+    for name, axis_type, axis_scale in zip(
+        axis_names, axis_types, tracks.scale, strict=True
+    ):
+        if axis_type != "space" or axis_scale == 1:
+            continue
+        pos_pixel_coords[name] = node_df[name].to_numpy() / axis_scale
+    if pos_pixel_coords:
+        graph.update_node_attrs(attrs=pos_pixel_coords, node_ids=graph.node_ids())
+
     extra: dict = {}
     if include_features:
         extra["funtracks"] = {"features": tracks.features.dump_json()}
@@ -253,7 +270,11 @@ def split_position_attr(tracks: Tracks) -> tuple[td.graph.GraphView, list[str] |
         new_graph.remove_node_attr_key(pos_key)
         return new_graph, new_keys
     elif pos_key is not None:
-        # Position is already split into separate attributes
-        return tracks.graph, list(pos_key)
+        # Position is already split into separate attributes. Detach so that
+        # downstream coordinate rescaling operates on an independent copy and
+        # does not mutate the live graph.
+        new_graph = tracks.graph.detach()
+        new_graph = new_graph.filter().subgraph()
+        return new_graph, list(pos_key)
     else:
         return tracks.graph, None
diff --git a/src/funtracks/import_export/geff/_import.py b/src/funtracks/import_export/geff/_import.py
index 060a481e..76483ea7 100644
--- a/src/funtracks/import_export/geff/_import.py
+++ b/src/funtracks/import_export/geff/_import.py
@@ -3,6 +3,7 @@
 import warnings
 from typing import TYPE_CHECKING
 
+import numpy as np
 import tracksdata as td
 from geff._typing import InMemoryGeff
 from geff.core_io._base_read import read_to_memory
@@ -166,7 +167,7 @@ def read_header(self, source_path: Path) -> None:
         self.importable_node_props = list(metadata.node_props_metadata.keys())
         self.importable_edge_props = list(metadata.edge_props_metadata.keys())
 
-        # Store axes metadata for use in infer_node_name_map
+        # Store axes metadata for use in infer_node_name_map and scale resolution
         self._geff_axes = metadata.axes or []
 
         # Read funtracks FeatureDict from GEFF extra metadata if present
@@ -313,6 +314,61 @@ def load_source(
         if self.ndim is None:
             self.ndim = ndim
 
+    def _scale_from_axes(self) -> list[float] | None:
+        """Return the per-axis scale stored in the geff metadata.
+
+        The scale is ordered as ``[time, *space]`` (mirroring
+        :meth:`infer_node_name_map`), matching the funtracks scale convention
+        ([time, z, y, x]) and the position column order. Returns None when the
+        metadata has no scaled axes (e.g. external geffs without axis scales).
+        """
+        axes = getattr(self, "_geff_axes", [])
+        time_axes = [ax for ax in axes if ax.type == "time"]
+        space_axes = [ax for ax in axes if ax.type == "space"]
+        ordered_axes = time_axes + space_axes
+        if ordered_axes and all(ax.scale is not None for ax in ordered_axes):
+            return [float(ax.scale) for ax in ordered_axes]
+        return None
+
+    def _resolve_import_scale(self, scale: list[float] | None) -> list[float] | None:
+        """For GEFF, the per-axis scale stored in the metadata is used when present, the
+        caller-provided scale is only used when the metadata has no scaled axes.
+        """
+        geff_scale = self._scale_from_axes()
+        return geff_scale if geff_scale is not None else scale
+
+    def _scale_to_world_coords(self, scale: list[float] | None) -> None:
+        """Convert stored pixel coordinates to world coordinates.
+
+        geff stores spatial coordinates in pixels together with a per-axis scale;
+        funtracks keeps positions in world coordinates, so each spatial position
+        column is multiplied by its scale (from `builder._resolve_import_scale`), so
+        the conversion stays consistent with the scale attached to the Tracks.
+        """
+        if scale is None or self.in_memory_geff is None:
+            return
+
+        spatial_scale = list(scale)[1:]  # drop the time axis
+        pos_cols = self.position_attr
+        node_props = self.in_memory_geff["node_props"]
+
+        # Composite position stored as separate columns (e.g. "z", "y", "x").
+        if (
+            pos_cols
+            and len(pos_cols) == len(spatial_scale)
+            and all(col in node_props for col in pos_cols)
+        ):
+            for col, s in zip(pos_cols, spatial_scale, strict=True):
+                if s != 1:
+                    node_props[col]["values"] = node_props[col]["values"] * s
+            return
+
+        # Position stored as a single vector column named "pos".
+        if "pos" in node_props:
+            values = node_props["pos"]["values"]
+            if values.ndim == 2 and values.shape[1] == len(spatial_scale):
+                node_props["pos"]["values"] = values * np.asarray(spatial_scale)
+
 
 def import_from_geff(
     directory: Path,
@@ -334,7 +390,10 @@ def import_from_geff(
             - For multi-value features like position, use a list: {"pos": ["y", "x"]}
             If None, property names are auto-inferred using fuzzy matching.
         segmentation_path: Optional path to segmentation data
-        scale: Optional spatial scale
+        scale: Optional scale ([time, z, y, x]). For a GEFF import the per-axis
+            scale stored in the metadata is authoritative and is always used when
+            present; this argument is only a fallback for geffs whose axes carry
+            no scale.
         edge_name_map: Optional mapping from standard funtracks keys to GEFF
             edge property names. Example: {"iou": "overlap"}
         database: Optional path to a SQLite database file for backing storage.
diff --git a/tests/import_export/test_import_from_geff.py b/tests/import_export/test_import_from_geff.py
index 5d6fdbd0..60192d21 100644
--- a/tests/import_export/test_import_from_geff.py
+++ b/tests/import_export/test_import_from_geff.py
@@ -66,12 +66,15 @@ def valid_segmentation():
     times = [1, 2, 3, 4, 5]
     x = [1.0, 0.775, 0.55, 0.325, 0.1]
     y = [100, 200, 300, 400, 500]
-    scale = [1, 1, 100]
     seg_ids = np.array([10, 20, 30, 40, 50])
 
+    # The geff stores pixel coordinates. On import each spatial coordinate is
+    # multiplied by its scale to get world coordinates, and the validation
+    # converts back to pixels via ``pixel = world / scale``, recovering the
+    # original stored coordinate. So place each seg id at the pixel equal to the
+    # stored geff coordinate.
     for t, y_val, x_f, seg_id in zip(times, y, x, seg_ids, strict=False):
-        x = int(x_f * scale[2])
-        seg[t, y_val, x] = seg_id
+        seg[t, int(y_val), int(x_f)] = seg_id
     return seg
 
 

From 3d51c4236dd9d876958b2c3b4a32befbe8f74752 Mon Sep 17 00:00:00 2001
From: AnniekStok <anniek.stokkermans@gmail.com>
Date: Fri, 21 Aug 2026 18:24:49 +0200
Subject: [PATCH 2/5] keep coordinates in pixels, export scaling in csv

---
 .../annotators/_regionprops_annotator.py      |  22 ++--
 .../annotators/_regionprops_extended.py       |  12 +-
 .../candidate_graph/compute_graph.py          |  34 ++---
 src/funtracks/candidate_graph/utils.py        |  62 ++++++----
 src/funtracks/data_model/tracks.py            |  21 +++-
 src/funtracks/features/_node_features.py      |   4 +
 .../import_export/_tracks_builder.py          |  27 ++--
 src/funtracks/import_export/_validation.py    |  15 +--
 src/funtracks/import_export/csv/_export.py    |  13 ++
 src/funtracks/import_export/geff/_export.py   |  25 +---
 src/funtracks/import_export/geff/_import.py   |  36 +-----
 tests/annotators/test_annotator_registry.py   |  37 ++++--
 .../annotators/test_regionprops_annotator.py  |  31 +++--
 tests/candidate_graph/test_compute_graph.py   |  15 ++-
 tests/candidate_graph/test_graph_utils.py     |  40 ++++--
 tests/data_model/test_solution_tracks.py      |  18 ++-
 tests/import_export/test_csv_export.py        | 117 +++++++++++++++++-
 tests/import_export/test_export_to_geff.py    |  37 ++++++
 tests/import_export/test_import_from_geff.py  |  11 +-
 19 files changed, 374 insertions(+), 203 deletions(-)

diff --git a/src/funtracks/annotators/_regionprops_annotator.py b/src/funtracks/annotators/_regionprops_annotator.py
index ef07befc..e567c266 100644
--- a/src/funtracks/annotators/_regionprops_annotator.py
+++ b/src/funtracks/annotators/_regionprops_annotator.py
@@ -31,16 +31,15 @@
 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.
@@ -48,10 +47,7 @@ def _centroid(mask: Mask, spacing: tuple[float, ...] | None) -> list[float]:
     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:
diff --git a/src/funtracks/annotators/_regionprops_extended.py b/src/funtracks/annotators/_regionprops_extended.py
index d781b727..7420f7ed 100644
--- a/src/funtracks/annotators/_regionprops_extended.py
+++ b/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,
diff --git a/src/funtracks/candidate_graph/compute_graph.py b/src/funtracks/candidate_graph/compute_graph.py
index 390a57e1..be950bed 100644
--- a/src/funtracks/candidate_graph/compute_graph.py
+++ b/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
diff --git a/src/funtracks/candidate_graph/utils.py b/src/funtracks/candidate_graph/utils.py
index f1fee38c..7d07f2e1 100644
--- a/src/funtracks/candidate_graph/utils.py
+++ b/src/funtracks/candidate_graph/utils.py
@@ -26,8 +26,7 @@ def nodes_from_segmentation(
 
     Each node will have the following attributes:
         - t
-        - pos
-        - area
+        - pos (in pixel coordinates)
         - mask (if mask=True): cropped boolean mask (tracksdata Mask object)
         - bbox (if mask=True): bounding box as [min_0, ..., max_0, ...] int array
 
@@ -38,8 +37,9 @@ 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.
-            Defaults to None, which implies the data is isotropic.
+            Positions and masks stay in pixel coordinates, so this is currently only
+            validated for length; measurements in world units are computed later, by
+            the annotators. 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
             cost) and regionprop.bbox. Including them here avoids a separate write
@@ -65,7 +65,10 @@ def nodes_from_segmentation(
             f"Scale {scale} should have {segmentation.ndim} dims"
         )
 
-    node_attributes = ["pos", "area", "mask", "bbox"] if mask else ["pos", "area"]
+    # No "area" column: it is a regionprops feature now, computed from the mask when
+    # the Tracks is built. Declaring it here would leave a zero-filled column behind,
+    # which Tracks would then activate as-is instead of computing.
+    node_attributes = ["pos", "mask", "bbox"] if mask else ["pos"]
     cand_graph = create_empty_graph(
         node_attributes=node_attributes,
         position_attrs=["pos"],
@@ -80,7 +83,7 @@ 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:]))
+        props = regionprops(segs)
         for regionprop in props:
             node_id = regionprop.label
             if node_id in seen_ids:
@@ -89,7 +92,6 @@ def nodes_from_segmentation(
             attrs: dict[str, Any] = {
                 "t": t + t_start,
                 "pos": list(regionprop.centroid),
-                "area": float(regionprop.area),
             }
             if mask:
                 # regionprop.image is the cropped boolean mask within the bbox,
@@ -114,20 +116,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 +138,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 +182,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 +196,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 +210,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 +221,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 +230,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 +240,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 +252,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)
 
diff --git a/src/funtracks/data_model/tracks.py b/src/funtracks/data_model/tracks.py
index 4b6fb278..0f4f0bf3 100644
--- a/src/funtracks/data_model/tracks.py
+++ b/src/funtracks/data_model/tracks.py
@@ -52,6 +52,13 @@ class Tracks:
     The graph nodes represent detections and must have a time attribute and
     position attribute. Edges in the graph represent links across time.
 
+    Coordinate convention: node positions (and segmentation masks/bounding boxes)
+    are stored in pixel coordinates. ``scale`` holds the size of one pixel per
+    dimension, and is what converts those coordinates to world units. Derived
+    measurements that are physically meaningful (area/volume, perimeter, ellipsoid
+    axes) are computed with the scale applied and are therefore already in world
+    units.
+
     Attributes:
         graph_full (td.graph.BaseGraph): The full graph (first-class): every node/edge
             ever known, including soft-deleted (solution=False) candidates. Nodes
@@ -60,7 +67,8 @@ class Tracks:
             graph_full; the user-visible tracking solution.
         features (FeatureDict): Dictionary of features tracked on graph nodes/edges.
         annotators (AnnotatorRegistry): List of annotators that compute features.
-        scale (list[float] | None): How much to scale each dimension by, including time.
+        scale (list[float] | None): The size of one pixel in each dimension,
+            including time. Converts the stored pixel coordinates to world units.
         ndim (int): Number of dimensions (3 for 2D+time, 4 for 3D+time).
     """
 
@@ -381,11 +389,13 @@ def _check_existing_feature(self, key: str) -> bool:
         return key in node_attrs
 
     def _setup_core_computed_features(self) -> None:
-        """Sets up core computed position features.
+        """Sets up core computed position and size features.
+
+        Registers the position and area features from the RegionpropsAnnotator into
+        the FeatureDict, activating each if it already exists on the graph or
+        computing it otherwise. Track-id features are handled separately by
+        _ensure_track_features.
 
-        Registers the position feature from the RegionpropsAnnotator into the
-        FeatureDict, activating it if it already exists on the graph or computing it
-        otherwise. Track-id features are handled separately by _ensure_track_features.
         """
         # Import here to avoid circular dependency
         from funtracks.annotators import RegionpropsAnnotator
@@ -397,6 +407,7 @@ def _setup_core_computed_features(self) -> None:
                 if self.features.position_key is None:
                     self.features.position_key = pos_key
                 core_features.append(pos_key)
+                core_features.append(annotator.area_key)
         self._register_core_features(core_features)
 
     def _register_core_features(self, keys: list[str]) -> None:
diff --git a/src/funtracks/features/_node_features.py b/src/funtracks/features/_node_features.py
index c78c7e46..211f318f 100644
--- a/src/funtracks/features/_node_features.py
+++ b/src/funtracks/features/_node_features.py
@@ -88,6 +88,10 @@ def SegBbox(ndim: int) -> Feature:
 def Position(axes: Sequence[str]) -> Feature:
     """A feature to hold the position of a node (time not included).
 
+    Positions are always expressed in pixel coordinates, matching how the
+    segmentation is indexed and how coordinates are stored in geff. Multiply by
+    the spatial part of ``Tracks.scale`` to get world units.
+
     Args:
         axes (Sequence[str]): A sequence of the axis names. Used to infer the number
             of dimensions of the position values.
diff --git a/src/funtracks/import_export/_tracks_builder.py b/src/funtracks/import_export/_tracks_builder.py
index 30bd66e8..d9f208fd 100644
--- a/src/funtracks/import_export/_tracks_builder.py
+++ b/src/funtracks/import_export/_tracks_builder.py
@@ -651,7 +651,7 @@ def handle_segmentation(
         # sample_node = next(iter(graph.node_ids()))
         has_position = "pos" in graph.node_attr_keys()
         if has_position:
-            validate_graph_seg_match(graph, seg_array, scale, self.axis_names)
+            validate_graph_seg_match(graph, seg_array, self.axis_names)
 
         # Check if relabeling is needed (seg_id != node_id)
         node_props = self.in_memory_geff["node_props"]
@@ -755,15 +755,6 @@ def _resolve_import_scale(self, scale: list[float] | None) -> list[float] | None
         """
         return scale
 
-    def _scale_to_world_coords(self, scale: list[float] | None) -> None:
-        """Rescale the loaded position columns in ``self.in_memory_geff``.
-
-        No-op by default. Formats that store positions in a different coordinate
-        space than funtracks' world coordinates (e.g. GEFF stores pixels) override
-        this to convert using the resolved import ``scale``.
-        """
-        return
-
     def build(
         self,
         source: Path | pd.DataFrame,
@@ -826,21 +817,17 @@ def build(
         # Validate node_name_map is complete and valid
         self.validate_name_map(has_segmentation=segmentation is not None)
 
-        # Resolve the scale to use for this import. By default the caller-provided
-        # scale is used (e.g. from a scale widget for CSV imports); format-specific
-        # builders may override this (the GEFF builder reads the authoritative
-        # scale from the geff metadata).
-        scale = self._resolve_import_scale(scale)
-
         # 1. Load source data to InMemoryGeff
         self.load_source(source, self.node_name_map)
         if self.in_memory_geff is None:
             raise ValueError("load_source() must populate self.in_memory_geff")
 
-        # Rescale loaded positions into world coordinates if the format requires
-        # it (GEFF stores pixel coordinates). Uses the same resolved scale that is
-        # attached to the resulting Tracks, so the two stay consistent.
-        self._scale_to_world_coords(scale)
+        # Resolve the scale to use for this import. By default the caller-provided
+        # scale is used (e.g. from a scale widget); format-specific builders may
+        # override this to read a scale out of the source (the GEFF builder from the
+        # axes metadata). Runs after load_source so those builders can use anything the
+        # load turned up.
+        scale = self._resolve_import_scale(scale)
 
         # 2. Combine multi-value feature columns
         self._combine_multi_value_props(
diff --git a/src/funtracks/import_export/_validation.py b/src/funtracks/import_export/_validation.py
index 27cc028c..bd65929a 100644
--- a/src/funtracks/import_export/_validation.py
+++ b/src/funtracks/import_export/_validation.py
@@ -24,18 +24,17 @@
 def validate_graph_seg_match(
     graph: td.graph.BaseGraph,
     segmentation: da.Array,
-    scale: list[float],
     position_attr: list[str],
 ) -> None:
     """Validate if the graph matches the provided segmentation data.
 
     Checks if the seg_id value of the last node matches the pixel value at the
-    (scaled) node coordinates.
+    node coordinates. Node coordinates and the segmentation are both in pixel
+    space, so no scaling is applied.
 
     Args:
         graph: tracksdata graph with standard keys
         segmentation: Segmentation data (dask array)
-        scale: Scaling information (pixel to world coordinates)
         position_attr: Position keys (e.g., ["y", "x"] or ["z", "y", "x"])
 
     Raises:
@@ -72,17 +71,13 @@ def validate_graph_seg_match(
 
     # Check bounds
     for i, (c, s) in enumerate(zip(coord, segmentation.shape, strict=False)):
-        pixel_coord = int(c / scale[i])
-        if not (0 <= pixel_coord < s):
+        if not (0 <= int(c) < s):
             raise ValueError(
-                f"Coordinate {i} ({c}) is out of bounds for segmentation shape {s} "
-                f"(pixel coord: {pixel_coord})"
+                f"Coordinate {i} ({c}) is out of bounds for segmentation shape {s}"
             )
 
     # Check if the segmentation pixel value at the coordinates matches the seg id
-    seg_id_at_coord, errors = has_seg_ids_at_coords(
-        segmentation, [coord], [seg_id], tuple(1 / s for s in scale)
-    )
+    seg_id_at_coord, errors = has_seg_ids_at_coords(segmentation, [coord], [seg_id])
     if not seg_id_at_coord:
         error_msg = "Error testing seg id:\n" + "\n".join(f"- {e}" for e in errors)
         raise ValueError(error_msg)
diff --git a/src/funtracks/import_export/csv/_export.py b/src/funtracks/import_export/csv/_export.py
index f4d42c35..0eb9c7c1 100644
--- a/src/funtracks/import_export/csv/_export.py
+++ b/src/funtracks/import_export/csv/_export.py
@@ -36,6 +36,10 @@ def export_to_csv(
     and all registered features. Optionally also exports the segmentation as zarr or
     tiff. If a color dictionary is provided, it will also export the tracklet colors.
 
+    Coordinates are written as pixel coordinates. The spatial scale is written as one
+    constant column per axis (``z_scale``/``y_scale``/``x_scale``, defaulting to 1.0),
+    so a reader can convert to world units.
+
     Args:
         tracks: Tracks object containing the tracking data to export
         outfile: Path to output CSV file
@@ -164,6 +168,14 @@ def convert_numpy_to_python(value):
 
             column_map[feature_name] = names
 
+    # Exported coordinates are pixel coordinates, so carry the voxel size along in
+    # one constant column per spatial axis - a CSV has nowhere else to put it.
+    axes = ["z", "y", "x"] if tracks.ndim == 4 else ["y", "x"]
+    scale_columns = [f"{axis}_scale" for axis in axes]
+    # tracks.scale is optional; an absent scale means one world unit per pixel
+    scale_values = [1.0] * len(axes) if tracks.scale is None else list(tracks.scale[1:])
+    header.extend(scale_columns)
+
     # Determine which nodes to export
     if node_ids is None:
         nodes_to_keep = graph.node_ids()
@@ -206,6 +218,7 @@ def convert_numpy_to_python(value):
 
             row[cast(str, column_map["tracklet_id"])] = tracks.get_track_id(node_id)
 
+        row.update(zip(scale_columns, scale_values, strict=True))
         rows.append(row)
 
     df = pd.DataFrame(rows)
diff --git a/src/funtracks/import_export/geff/_export.py b/src/funtracks/import_export/geff/_export.py
index ab45c4ba..e02b5545 100644
--- a/src/funtracks/import_export/geff/_export.py
+++ b/src/funtracks/import_export/geff/_export.py
@@ -6,7 +6,6 @@
 )
 
 import geff_spec
-import numpy as np
 import polars as pl
 import tracksdata as td
 from geff_spec import GeffMetadata
@@ -178,22 +177,6 @@ def _build_geff_metadata(
             }
         )
 
-    # Convert world coordinates to pixel coordinates by dividing each spatial
-    # position attribute by its axis scale, to align with the geff format, which
-    # stores pixel coordinates together with a per-axis scale in the metadata.
-    # `graph` is a detached copy (see split_position_attr), so this does not
-    # mutate the caller's live graph.
-    pos_pixel_coords: dict[str, np.ndarray] = {}
-    node_df = graph.node_attrs()
-    for name, axis_type, axis_scale in zip(
-        axis_names, axis_types, tracks.scale, strict=True
-    ):
-        if axis_type != "space" or axis_scale == 1:
-            continue
-        pos_pixel_coords[name] = node_df[name].to_numpy() / axis_scale
-    if pos_pixel_coords:
-        graph.update_node_attrs(attrs=pos_pixel_coords, node_ids=graph.node_ids())
-
     extra: dict = {}
     if include_features:
         extra["funtracks"] = {"features": tracks.features.dump_json()}
@@ -277,11 +260,7 @@ def split_position_attr(tracks: Tracks) -> tuple[td.graph.GraphView, list[str] |
         new_graph.remove_node_attr_key(pos_key)
         return new_graph, new_keys
     elif pos_key is not None:
-        # Position is already split into separate attributes. Detach so that
-        # downstream coordinate rescaling operates on an independent copy and
-        # does not mutate the live graph.
-        new_graph = tracks.graph.detach()
-        new_graph = new_graph.filter().subgraph()
-        return new_graph, list(pos_key)
+        # Position is already split into separate attributes
+        return tracks.graph_solution, list(pos_key)
     else:
         return tracks.graph_solution, None
diff --git a/src/funtracks/import_export/geff/_import.py b/src/funtracks/import_export/geff/_import.py
index 3078131d..93487bff 100644
--- a/src/funtracks/import_export/geff/_import.py
+++ b/src/funtracks/import_export/geff/_import.py
@@ -3,7 +3,6 @@
 import warnings
 from typing import TYPE_CHECKING
 
-import numpy as np
 import tracksdata as td
 import zarr
 from geff._typing import InMemoryGeff
@@ -373,43 +372,12 @@ def _scale_from_axes(self) -> list[float] | None:
 
     def _resolve_import_scale(self, scale: list[float] | None) -> list[float] | None:
         """For GEFF, the per-axis scale stored in the metadata is used when present, the
-        caller-provided scale is only used when the metadata has no scaled axes.
+        caller-provided scale is only used when the metadata has no scaled axes. The
+        coordinates are imported as is, since they are expected to be pixel coordinates.
         """
         geff_scale = self._scale_from_axes()
         return geff_scale if geff_scale is not None else scale
 
-    def _scale_to_world_coords(self, scale: list[float] | None) -> None:
-        """Convert stored pixel coordinates to world coordinates.
-
-        geff stores spatial coordinates in pixels together with a per-axis scale;
-        funtracks keeps positions in world coordinates, so each spatial position
-        column is multiplied by its scale (from `builder._resolve_import_scale`), so
-        the conversion stays consistent with the scale attached to the Tracks.
-        """
-        if scale is None or self.in_memory_geff is None:
-            return
-
-        spatial_scale = list(scale)[1:]  # drop the time axis
-        pos_cols = self.position_attr
-        node_props = self.in_memory_geff["node_props"]
-
-        # Composite position stored as separate columns (e.g. "z", "y", "x").
-        if (
-            pos_cols
-            and len(pos_cols) == len(spatial_scale)
-            and all(col in node_props for col in pos_cols)
-        ):
-            for col, s in zip(pos_cols, spatial_scale, strict=True):
-                if s != 1:
-                    node_props[col]["values"] = node_props[col]["values"] * s
-            return
-
-        # Position stored as a single vector column named "pos".
-        if "pos" in node_props:
-            values = node_props["pos"]["values"]
-            if values.ndim == 2 and values.shape[1] == len(spatial_scale):
-                node_props["pos"]["values"] = values * np.asarray(spatial_scale)
-
 
 def import_from_geff(
     directory: Path,
diff --git a/tests/annotators/test_annotator_registry.py b/tests/annotators/test_annotator_registry.py
index 30996ff3..f4462f56 100644
--- a/tests/annotators/test_annotator_registry.py
+++ b/tests/annotators/test_annotator_registry.py
@@ -69,13 +69,13 @@ def test_enable_disable_features(graph_2d_with_segmentation):
     assert "t" in tracks.features
     assert tracks.graph_solution.nodes[nodes[0]]["pos"] is not None
 
-    # area and other features should NOT be in tracks.features initially
-    assert "area" not in tracks.features
+    # area is a core feature and is on from the start; the rest are opt-in
+    assert "area" in tracks.features
     assert "iou" not in tracks.features
     assert "circularity" not in tracks.features
 
     # Enable multiple features at once
-    tracks.enable_features(["area", "iou", "circularity"])
+    tracks.enable_features(["iou", "circularity"])
 
     # Features should now be in FeatureDict
     assert "iou" in tracks.features
@@ -105,23 +105,38 @@ def test_enable_disable_features(graph_2d_with_segmentation):
     assert "circularity" not in tracks.features
 
 
-def test_area_on_graph_not_auto_activated(graph_2d_with_segmentation):
-    """Area pre-populated on a raw graph must not be auto-activated.
+def test_area_activated_by_default(graph_2d_with_segmentation):
+    """Area is a core computed feature: on by default wherever there is a mask.
 
-    Lock in the behavior introduced when area was removed from the core
-    auto-detected features: even though the graph already carries area values,
-    Tracks(graph) without an explicit FeatureDict must leave area out of
-    tracks.features. Callers opt in via enable_features(["area"]).
+    It comes out of the same regionprops pass as the centroid, so it is nearly
+    free, and it is the measurement users reach for first. Values already on the
+    graph are reused rather than recomputed.
     """
     assert "area" in graph_2d_with_segmentation.node_attr_keys()
+    node = next(iter(graph_2d_with_segmentation.node_ids()))
+    existing = graph_2d_with_segmentation.nodes[node]["area"]
 
     tracks = Tracks(graph_2d_with_segmentation, ndim=3, **track_attrs)
 
+    assert "area" in tracks.features
+    assert tracks.graph_solution.nodes[node]["area"] == existing
+
+    # and it can still be turned off
+    tracks.disable_features(["area"])
     assert "area" not in tracks.features
-    assert "area" in tracks.annotators.all_features
 
-    tracks.enable_features(["area"])
+
+def test_area_computed_when_absent_from_graph(get_graph):
+    """With no area column on the graph, the default activation computes it."""
+    graph = get_graph(3, with_seg=True)
+    if "area" in graph.node_attr_keys():
+        graph.remove_node_attr_key("area")
+
+    tracks = Tracks(graph, ndim=3, **track_attrs)
+
     assert "area" in tracks.features
+    node = next(iter(tracks.graph_solution.node_ids()))
+    assert tracks.graph_solution.nodes[node]["area"] > 0
 
 
 def test_get_available_features(graph_2d_with_segmentation):
diff --git a/tests/annotators/test_regionprops_annotator.py b/tests/annotators/test_regionprops_annotator.py
index 7ee98ee3..16f5f98f 100644
--- a/tests/annotators/test_regionprops_annotator.py
+++ b/tests/annotators/test_regionprops_annotator.py
@@ -133,15 +133,17 @@ def test_missing_seg(self, get_graph, ndim):
         # Should not raise an error, just return silently
         rp_ann.compute()  # No error expected
 
-    def test_centroid_world_coords_with_scale(self, get_graph, ndim):
-        """Centroid in 'pos' must be pixel_centroid * scale (world units).
-
-        Without the fix, skimage returns local_centroid * spacing + bbox_min_pixel
-        (mixed units). The correct formula is (local_centroid + bbox_min_pixel) *
-        spacing = pixel_centroid * spacing.
+    @pytest.mark.parametrize("only_pos", [True, False])
+    def test_centroid_pixel_coords_with_scale(self, get_graph, ndim, only_pos):
+        """Centroid in 'pos' must stay in pixel coordinates, whatever the scale.
 
         Node 6 has a cube/square at corner (96, 96, ...) with width 4,
-        so pixel centroid = 97.5 in each spatial axis.
+        so the pixel centroid is 97.5 in each spatial axis.
+
+        Both computation paths are covered: the fast mask-only path taken when
+        'pos' is the sole requested feature, and the skimage regionprops path
+        (which runs with the scale as spacing) taken when other features are
+        requested alongside it.
         """
         graph = get_graph(ndim, with_seg=True)
         if ndim == 3:
@@ -153,22 +155,19 @@ def test_centroid_world_coords_with_scale(self, get_graph, ndim):
 
         tracks = Tracks(graph, ndim=ndim, scale=scale, **track_attrs)
         # Force recomputation so regionprops runs with the given scale as spacing
-        tracks.enable_features(["pos"])
+        tracks.enable_features(["pos"] if only_pos else ["pos", "area"])
 
         pos = np.array(tracks.graph_solution.nodes[6]["pos"])
-        expected = pixel_centroid * np.array(scale[1:])
 
-        bug_value = np.array([1.5] * len(pixel_centroid)) * np.array(
-            scale[1:]
-        ) + np.array([96.0] * len(pixel_centroid))
+        world_value = pixel_centroid * np.array(scale[1:])
         np.testing.assert_allclose(
             pos,
-            expected,
+            pixel_centroid,
             atol=0.1,
             err_msg=(
-                f"World centroid must be pixel_centroid * scale. "
-                f"Got {pos}, expected {expected}. "
-                f"Bug value would be local_centroid * scale + bbox_min = {bug_value}"
+                f"Centroid must be stored in pixel coordinates. "
+                f"Got {pos}, expected {pixel_centroid}. "
+                f"World coordinates would be {world_value}"
             ),
         )
 
diff --git a/tests/candidate_graph/test_compute_graph.py b/tests/candidate_graph/test_compute_graph.py
index ace7e57f..6bf97e30 100644
--- a/tests/candidate_graph/test_compute_graph.py
+++ b/tests/candidate_graph/test_compute_graph.py
@@ -20,9 +20,10 @@ def test_graph_from_segmentation_2d(get_tracks):
     # Same node IDs as the segmentation labels
     assert set(cand_graph.node_ids()) == set(tracks.graph_solution.node_ids())
 
-    # t, pos, area must match the source graph for every node
+    # t and pos must match the source graph for every node; measurements are
+    # not computed on the candidate graph
     for node in cand_graph.node_ids():
-        for key in ["t", "pos", "area"]:
+        for key in ["t", "pos"]:
             assert np.array(cand_graph.nodes[node][key]) == pytest.approx(
                 np.array(tracks.graph_solution.nodes[node][key]), abs=0.01
             )
@@ -70,7 +71,7 @@ def test_graph_from_segmentation_3d(get_tracks):
     assert set(cand_graph.node_ids()) == set(tracks.graph_solution.node_ids())
 
     for node in cand_graph.node_ids():
-        for key in ["t", "pos", "area"]:
+        for key in ["t", "pos"]:
             assert np.array(cand_graph.nodes[node][key]) == pytest.approx(
                 np.array(tracks.graph_solution.nodes[node][key]), abs=0.01
             )
@@ -156,7 +157,7 @@ def test_graph_from_segmentation_t_start_zero_matches_default(get_tracks):
     assert set(explicit_graph.node_ids()) == set(default_graph.node_ids())
 
     for node in default_graph.node_ids():
-        for key in ["t", "pos", "area"]:
+        for key in ["t", "pos"]:
             assert np.array(explicit_graph.nodes[node][key]) == pytest.approx(
                 np.array(default_graph.nodes[node][key]), abs=0.01
             )
@@ -188,11 +189,13 @@ def test_graph_from_points_list():
     assert cand_graph.num_edges() == 3
     assert len(list(cand_graph.predecessors(3))) == 0
 
-    # test scale
+    # test scale: distances are measured in world units, so stretching x by 5
+    # pushes every candidate beyond max_edge_distance. The stored positions stay
+    # in pixel coordinates.
     cand_graph = compute_graph_from_points_list(
         points_list, max_edge_distance=3, scale=[1, 1, 1, 5]
     )
     assert cand_graph.num_edges() == 0
     assert len(list(cand_graph.predecessors(3))) == 0
-    assert np.array(cand_graph.nodes[0]["pos"]) == pytest.approx([1, 1, 5])
+    assert np.array(cand_graph.nodes[0]["pos"]) == pytest.approx([1, 1, 1])
     assert cand_graph.nodes[0]["t"] == 0
diff --git a/tests/candidate_graph/test_graph_utils.py b/tests/candidate_graph/test_graph_utils.py
index dd464ade..24d8ae2e 100644
--- a/tests/candidate_graph/test_graph_utils.py
+++ b/tests/candidate_graph/test_graph_utils.py
@@ -29,20 +29,20 @@ def test_nodes_from_segmentation_2d(get_tracks):
     )
     assert sorted(node_graph.node_ids()) == [1, 2, 3, 4, 5, 6]
     assert node_graph.nodes[2]["t"] == 1
-    assert node_graph.nodes[2]["area"] == 305
     assert np.array_equal(node_graph.nodes[2]["pos"], np.array([20, 80]))
+    # measurements are not computed here; they are regionprops features on Tracks
+    assert "area" not in node_graph.node_attr_keys()
 
     assert node_frame_dict[0] == [1]
     assert Counter(node_frame_dict[1]) == Counter([2, 3])
 
-    # test with scaling
+    # a scale does not move the positions: they stay in pixel coordinates
     node_graph, node_frame_dict = nodes_from_segmentation(
         segmentation=segmentation_2d, scale=[1, 1, 2]
     )
     assert sorted(node_graph.node_ids()) == [1, 2, 3, 4, 5, 6]
     assert node_graph.nodes[2]["t"] == 1
-    assert node_graph.nodes[2]["area"] == 610
-    assert np.array_equal(node_graph.nodes[2]["pos"], np.array([20, 160]))
+    assert np.array_equal(node_graph.nodes[2]["pos"], np.array([20, 80]))
 
     assert node_frame_dict[0] == [1]
     assert Counter(node_frame_dict[1]) == Counter([2, 3])
@@ -58,20 +58,20 @@ def test_nodes_from_segmentation_3d(get_tracks):
     )
     assert sorted(node_graph.node_ids()) == [1, 2, 3, 4, 5, 6]
     assert node_graph.nodes[2]["t"] == 1
-    assert node_graph.nodes[2]["area"] == 4169
     assert np.array_equal(node_graph.nodes[2]["pos"], np.array([20, 50, 80]))
+    # measurements are not computed here; they are regionprops features on Tracks
+    assert "area" not in node_graph.node_attr_keys()
 
     assert node_frame_dict[0] == [1]
     assert Counter(node_frame_dict[1]) == Counter([2, 3])
 
-    # test with scaling
+    # a scale does not move the positions: they stay in pixel coordinates
     node_graph, node_frame_dict = nodes_from_segmentation(
         segmentation=segmentation_3d, scale=[1, 1, 4.5, 1]
     )
     assert sorted(node_graph.node_ids()) == [1, 2, 3, 4, 5, 6]
     assert node_graph.nodes[2]["t"] == 1
-    assert node_graph.nodes[2]["area"] == 4169 * 4.5
-    assert np.array_equal(node_graph.nodes[2]["pos"], np.array([20.0, 225.0, 80.0]))
+    assert np.array_equal(node_graph.nodes[2]["pos"], np.array([20.0, 50.0, 80.0]))
 
     assert node_frame_dict[0] == [1]
     assert Counter(node_frame_dict[1]) == Counter([2, 3])
@@ -96,6 +96,30 @@ def test_add_cand_edges_3d(get_tracks):
     assert sorted(node_graph.edge_list()) == [[1, 3]]
 
 
+def test_add_cand_edges_with_scale(get_tracks):
+    """max_edge_distance is in world units, so anisotropy changes which edges fit.
+
+    Positions are in pixels: dist(1→3) ≈ 11.2 pixels, made up of a 10 pixel step
+    along y and a 5 pixel step along x. Stretching x by 4 makes that ≈ 22.4, which
+    no longer fits in a 15 unit budget.
+    """
+    tracks = get_tracks(ndim=3, with_seg=True)
+    segmentation_2d = np.asarray(tracks.segmentation)
+    node_graph, node_frame_dict = nodes_from_segmentation(segmentation_2d)
+
+    add_cand_edges(node_graph, max_edge_distance=15, node_frame_dict=node_frame_dict)
+    assert sorted(node_graph.edge_list()) == [[1, 3]]
+
+    node_graph, node_frame_dict = nodes_from_segmentation(segmentation_2d)
+    add_cand_edges(
+        node_graph,
+        max_edge_distance=15,
+        node_frame_dict=node_frame_dict,
+        scale=[1, 1, 4],
+    )
+    assert sorted(node_graph.edge_list()) == []
+
+
 def test_compute_node_frame_dict(get_tracks):
     tracks = get_tracks(ndim=3, with_seg=True)
     segmentation_2d = np.asarray(tracks.segmentation)
diff --git a/tests/data_model/test_solution_tracks.py b/tests/data_model/test_solution_tracks.py
index 4f1afa81..f3a7b20c 100644
--- a/tests/data_model/test_solution_tracks.py
+++ b/tests/data_model/test_solution_tracks.py
@@ -85,17 +85,20 @@ def test_export_to_csv_with_display_names(
 
     assert len(lines) == tracks.graph_solution.num_nodes() + 1  # add header
 
-    # With display names: ID, Parent ID, Time, y, x, Tracklet ID,
-    # Lineage ID, Area
+    # With display names: ID, Parent ID, Time, y, x, Area, Tracklet ID,
+    # Lineage ID, and the per-axis scale columns. Area sits with the other
+    # regionprops features because it is registered at construction.
     header = [
         "ID",
         "Parent ID",
         "Time",
         "y",
         "x",
+        "Area",
         "Tracklet ID",
         "Lineage ID",
-        "Area",
+        "y_scale",
+        "x_scale",
     ]
     assert lines[0].strip().split(",") == header
 
@@ -109,8 +112,8 @@ def test_export_to_csv_with_display_names(
 
     assert len(lines) == tracks.graph_solution.num_nodes() + 1  # add header
 
-    # With display names: ID, Parent ID, Time, z, y, x,
-    # Tracklet ID, Lineage ID, Volume
+    # With display names: ID, Parent ID, Time, z, y, x, Volume,
+    # Tracklet ID, Lineage ID, and the per-axis scale columns
     header = [
         "ID",
         "Parent ID",
@@ -118,9 +121,12 @@ def test_export_to_csv_with_display_names(
         "z",
         "y",
         "x",
+        "Volume",
         "Tracklet ID",
         "Lineage ID",
-        "Volume",
+        "z_scale",
+        "y_scale",
+        "x_scale",
     ]
     assert lines[0].strip().split(",") == header
 
diff --git a/tests/import_export/test_csv_export.py b/tests/import_export/test_csv_export.py
index ae13a721..5083c99c 100644
--- a/tests/import_export/test_csv_export.py
+++ b/tests/import_export/test_csv_export.py
@@ -9,8 +9,25 @@
 @pytest.mark.parametrize(
     ("ndim", "expected_header"),
     [
-        (3, ["t", "y", "x", "id", "parent_id", "track_id"]),
-        (4, ["t", "z", "y", "x", "id", "parent_id", "track_id"]),
+        (
+            3,
+            ["t", "y", "x", "id", "parent_id", "track_id", "y_scale", "x_scale"],
+        ),
+        (
+            4,
+            [
+                "t",
+                "z",
+                "y",
+                "x",
+                "id",
+                "parent_id",
+                "track_id",
+                "z_scale",
+                "y_scale",
+                "x_scale",
+            ],
+        ),
     ],
     ids=["2d", "3d"],
 )
@@ -27,18 +44,36 @@ def test_export_solution_to_csv(get_tracks, tmp_path, ndim, expected_header):
     assert lines[0].strip().split(",") == expected_header
 
     # Check first data line (node 1: t=0, pos=[50, 50] or [50, 50, 50], track_id=1)
+    # trailing values are the per-axis scale columns, 1.0 when tracks has no scale
     if ndim == 3:
-        expected_line1 = ["0", "50.0", "50.0", "1", "", "1"]
+        expected_line1 = ["0", "50.0", "50.0", "1", "", "1", "1.0", "1.0"]
     else:
-        expected_line1 = ["0", "50.0", "50.0", "50.0", "1", "", "1"]
+        expected_line1 = ["0", "50.0", "50.0", "50.0", "1", "", "1", "1.0", "1.0", "1.0"]
     assert lines[1].strip().split(",") == expected_line1
 
 
 @pytest.mark.parametrize(
     ("ndim", "expected_header"),
     [
-        (3, ["t", "y", "x", "id", "parent_id", "track_id"]),
-        (4, ["t", "z", "y", "x", "id", "parent_id", "track_id"]),
+        (
+            3,
+            ["t", "y", "x", "id", "parent_id", "track_id", "y_scale", "x_scale"],
+        ),
+        (
+            4,
+            [
+                "t",
+                "z",
+                "y",
+                "x",
+                "id",
+                "parent_id",
+                "track_id",
+                "z_scale",
+                "y_scale",
+                "x_scale",
+            ],
+        ),
     ],
     ids=["2d", "3d"],
 )
@@ -322,3 +357,73 @@ def test_export_full_vs_solution(get_tracks, tmp_path):
     assert "Solution" in disp.columns
     assert not bool(disp.loc[disp["ID"] == 5, "Solution"].iloc[0])
     assert bool(disp.loc[disp["ID"] == 1, "Solution"].iloc[0])
+
+
+@pytest.mark.parametrize("ndim", [3, 4])
+def test_scale_columns_carry_the_voxel_size(get_tracks, tmp_path, ndim):
+    """Coordinates are exported in pixels, with the voxel size in its own columns.
+
+    A CSV has no metadata block, so the scale rides along as one constant column
+    per spatial axis. A reader (or the user, via the import dialog) can then
+    convert the pixel coordinates to world units.
+    """
+    import pandas as pd
+
+    scale = [1.0, 2.0, 0.416] if ndim == 3 else [1.0, 2.0, 0.416, 0.416]
+    tracks = get_tracks(ndim=ndim, with_seg=False, prefill_track_ids=True)
+    tracks.scale = scale
+
+    temp_file = tmp_path / "test_scale.csv"
+    export_to_csv(tracks, temp_file)
+
+    df = pd.read_csv(temp_file)
+    axes = ["y", "x"] if ndim == 3 else ["z", "y", "x"]
+    for axis, value in zip(axes, scale[1:], strict=True):
+        assert (df[f"{axis}_scale"] == value).all()
+
+    # Positions are pixel coordinates, written straight from the graph.
+    node_ids = sorted(tracks.graph_solution.node_ids())
+    positions = tracks.get_positions(node_ids)
+    exported = df.set_index("id").loc[node_ids, axes].to_numpy()
+    np.testing.assert_allclose(exported, positions)
+
+
+@pytest.mark.parametrize("ndim", [3, 4])
+def test_scale_columns_default_to_one_without_a_scale(get_tracks, tmp_path, ndim):
+    """tracks.scale is optional; an absent scale means one world unit per pixel."""
+    import pandas as pd
+
+    tracks = get_tracks(ndim=ndim, with_seg=False, prefill_track_ids=True)
+    assert tracks.scale is None
+
+    temp_file = tmp_path / "test_no_scale.csv"
+    export_to_csv(tracks, temp_file)
+
+    df = pd.read_csv(temp_file)
+    axes = ["y", "x"] if ndim == 3 else ["z", "y", "x"]
+    for axis in axes:
+        assert (df[f"{axis}_scale"] == 1.0).all()
+
+
+def test_scale_columns_are_not_imported_as_features(get_tracks, tmp_path):
+    """A round trip must not turn the scale columns into per-node features."""
+    import pandas as pd
+
+    from funtracks.import_export import tracks_from_df
+
+    tracks = get_tracks(ndim=3, with_seg=False, prefill_track_ids=True)
+    tracks.scale = [1.0, 2.0, 0.416]
+    temp_file = tmp_path / "test_roundtrip.csv"
+    export_to_csv(tracks, temp_file)
+
+    loaded = tracks_from_df(pd.read_csv(temp_file))
+
+    for column in ("y_scale", "x_scale"):
+        assert column not in loaded.features
+        assert column not in loaded.graph_solution.node_attr_keys()
+
+    # the coordinates themselves come back untouched
+    node_ids = sorted(tracks.graph_solution.node_ids())
+    np.testing.assert_allclose(
+        loaded.get_positions(node_ids), tracks.get_positions(node_ids)
+    )
diff --git a/tests/import_export/test_export_to_geff.py b/tests/import_export/test_export_to_geff.py
index 573c7030..427710d7 100644
--- a/tests/import_export/test_export_to_geff.py
+++ b/tests/import_export/test_export_to_geff.py
@@ -398,6 +398,43 @@ def test_write_to_geff_roundtrip(get_tracks, ndim, with_seg, tmp_path):
         )
 
 
+@pytest.mark.parametrize("ndim", [3, 4])
+def test_geff_roundtrip_keeps_pixel_coordinates(get_tracks, ndim, tmp_path):
+    """Positions are pixel coordinates on both sides of a geff round trip.
+
+    geff stores pixel coordinates plus a per-axis scale, which is exactly how
+    funtracks stores them, so an anisotropic scale must not move any coordinate:
+    not the ones written into the store, and not the ones read back.
+    """
+    scale = [1.0, 2.0, 3.0] if ndim == 3 else [1.0, 2.0, 3.0, 4.0]
+    axis_names = ["y", "x"] if ndim == 3 else ["z", "y", "x"]
+
+    tracks = get_tracks(ndim=ndim, with_seg=False, prefill_track_ids=True)
+    tracks.scale = scale
+    node_ids = sorted(tracks.graph_solution.node_ids())
+    original_pos = tracks.get_positions(node_ids)
+
+    geff_path = tmp_path / "my_tracks.geff"
+    write_to_geff(tracks, geff_path)
+
+    # The coordinates in the store are the unchanged pixel coordinates, and the
+    # scale is carried alongside them in the axes metadata.
+    z = zarr.open(str(geff_path), mode="r")
+    stored = np.stack(
+        [np.asarray(z[f"nodes/props/{name}/values"][:]) for name in axis_names], axis=1
+    )
+    np.testing.assert_allclose(stored, original_pos)
+    axes = {ax["name"]: ax["scale"] for ax in dict(z.attrs)["geff"]["axes"]}
+    assert [axes[name] for name in axis_names] == scale[1:]
+
+    # Exporting must not have touched the live graph either.
+    np.testing.assert_allclose(tracks.get_positions(node_ids), original_pos)
+
+    loaded = import_from_geff(geff_path)
+    np.testing.assert_allclose(loaded.get_positions(node_ids), original_pos)
+    assert list(loaded.scale) == scale
+
+
 def test_write_to_geff_no_parent_container(get_tracks, tmp_path):
     """write_to_geff writes directly to the path — no parent .zgroup or
     tracks.geff subfolder."""
diff --git a/tests/import_export/test_import_from_geff.py b/tests/import_export/test_import_from_geff.py
index f514269b..1dac339d 100644
--- a/tests/import_export/test_import_from_geff.py
+++ b/tests/import_export/test_import_from_geff.py
@@ -73,11 +73,8 @@ def valid_segmentation():
     y = [100, 200, 300, 400, 500]
     seg_ids = np.array([10, 20, 30, 40, 50])
 
-    # The geff stores pixel coordinates. On import each spatial coordinate is
-    # multiplied by its scale to get world coordinates, and the validation
-    # converts back to pixels via ``pixel = world / scale``, recovering the
-    # original stored coordinate. So place each seg id at the pixel equal to the
-    # stored geff coordinate.
+    # The geff stores pixel coordinates and funtracks imports them unchanged, so
+    # each seg id sits at the pixel equal to the stored geff coordinate.
     for t, y_val, x_f, seg_id in zip(times, y, x, seg_ids, strict=False):
         seg[t, int(y_val), int(x_f)] = seg_id
     return seg
@@ -332,7 +329,9 @@ def test_tracks_with_segmentation(valid_geff, invalid_geff, valid_segmentation,
         pos[0],  # y
         pos[1],  # x
     ]
-    coords = tuple(int(c * 1 / s) for c, s in zip(coords, scale, strict=True))
+    # Positions are stored in pixel coordinates, so they index the segmentation
+    # directly - no conversion by the scale needed.
+    coords = tuple(int(c) for c in coords)
     assert (
         valid_segmentation[tuple(coords)] == 50
     )  # in original segmentation, the pixel value is equal to seg_id

From 36643bebdb276efbde8bc162eeec2167f00c2571 Mon Sep 17 00:00:00 2001
From: AnniekStok <anniek.stokkermans@gmail.com>
Date: Fri, 21 Aug 2026 18:57:29 +0200
Subject: [PATCH 3/5] put area computation back on compute_graph_from_seg, but
 scale it

---
 src/funtracks/candidate_graph/utils.py      | 18 ++++++----
 src/funtracks/data_model/tracks.py          | 11 +++---
 tests/annotators/test_annotator_registry.py | 37 ++++++---------------
 tests/candidate_graph/test_compute_graph.py |  9 +++--
 tests/candidate_graph/test_graph_utils.py   | 12 +++----
 tests/data_model/test_solution_tracks.py    | 13 ++++----
 tests/import_export/test_csv_export.py      | 24 -------------
 7 files changed, 42 insertions(+), 82 deletions(-)

diff --git a/src/funtracks/candidate_graph/utils.py b/src/funtracks/candidate_graph/utils.py
index 7d07f2e1..b448719e 100644
--- a/src/funtracks/candidate_graph/utils.py
+++ b/src/funtracks/candidate_graph/utils.py
@@ -27,6 +27,7 @@ def nodes_from_segmentation(
     Each node will have the following attributes:
         - t
         - 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
 
@@ -37,9 +38,9 @@ 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).
-            Positions and masks stay in pixel coordinates, so this is currently only
-            validated for length; measurements in world units are computed later, by
-            the annotators. Defaults to None, which implies the data is isotropic.
+            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
             cost) and regionprop.bbox. Including them here avoids a separate write
@@ -65,10 +66,10 @@ def nodes_from_segmentation(
             f"Scale {scale} should have {segmentation.ndim} dims"
         )
 
-    # No "area" column: it is a regionprops feature now, computed from the mask when
-    # the Tracks is built. Declaring it here would leave a zero-filled column behind,
-    # which Tracks would then activate as-is instead of computing.
-    node_attributes = ["pos", "mask", "bbox"] if mask else ["pos"]
+    # 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,
         position_attrs=["pos"],
@@ -83,6 +84,8 @@ def nodes_from_segmentation(
     for t in tqdm(range(len(segmentation)), desc="Extracting nodes from segmentation"):
         segs = segmentation[t]
         nodes_in_frame: list[int] = []
+        # 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
@@ -92,6 +95,7 @@ def nodes_from_segmentation(
             attrs: dict[str, Any] = {
                 "t": t + t_start,
                 "pos": list(regionprop.centroid),
+                "area": float(regionprop.area) * voxel_volume,
             }
             if mask:
                 # regionprop.image is the cropped boolean mask within the bbox,
diff --git a/src/funtracks/data_model/tracks.py b/src/funtracks/data_model/tracks.py
index 0f4f0bf3..81e8c338 100644
--- a/src/funtracks/data_model/tracks.py
+++ b/src/funtracks/data_model/tracks.py
@@ -389,13 +389,11 @@ def _check_existing_feature(self, key: str) -> bool:
         return key in node_attrs
 
     def _setup_core_computed_features(self) -> None:
-        """Sets up core computed position and size features.
-
-        Registers the position and area features from the RegionpropsAnnotator into
-        the FeatureDict, activating each if it already exists on the graph or
-        computing it otherwise. Track-id features are handled separately by
-        _ensure_track_features.
+        """Sets up core computed position features.
 
+        Registers the position feature from the RegionpropsAnnotator into the
+        FeatureDict, activating it if it already exists on the graph or computing it
+        otherwise. Track-id features are handled separately by _ensure_track_features.
         """
         # Import here to avoid circular dependency
         from funtracks.annotators import RegionpropsAnnotator
@@ -407,7 +405,6 @@ def _setup_core_computed_features(self) -> None:
                 if self.features.position_key is None:
                     self.features.position_key = pos_key
                 core_features.append(pos_key)
-                core_features.append(annotator.area_key)
         self._register_core_features(core_features)
 
     def _register_core_features(self, keys: list[str]) -> None:
diff --git a/tests/annotators/test_annotator_registry.py b/tests/annotators/test_annotator_registry.py
index f4462f56..30996ff3 100644
--- a/tests/annotators/test_annotator_registry.py
+++ b/tests/annotators/test_annotator_registry.py
@@ -69,13 +69,13 @@ def test_enable_disable_features(graph_2d_with_segmentation):
     assert "t" in tracks.features
     assert tracks.graph_solution.nodes[nodes[0]]["pos"] is not None
 
-    # area is a core feature and is on from the start; the rest are opt-in
-    assert "area" in tracks.features
+    # area and other features should NOT be in tracks.features initially
+    assert "area" not in tracks.features
     assert "iou" not in tracks.features
     assert "circularity" not in tracks.features
 
     # Enable multiple features at once
-    tracks.enable_features(["iou", "circularity"])
+    tracks.enable_features(["area", "iou", "circularity"])
 
     # Features should now be in FeatureDict
     assert "iou" in tracks.features
@@ -105,38 +105,23 @@ def test_enable_disable_features(graph_2d_with_segmentation):
     assert "circularity" not in tracks.features
 
 
-def test_area_activated_by_default(graph_2d_with_segmentation):
-    """Area is a core computed feature: on by default wherever there is a mask.
+def test_area_on_graph_not_auto_activated(graph_2d_with_segmentation):
+    """Area pre-populated on a raw graph must not be auto-activated.
 
-    It comes out of the same regionprops pass as the centroid, so it is nearly
-    free, and it is the measurement users reach for first. Values already on the
-    graph are reused rather than recomputed.
+    Lock in the behavior introduced when area was removed from the core
+    auto-detected features: even though the graph already carries area values,
+    Tracks(graph) without an explicit FeatureDict must leave area out of
+    tracks.features. Callers opt in via enable_features(["area"]).
     """
     assert "area" in graph_2d_with_segmentation.node_attr_keys()
-    node = next(iter(graph_2d_with_segmentation.node_ids()))
-    existing = graph_2d_with_segmentation.nodes[node]["area"]
 
     tracks = Tracks(graph_2d_with_segmentation, ndim=3, **track_attrs)
 
-    assert "area" in tracks.features
-    assert tracks.graph_solution.nodes[node]["area"] == existing
-
-    # and it can still be turned off
-    tracks.disable_features(["area"])
     assert "area" not in tracks.features
+    assert "area" in tracks.annotators.all_features
 
-
-def test_area_computed_when_absent_from_graph(get_graph):
-    """With no area column on the graph, the default activation computes it."""
-    graph = get_graph(3, with_seg=True)
-    if "area" in graph.node_attr_keys():
-        graph.remove_node_attr_key("area")
-
-    tracks = Tracks(graph, ndim=3, **track_attrs)
-
+    tracks.enable_features(["area"])
     assert "area" in tracks.features
-    node = next(iter(tracks.graph_solution.node_ids()))
-    assert tracks.graph_solution.nodes[node]["area"] > 0
 
 
 def test_get_available_features(graph_2d_with_segmentation):
diff --git a/tests/candidate_graph/test_compute_graph.py b/tests/candidate_graph/test_compute_graph.py
index 6bf97e30..c0993a7f 100644
--- a/tests/candidate_graph/test_compute_graph.py
+++ b/tests/candidate_graph/test_compute_graph.py
@@ -20,10 +20,9 @@ def test_graph_from_segmentation_2d(get_tracks):
     # Same node IDs as the segmentation labels
     assert set(cand_graph.node_ids()) == set(tracks.graph_solution.node_ids())
 
-    # t and pos must match the source graph for every node; measurements are
-    # not computed on the candidate graph
+    # t, pos, area must match the source graph for every node
     for node in cand_graph.node_ids():
-        for key in ["t", "pos"]:
+        for key in ["t", "pos", "area"]:
             assert np.array(cand_graph.nodes[node][key]) == pytest.approx(
                 np.array(tracks.graph_solution.nodes[node][key]), abs=0.01
             )
@@ -71,7 +70,7 @@ def test_graph_from_segmentation_3d(get_tracks):
     assert set(cand_graph.node_ids()) == set(tracks.graph_solution.node_ids())
 
     for node in cand_graph.node_ids():
-        for key in ["t", "pos"]:
+        for key in ["t", "pos", "area"]:
             assert np.array(cand_graph.nodes[node][key]) == pytest.approx(
                 np.array(tracks.graph_solution.nodes[node][key]), abs=0.01
             )
@@ -157,7 +156,7 @@ def test_graph_from_segmentation_t_start_zero_matches_default(get_tracks):
     assert set(explicit_graph.node_ids()) == set(default_graph.node_ids())
 
     for node in default_graph.node_ids():
-        for key in ["t", "pos"]:
+        for key in ["t", "pos", "area"]:
             assert np.array(explicit_graph.nodes[node][key]) == pytest.approx(
                 np.array(default_graph.nodes[node][key]), abs=0.01
             )
diff --git a/tests/candidate_graph/test_graph_utils.py b/tests/candidate_graph/test_graph_utils.py
index 24d8ae2e..af35f64e 100644
--- a/tests/candidate_graph/test_graph_utils.py
+++ b/tests/candidate_graph/test_graph_utils.py
@@ -29,19 +29,19 @@ def test_nodes_from_segmentation_2d(get_tracks):
     )
     assert sorted(node_graph.node_ids()) == [1, 2, 3, 4, 5, 6]
     assert node_graph.nodes[2]["t"] == 1
+    assert node_graph.nodes[2]["area"] == 305
     assert np.array_equal(node_graph.nodes[2]["pos"], np.array([20, 80]))
-    # measurements are not computed here; they are regionprops features on Tracks
-    assert "area" not in node_graph.node_attr_keys()
 
     assert node_frame_dict[0] == [1]
     assert Counter(node_frame_dict[1]) == Counter([2, 3])
 
-    # a scale does not move the positions: they stay in pixel coordinates
+    # with a scale the area is in world units, while the position stays in pixels
     node_graph, node_frame_dict = nodes_from_segmentation(
         segmentation=segmentation_2d, scale=[1, 1, 2]
     )
     assert sorted(node_graph.node_ids()) == [1, 2, 3, 4, 5, 6]
     assert node_graph.nodes[2]["t"] == 1
+    assert node_graph.nodes[2]["area"] == 610
     assert np.array_equal(node_graph.nodes[2]["pos"], np.array([20, 80]))
 
     assert node_frame_dict[0] == [1]
@@ -58,19 +58,19 @@ def test_nodes_from_segmentation_3d(get_tracks):
     )
     assert sorted(node_graph.node_ids()) == [1, 2, 3, 4, 5, 6]
     assert node_graph.nodes[2]["t"] == 1
+    assert node_graph.nodes[2]["area"] == 4169
     assert np.array_equal(node_graph.nodes[2]["pos"], np.array([20, 50, 80]))
-    # measurements are not computed here; they are regionprops features on Tracks
-    assert "area" not in node_graph.node_attr_keys()
 
     assert node_frame_dict[0] == [1]
     assert Counter(node_frame_dict[1]) == Counter([2, 3])
 
-    # a scale does not move the positions: they stay in pixel coordinates
+    # with a scale the area is in world units, while the position stays in pixels
     node_graph, node_frame_dict = nodes_from_segmentation(
         segmentation=segmentation_3d, scale=[1, 1, 4.5, 1]
     )
     assert sorted(node_graph.node_ids()) == [1, 2, 3, 4, 5, 6]
     assert node_graph.nodes[2]["t"] == 1
+    assert node_graph.nodes[2]["area"] == 4169 * 4.5
     assert np.array_equal(node_graph.nodes[2]["pos"], np.array([20.0, 50.0, 80.0]))
 
     assert node_frame_dict[0] == [1]
diff --git a/tests/data_model/test_solution_tracks.py b/tests/data_model/test_solution_tracks.py
index f3a7b20c..07e04ae5 100644
--- a/tests/data_model/test_solution_tracks.py
+++ b/tests/data_model/test_solution_tracks.py
@@ -85,18 +85,17 @@ def test_export_to_csv_with_display_names(
 
     assert len(lines) == tracks.graph_solution.num_nodes() + 1  # add header
 
-    # With display names: ID, Parent ID, Time, y, x, Area, Tracklet ID,
-    # Lineage ID, and the per-axis scale columns. Area sits with the other
-    # regionprops features because it is registered at construction.
+    # With display names: ID, Parent ID, Time, y, x, Tracklet ID,
+    # Lineage ID, Area, and the per-axis scale columns
     header = [
         "ID",
         "Parent ID",
         "Time",
         "y",
         "x",
-        "Area",
         "Tracklet ID",
         "Lineage ID",
+        "Area",
         "y_scale",
         "x_scale",
     ]
@@ -112,8 +111,8 @@ def test_export_to_csv_with_display_names(
 
     assert len(lines) == tracks.graph_solution.num_nodes() + 1  # add header
 
-    # With display names: ID, Parent ID, Time, z, y, x, Volume,
-    # Tracklet ID, Lineage ID, and the per-axis scale columns
+    # With display names: ID, Parent ID, Time, z, y, x,
+    # Tracklet ID, Lineage ID, Volume, and the per-axis scale columns
     header = [
         "ID",
         "Parent ID",
@@ -121,9 +120,9 @@ def test_export_to_csv_with_display_names(
         "z",
         "y",
         "x",
-        "Volume",
         "Tracklet ID",
         "Lineage ID",
+        "Volume",
         "z_scale",
         "y_scale",
         "x_scale",
diff --git a/tests/import_export/test_csv_export.py b/tests/import_export/test_csv_export.py
index 5083c99c..9ce488ec 100644
--- a/tests/import_export/test_csv_export.py
+++ b/tests/import_export/test_csv_export.py
@@ -403,27 +403,3 @@ def test_scale_columns_default_to_one_without_a_scale(get_tracks, tmp_path, ndim
     axes = ["y", "x"] if ndim == 3 else ["z", "y", "x"]
     for axis in axes:
         assert (df[f"{axis}_scale"] == 1.0).all()
-
-
-def test_scale_columns_are_not_imported_as_features(get_tracks, tmp_path):
-    """A round trip must not turn the scale columns into per-node features."""
-    import pandas as pd
-
-    from funtracks.import_export import tracks_from_df
-
-    tracks = get_tracks(ndim=3, with_seg=False, prefill_track_ids=True)
-    tracks.scale = [1.0, 2.0, 0.416]
-    temp_file = tmp_path / "test_roundtrip.csv"
-    export_to_csv(tracks, temp_file)
-
-    loaded = tracks_from_df(pd.read_csv(temp_file))
-
-    for column in ("y_scale", "x_scale"):
-        assert column not in loaded.features
-        assert column not in loaded.graph_solution.node_attr_keys()
-
-    # the coordinates themselves come back untouched
-    node_ids = sorted(tracks.graph_solution.node_ids())
-    np.testing.assert_allclose(
-        loaded.get_positions(node_ids), tracks.get_positions(node_ids)
-    )

From 512d297f8d363d0645160fb747048bb10d756ad0 Mon Sep 17 00:00:00 2001
From: AnniekStok <anniek.stokkermans@gmail.com>
Date: Mon, 24 Aug 2026 18:30:53 +0200
Subject: [PATCH 4/5] add scale_dependent field to features, to update features
 when the scaling is changed

---
 src/funtracks/data_model/tracks.py            |  42 +++++++
 src/funtracks/features/_feature.py            |   6 +
 .../features/_regionprops_features.py         |   5 +
 tests/data_model/test_tracks.py               | 103 ++++++++++++++++++
 tests/import_export/test_internal_format.py   |   3 +
 5 files changed, 159 insertions(+)

diff --git a/src/funtracks/data_model/tracks.py b/src/funtracks/data_model/tracks.py
index 81e8c338..788b42ef 100644
--- a/src/funtracks/data_model/tracks.py
+++ b/src/funtracks/data_model/tracks.py
@@ -228,6 +228,48 @@ def __init__(
         # and computed. A provided FeatureDict that omitted them is completed here.
         self._ensure_track_features()
 
+    def update_scale(self, scale: list[float] | None) -> None:
+        """Set a new voxel size, recompute the features derived from it, and refresh.
+
+        Recompute regionprops features that use the scale information.
+
+        Args:
+            scale (list[float] | None): The size of one pixel in each dimension,
+                including time. None means all dimensions are unscaled.
+
+        Raises:
+            ValueError: If the scale does not hold exactly one value per dimension.
+        """
+        if scale is not None:
+            if len(scale) != self.ndim:
+                raise ValueError(
+                    f"Scale {scale} must have one value per dimension ({self.ndim})"
+                )
+            scale = [float(s) for s in scale]
+
+        if scale == self.scale:
+            return
+
+        # Update scale dependent features.
+        keys = [
+            key
+            for key, (feature, active) in self.annotators.all_features.items()
+            if active and feature.get("scale_dependent", False)
+        ]
+
+        previous = self.scale
+        self.scale = scale
+        try:
+            if keys:
+                self.annotators.compute(keys)
+        except Exception:
+            # Put the old scale back in case of an error (e.g. anisotropic xy scale while
+            # features like perimeter require isotropic xy scaling).
+            self.scale = previous
+            raise
+
+        self.refresh.emit()
+
     def _get_feature_set(
         self,
         time_attr: str | None,
diff --git a/src/funtracks/features/_feature.py b/src/funtracks/features/_feature.py
index 8bf982fc..5891be5f 100644
--- a/src/funtracks/features/_feature.py
+++ b/src/funtracks/features/_feature.py
@@ -40,6 +40,11 @@ class Feature(TypedDict):
             from this feature and should be cascade-deleted when this feature
             is removed. For example, SegMask lists the bbox key here so that
             deleting the mask also deletes the bounding box.
+        scale_dependent (bool): Optional. True if the value is computed with
+            ``Tracks.scale`` and therefore has to be recomputed when the scale
+            changes (see `Tracks.update_scale`). Features measured in pixel space
+            (`Position`) or derived from graph structure (`TrackletID`, `IoU`)
+            leave this unset.
     """
 
     feature_type: FeatureType
@@ -50,3 +55,4 @@ class Feature(TypedDict):
     default_value: Any
     spatial_dims: NotRequired[bool]
     derived_features: NotRequired[list[str]]
+    scale_dependent: NotRequired[bool]
diff --git a/src/funtracks/features/_regionprops_features.py b/src/funtracks/features/_regionprops_features.py
index 7f51fc8e..ca548855 100644
--- a/src/funtracks/features/_regionprops_features.py
+++ b/src/funtracks/features/_regionprops_features.py
@@ -24,6 +24,7 @@ def Area(ndim: int = 3) -> Feature:
         "num_values": 1,
         "display_name": "Area" if ndim == 3 else "Volume",
         "default_value": None,
+        "scale_dependent": True,
     }
 
 
@@ -71,6 +72,7 @@ def EllipsoidAxes(ndim: int | None = 4) -> Feature:
         "value_names": value_names,
         "default_value": None,
         "spatial_dims": True,
+        "scale_dependent": True,
     }
 
 
@@ -90,6 +92,8 @@ def Circularity(ndim: int = 3) -> Feature:
         "num_values": 1,
         "display_name": "Circularity" if ndim == 3 else "Sphericity",
         "default_value": None,
+        "scale_dependent": True,  # because area and perimeter scale differently with
+        # anisotropy
     }
 
 
@@ -109,4 +113,5 @@ def Perimeter(ndim: int = 3) -> Feature:
         "num_values": 1,
         "display_name": "Perimeter" if ndim == 3 else "Surface Area",
         "default_value": None,
+        "scale_dependent": True,
     }
diff --git a/tests/data_model/test_tracks.py b/tests/data_model/test_tracks.py
index 06f35b81..9fa889c5 100644
--- a/tests/data_model/test_tracks.py
+++ b/tests/data_model/test_tracks.py
@@ -350,3 +350,106 @@ def test_update_mask_syncs_bbox(graph_2d_with_segmentation):
 
     assert stored_mask is new_mask
     assert np.array_equal(stored_bbox, new_mask.bbox)
+
+
+@pytest.fixture
+def scaled_tracks(get_graph) -> Tracks:
+    """2D+time tracks with segmentation, an unscaled scale, and size features on."""
+    tracks = Tracks(
+        get_graph(3, with_seg=True), ndim=3, scale=[1.0, 1.0, 1.0], **track_attrs
+    )
+    tracks.enable_features(["pos", "area", "ellipse_axis_radii"])
+    return tracks
+
+
+def _node_attrs(tracks: Tracks, key: str) -> dict[int, object]:
+    """Every node's value for one feature, as plain floats/lists for comparison."""
+    return {
+        node: np.asarray(tracks.get_node_attr(node, key)).tolist()
+        for node in tracks.graph_full.node_ids()
+    }
+
+
+def test_update_scale(scaled_tracks):
+    areas_before = _node_attrs(scaled_tracks, "area")
+    positions_before = _node_attrs(scaled_tracks, "pos")
+    refreshes = []
+    scaled_tracks.refresh.connect(lambda *args: refreshes.append(args))
+
+    scaled_tracks.update_scale([1, 2, 3])
+
+    assert scaled_tracks.scale == [1.0, 2.0, 3.0]
+    assert len(refreshes) == 1
+
+    # areas are in world units, so they follow the new voxel size
+    areas_after = _node_attrs(scaled_tracks, "area")
+    for node, before in areas_before.items():
+        assert areas_after[node] == pytest.approx(before * 6.0)
+
+    # positions stay in pixel coordinates
+    assert _node_attrs(scaled_tracks, "pos") == positions_before
+
+
+def test_update_scale_unchanged_is_noop(scaled_tracks):
+    refreshes = []
+    scaled_tracks.refresh.connect(lambda *args: refreshes.append(args))
+    scaled_tracks.update_scale([1, 1, 1])
+    assert refreshes == []
+
+
+def test_update_scale_skips_disabled_features(get_graph):
+    tracks = Tracks(
+        get_graph(3, with_seg=True), ndim=3, scale=[1.0, 1.0, 1.0], **track_attrs
+    )
+    tracks.enable_features(["pos"])
+
+    tracks.update_scale([1.0, 2.0, 2.0])
+
+    assert tracks.scale == [1.0, 2.0, 2.0]
+    assert "area" not in tracks.features
+
+
+def test_update_scale_wrong_ndim(scaled_tracks):
+    with pytest.raises(ValueError, match="one value per dimension"):
+        scaled_tracks.update_scale([1.0, 1.0])
+
+
+def test_update_scale_restores_scale_when_recompute_fails(get_graph):
+    """An anisotropic scale is rejected by skimage's perimeter, so nothing changes.
+
+    The tracks must not be left holding the new scale with values measured under
+    the old one.
+    """
+    tracks = Tracks(
+        get_graph(3, with_seg=True), ndim=3, scale=[1.0, 1.0, 1.0], **track_attrs
+    )
+    tracks.enable_features(["area", "perimeter"])
+    areas_before = _node_attrs(tracks, "area")
+    refreshes = []
+    tracks.refresh.connect(lambda *args: refreshes.append(args))
+
+    with pytest.raises(NotImplementedError):
+        tracks.update_scale([1.0, 2.0, 3.0])
+
+    assert tracks.scale == [1.0, 1.0, 1.0]
+    assert _node_attrs(tracks, "area") == areas_before
+    assert refreshes == []
+
+
+def test_update_scale_leaves_pixel_space_features_alone(get_graph, monkeypatch):
+    """Position is `centroid_pixel`, so it is not even handed to the annotator."""
+    tracks = Tracks(
+        get_graph(3, with_seg=True), ndim=3, scale=[1.0, 1.0, 1.0], **track_attrs
+    )
+    tracks.enable_features(["pos", "area"])
+
+    computed = []
+    annotator = tracks.annotators[0]
+    original = annotator.compute
+    monkeypatch.setattr(
+        annotator, "compute", lambda keys=None: computed.append(keys) or original(keys)
+    )
+
+    tracks.update_scale([1.0, 2.0, 3.0])
+
+    assert computed == [["area"]]
diff --git a/tests/import_export/test_internal_format.py b/tests/import_export/test_internal_format.py
index 8e9871b3..56efdf0f 100644
--- a/tests/import_export/test_internal_format.py
+++ b/tests/import_export/test_internal_format.py
@@ -42,6 +42,9 @@ def test_save_load(
         tracks_feature = tracks.features[key]
 
         for attr_name, attr_value in tracks_feature.items():
+            if attr_name not in loaded_feature:
+                # skip fields that are added after the snapshot was made.
+                continue
             loaded_attr_value = loaded_feature[attr_name]
 
             # For sequence attributes, cast to list to compare (handles tuple vs list)

From a43a6d48220f46f355791bb2d40d7b085d8d823f Mon Sep 17 00:00:00 2001
From: Caroline Malin-Mayor <malinmayorc@janelia.hhmi.org>
Date: Thu, 24 Sep 2026 20:55:30 -0400
Subject: [PATCH 5/5] Add property for pixel or world unit positions

---
 .../annotators/_regionprops_annotator.py      | 47 ++++++++---
 src/funtracks/data_model/tracks.py            | 51 ++++++++++--
 src/funtracks/features/_feature.py            |  6 +-
 src/funtracks/features/_node_features.py      |  7 +-
 .../annotators/test_regionprops_annotator.py  | 79 +++++++++++++++++++
 tests/data_model/test_tracks.py               | 26 ++++++
 6 files changed, 193 insertions(+), 23 deletions(-)

diff --git a/src/funtracks/annotators/_regionprops_annotator.py b/src/funtracks/annotators/_regionprops_annotator.py
index e567c266..bee9da68 100644
--- a/src/funtracks/annotators/_regionprops_annotator.py
+++ b/src/funtracks/annotators/_regionprops_annotator.py
@@ -23,6 +23,7 @@
 if TYPE_CHECKING:
     from funtracks.actions import BasicAction
     from funtracks.data_model import Tracks
+    from funtracks.data_model.tracks import PositionUnits
 
 DEFAULT_POS_KEY = "pos"
 DEFAULT_AREA_KEY = "area"
@@ -31,15 +32,17 @@
 DEFAULT_PERIMETER_KEY = "perimeter"
 
 
-def _centroid(mask: Mask) -> list[float]:
-    """Centroid in pixel coordinates, read directly from the mask array.
+def _centroid(mask: Mask, spacing: tuple[float, ...] | None) -> list[float]:
+    """Centroid read directly from the mask array, in pixel or world coordinates.
 
-    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.
+    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 to scale the centroid by, or
+            None to return it in pixel coordinates.
 
     Returns:
         The centroid coordinates, one float per spatial dimension.
@@ -47,7 +50,10 @@ def _centroid(mask: Mask) -> list[float]:
     arr = mask.mask
     bbox_min = mask.bbox[: arr.ndim]
     local = np.array([idx.mean() for idx in np.nonzero(arr)])
-    return [float(v) for v in local + bbox_min]
+    pixel = local + bbox_min
+    if spacing is not None:
+        pixel = pixel * np.asarray(spacing)
+    return [float(v) for v in pixel]
 
 
 class FeatureSpec(NamedTuple):
@@ -74,9 +80,11 @@ 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.
+    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.
@@ -109,6 +117,7 @@ def __init__(
 
         specs = RegionpropsAnnotator._define_features(
             tracks.ndim,
+            position_units=tracks.position_units,
         )
         # update position key in spec
         if self.pos_key != DEFAULT_POS_KEY:
@@ -130,6 +139,7 @@ def __init__(
     def _define_features(
         cls,
         ndim: int,
+        position_units: PositionUnits = "pixel",
     ) -> list[FeatureSpec]:
         """Define all supported regionprops features along with keys and function names.
 
@@ -138,6 +148,10 @@ def _define_features(
 
         Args:
             ndim: Total number of dimensions including time (3 or 4)
+            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,
@@ -147,8 +161,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_pixel"),
+            FeatureSpec(
+                DEFAULT_POS_KEY,
+                position_feature,
+                "centroid" if position_units == "world" else "centroid_pixel",
+            ),
             FeatureSpec(DEFAULT_AREA_KEY, Area(ndim=ndim), "area"),
             # TODO: Add in intensity when image is passed
             # FeatureSpec("intensity", Intensity(ndim=ndim), "intensity"),
@@ -200,6 +222,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 = keys_to_compute == [self.pos_key]
+        # 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
 
         for node_id in self.graph.node_ids():
             if not self.graph.has_node(node_id):
@@ -207,7 +232,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))
+                all_values[self.pos_key].append(_centroid(mask, pos_spacing))
                 continue
             (region,) = regionprops_extended(mask, spacing=spacing)
             for key in keys_to_compute:
diff --git a/src/funtracks/data_model/tracks.py b/src/funtracks/data_model/tracks.py
index 788b42ef..8a6a2b3c 100644
--- a/src/funtracks/data_model/tracks.py
+++ b/src/funtracks/data_model/tracks.py
@@ -5,6 +5,7 @@
 from typing import (
     TYPE_CHECKING,
     Any,
+    Literal,
     TypeAlias,
 )
 from warnings import warn
@@ -42,6 +43,7 @@
 Edge: TypeAlias = tuple[Node, Node]
 AttrValues: TypeAlias = list[AttrValue]
 Attrs: TypeAlias = dict[str, AttrValues]
+PositionUnits: TypeAlias = Literal["pixel", "world"]
 
 logger = logging.getLogger(__name__)
 
@@ -52,12 +54,18 @@ class Tracks:
     The graph nodes represent detections and must have a time attribute and
     position attribute. Edges in the graph represent links across time.
 
-    Coordinate convention: node positions (and segmentation masks/bounding boxes)
-    are stored in pixel coordinates. ``scale`` holds the size of one pixel per
-    dimension, and is what converts those coordinates to world units. Derived
-    measurements that are physically meaningful (area/volume, perimeter, ellipsoid
-    axes) are computed with the scale applied and are therefore already in world
-    units.
+    Coordinate convention: segmentation masks and bounding boxes are always stored
+    in pixel coordinates, since they index directly into the segmentation array.
+    ``scale`` holds the size of one pixel per dimension and is what converts pixel
+    coordinates to world units; it describes the segmentation regardless of which
+    units node positions happen to be in. Derived measurements that are physically
+    meaningful (area/volume, perimeter, ellipsoid axes) are always computed with
+    the scale applied and are therefore in world units.
+
+    Node positions themselves can be stored either way - see ``position_units``.
+    This matters for data that already lives in world units (e.g. a large existing
+    database of point detections): funtracks does not rewrite positions in place
+    just to normalize them to pixel coordinates.
 
     Attributes:
         graph_full (td.graph.BaseGraph): The full graph (first-class): every node/edge
@@ -68,7 +76,12 @@ class Tracks:
         features (FeatureDict): Dictionary of features tracked on graph nodes/edges.
         annotators (AnnotatorRegistry): List of annotators that compute features.
         scale (list[float] | None): The size of one pixel in each dimension,
-            including time. Converts the stored pixel coordinates to world units.
+            including time. Describes the segmentation, independent of
+            ``position_units``; see ``position_scale`` for the scale to apply to
+            node positions specifically.
+        position_units (Literal["pixel", "world"]): Whether node positions (as
+            returned by get_position/get_positions) are stored in pixel or world
+            coordinates. Defaults to "pixel".
         ndim (int): Number of dimensions (3 for 2D+time, 4 for 3D+time).
     """
 
@@ -85,6 +98,7 @@ def __init__(
         scale: list[float] | None = None,
         ndim: int | None = None,
         features: FeatureDict | None = None,
+        position_units: PositionUnits = "pixel",
         _segmentation: GraphArrayView | None = None,
     ):
         """Initialize a Tracks object.
@@ -115,6 +129,11 @@ def __init__(
                 Assumes that all features in the dict already exist on the graph (will
                 be activated but not recomputed). If None, core computed features (pos,
                 tracklet_id) are auto-detected by checking if they exist on the graph.
+            position_units (Literal["pixel", "world"]): Whether the graph's node
+                positions are already in pixel or in world coordinates. Defaults to
+                "pixel". Positions are never rewritten to match this - it only tells
+                funtracks (and a RegionpropsAnnotator, if segmentation is present)
+                which unit the stored/computed positions are in.
             _segmentation (GraphArrayView | None): Internal parameter for reusing an
                 existing GraphArrayView instance. Not intended for public use.
         """
@@ -182,6 +201,7 @@ def __init__(
         else:
             self.segmentation = None
         self.scale = scale
+        self.position_units: PositionUnits = position_units
         self.ndim = self._compute_ndim(
             self.segmentation.shape if self.segmentation is not None else None,
             scale,
@@ -270,6 +290,23 @@ def update_scale(self, scale: list[float] | None) -> None:
 
         self.refresh.emit()
 
+    @property
+    def position_scale(self) -> list[float] | None:
+        """The scale to apply to get_position()/get_positions() output for world units.
+
+        ``scale`` always describes the segmentation (pixel size), regardless of
+        ``position_units``. This is the scale a position consumer should use
+        instead: identity when positions are already in world units (nothing left
+        to convert), or ``scale`` itself when they still need converting.
+
+        A napari points/tracks layer, for instance, should be constructed with
+        ``data=tracks.get_position(...)`` and ``scale=tracks.position_scale`` - never
+        ``tracks.scale`` directly, which would double-scale a world-unit position.
+        """
+        if self.position_units == "world":
+            return None
+        return self.scale
+
     def _get_feature_set(
         self,
         time_attr: str | None,
diff --git a/src/funtracks/features/_feature.py b/src/funtracks/features/_feature.py
index 5891be5f..507ef236 100644
--- a/src/funtracks/features/_feature.py
+++ b/src/funtracks/features/_feature.py
@@ -43,8 +43,10 @@ class Feature(TypedDict):
         scale_dependent (bool): Optional. True if the value is computed with
             ``Tracks.scale`` and therefore has to be recomputed when the scale
             changes (see `Tracks.update_scale`). Features measured in pixel space
-            (`Position`) or derived from graph structure (`TrackletID`, `IoU`)
-            leave this unset.
+            or derived from graph structure (`TrackletID`, `IoU`) leave this unset.
+            `Position` is scale-dependent only when `Tracks.position_units` is
+            "world" (a computed world-unit centroid moves when the scale changes);
+            when it is "pixel" (the default), `Position` leaves this unset too.
     """
 
     feature_type: FeatureType
diff --git a/src/funtracks/features/_node_features.py b/src/funtracks/features/_node_features.py
index 211f318f..ceccdbc1 100644
--- a/src/funtracks/features/_node_features.py
+++ b/src/funtracks/features/_node_features.py
@@ -88,9 +88,10 @@ def SegBbox(ndim: int) -> Feature:
 def Position(axes: Sequence[str]) -> Feature:
     """A feature to hold the position of a node (time not included).
 
-    Positions are always expressed in pixel coordinates, matching how the
-    segmentation is indexed and how coordinates are stored in geff. Multiply by
-    the spatial part of ``Tracks.scale`` to get world units.
+    Positions are expressed in pixel or world coordinates according to
+    ``Tracks.position_units`` (pixel by default). Use ``Tracks.position_scale``
+    (not ``Tracks.scale`` directly) to convert to world units - it is already
+    identity when positions are stored in world units.
 
     Args:
         axes (Sequence[str]): A sequence of the axis names. Used to infer the number
diff --git a/tests/annotators/test_regionprops_annotator.py b/tests/annotators/test_regionprops_annotator.py
index 16f5f98f..eb1fa623 100644
--- a/tests/annotators/test_regionprops_annotator.py
+++ b/tests/annotators/test_regionprops_annotator.py
@@ -171,6 +171,85 @@ def test_centroid_pixel_coords_with_scale(self, get_graph, ndim, only_pos):
             ),
         )
 
+    @pytest.mark.parametrize("only_pos", [True, False])
+    def test_centroid_world_coords_with_position_units_world(
+        self, get_graph, ndim, only_pos
+    ):
+        """With position_units="world", 'pos' must be in world coordinates.
+
+        Same node 6 cube as test_centroid_pixel_coords_with_scale, but with
+        position_units="world": the centroid must come out scaled, through both
+        the fast mask-only path and the skimage regionprops path.
+        """
+        graph = get_graph(ndim, with_seg=True)
+        if ndim == 3:
+            scale = [1.0, 2.0, 3.0]
+            pixel_centroid = np.array([97.5, 97.5])
+        else:
+            scale = [1.0, 2.0, 3.0, 4.0]
+            pixel_centroid = np.array([97.5, 97.5, 97.5])
+
+        tracks = Tracks(
+            graph,
+            ndim=ndim,
+            scale=scale,
+            position_units="world",
+            **track_attrs,
+        )
+        tracks.enable_features(["pos"] if only_pos else ["pos", "area"])
+
+        pos = np.array(tracks.graph_solution.nodes[6]["pos"])
+        world_centroid = pixel_centroid * np.array(scale[1:])
+
+        np.testing.assert_allclose(
+            pos,
+            world_centroid,
+            atol=0.1,
+            err_msg=(
+                f"Centroid must be scaled to world coordinates when "
+                f"position_units='world'. Got {pos}, expected {world_centroid}."
+            ),
+        )
+
+    def test_position_scale_dependent_only_in_world_mode(self, get_graph, ndim):
+        """The 'pos' feature is scale_dependent iff position_units is 'world'.
+
+        area etc. are always scale_dependent; 'pos' should only join them when
+        it is actually derived using the scale (world mode).
+        """
+        graph = get_graph(ndim, with_seg=True)
+
+        pixel_tracks = Tracks(graph, ndim=ndim, **track_attrs)
+        assert pixel_tracks.features["pos"].get("scale_dependent", False) is False
+
+        world_graph = get_graph(ndim, with_seg=True)
+        world_tracks = Tracks(
+            world_graph, ndim=ndim, position_units="world", **track_attrs
+        )
+        assert world_tracks.features["pos"].get("scale_dependent", False) is True
+
+    def test_update_scale_recomputes_position_in_world_mode(self, get_graph, ndim):
+        """Changing tracks.scale must recompute 'pos' when it is in world units."""
+        graph = get_graph(ndim, with_seg=True)
+        scale_a = [1.0] * ndim
+        scale_b = [1.0, 4.0, 4.0] if ndim == 3 else [1.0, 4.0, 4.0, 4.0]
+
+        tracks = Tracks(
+            graph,
+            ndim=ndim,
+            scale=scale_a,
+            position_units="world",
+            **track_attrs,
+        )
+        tracks.enable_features(["pos"])
+        before = np.array(tracks.get_position(6))
+
+        tracks.update_scale(scale_b)
+
+        after = np.array(tracks.get_position(6))
+        ratio = np.array(scale_b[1:]) / np.array(scale_a[1:])
+        np.testing.assert_allclose(after, before * ratio, atol=0.1)
+
     def test_ignores_irrelevant_actions(self, get_graph, ndim):
         """Test that RegionpropsAnnotator ignores actions that don't affect
         segmentation.
diff --git a/tests/data_model/test_tracks.py b/tests/data_model/test_tracks.py
index 0b63e7f0..661406a0 100644
--- a/tests/data_model/test_tracks.py
+++ b/tests/data_model/test_tracks.py
@@ -180,6 +180,32 @@ def test_set_positions_list(graph_2d_list):
     )
 
 
+def test_position_scale_defaults_to_pixel(graph_2d_with_segmentation):
+    """position_units defaults to "pixel", so position_scale is just tracks.scale."""
+    scale = [1.0, 2.0, 3.0]
+    tracks = Tracks(graph_2d_with_segmentation, ndim=3, scale=scale, **track_attrs)
+    assert tracks.position_units == "pixel"
+    assert tracks.position_scale == scale
+    # tracks.scale itself must stay the real segmentation scale either way
+    assert tracks.scale == scale
+
+
+def test_position_scale_is_identity_in_world_mode(graph_2d_with_segmentation):
+    """In world mode positions need no further scaling, but tracks.scale must not
+    change: the segmentation is still pixel-indexed and still needs it."""
+    scale = [1.0, 2.0, 3.0]
+    tracks = Tracks(
+        graph_2d_with_segmentation,
+        ndim=3,
+        scale=scale,
+        position_units="world",
+        **track_attrs,
+    )
+    assert tracks.position_units == "world"
+    assert tracks.position_scale is None
+    assert tracks.scale == scale
+
+
 def test_get_mask_none(graph_2d_with_track_id):
     tracks = Tracks(graph_2d_with_track_id, ndim=3, **track_attrs)
     assert tracks.get_mask(1) is None