Physical scale: detect the audited shape (no unit / [1,1,1]) and use the manifest pitch

lasagna/scripts/download_omezarr.py

@@ -446,7 +446,67 @@ def _discover_levels(bucket: str, prefix: str, anon: bool) -> list[int]:
     return sorted(levels)
 
 
-def _parse_multiscales(zattrs: dict) -> dict[int, list[float]]:
+def _axes_have_units(zattrs: dict) -> bool:
+    """True when every spatial axis in the OME metadata declares a unit.
+
+    This is the signal the catalog audit found: the affected stores publish an explicit scale of
+    [1, 1, 1] and omit `unit` on every axis, so a check that only looks for a *missing* scale
+    transform never fires on them (review on ScrollPrize/villa#2137).
+    """
+    ms_list = zattrs.get("multiscales") or []
+    if not ms_list:
+        return False
+    axes = [a for a in (ms_list[0].get("axes") or []) if isinstance(a, dict)]
+    spatial = [a for a in axes if a.get("type") not in ("channel", "time")]
+    return bool(spatial) and all(a.get("unit") for a in spatial)
+
+
+_MANIFEST_PITCH: dict[str, float] = {}
+
+
+def _pitch_from_manifest(prefix: str) -> float | None:
+    """Read properties.pixel_size_um for this store out of the catalog manifest, metadata only.
+
+    The catalog already records the pitch for most of the affected stores, so the repair is a lookup
+    rather than a re-render. Returns None when the manifest has no pitch for this store either.
+    """
+    if not _MANIFEST_PITCH:
+        try:
+            import gzip
+            import json as _json
+            import urllib.request
+
+            url = "https://vesuvius-challenge-open-data.s3.amazonaws.com/metadata.min.json"
+            request = urllib.request.Request(url, headers={"User-Agent": "lasagna-download-omezarr"})
+            with urllib.request.urlopen(request, timeout=60) as response:
+                raw = response.read()
+            if raw[:2] == b"\x1f\x8b":
+                raw = gzip.decompress(raw)
+            doc = _json.loads(raw.decode("utf-8"))
+            for info in (doc.get("samples") or {}).values():
+                for section in ("volumes", "segments"):
+                    for item in (info.get(section) or {}).values():
+                        pitch = (item.get("properties") or {}).get("pixel_size_um")
+                        if not pitch:
+                            continue
+                        for datum in item.get("data") or []:
+                            for origin in datum.get("origins") or []:
+                                path = (origin.get("path") or "").strip("/").rstrip("/")
+                                if path:
+                                    _MANIFEST_PITCH[path] = float(pitch)
+        except Exception:
+            return None
+    return _MANIFEST_PITCH.get(prefix.strip("/").rstrip("/"))
+
+
+def _parse_multiscales(zattrs: dict, missing: set[int] | None = None) -> dict[int, list[float]]:
+    """Map level -> ZYX scale from OME metadata.
+
+    A level whose dataset has no `scale` coordinate transformation is recorded as
+    1.0 per voxel (the OME default for a missing transform). Such levels are also
+    reported through `missing`, so callers can warn rather than silently proceed
+    with a unitless scale.
+    """
     ms_list = zattrs.get("multiscales", [])
     if not ms_list:
         return {}
@@ -461,11 +521,15 @@ def _parse_multiscales(zattrs: dict) -> dict[int, list[float]]:
             continue
         transforms = ds.get("coordinateTransformations", [])
         scale = [1.0, 1.0, 1.0]
+        has_scale = False
         for t in transforms:
             if t.get("type") == "scale":
                 s = t["scale"]
                 scale = [float(v) for v in s[-n_spatial:]]
+                has_scale = True
                 break
+        if not has_scale and missing is not None:
+            missing.add(lvl)
         result[lvl] = scale
     return result
 
@@ -1114,6 +1178,11 @@ def parse_args(argv: list[str] | None = None) -> argparse.Namespace:
                    help="Z slice range at scale 0: z0,z1 (downloads all X/Y for that Z range)")
     p.add_argument("--workers", type=_positive_int, default=64,
                    help="Parallel download threads (default: 64)")
+    p.add_argument(
+        "--require-scale",
+        action="store_true",
+        help="Fail instead of proceeding when OME metadata carries no physical scale",
+    )
     p.add_argument("--order", choices=("z", "morton", "random"), default="z",
                    help="Chunk queue order: z is z-major z/y/x, morton is a Z-order curve (default: z)")
     p.add_argument("--no-remote-inventory", action="store_true",
@@ -1181,7 +1250,8 @@ def main(argv: list[str] | None = None) -> int:
     except (botocore.exceptions.ClientError, FileNotFoundError):
         zattrs = {}
 
-    multiscales = _parse_multiscales(zattrs)
+    missing_scale: set[int] = set()
+    multiscales = _parse_multiscales(zattrs, missing_scale)
     available_levels = _discover_levels(bucket, prefix, anon)
     if not available_levels:
         print("ERROR: no OME-Zarr levels found", file=sys.stderr)
@@ -1228,6 +1298,36 @@ def main(argv: list[str] | None = None) -> int:
     if not multiscales:
         for lvl in levels:
             multiscales[lvl] = [1.0, 1.0, 1.0]
+        missing_scale.update(levels)
+
+    # The defect the catalog audit reported is not a missing transform: the affected stores publish
+    # an explicit [1, 1, 1] scale and no unit on any axis. Check the three conditions separately and,
+    # when the catalog records the pitch, report it - that value is the repair.
+    unitless = not _axes_have_units(zattrs)
+    placeholder = sorted(lvl for lvl, scale in multiscales.items() if scale == [1.0, 1.0, 1.0])
+    if missing_scale or unitless or placeholder:
+        reasons = []
+        if missing_scale:
+            reasons.append(f"level(s) {sorted(missing_scale)} carry no scale transform")
+        if unitless:
+            reasons.append("no axes[].unit is declared")
+        if placeholder:
+            reasons.append(f"level(s) {placeholder} carry the [1, 1, 1] placeholder scale")
+        declared = _pitch_from_manifest(prefix)
+        advice = (
+            f"the catalog manifest records pixel_size_um={declared}; use that value"
+            if declared is not None
+            else "the catalog manifest records no pitch for this store either"
+        )
+        print(
+            f"WARN: {prefix} does not describe its physical scale: " + "; ".join(reasons)
+            + f". {advice}. Offsets and distances derived from 1.0 per voxel will be wrong.",
+            file=sys.stderr,
+        )
+        if args.require_scale:
+            print("ERROR: --require-scale was set and the physical scale is not usable",
+                  file=sys.stderr)
+            return 1
 
     # --- Run ---
     local_chunk_keys: dict[int, set[str]] = {}

lasagna/scripts/zarr_threshold_masks.py

@@ -7,6 +7,7 @@
 import json
 import multiprocessing
 import os
+import re
 import shutil
 import sys
 import time
@@ -131,7 +132,23 @@ def _infer_zyx_scale(arr, array_path: Path) -> tuple[float, float, float]:
         return scale
     scale = _scale_from_ome_parent(array_path)
     if scale is not None:
+        # A [1, 1, 1] parent scale on a store whose name states a pitch is the placeholder case the
+        # catalog audit found, not a measured value - warn instead of returning it silently.
+        stated = re.search(r"(\d+(?:\.\d+)?)um", str(array_path), re.IGNORECASE)
+        if stated and tuple(float(v) for v in scale) == (1.0, 1.0, 1.0):
+            print(
+                f"WARN: parent OME-Zarr metadata for {array_path} declares a [1, 1, 1] scale with no "
+                f"axes[].unit, while the path states {stated.group(1)}um; distances derived from 1.0 "
+                "per voxel will be wrong unless the voxel size really is 1.0.",
+                file=sys.stderr,
+            )
         return scale
+    print(
+        f"WARN: no physical scale found for {array_path} (neither preprocess params nor "
+        "parent OME-Zarr metadata carry one); assuming 1.0 per voxel. Physical distances "
+        "derived from this scale will be wrong unless the voxel size is known to be 1.0.",
+        file=sys.stderr,
+    )
     return 1.0, 1.0, 1.0
 
 

lasagna/tests/test_download_omezarr.py

@@ -368,3 +368,106 @@ def fake_list(_bucket, prefix, _anon):
         "3.0.0",
     }
     assert calls == ["root/0/1.", "root/0/3."]
+
+
+def test_parse_multiscales_reports_levels_without_a_scale_transform():
+    zattrs = {
+        "multiscales": [
+            {
+                "datasets": [
+                    {
+                        "path": "0",
+                        "coordinateTransformations": [
+                            {"type": "scale", "scale": [8.64, 8.64, 8.64]}
+                        ],
+                    },
+                    {"path": "1"},
+                    {
+                        "path": "2",
+                        "coordinateTransformations": [
+                            {"type": "translation", "translation": [0.0, 0.0, 0.0]}
+                        ],
+                    },
+                ]
+            }
+        ]
+    }
+    missing: set[int] = set()
+
+    assert dl._parse_multiscales(zattrs, missing) == {
+        0: [8.64, 8.64, 8.64],
+        1: [1.0, 1.0, 1.0],
+        2: [1.0, 1.0, 1.0],
+    }
+    assert missing == {1, 2}
+
+
+def test_parse_multiscales_without_missing_sink_keeps_previous_behaviour():
+    assert dl._parse_multiscales({}) == {}
+    assert dl._parse_multiscales({"multiscales": [{"datasets": [{"path": "0"}]}]}) == {
+        0: [1.0, 1.0, 1.0]
+    }
+
+
+# ---------------------------------------------------------------------------------------------
+# The defect the catalog audit reported is NOT a missing scale transform: the affected stores
+# publish an explicit [1, 1, 1] scale and no unit on any axis. These tests pin the check to that
+# condition, because the first version of this patch only fired on the missing-transform case and
+# never warned on the stores it was written for (review on ScrollPrize/villa#2137).
+
+
+def test_axes_have_units_is_false_for_the_audited_shape():
+    zattrs = {
+        "multiscales": [
+            {
+                "axes": [{"name": "z"}, {"name": "y"}, {"name": "x"}],
+                "datasets": [
+                    {"path": "0", "coordinateTransformations": [{"type": "scale", "scale": [1.0, 1.0, 1.0]}]}
+                ],
+            }
+        ]
+    }
+    assert download_omezarr._axes_have_units(zattrs) is False
+
+
+def test_axes_have_units_is_true_when_every_axis_declares_a_unit():
+    zattrs = {
+        "multiscales": [
+            {
+                "axes": [
+                    {"name": "z", "unit": "micrometer"},
+                    {"name": "y", "unit": "micrometer"},
+                    {"name": "x", "unit": "micrometer"},
+                ]
+            }
+        ]
+    }
+    assert download_omezarr._axes_have_units(zattrs) is True
+
+
+def test_axes_have_units_ignores_a_channel_axis_without_a_unit():
+    zattrs = {
+        "multiscales": [
+            {
+                "axes": [
+                    {"name": "c", "type": "channel"},
+                    {"name": "z", "unit": "micrometer"},
+                    {"name": "y", "unit": "micrometer"},
+                    {"name": "x", "unit": "micrometer"},
+                ]
+            }
+        ]
+    }
+    assert download_omezarr._axes_have_units(zattrs) is True
+
+
+def test_pitch_from_manifest_returns_the_recorded_pitch(monkeypatch):
+    monkeypatch.setattr(
+        download_omezarr, "_MANIFEST_PITCH", {"PHerc0009B/volumes/x.zarr": 8.64}
+    )
+    assert download_omezarr._pitch_from_manifest("PHerc0009B/volumes/x.zarr/") == 8.64
+
+
+def test_pitch_from_manifest_returns_none_when_the_catalog_has_no_pitch(monkeypatch):
+    monkeypatch.setattr(download_omezarr, "_MANIFEST_PITCH", {"other/store.zarr": 2.4})
+    assert download_omezarr._pitch_from_manifest("PHercParis2Fr47/legacy/x.zarr") is None

lasagna/tests/test_zarr_threshold_masks_scale.py

@@ -0,0 +1,67 @@
+"""A missing physical scale must be visible, not silently replaced by 1.0."""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+import pytest
+
+zarr_threshold_masks = pytest.importorskip("lasagna.scripts.zarr_threshold_masks")
+
+
+def test_infer_zyx_scale_warns_when_no_scale_is_available(monkeypatch, capsys):
+    monkeypatch.setattr(zarr_threshold_masks, "_scale_from_preprocess_params", lambda arr: None)
+    monkeypatch.setattr(zarr_threshold_masks, "_scale_from_ome_parent", lambda path: None)
+
+    scale = zarr_threshold_masks._infer_zyx_scale(None, Path("/nonexistent/level0"))
+
+    assert scale == (1.0, 1.0, 1.0)
+    err = capsys.readouterr().err
+    assert "no physical scale found" in err
+    assert "assuming 1.0 per voxel" in err
+
+
+def test_infer_zyx_scale_prefers_metadata_without_warning(monkeypatch, capsys):
+    monkeypatch.setattr(zarr_threshold_masks, "_scale_from_preprocess_params", lambda arr: None)
+    monkeypatch.setattr(zarr_threshold_masks, "_scale_from_ome_parent", lambda path: (8.64, 8.64, 8.64))
+
+    assert zarr_threshold_masks._infer_zyx_scale(None, Path("/nonexistent/level0")) == (
+        8.64,
+        8.64,
+        8.64,
+    )
+    assert "no physical scale found" not in capsys.readouterr().err
+
+
+def test_a_111_parent_scale_on_a_named_pitch_warns(monkeypatch, capsys):
+    """The audited shape: parent OME metadata returns [1,1,1] while the path states 8.64um."""
+    monkeypatch.setattr(zarr_threshold_masks, "_scale_from_preprocess_params", lambda arr: None)
+    monkeypatch.setattr(
+        zarr_threshold_masks, "_scale_from_ome_parent", lambda path: (1.0, 1.0, 1.0)
+    )
+    scale = zarr_threshold_masks._infer_zyx_scale(
+        None, Path("/data/PHerc0009B/volumes/20250521125136-8.640um-masked.zarr/0")
+    )
+    assert scale == (1.0, 1.0, 1.0)
+    err = capsys.readouterr().err
+    assert "8.640um" in err
+    assert "[1, 1, 1]" in err
+
+
+def test_a_111_parent_scale_without_a_stated_pitch_stays_quiet(monkeypatch, capsys):
+    monkeypatch.setattr(zarr_threshold_masks, "_scale_from_preprocess_params", lambda arr: None)
+    monkeypatch.setattr(
+        zarr_threshold_masks, "_scale_from_ome_parent", lambda path: (1.0, 1.0, 1.0)
+    )
+    zarr_threshold_masks._infer_zyx_scale(None, Path("/data/store/zarr/0"))
+    assert capsys.readouterr().err == ""
+
+
+def test_a_real_scale_is_returned_without_a_warning(monkeypatch, capsys):
+    monkeypatch.setattr(zarr_threshold_masks, "_scale_from_preprocess_params", lambda arr: None)
+    monkeypatch.setattr(
+        zarr_threshold_masks, "_scale_from_ome_parent", lambda path: (2.4, 2.4, 2.4)
+    )
+    scale = zarr_threshold_masks._infer_zyx_scale(None, Path("/data/x-8.64um.zarr/0"))
+    assert scale == (2.4, 2.4, 2.4)
+    assert capsys.readouterr().err == ""

spiral-fitting/tests/test_tracks_to_ome_zarr.py

@@ -182,6 +182,28 @@ def test_pipeline_schedule_does_not_change_output(self):
                 zarr.open_group(pipelined_path, mode='r')['0'][:],
             )
 
+    def test_unitless_scale_is_reported_when_no_voxel_size_is_available(self):
+        """A physical scale that had to be assumed must be visible, not silent."""
+        with tempfile.TemporaryDirectory() as temporary:
+            temporary = Path(temporary)
+            database_path = temporary / 'tracks.dbm'
+            output_path = temporary / 'unitless.ome.zarr'
+            with dbm.open(str(database_path), 'c') as database:
+                database[b'a'] = pickle.dumps([
+                    np.asarray([[1, 1, 1], [1, 1, 8]], dtype=np.int32),
+                ])
+
+            result = CliRunner().invoke(main, [
+                str(database_path), '--shape', '16,16,16', '--chunk', '16',
+                '--out', str(output_path), '--workers', '1',
+            ])
+            if result.exception:
+                raise result.exception
+
+            # The warning names the option to pass; the written value is unchanged.
+            self.assertIn('no --voxel-size', result.output)
+
+
 
 if __name__ == '__main__':
     unittest.main()

spiral-fitting/tracks_to_ome_zarr.py

@@ -611,6 +611,12 @@ def main(track_dbms, output_path, shape_text, reference_path, like_group, voxel_
         scale = _parse_zyx(voxel_size, '--voxel-size', float)
     else:
         scale = reference_scale or (1.0, 1.0, 1.0)
+        if not reference_scale:
+            click.echo(
+                'WARN: no --voxel-size and no --like reference provide a physical scale; '
+                'writing 1.0 per voxel into the output OME-Zarr metadata. Pass '
+                '--voxel-size Z,Y,X if the voxel size is known.'
+            )
     unit = unit or reference_unit
 
     parameters = {