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..b448719e 100644
--- a/src/funtracks/candidate_graph/utils.py
+++ b/src/funtracks/candidate_graph/utils.py
@@ -26,8 +26,8 @@ def nodes_from_segmentation(
 
     Each node will have the following attributes:
         - t
-        - pos
-        - area
+        - pos (in pixel coordinates)
+        - area (in world units, i.e. scaled by the voxel size)
         - mask (if mask=True): cropped boolean mask (tracksdata Mask object)
         - bbox (if mask=True): bounding box as [min_0, ..., max_0, ...] int array
 
@@ -38,7 +38,8 @@ def nodes_from_segmentation(
             with funtracks.utils.ensure_unique_labels before calling this function.
         scale (list[float] | None, optional): The scale of the segmentation data in all
             dimensions (including time, which should have a dummy 1 value).
-            Will be used to rescale the point locations and attribute computations.
+            Positions and masks stay in pixel coordinates; the scale is only applied
+            to measurements that are expressed in world units (here: area).
             Defaults to None, which implies the data is isotropic.
         mask (bool, optional): Whether to include mask and bbox attributes for each
             node. Uses regionprop.image (already computed by regionprops at no extra
@@ -65,6 +66,9 @@ def nodes_from_segmentation(
             f"Scale {scale} should have {segmentation.ndim} dims"
         )
 
+    # Voxel volume, used to express the pixel-count area in world units.
+    voxel_volume = float(np.prod(scale[1:]))
+
     node_attributes = ["pos", "area", "mask", "bbox"] if mask else ["pos", "area"]
     cand_graph = create_empty_graph(
         node_attributes=node_attributes,
@@ -80,7 +84,9 @@ def nodes_from_segmentation(
     for t in tqdm(range(len(segmentation)), desc="Extracting nodes from segmentation"):
         segs = segmentation[t]
         nodes_in_frame: list[int] = []
-        props = regionprops(segs, spacing=tuple(scale[1:]))
+        # No spacing: the centroid must stay in pixel coordinates. The area is converted
+        # to world coordinates
+        props = regionprops(segs)
         for regionprop in props:
             node_id = regionprop.label
             if node_id in seen_ids:
@@ -89,7 +95,7 @@ def nodes_from_segmentation(
             attrs: dict[str, Any] = {
                 "t": t + t_start,
                 "pos": list(regionprop.centroid),
-                "area": float(regionprop.area),
+                "area": float(regionprop.area) * voxel_volume,
             }
             if mask:
                 # regionprop.image is the cropped boolean mask within the bbox,
@@ -114,20 +120,19 @@ def nodes_from_segmentation(
 
 def nodes_from_points_list(
     points_list: np.ndarray,
-    scale: list[float] | None = None,
 ) -> tuple[td.graph.BaseGraph, dict[int, list[Any]]]:
     """Extract candidate nodes from a list of points. Uses the index of the
     point in the list as its unique id.
     Returns a tracksdata graph with only nodes, and also a dictionary from frames to
     node_ids for efficient edge adding.
 
+    The points are stored as provided: funtracks keeps positions in pixel
+    coordinates, so the caller should pass pixel coordinates here.
+
     Args:
         points_list (np.ndarray): An NxD numpy array with N points and D
-            (3 or 4) dimensions. Dimensions should be in order (t, [z], y, x).
-        scale (list[float] | None, optional): Amount to scale the points in each
-            dimension (including time). Only needed if the provided points are in
-            "voxel" coordinates instead of world coordinates. Defaults to None, which
-            implies the data is isotropic.
+            (3 or 4) dimensions, in pixel coordinates. Dimensions should be in
+            order (t, [z], y, x).
 
     Returns:
         tuple[td.graph.BaseGraph, dict[int, list[Any]]]: A candidate graph with only
@@ -137,12 +142,6 @@ def nodes_from_points_list(
 
     ndim = points_list.shape[1]
 
-    if scale is not None:
-        assert len(scale) == ndim, (
-            f"Cannot scale points with {ndim} dims by factor {scale}"
-        )
-        points_list = points_list * np.array(scale)
-
     cand_graph = create_empty_graph(
         node_attributes=["pos"],
         position_attrs=["pos"],
@@ -187,7 +186,11 @@ def _compute_node_frame_dict(cand_graph: td.graph.BaseGraph) -> dict[int, list[A
     return node_frame_dict
 
 
-def create_kdtree(cand_graph: td.graph.BaseGraph, node_ids: list[Any]) -> KDTree:
+def create_kdtree(
+    cand_graph: td.graph.BaseGraph,
+    node_ids: list[Any],
+    spatial_scale: np.ndarray | None = None,
+) -> KDTree:
     """Create a kdtree with the given nodes from the candidate graph.
     Will fail if provided node ids are not in the candidate graph.
 
@@ -197,6 +200,10 @@ def create_kdtree(cand_graph: td.graph.BaseGraph, node_ids: list[Any]) -> KDTree
             include in the KDTree. Useful for limiting to one time frame.
             Must be a list (not a generic iterable) to preserve order for
             correct mapping of query_ball_tree results back to node ids.
+        spatial_scale (np.ndarray | None, optional): Per-spatial-dimension voxel
+            size. Positions are stored in pixel coordinates; multiplying by this
+            puts the tree in world units so distance queries are in world units
+            too. Defaults to None (positions used as-is).
 
     Returns:
         KDTree: A KDTree containing the positions of the given nodes,
@@ -207,7 +214,9 @@ def create_kdtree(cand_graph: td.graph.BaseGraph, node_ids: list[Any]) -> KDTree
     # graph attribute schema on every call, which dominates runtime at scale.
     df = cand_graph.filter(node_ids=node_ids).node_attrs(attr_keys=["node_id", "pos"])
     pos_by_id = dict(zip(df["node_id"].to_list(), df["pos"].to_list(), strict=True))
-    positions = [pos_by_id[nid] for nid in node_ids]
+    positions = np.asarray([pos_by_id[nid] for nid in node_ids], dtype=float)
+    if spatial_scale is not None:
+        positions = positions * spatial_scale
     return KDTree(positions)
 
 
@@ -216,6 +225,7 @@ def add_cand_edges(
     max_edge_distance: float,
     node_frame_dict: None | dict[int, list[Any]] = None,
     iou_dict: dict[int, dict[int, float]] | None = None,
+    scale: list[float] | None = None,
 ) -> None:
     """Add candidate edges to a candidate graph by connecting all nodes in adjacent
     frames that are closer than max_edge_distance.
@@ -224,8 +234,8 @@ def add_cand_edges(
         cand_graph (td.graph.BaseGraph): Candidate graph with only nodes populated.
             Will be modified in-place to add edges.
         max_edge_distance (float): Maximum distance that objects can travel between
-            frames. All nodes within this distance in adjacent frames will by connected
-            with a candidate edge.
+            frames, in world units. All nodes within this distance in adjacent frames
+            will by connected with a candidate edge.
         node_frame_dict (dict[int, list[Any]] | None, optional): A mapping from frames
             to node ids. If not provided, it will be computed from cand_graph. Defaults
             to None.
@@ -234,6 +244,10 @@ def add_cand_edges(
             provided, the iou value is included directly in each edge dict passed to
             bulk_add_edges, avoiding a separate per-edge update pass after insertion.
             Defaults to None (no iou attribute added).
+        scale (list[float] | None, optional): The scale of the data in all dimensions
+            (including time). Node positions are in pixel coordinates, so they are
+            multiplied by the spatial part of the scale before measuring distances,
+            keeping max_edge_distance in world units. Defaults to None (isotropic).
     """
     logger.info("Extracting candidate edges")
     if not node_frame_dict:
@@ -242,14 +256,16 @@ def add_cand_edges(
     if iou_dict is not None:
         cand_graph.add_edge_attr_key("iou", default_value=0.0, dtype=pl.Float64)
 
+    spatial_scale = None if scale is None else np.asarray(scale[1:], dtype=float)
+
     frames = sorted(node_frame_dict.keys())
     prev_node_ids = node_frame_dict[frames[0]]
-    prev_kdtree = create_kdtree(cand_graph, prev_node_ids)
+    prev_kdtree = create_kdtree(cand_graph, prev_node_ids, spatial_scale)
     for frame in tqdm(frames, desc="Adding candidate edges"):
         if frame + 1 not in node_frame_dict:
             continue
         next_node_ids = node_frame_dict[frame + 1]
-        next_kdtree = create_kdtree(cand_graph, next_node_ids)
+        next_kdtree = create_kdtree(cand_graph, next_node_ids, spatial_scale)
 
         matched_indices = prev_kdtree.query_ball_tree(next_kdtree, max_edge_distance)
 
diff --git a/src/funtracks/data_model/tracks.py b/src/funtracks/data_model/tracks.py
index 4b6fb278..788b42ef 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).
     """
 
@@ -220,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/_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/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/src/funtracks/import_export/_tracks_builder.py b/src/funtracks/import_export/_tracks_builder.py
index 65681688..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"]
@@ -746,6 +746,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 build(
         self,
         source: Path | pd.DataFrame,
@@ -813,6 +822,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
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/_import.py b/src/funtracks/import_export/geff/_import.py
index 923b1248..93487bff 100644
--- a/src/funtracks/import_export/geff/_import.py
+++ b/src/funtracks/import_export/geff/_import.py
@@ -214,7 +214,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
@@ -354,6 +354,30 @@ 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. 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 import_from_geff(
     directory: Path,
@@ -375,7 +399,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/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 9dcb7a03..bf66254f 100644
--- a/tests/candidate_graph/test_compute_graph.py
+++ b/tests/candidate_graph/test_compute_graph.py
@@ -192,11 +192,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..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 a277eedc..1831c0fe 100644
--- a/tests/data_model/test_solution_tracks.py
+++ b/tests/data_model/test_solution_tracks.py
@@ -154,7 +154,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",
@@ -164,6 +164,8 @@ def test_export_to_csv_with_display_names(
         "Tracklet ID",
         "Lineage ID",
         "Area",
+        "y_scale",
+        "x_scale",
     ]
     assert lines[0].strip().split(",") == header
 
@@ -178,7 +180,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",
@@ -189,6 +191,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 ca7bbf73..0b63e7f0 100644
--- a/tests/data_model/test_tracks.py
+++ b/tests/data_model/test_tracks.py
@@ -359,3 +359,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_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 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 12bbeac1..1dac339d 100644
--- a/tests/import_export/test_import_from_geff.py
+++ b/tests/import_export/test_import_from_geff.py
@@ -71,12 +71,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
 
 
@@ -329,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
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)