diff --git a/src/cedalion/geometry/landmarks.py b/src/cedalion/geometry/landmarks.py
index 038a577b..37e25ba7 100644
--- a/src/cedalion/geometry/landmarks.py
+++ b/src/cedalion/geometry/landmarks.py
@@ -8,6 +8,7 @@
 import xarray as xr
 import pyvista as pv
 
+import cedalion
 import cedalion.vis.blocks as vbx
 from cedalion import cite
 
@@ -198,6 +199,11 @@ def __init__(self, scalp_surface: Surface, landmarks: LabeledPoints):
             landmarks (LabeledPoints): positions of "Nz", "Iz", "LPA", "RPA"
         """
         cite("Oostenveld2001")
+        if scalp_surface.units != cedalion.units.mm:
+            raise NotImplementedError(
+                "LandmarksBuilder1010 currently requires a scalp surface in "
+                f"mm; got units={scalp_surface.units}."
+            )
         if isinstance(scalp_surface, TrimeshSurface):
             scalp_surface = VTKSurface.from_trimeshsurface(scalp_surface)
 
diff --git a/src/cedalion/geometry/photogrammetry/processors.py b/src/cedalion/geometry/photogrammetry/processors.py
index 469dbe30..6bd2c0c1 100644
--- a/src/cedalion/geometry/photogrammetry/processors.py
+++ b/src/cedalion/geometry/photogrammetry/processors.py
@@ -194,7 +194,11 @@ def process(
             surface normal vectors
             aux. detail object if detail == True
         """
-        assert surface.units == units.mm  # FIXME Einstar yields mm. Allow other units.
+        if surface.units != units.mm:
+            raise NotImplementedError(
+                "ColoredStickerProcessor currently supports mm-scale "
+                f"surfaces only; got units={surface.units}."
+            )
 
         vertex_colors, h, s, v = self._extract_vertex_color_data(surface)
 
diff --git a/src/cedalion/io/photogrammetry.py b/src/cedalion/io/photogrammetry.py
index ee0e944f..9609ae52 100644
--- a/src/cedalion/io/photogrammetry.py
+++ b/src/cedalion/io/photogrammetry.py
@@ -1,11 +1,12 @@
 """Module for reading photogrammetry output file formats."""
 
+import cedalion
 import cedalion.dataclasses as cdc
 import numpy as np
 from collections import OrderedDict
 
 
-def read_photogrammetry_einstar(fn):
+def read_photogrammetry_einstar(fn, units=cedalion.units.mm):
     """Read optodes and fiducials from photogrammetry pipeline.
 
     This method reads the output file as returned by the
@@ -13,6 +14,9 @@ def read_photogrammetry_einstar(fn):
 
     Args:
         fn (str): The filename of the einstar photogrammetry output file.
+        units: Units of the coordinates in the file. Defaults to millimetres
+            (Einstar convention). Pass e.g. ``cedalion.units.m`` for exports
+            in metres.
 
     Returns:
         tuple: A tuple containing:
@@ -24,7 +28,7 @@ def read_photogrammetry_einstar(fn):
     fiducials, optodes = read_einstar(fn)
     fiducials = OrderedDict((k, v) for k, v in fiducials.items() if v != [])
     optodes = OrderedDict((k, v) for k, v in optodes.items() if v != [])
-    fiducials, optodes = opt_fid_to_xr(fiducials, optodes)
+    fiducials, optodes = opt_fid_to_xr(fiducials, optodes, units=units)
     return fiducials, optodes
 
 
@@ -53,12 +57,13 @@ def read_einstar(fn):
     return fiducials, optodes
 
 
-def opt_fid_to_xr(fiducials, optodes):
+def opt_fid_to_xr(fiducials, optodes, units=cedalion.units.mm):
     """Convert OrderedDicts fiducials and optodes to cedalion LabeledPoints objects.
 
     Args:
         fiducials (OrderedDict): The fiducials as an OrderedDict.
         optodes (OrderedDict): The optodes as an OrderedDict.
+        units: Units to attach to the returned LabeledPoints. Defaults to mm.
 
     Returns:
         tuple: A tuple containing:
@@ -81,8 +86,8 @@ def opt_fid_to_xr(fiducials, optodes):
         opt_coords = np.array(list(optodes.values()))
 
     fiducials = cdc.build_labeled_points(
-        fidu_coords, labels=list(fiducials.keys()), crs=CRS
-    )  # , units="mm")
+        fidu_coords, labels=list(fiducials.keys()), crs=CRS, units=units,
+    )
 
     types = []
     for lab in list(optodes.keys()):
@@ -93,6 +98,7 @@ def opt_fid_to_xr(fiducials, optodes):
         else:
             types.append(cdc.PointType(0))
     optodes = cdc.build_labeled_points(
-        opt_coords, labels=list(optodes.keys()), crs=CRS, types=types
+        opt_coords, labels=list(optodes.keys()), crs=CRS, types=types,
+        units=units,
     )
     return fiducials, optodes
diff --git a/src/cedalion/io/probe_geometry.py b/src/cedalion/io/probe_geometry.py
index 561e43d2..19e650ed 100644
--- a/src/cedalion/io/probe_geometry.py
+++ b/src/cedalion/io/probe_geometry.py
@@ -249,17 +249,25 @@ def read_digpts(fname: str, units: str="mm") -> xr.DataArray:
     return result
 
 
-def read_einstar_obj(fname: str) -> TrimeshSurface:
-    """Read a textured triangle mesh generated by Einstar 3D scanners.
+def read_einstar_obj(
+    fname: str,
+    crs: str = "digitized",
+    units=cedalion.units.mm,
+) -> TrimeshSurface:
+    """Read a textured triangle mesh (Einstar, Scaniverse, ...).
 
     Args:
         fname: Path to the ``.obj`` mesh file.
+        crs: Coordinate reference system tag. Defaults to ``"digitized"``.
+        units: Units of the vertex coordinates in the file. Defaults to
+            millimetres (Einstar convention). Pass e.g. ``cedalion.units.m``
+            for Scaniverse exports.
 
     Returns:
-        TrimeshSurface in the ``"digitized"`` CRS with millimetre units.
+        TrimeshSurface in the given CRS and units.
     """
     mesh = trimesh.load(fname)
-    return TrimeshSurface(mesh, crs="digitized", units=cedalion.units.mm)
+    return TrimeshSurface(mesh, crs=crs, units=units)
 
 
 def read_fieldtrip_elc(fname : Path | str) -> cdt.LabeledPoints:
diff --git a/src/cedalion/vis/anatomy/optode_selector.py b/src/cedalion/vis/anatomy/optode_selector.py
index b382211a..db491148 100644
--- a/src/cedalion/vis/anatomy/optode_selector.py
+++ b/src/cedalion/vis/anatomy/optode_selector.py
@@ -1,4 +1,5 @@
 import pyvista as pv
+import cedalion
 from cedalion.dataclasses import PointType
 import cedalion.dataclasses as cdc
 import numpy as np
@@ -35,6 +36,20 @@ class OptodeSelector:
         - Masha Iudina | mashayudi@gmail.com | 2024
     """
     def __init__(self, surface, points, normals=None, plotter=None, labels = None):
+        if surface.units != cedalion.units.mm:
+            raise NotImplementedError(
+                "OptodeSelector currently requires a surface in mm; got "
+                f"units={surface.units}. Convert the surface to mm before "
+                "picking (e.g. `surface.mesh.apply_scale(factor); "
+                "surface.units = cedalion.units.mm`)."
+            )
+        if points.pint.units != cedalion.units.mm:
+            raise NotImplementedError(
+                "OptodeSelector currently requires points in mm; got "
+                f"units={points.pint.units}. Convert with "
+                "`points = points.pint.to('mm')` first."
+            )
+
         self.points = points
         self.normals = normals
         self.surface = surface
@@ -65,8 +80,6 @@ def plot(self):
             PointType.ELECTRODE: 3,
         }
 
-        # points = points.pint.to("mm").pint.dequantify()  # FIXME unit handling
-        # points = points.pint.dequantify()  # FIXME unit handling
         for type, x in points.groupby("type"):
             for i_point in range(len(x)):
 
diff --git a/src/cedalion/vis/blocks.py b/src/cedalion/vis/blocks.py
index 0f5eb88b..3875676e 100644
--- a/src/cedalion/vis/blocks.py
+++ b/src/cedalion/vis/blocks.py
@@ -14,6 +14,7 @@
 from vtk.util.numpy_support import numpy_to_vtk
 from numpy.typing import ArrayLike
 
+import cedalion
 import cedalion.dataclasses as cdc
 import cedalion.typing as cdt
 from cedalion.dataclasses import PointType
@@ -162,6 +163,14 @@ def plot_surface(
         - Masha Iudina | mashayudi@gmail.com | 2024
     """
 
+    if pick_landmarks and surface.units != cedalion.units.mm:
+        raise NotImplementedError(
+            "plot_surface(pick_landmarks=...) currently requires a surface "
+            f"in mm; got units={surface.units}. Convert the surface to mm "
+            "before picking (e.g. `surface.mesh.apply_scale(factor); "
+            "surface.units = cedalion.units.mm`)."
+        )
+
     if isinstance(surface, cdc.VTKSurface):
         mesh = surface.mesh
     elif isinstance(surface, cdc.TrimeshSurface):
diff --git a/tests/test_landmarks.py b/tests/test_landmarks.py
index 5f7e8c76..c9867882 100644
--- a/tests/test_landmarks.py
+++ b/tests/test_landmarks.py
@@ -7,6 +7,7 @@
 import trimesh
 import xarray as xr
 
+import cedalion
 import cedalion.dataclasses as cdc
 from cedalion.geometry.landmarks import (
     LandmarksBuilder1010,
@@ -262,7 +263,7 @@ def _icosphere_scalp():
     sphere = trimesh.creation.icosphere(subdivisions=4, radius=90.0)
     # squash into a vaguely head-shaped ellipsoid (a, b, c)
     sphere.vertices = sphere.vertices * np.array([85.0, 100.0, 90.0]) / 90.0
-    surf = cdc.TrimeshSurface(sphere, crs="ras", units="mm")
+    surf = cdc.TrimeshSurface(sphere, crs="ras", units=cedalion.units.mm)
 
     landmarks = cdc.build_labeled_points(
         [