diff --git a/src/funtracks/annotators/_regionprops_annotator.py b/src/funtracks/annotators/_regionprops_annotator.py
index bf1ca1ac..e27cbfa4 100644
--- a/src/funtracks/annotators/_regionprops_annotator.py
+++ b/src/funtracks/annotators/_regionprops_annotator.py
@@ -28,6 +28,7 @@
 
     from funtracks.actions import BasicAction
     from funtracks.data_model import Tracks
+    from funtracks.data_model.tracks import PositionUnits
 
     IntensityImage: TypeAlias = np.ndarray | da.Array
 
@@ -40,15 +41,16 @@
 
 
 def _centroid(mask: Mask, spacing: tuple[float, ...] | None) -> list[float]:
-    """Centroid in world units, read directly from the mask array.
+    """Centroid read directly from the mask array, in pixel or world coordinates.
 
-    Equivalent to ``ExtendedRegionProperties.centroid`` (``(local_centroid + bbox_min)
-    * spacing``) but skips the skimage regionprops machinery (find_objects, region
-    caching), which is wasted overhead when only the centroid is needed.
+    Equivalent to ``ExtendedRegionProperties.centroid_pixel``/``.centroid`` but
+    skips the skimage regionprops machinery (find_objects, region caching), which
+    is wasted overhead when only the centroid is needed.
 
     Args:
         mask: A Mask object representing one detection.
-        spacing: Voxel spacing per spatial dimension, or None for unit spacing.
+        spacing: Voxel spacing per spatial dimension to scale the centroid by, or
+            None to return it in pixel coordinates.
 
     Returns:
         The centroid coordinates, one float per spatial dimension.
@@ -56,10 +58,10 @@ 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
+    pixel = local + bbox_min
     if spacing is not None:
-        world = world * np.asarray(spacing)
-    return [float(v) for v in world]
+        pixel = pixel * np.asarray(spacing)
+    return [float(v) for v in pixel]
 
 
 def _bbox_slicing(mask: Mask) -> tuple[slice, ...]:
@@ -146,6 +148,12 @@ class RegionpropsAnnotator(GraphAnnotator):
     - perimeter/surface area
     - mean intensity (requires an intensity image, see ``set_intensity_images``)
 
+    The centroid is stored in pixel or world coordinates according to
+    ``tracks.position_units`` (pixel by default, matching how the segmentation
+    itself is indexed). All size/shape measurements are always computed with
+    ``tracks.scale`` as the voxel spacing and are therefore in world units,
+    regardless of ``position_units``.
+
     Defaults to computing all features, but individual ones can be turned off by changing
     the self.include value at the corresponding index to the feature in self.features.
     """
@@ -211,6 +219,7 @@ def __init__(
         specs = RegionpropsAnnotator._define_features(
             tracks.ndim,
             self.channel_names,
+            position_units=tracks.position_units,
         )
         # Replace the default position spec with the requested key(s): one spec
         # holding the whole coordinate, or one scalar spec per axis. Keyed off the
@@ -355,6 +364,7 @@ def _define_features(
         cls,
         ndim: int,
         channel_names: Sequence[str] | None = None,
+        position_units: PositionUnits = "pixel",
     ) -> list[FeatureSpec]:
         """Define all supported regionprops features along with keys and function names.
 
@@ -365,6 +375,10 @@ def _define_features(
             ndim: Total number of dimensions including time (3 or 4)
             channel_names: Display names for the intensity channels, one per channel.
                 Controls how many values the intensity feature holds.
+            position_units: Whether the position feature should be computed (and
+                recomputed) in pixel or world coordinates. World positions are
+                ``scale_dependent``: they need re-deriving whenever ``tracks.scale``
+                changes, the same as area/perimeter/etc. Pixel positions do not.
 
         Returns:
             list[FeatureSpec]: List of feature specifications with key, feature,
@@ -374,8 +388,16 @@ def _define_features(
         # Default to 3D when ndim is None to enable matching all position columns
         axis_names = ["z", "y", "x"] if ndim is None or ndim == 4 else ["y", "x"]
 
+        position_feature = Position(axes=axis_names)
+        if position_units == "world":
+            position_feature["scale_dependent"] = True
+
         return [
-            FeatureSpec(DEFAULT_POS_KEY, Position(axes=axis_names), "centroid"),
+            FeatureSpec(
+                DEFAULT_POS_KEY,
+                position_feature,
+                "centroid" if position_units == "world" else "centroid_pixel",
+            ),
             FeatureSpec(DEFAULT_AREA_KEY, Area(ndim=ndim), "area"),
             FeatureSpec(
                 DEFAULT_INTENSITY_KEY, Intensity(channel_names), "intensity_mean"
@@ -468,6 +490,9 @@ def compute(self, feature_keys: list[str] | None = None) -> None:
         # regionprops pass anyway and its centroid comes for free, so pos goes through
         # the normal path with everything else.
         fast_pos = set(keys_to_compute) == set(self.pos_keys)
+        # Only scale the fast-path centroid when positions are stored in world units;
+        # centroid_pixel is the pixel-coordinate path and must stay unscaled.
+        pos_spacing = spacing if self.tracks.position_units == "world" else None
 
         node_ids = [
             node_id for node_id in self.graph.node_ids() if self.graph.has_node(node_id)
@@ -491,7 +516,7 @@ def compute(self, feature_keys: list[str] | None = None) -> None:
             mask = self.graph.nodes[node_id]["mask"]
             all_node_ids.append(node_id)
             if fast_pos:
-                centroid = _centroid(mask, spacing)
+                centroid = _centroid(mask, pos_spacing)
                 for key in keys_to_compute:
                     all_values[key].append(self._pos_value(key, centroid))
                 continue
diff --git a/src/funtracks/annotators/_regionprops_extended.py b/src/funtracks/annotators/_regionprops_extended.py
index dce95391..5b289fd2 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):
         """
@@ -233,7 +243,7 @@ def regionprops_extended(
 
     region = mask.regionprops(spacing=spacing, intensity_image=intensity_image)
 
-    bbox_min = np.array(mask.bbox[: mask._mask.ndim]) if spacing is not None else None
+    bbox_min = np.array(mask.bbox[: mask._mask.ndim])
 
     extended_region = ExtendedRegionProperties(
         slice=region.slice,
diff --git a/src/funtracks/candidate_graph/compute_graph.py b/src/funtracks/candidate_graph/compute_graph.py
index 615adae5..818d9d05 100644
--- a/src/funtracks/candidate_graph/compute_graph.py
+++ b/src/funtracks/candidate_graph/compute_graph.py
@@ -18,19 +18,20 @@ def compute_graph_from_seg(
     backend: str = "memory",
 ) -> td.graph.BaseGraph:
     """Construct a candidate graph from a segmentation array. Nodes are placed at the
-    centroid of each segmentation and edges are added for all nodes in adjacent frames
-    within max_edge_distance.
+    centroid of each segmentation (in pixel coordinates) and edges are added for all
+    nodes in adjacent frames within max_edge_distance.
 
     Args:
         segmentation (np.ndarray): A numpy array with integer labels and dimensions
             (t, [z], y, x).
-        max_edge_distance (float): Maximum distance that objects can travel between
-            frames. All nodes with centroids within this distance in adjacent frames
-            will by connected with a candidate edge.
+        max_edge_distance (float): Maximum distance in world units that objects can
+            travel between frames. All nodes with centroids within this distance in
+            adjacent frames will by connected with a candidate edge.
         iou (bool, optional): Whether to include IOU on the candidate graph.
             Defaults to False.
         scale (list[float] | None, optional): The scale of the segmentation data.
-            Will be used to rescale the point locations and attribute computations.
+            Node positions stay in pixel coordinates; the scale is used for
+            measurements in world units (area) and for the distance computations.
             Defaults to None, which implies the data is isotropic.
         t_start (int, optional): The time value to assign to the first frame of the
             segmentation. Frame i will get t = t_start + i. Useful when the
@@ -61,6 +62,7 @@ def compute_graph_from_seg(
         max_edge_distance=max_edge_distance,
         node_frame_dict=node_frame_dict,
         iou_dict=iou_dict,
+        scale=scale,
     )
 
     logger.info("Candidate edges: %d", cand_graph.num_edges())
@@ -83,28 +85,28 @@ def compute_graph_from_points_list(
 
     Args:
         points_list (np.ndarray): An NxD numpy array with N points and D
-            (3 or 4) dimensions. Dimensions should be in order  (t, [z], y, x).
-        max_edge_distance (float): Maximum distance that objects can travel between
-            frames. All nodes with centroids within this distance in adjacent frames
-            will by connected with a candidate edge.
-        scale (list[float] | None, optional): Amount to scale the points in each
-            dimension. Only needed if the provided points are in "voxel" coordinates
-            instead of world coordinates. Defaults to None, which implies the data is
-            isotropic.
+            (3 or 4) dimensions, in pixel coordinates. Dimensions should be in
+            order (t, [z], y, x).
+        max_edge_distance (float): Maximum distance in world units that objects can
+            travel between frames. All nodes with centroids within this distance in
+            adjacent frames will by connected with a candidate edge.
+        scale (list[float] | None, optional): The scale of the data in each dimension
+            (including time). The points are stored unchanged (in pixel coordinates);
+            the scale is only used to measure distances in world units. Defaults to
+            None, which implies the data is isotropic.
         backend (str, optional): Graph backend, "memory" or "sql". Defaults to "memory".
 
     Returns:
         td.graph.BaseGraph: A candidate graph that can be passed to the motile solver.
     """
     # add nodes
-    cand_graph, node_frame_dict = nodes_from_points_list(
-        points_list, scale=scale, backend=backend
-    )
+    cand_graph, node_frame_dict = nodes_from_points_list(points_list, backend=backend)
     logger.info("Candidate nodes: %d", cand_graph.num_nodes())
     # add edges
     add_cand_edges(
         cand_graph,
         max_edge_distance=max_edge_distance,
         node_frame_dict=node_frame_dict,
+        scale=scale,
     )
     return cand_graph
diff --git a/src/funtracks/candidate_graph/utils.py b/src/funtracks/candidate_graph/utils.py
index e000f5a1..77084137 100644
--- a/src/funtracks/candidate_graph/utils.py
+++ b/src/funtracks/candidate_graph/utils.py
@@ -27,8 +27,8 @@ def nodes_from_segmentation(
 
     Each node will have the following attributes:
         - t
-        - pos
-        - area
+        - pos (in pixel coordinates)
+        - area (in world units, i.e. scaled by the voxel size)
         - mask (if mask=True): cropped boolean mask (tracksdata Mask object)
         - bbox (if mask=True): bounding box as [min_0, ..., max_0, ...] int array
 
@@ -39,7 +39,8 @@ def nodes_from_segmentation(
             with funtracks.utils.ensure_unique_labels before calling this function.
         scale (list[float] | None, optional): The scale of the segmentation data in all
             dimensions (including time, which should have a dummy 1 value).
-            Will be used to rescale the point locations and attribute computations.
+            Positions and masks stay in pixel coordinates; the scale is only applied
+            to measurements that are expressed in world units (here: area).
             Defaults to None, which implies the data is isotropic.
         mask (bool, optional): Whether to include mask and bbox attributes for each
             node. Uses regionprop.image (already computed by regionprops at no extra
@@ -67,6 +68,9 @@ def nodes_from_segmentation(
             f"Scale {scale} should have {segmentation.ndim} dims"
         )
 
+    # Voxel volume, used to express the pixel-count area in world units.
+    voxel_volume = float(np.prod(scale[1:]))
+
     node_attributes = ["pos", "area", "mask", "bbox"] if mask else ["pos", "area"]
     cand_graph = create_empty_graph(
         node_attributes=node_attributes,
@@ -83,7 +87,9 @@ def nodes_from_segmentation(
     for t in tqdm(range(len(segmentation)), desc="Extracting nodes from segmentation"):
         segs = segmentation[t]
         nodes_in_frame: list[int] = []
-        props = regionprops(segs, spacing=tuple(scale[1:]))
+        # No spacing: the centroid must stay in pixel coordinates. The area is converted
+        # to world coordinates
+        props = regionprops(segs)
         for regionprop in props:
             node_id = regionprop.label
             if node_id in seen_ids:
@@ -92,7 +98,7 @@ def nodes_from_segmentation(
             attrs: dict[str, Any] = {
                 "t": t + t_start,
                 "pos": list(regionprop.centroid),
-                "area": float(regionprop.area),
+                "area": float(regionprop.area) * voxel_volume,
             }
             if mask:
                 # regionprop.image is the cropped boolean mask within the bbox,
@@ -117,7 +123,6 @@ def nodes_from_segmentation(
 
 def nodes_from_points_list(
     points_list: np.ndarray,
-    scale: list[float] | None = None,
     backend: str = "memory",
 ) -> tuple[td.graph.BaseGraph, dict[int, list[Any]]]:
     """Extract candidate nodes from a list of points. Uses the index of the
@@ -125,13 +130,13 @@ def nodes_from_points_list(
     Returns a tracksdata graph with only nodes, and also a dictionary from frames to
     node_ids for efficient edge adding.
 
+    The points are stored as provided: funtracks keeps positions in pixel
+    coordinates, so the caller should pass pixel coordinates here.
+
     Args:
         points_list (np.ndarray): An NxD numpy array with N points and D
-            (3 or 4) dimensions. Dimensions should be in order (t, [z], y, x).
-        scale (list[float] | None, optional): Amount to scale the points in each
-            dimension (including time). Only needed if the provided points are in
-            "voxel" coordinates instead of world coordinates. Defaults to None, which
-            implies the data is isotropic.
+            (3 or 4) dimensions, in pixel coordinates. Dimensions should be in
+            order (t, [z], y, x).
         backend (str, optional): Graph backend, "memory" or "sql". Defaults to "memory".
 
     Returns:
@@ -142,12 +147,6 @@ def nodes_from_points_list(
 
     ndim = points_list.shape[1]
 
-    if scale is not None:
-        assert len(scale) == ndim, (
-            f"Cannot scale points with {ndim} dims by factor {scale}"
-        )
-        points_list = points_list * np.array(scale)
-
     cand_graph = create_empty_graph(
         node_attributes=["pos"],
         position_attrs=["pos"],
@@ -193,7 +192,11 @@ def _compute_node_frame_dict(cand_graph: td.graph.BaseGraph) -> dict[int, list[A
     return node_frame_dict
 
 
-def create_kdtree(cand_graph: td.graph.BaseGraph, node_ids: list[Any]) -> KDTree:
+def create_kdtree(
+    cand_graph: td.graph.BaseGraph,
+    node_ids: list[Any],
+    spatial_scale: np.ndarray | None = None,
+) -> KDTree:
     """Create a kdtree with the given nodes from the candidate graph.
     Will fail if provided node ids are not in the candidate graph.
 
@@ -203,6 +206,10 @@ def create_kdtree(cand_graph: td.graph.BaseGraph, node_ids: list[Any]) -> KDTree
             include in the KDTree. Useful for limiting to one time frame.
             Must be a list (not a generic iterable) to preserve order for
             correct mapping of query_ball_tree results back to node ids.
+        spatial_scale (np.ndarray | None, optional): Per-spatial-dimension voxel
+            size. Positions are stored in pixel coordinates; multiplying by this
+            puts the tree in world units so distance queries are in world units
+            too. Defaults to None (positions used as-is).
 
     Returns:
         KDTree: A KDTree containing the positions of the given nodes,
@@ -213,7 +220,9 @@ def create_kdtree(cand_graph: td.graph.BaseGraph, node_ids: list[Any]) -> KDTree
     # graph attribute schema on every call, which dominates runtime at scale.
     df = cand_graph.filter(node_ids=node_ids).node_attrs(attr_keys=["node_id", "pos"])
     pos_by_id = dict(zip(df["node_id"].to_list(), df["pos"].to_list(), strict=True))
-    positions = [pos_by_id[nid] for nid in node_ids]
+    positions = np.asarray([pos_by_id[nid] for nid in node_ids], dtype=float)
+    if spatial_scale is not None:
+        positions = positions * spatial_scale
     return KDTree(positions)
 
 
@@ -222,6 +231,7 @@ def add_cand_edges(
     max_edge_distance: float,
     node_frame_dict: None | dict[int, list[Any]] = None,
     iou_dict: dict[int, dict[int, float]] | None = None,
+    scale: list[float] | None = None,
 ) -> None:
     """Add candidate edges to a candidate graph by connecting all nodes in adjacent
     frames that are closer than max_edge_distance.
@@ -230,8 +240,8 @@ def add_cand_edges(
         cand_graph (td.graph.BaseGraph): Candidate graph with only nodes populated.
             Will be modified in-place to add edges.
         max_edge_distance (float): Maximum distance that objects can travel between
-            frames. All nodes within this distance in adjacent frames will by connected
-            with a candidate edge.
+            frames, in world units. All nodes within this distance in adjacent frames
+            will by connected with a candidate edge.
         node_frame_dict (dict[int, list[Any]] | None, optional): A mapping from frames
             to node ids. If not provided, it will be computed from cand_graph. Defaults
             to None.
@@ -240,6 +250,10 @@ def add_cand_edges(
             provided, the iou value is included directly in each edge dict passed to
             bulk_add_edges, avoiding a separate per-edge update pass after insertion.
             Defaults to None (no iou attribute added).
+        scale (list[float] | None, optional): The scale of the data in all dimensions
+            (including time). Node positions are in pixel coordinates, so they are
+            multiplied by the spatial part of the scale before measuring distances,
+            keeping max_edge_distance in world units. Defaults to None (isotropic).
     """
     logger.info("Extracting candidate edges")
     if not node_frame_dict:
@@ -248,14 +262,16 @@ def add_cand_edges(
     if iou_dict is not None:
         cand_graph.add_edge_attr_key("iou", default_value=0.0, dtype=pl.Float64)
 
+    spatial_scale = None if scale is None else np.asarray(scale[1:], dtype=float)
+
     frames = sorted(node_frame_dict.keys())
     prev_node_ids = node_frame_dict[frames[0]]
-    prev_kdtree = create_kdtree(cand_graph, prev_node_ids)
+    prev_kdtree = create_kdtree(cand_graph, prev_node_ids, spatial_scale)
     for frame in tqdm(frames, desc="Adding candidate edges"):
         if frame + 1 not in node_frame_dict:
             continue
         next_node_ids = node_frame_dict[frame + 1]
-        next_kdtree = create_kdtree(cand_graph, next_node_ids)
+        next_kdtree = create_kdtree(cand_graph, next_node_ids, spatial_scale)
 
         matched_indices = prev_kdtree.query_ball_tree(next_kdtree, max_edge_distance)
 
diff --git a/src/funtracks/data_model/tracks.py b/src/funtracks/data_model/tracks.py
index 8c1666fc..1b6dc7dd 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,
     cast,
 )
@@ -50,6 +51,7 @@
 Node: TypeAlias = int
 Edge: TypeAlias = tuple[Node, Node]
 Attrs: TypeAlias = dict[str, list[AttrValue]]
+PositionUnits: TypeAlias = Literal["pixel", "world"]
 
 logger = logging.getLogger(__name__)
 
@@ -60,6 +62,19 @@ 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: 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
             ever known, including soft-deleted (solution=False) candidates. Nodes
@@ -68,7 +83,13 @@ 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. 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 +106,7 @@ def __init__(
         scale: list[float] | None = None,
         ndim: int | None = None,
         features: FeatureDict | None = None,
+        position_units: PositionUnits = "pixel",
         intensity_images: Sequence[IntensityImage] | None = None,
         channel_names: Sequence[str] | None = None,
         _segmentation: GraphArrayView | None = None,
@@ -117,6 +139,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.
             intensity_images (Sequence[IntensityImage] | None): Raw images to measure
                 intensity in, one per channel, each shaped like the segmentation
                 (t, [z], y, x). May be lazy (e.g. dask arrays); each node's bounding
@@ -200,6 +227,7 @@ def __init__(
             self.segmentation = None
         if scale is not 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,
@@ -256,7 +284,8 @@ def scale(self) -> list[float] | None:
 
         Backed by ``graph_full.metadata["scale"]`` (one copy, no staleness), which
         is tracksdata's own convention and is spatial-only - this property adds/
-        strips the dummy time entry at the boundary.
+        strips the dummy time entry at the boundary. Always describes the
+        segmentation's pixel size, independent of ``position_units``.
         """
         spatial_scale = self.graph_full.metadata.get("scale")
         if spatial_scale is None:
@@ -270,6 +299,65 @@ def scale(self, value: list[float] | None) -> None:
         else:
             self.graph_full.metadata["scale"] = list(value[1:])
 
+    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()
+
+    @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 a97323b5..0a5e9bb7 100644
--- a/src/funtracks/features/_feature.py
+++ b/src/funtracks/features/_feature.py
@@ -38,6 +38,13 @@ 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
+            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
@@ -48,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/_node_features.py b/src/funtracks/features/_node_features.py
index 53b1a080..8ec71fb0 100644
--- a/src/funtracks/features/_node_features.py
+++ b/src/funtracks/features/_node_features.py
@@ -88,6 +88,11 @@ 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 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
             of dimensions of the position values.
diff --git a/src/funtracks/features/_regionprops_features.py b/src/funtracks/features/_regionprops_features.py
index 953fbe66..eebdf79f 100644
--- a/src/funtracks/features/_regionprops_features.py
+++ b/src/funtracks/features/_regionprops_features.py
@@ -21,6 +21,7 @@ def Area(ndim: int = 3) -> Feature:
         "num_values": 1,
         "display_name": "Area" if ndim == 3 else "Volume",
         "default_value": None,
+        "scale_dependent": True,
     }
 
 
@@ -80,6 +81,7 @@ def EllipsoidAxes(ndim: int | None = 4) -> Feature:
         "value_names": value_names,
         "default_value": None,
         "spatial_dims": True,
+        "scale_dependent": True,
     }
 
 
@@ -99,6 +101,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
     }
 
 
@@ -118,4 +122,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/src/funtracks/import_export/_tracks_builder.py b/src/funtracks/import_export/_tracks_builder.py
index 31814092..0cb20c39 100644
--- a/src/funtracks/import_export/_tracks_builder.py
+++ b/src/funtracks/import_export/_tracks_builder.py
@@ -16,7 +16,7 @@
 import tracksdata as td
 from geff._typing import InMemoryGeff
 
-from funtracks.data_model.tracks import Tracks
+from funtracks.data_model.tracks import PositionUnits, Tracks
 from funtracks.features import Feature, FeatureDict, Position
 from funtracks.import_export._import_segmentation import (
     load_segmentation,
@@ -117,6 +117,10 @@ def __init__(self) -> None:
         # A complete FeatureDict recovered from the source, when it carries one (geff
         # written by funtracks does). None means the features get inferred instead.
         self.features: FeatureDict | None = None
+        # Whether "pos" ends up in pixel or world coordinates after apply_points_scale.
+        # Pixel by default; a format that actually converts positions (e.g. the GEFF
+        # builder, when the file declares real-world axes) sets this to "world".
+        self.position_units: PositionUnits = "pixel"
 
         # Name maps: {standard_key -> source_property_name(s)}
         # Keys are standard funtracks attribute names (e.g., "time", "pos", "seg_id")
@@ -395,11 +399,12 @@ def apply_points_scale(self) -> None:  # noqa: B027
         """Scale combined ``pos`` values in place, if the source declares a scale.
 
         No-op by default (deliberately not ``@abstractmethod``: this default is a
-        valid, usable behavior, not a contract every subclass must fulfill).
-        Funtracks requires ``pos`` to be in world units; formats that can store
-        points in a different unit alongside an explicit scale (e.g. GEFF's
+        valid, usable behavior, not a contract every subclass must fulfill), which
+        leaves ``self.position_units`` at its default ``"pixel"``. Formats that can
+        store points in a different unit alongside an explicit scale (e.g. GEFF's
         ``axes.scale``) override this to multiply ``pos`` by that scale after
-        :meth:`_combine_multi_value_props` has assembled it.
+        :meth:`_combine_multi_value_props` has assembled it, and to set
+        ``self.position_units = "world"`` when they do.
         """
 
     def validate(self) -> None:
@@ -672,7 +677,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"]
@@ -767,6 +772,15 @@ 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 _retarget_position_key(self, graph: td.graph.BaseGraph) -> None:
         """Point a restored FeatureDict's position_key at the column the graph has.
 
@@ -866,6 +880,13 @@ def build(
         if self.in_memory_geff is None:
             raise ValueError("load_source() must populate self.in_memory_geff")
 
+        # 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(
             self.in_memory_geff["node_props"], self.node_name_map
@@ -909,6 +930,7 @@ def build(
                 ndim=self.ndim,
                 scale=scale,
                 features=self.features,
+                position_units=self.position_units,
             )
         else:
             # The builder always produces a solution, so declare tracklet/lineage
@@ -922,6 +944,7 @@ def build(
                 lineage_attr="lineage_id",
                 ndim=self.ndim,
                 scale=scale,
+                position_units=self.position_units,
             )
 
         # 8. Enable and register features from name maps
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 0662dfab..a9073f31 100644
--- a/src/funtracks/import_export/geff/_export.py
+++ b/src/funtracks/import_export/geff/_export.py
@@ -178,10 +178,23 @@ def _build_geff_metadata(
         else ["time", "space", "space", "space"]
     )
 
-    # Create axes metadata. Points are always in world units, so no "scale" is
-    # written to the geff axes.
+    # Create axes metadata. When position_units is "world" points are already in
+    # world units, so no "scale" is written to the geff axes (world is the
+    # default a foreign reader should assume for an axis with no scale at all).
+    # When "pixel", the segmentation scale is written into axes.scale too, so a
+    # round trip through geff (or a foreign geff reader) can tell these apart
+    # from a genuine unit conversion: apply_points_scale on import recognizes an
+    # axes scale that matches tracks.scale as "this describes the segmentation,
+    # not a conversion" and leaves pos alone.
+    axis_scale = None
+    if tracks.position_units == "pixel" and tracks.scale is not None:
+        axis_scale = dict(zip(axis_names, tracks.scale, strict=True))
     axes = [
-        {"name": name, "type": axis_type}
+        {
+            "name": name,
+            "type": axis_type,
+            **({"scale": axis_scale[name]} if axis_scale is not None else {}),
+        }
         for name, axis_type in zip(axis_names, axis_types, strict=True)
     ]
 
diff --git a/src/funtracks/import_export/geff/_import.py b/src/funtracks/import_export/geff/_import.py
index eb9ee3b5..a16b9d90 100644
--- a/src/funtracks/import_export/geff/_import.py
+++ b/src/funtracks/import_export/geff/_import.py
@@ -215,7 +215,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
@@ -339,16 +339,65 @@ def _axes_scale(self) -> list[float] | None:
             return None
         return [1.0 if s is None else float(s) for s in scales]
 
-    def apply_points_scale(self) -> None:
-        """Multiply ``pos`` by the geff axes scale, so it ends up in world units.
+    def _spatial_axes_are_pixel(self) -> bool | None:
+        """Whether the geff's spatial axes explicitly declare a "pixel" unit.
+
+        Per the geff/OME-Zarr spec, "pixel" is a valid space unit and is the
+        clearest signal available: a file that says so means what it says,
+        regardless of what its ``scale`` values happen to be. Returns None when
+        no spatial axis declares a unit at all, so the caller can fall back to
+        comparing scales instead.
 
-        No-op for legacy funtracks GEFFs (see ``read_header``): their axes.scale
-        is actually a mislabeled segmentation scale, and their ``pos`` is already
-        in world units, so applying it here would double-scale positions.
+        Returns:
+            True if every spatial axis with a declared unit says "pixel", False
+            if any says something else, None if none declare a unit.
+        """
+        geff_axes = getattr(self, "_geff_axes", [])
+        space_axes = [ax for ax in geff_axes if ax.type == "space"]
+        units = [ax.unit for ax in space_axes if ax.unit is not None]
+        if not units:
+            return None
+        return all(unit == "pixel" for unit in units)
+
+    def apply_points_scale(self) -> None:
+        """Decide whether ``pos`` is in pixel or world coordinates, and convert if so.
+
+        No-op for legacy funtracks GEFFs (see ``read_header``): their axes.scale is
+        actually a mislabeled segmentation scale, and their ``pos`` is already in
+        world units (the scale was baked into the stored values long ago), so this
+        marks them "world" without multiplying again.
+
+        Otherwise, "pixel" vs. "world" is decided by, in order:
+        1. The spatial axes' declared ``unit`` (see ``_spatial_axes_are_pixel``),
+           when present - the file says what it means, explicitly.
+        2. Whether the axes scale (``_axes_scale``) matches the segmentation scale
+           (``infer_segmentation_scale``). A geff whose axes scale is exactly the
+           segmentation's voxel size is describing the segmentation, not
+           converting ``pos`` - so ``pos`` is pixel and left unscaled. A axes
+           scale that differs from the segmentation scale (including an axes
+           scale with no segmentation to compare against) is a real unit
+           conversion, and is applied to ``pos``.
+        3. No declared unit and no declared scale at all: nothing suggests these
+           are anything other than pixel coordinates, so "pixel" is left as is.
         """
         if self._is_legacy_funtracks_geff:
+            self.position_units = "world"
             return
+
+        is_pixel = self._spatial_axes_are_pixel()
         scale = self._axes_scale()
+
+        if is_pixel is None:
+            if scale is None:
+                # No unit, no scale declared: nothing suggests these are anything
+                # but pixel coordinates.
+                return
+            seg_scale = self.infer_segmentation_scale()
+            is_pixel = scale == seg_scale
+
+        if is_pixel:
+            return
+
         if scale is None:
             return
         if self.in_memory_geff is None:
@@ -359,6 +408,7 @@ def apply_points_scale(self) -> None:
         # scale is [time, *spatial]; pos only holds the spatial dims. Multiply in
         # float64 to avoid int truncation with scale < 1
         pos["values"] = pos["values"] * np.asarray(scale[1:], dtype=np.float64)
+        self.position_units = "world"
 
     def infer_segmentation_scale(self) -> list[float] | None:
         """Determine ``Tracks.scale`` (segmentation voxel spacing) for this file.
@@ -476,15 +526,17 @@ 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 segmentation voxel scale (``Tracks.scale``) -- this
-            only ever sets the segmentation scale, never the points' scale. If
-            None, defaults to the scale recorded in the geff's tracksdata graph
-            metadata (see :meth:`GeffTracksBuilder.infer_segmentation_scale`),
-            and stays None when the file does not declare one. This is
-            independent of the geff axes' own ``scale``, which (per the geff
-            spec) describes how to convert stored positions to world units and
-            is always applied to ``pos`` on import -- regardless of what is
-            passed here -- so that funtracks' world-units invariant holds.
+        scale: Optional segmentation voxel scale (``Tracks.scale``) -- this only
+            ever sets the segmentation scale, never the points' scale. If None,
+            defaults to the scale recorded in the geff's tracksdata graph metadata
+            (see :meth:`GeffTracksBuilder.infer_segmentation_scale`), and stays
+            None when the file does not declare one. This is independent of the
+            geff axes' own ``scale``: per the geff spec, that describes how to
+            convert stored positions to world units, and is applied to ``pos`` on
+            import when it actually differs from the segmentation scale (see
+            :meth:`GeffTracksBuilder.apply_points_scale`) -- resulting
+            ``Tracks.position_units`` reflects whether that happened ("world") or
+            not ("pixel", the default when nothing on the file says otherwise).
         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/annotators/test_regionprops_annotator.py b/tests/annotators/test_regionprops_annotator.py
index 3a82cfe0..0fde70ed 100644
--- a/tests/annotators/test_regionprops_annotator.py
+++ b/tests/annotators/test_regionprops_annotator.py
@@ -416,15 +416,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:
@@ -436,25 +438,101 @@ 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"])
+
+        world_value = pixel_centroid * np.array(scale[1:])
+        np.testing.assert_allclose(
+            pos,
+            pixel_centroid,
+            atol=0.1,
+            err_msg=(
+                f"Centroid must be stored in pixel coordinates. "
+                f"Got {pos}, expected {pixel_centroid}. "
+                f"World coordinates would be {world_value}"
+            ),
+        )
+
+    @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"])
-        expected = pixel_centroid * np.array(scale[1:])
+        world_centroid = 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))
         np.testing.assert_allclose(
             pos,
-            expected,
+            world_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 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)
+
     @pytest.mark.parametrize("split_pos", [False, True])
     def test_position_follows_segmentation_edit(self, get_graph, ndim, split_pos):
         """Editing a mask must move the node's position, in whichever columns
diff --git a/tests/candidate_graph/test_compute_graph.py b/tests/candidate_graph/test_compute_graph.py
index 21fa0fca..cd09b94b 100644
--- a/tests/candidate_graph/test_compute_graph.py
+++ b/tests/candidate_graph/test_compute_graph.py
@@ -203,11 +203,13 @@ def test_graph_from_points_list(backend):
     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..af35f64e 100644
--- a/tests/candidate_graph/test_graph_utils.py
+++ b/tests/candidate_graph/test_graph_utils.py
@@ -35,14 +35,14 @@ def test_nodes_from_segmentation_2d(get_tracks):
     assert node_frame_dict[0] == [1]
     assert Counter(node_frame_dict[1]) == Counter([2, 3])
 
-    # test with scaling
+    # 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, 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])
@@ -64,14 +64,14 @@ def test_nodes_from_segmentation_3d(get_tracks):
     assert node_frame_dict[0] == [1]
     assert Counter(node_frame_dict[1]) == Counter([2, 3])
 
-    # test with scaling
+    # 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, 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 d458e114..8a57a3ee 100644
--- a/tests/data_model/test_solution_tracks.py
+++ b/tests/data_model/test_solution_tracks.py
@@ -173,7 +173,7 @@ 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
+    # Lineage ID, Area, and the per-axis scale columns
     header = [
         "ID",
         "Parent ID",
@@ -183,6 +183,8 @@ def test_export_to_csv_with_display_names(
         "Tracklet ID",
         "Lineage ID",
         "Area",
+        "y_scale",
+        "x_scale",
     ]
     assert lines[0].strip().split(",") == header
 
@@ -197,7 +199,7 @@ 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
+    # Tracklet ID, Lineage ID, Volume, and the per-axis scale columns
     header = [
         "ID",
         "Parent ID",
@@ -208,6 +210,9 @@ def test_export_to_csv_with_display_names(
         "Tracklet ID",
         "Lineage ID",
         "Volume",
+        "z_scale",
+        "y_scale",
+        "x_scale",
     ]
     assert lines[0].strip().split(",") == header
 
diff --git a/tests/data_model/test_tracks.py b/tests/data_model/test_tracks.py
index 73e4e10f..ae65802e 100644
--- a/tests/data_model/test_tracks.py
+++ b/tests/data_model/test_tracks.py
@@ -195,6 +195,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
@@ -378,6 +404,109 @@ def test_update_mask_syncs_bbox(graph_2d_with_segmentation):
     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"]]
+
+
 def test_has_track_id_at_time(get_tracks):
     """Track 3 spans t=1, 2 and 4, but not t=0 or t=3."""
 
diff --git a/tests/import_export/test_csv_export.py b/tests/import_export/test_csv_export.py
index ae13a721..9ce488ec 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,49 @@ 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()
diff --git a/tests/import_export/test_export_to_geff.py b/tests/import_export/test_export_to_geff.py
index b36353a6..57b0f86d 100644
--- a/tests/import_export/test_export_to_geff.py
+++ b/tests/import_export/test_export_to_geff.py
@@ -47,11 +47,6 @@ def _assert_valid_geff_export(export_dir, expected_num_nodes=None):
     attrs = dict(z.attrs)
     assert "geff" in attrs
     assert "axes" in attrs["geff"]
-    # Tracks.scale (segmentation spacing) must never be written to axes.scale:
-    # per the geff spec that field means "multiply pos by this to get world
-    # units", which pos already is. It belongs in graph metadata (_seg_scale).
-    for ax in attrs["geff"]["axes"]:
-        assert ax["scale"] is None
 
     if expected_num_nodes is not None:
         assert len(z["nodes/ids"][:]) == expected_num_nodes
@@ -357,12 +352,19 @@ def test_export_metadata(get_tracks, ndim, with_seg, tmp_path):
         assert seg_shape is None
 
     # Tracks.scale (segmentation spacing) is written to graph metadata, mirroring
-    # shape, never to axes.scale. tracksdata's own "scale" metadata is
-    # spatial-only (no time), unlike Tracks.scale (time first, dummy 1.0).
+    # shape. tracksdata's own "scale" metadata is spatial-only (no time), unlike
+    # Tracks.scale (time first, dummy 1.0).
     seg_scale = td.io.read_graph_metadata(export_dir / "tracks.geff").get("scale")
     assert seg_scale is not None
     assert list(seg_scale) == list(tracks.scale)[1:]
 
+    # position_units defaults to "pixel", so the same scale is also written into
+    # axes.scale: a foreign reader needs it there to reconstruct world units, even
+    # though funtracks itself keeps positions in pixel coordinates.
+    assert tracks.position_units == "pixel"
+    space_axes = [ax for ax in axes if ax["type"] == "space"]
+    assert [ax["scale"] for ax in space_axes] == list(tracks.scale)[1:]
+
     funtracks_extra = _funtracks_extra(export_dir / "tracks.geff")
     assert funtracks_extra is not None
     assert "features" in funtracks_extra
@@ -430,6 +432,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_geff_scale_roundtrip.py b/tests/import_export/test_geff_scale_roundtrip.py
index 6de60826..5feab6dc 100644
--- a/tests/import_export/test_geff_scale_roundtrip.py
+++ b/tests/import_export/test_geff_scale_roundtrip.py
@@ -1,11 +1,16 @@
 """Regression tests for the geff scale round-trip.
 
 Two scales, two homes: ``Tracks.scale`` (segmentation voxel spacing) lives in
-``graph.metadata["scale"]`` and is never written to axes.scale; GEFF's
-per-axis ``scale`` means "multiply pos by this to get world units" and is
-always applied to ``pos`` on import, independent of ``Tracks.scale``. GEFFs
-written by pre-fix funtracks are the one exception -- see
-``test_legacy_funtracks_geff_scale_migrates_without_double_scaling_points``.
+``graph.metadata["scale"]``. GEFF's per-axis ``scale`` means "multiply pos by
+this to get world units" per the geff spec, but funtracks also writes it equal
+to ``Tracks.scale`` when ``position_units`` is "pixel" (the default), so a
+foreign reader has a scale to work with even though funtracks itself keeps
+positions in pixel coordinates. On import, an axes scale that matches the
+segmentation scale is recognized as describing the segmentation rather than
+converting points, and is left unapplied to ``pos`` (``position_units`` stays
+"pixel"); a genuinely different axes scale is applied and ``position_units``
+becomes "world". GEFFs written by pre-fix funtracks are a further exception --
+see ``test_legacy_funtracks_geff_scale_migrates_without_double_scaling_points``.
 """
 
 import numpy as np
@@ -36,22 +41,32 @@ def test_scale_survives_export_to_geff_roundtrip(get_tracks, tmp_path):
     assert loaded.scale == [1.0, 2.0, 0.25, 0.5]
 
 
-def test_scale_not_written_to_axes(get_tracks, tmp_path):
-    """Tracks.scale (segmentation spacing) must not leak into axes.scale,
-    which per the geff spec means something else: how to convert pos to
-    world units."""
+def test_pixel_scale_written_to_axes_and_not_double_applied(get_tracks, tmp_path):
+    """Tracks.scale (segmentation spacing) is also written to axes.scale when
+    position_units is "pixel" (the default), so a foreign reader has a scale to
+    work with. On import it must be recognized as describing the segmentation,
+    not a unit conversion, and left unapplied to pos - it round-trips as pixel,
+    not silently promoted to world."""
     tracks = get_tracks(ndim=3, with_seg=False, prefill_track_ids=True)
     tracks.scale = [1.0, 0.25, 0.5]
+    assert tracks.position_units == "pixel"
+    node_ids = sorted(tracks.graph_solution.node_ids())
+    original_pos = tracks.get_positions(node_ids)
 
     geff_path = tmp_path / "tracks.geff"
     write_to_geff(tracks, geff_path)
 
     meta = GeffMetadata.read(geff_path)
-    assert all(ax.scale is None for ax in meta.axes)
+    space_axes = [ax for ax in meta.axes if ax.type == "space"]
+    assert [ax.scale for ax in space_axes] == [0.25, 0.5]
     # tracksdata's own "scale" metadata is spatial-only (no time), unlike
     # Tracks.scale (time first, dummy 1.0).
     assert meta.extra["tracksdata"]["scale"] == [0.25, 0.5]
 
+    loaded = import_from_geff(geff_path)
+    assert loaded.position_units == "pixel"
+    np.testing.assert_allclose(loaded.get_positions(node_ids), original_pos)
+
 
 def test_explicit_scale_overrides_metadata(get_tracks, tmp_path):
     """A caller-supplied scale must still win over what the file says."""
diff --git a/tests/import_export/test_import_from_geff.py b/tests/import_export/test_import_from_geff.py
index bd5555f9..deac8b23 100644
--- a/tests/import_export/test_import_from_geff.py
+++ b/tests/import_export/test_import_from_geff.py
@@ -72,12 +72,12 @@ 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 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):
-        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
 
 
@@ -335,7 +335,9 @@ def test_tracks_with_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
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)