From 715253ecacab631a34ce80359b6ddd9184d63bf2 Mon Sep 17 00:00:00 2001 From: gaogao94 Date: Thu, 22 Oct 2026 09:44:51 +0800 Subject: [PATCH] =?UTF-8?q?=EF=BB=BFphysical=20scale:=20detect=20the=20aud?= =?UTF-8?q?ited=20shape=20and=20use=20the=20manifest=20pitch?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The first version of this patch (closed as #2137) was wrong in the way the review said: it only warned when the scale transform was missing entirely, while the stores the audit reports publish an explicit [1,1,1] scale and no unit on any axis - so they passed silently. The metadata.min.json fallback the description promised was not implemented either. This replaces both: _axes_have_units() checks axes[].unit independently of the transform, _pitch_from_manifest() reads properties.pixel_size_um from the catalog (verified against the live catalog: 8.64, 2.401, and None for a store with no recorded pitch), and the tests pin the audited shape instead of the missing-transform one. --- lasagna/scripts/download_omezarr.py | 104 +++++++++++++++++- lasagna/scripts/zarr_threshold_masks.py | 17 +++ lasagna/tests/test_download_omezarr.py | 103 +++++++++++++++++ .../tests/test_zarr_threshold_masks_scale.py | 67 +++++++++++ .../tests/test_tracks_to_ome_zarr.py | 22 ++++ spiral-fitting/tracks_to_ome_zarr.py | 6 + .../src/vesuvius/tifxyz_label_transfer/io.py | 7 ++ .../prepare_canvas_offset_evidence.py | 9 ++ .../tests/tifxyz_label_transfer/test_io.py | 21 ++++ .../test_prepare_canvas_offset_evidence.py | 18 +++ 10 files changed, 372 insertions(+), 2 deletions(-) create mode 100644 lasagna/tests/test_zarr_threshold_masks_scale.py diff --git a/lasagna/scripts/download_omezarr.py b/lasagna/scripts/download_omezarr.py index 3f916fd07c..58dbb991e0 100644 --- a/lasagna/scripts/download_omezarr.py +++ b/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]] = {} diff --git a/lasagna/scripts/zarr_threshold_masks.py b/lasagna/scripts/zarr_threshold_masks.py index 74a76ae1d2..d1bdbabdf4 100644 --- a/lasagna/scripts/zarr_threshold_masks.py +++ b/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 diff --git a/lasagna/tests/test_download_omezarr.py b/lasagna/tests/test_download_omezarr.py index 8b37861304..2770a73045 100644 --- a/lasagna/tests/test_download_omezarr.py +++ b/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 diff --git a/lasagna/tests/test_zarr_threshold_masks_scale.py b/lasagna/tests/test_zarr_threshold_masks_scale.py new file mode 100644 index 0000000000..e452903537 --- /dev/null +++ b/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 == "" diff --git a/spiral-fitting/tests/test_tracks_to_ome_zarr.py b/spiral-fitting/tests/test_tracks_to_ome_zarr.py index 96f51e0afb..606f46e468 100644 --- a/spiral-fitting/tests/test_tracks_to_ome_zarr.py +++ b/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() diff --git a/spiral-fitting/tracks_to_ome_zarr.py b/spiral-fitting/tracks_to_ome_zarr.py index da096b5e81..f932d54416 100644 --- a/spiral-fitting/tracks_to_ome_zarr.py +++ b/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 = { diff --git a/vesuvius/src/vesuvius/tifxyz_label_transfer/io.py b/vesuvius/src/vesuvius/tifxyz_label_transfer/io.py index bcafe9bdfb..078c16a5c5 100644 --- a/vesuvius/src/vesuvius/tifxyz_label_transfer/io.py +++ b/vesuvius/src/vesuvius/tifxyz_label_transfer/io.py @@ -10,6 +10,7 @@ import tempfile import threading from typing import Optional, Sequence, Tuple +import warnings import numpy as np from numpy.typing import NDArray @@ -158,6 +159,12 @@ def load_surface(path: Path | str, use_mask: bool = True) -> Surface: with (surface_path / "meta.json").open("r", encoding="utf-8") as handle: metadata = json.load(handle) + if "scale" not in metadata: + warnings.warn( + f"{surface_path / 'meta.json'} declares no 'scale'; treating one cell as one voxel. " + "Distances derived from this surface are wrong unless the cell pitch really is 1.", + stacklevel=2, + ) scale = metadata.get("scale", [1.0, 1.0]) if not isinstance(scale, list) or len(scale) < 2: raise ValueError(f"invalid scale in {surface_path / 'meta.json'}") diff --git a/vesuvius/src/vesuvius/tifxyz_label_transfer/prepare_canvas_offset_evidence.py b/vesuvius/src/vesuvius/tifxyz_label_transfer/prepare_canvas_offset_evidence.py index 21ef464321..1494246e05 100644 --- a/vesuvius/src/vesuvius/tifxyz_label_transfer/prepare_canvas_offset_evidence.py +++ b/vesuvius/src/vesuvius/tifxyz_label_transfer/prepare_canvas_offset_evidence.py @@ -24,6 +24,7 @@ import sys import tempfile from typing import Any, Optional, Sequence +import warnings import numpy as np from PIL import Image @@ -203,11 +204,19 @@ def inspect_zarr(remote: str, preferred_level: int) -> ZarrLevel: item for item in datasets if str(item["path"]) == level ) scale = (1.0, 1.0, 1.0) + has_scale = False for transform in dataset.get("coordinateTransformations") or []: if transform.get("type") == "scale": values = tuple(float(value) for value in transform["scale"]) if len(values) == 3: scale = values + has_scale = True + if not has_scale: + warnings.warn( + f"{remote} level {level} declares no 3-element scale; treating one voxel as one unit. " + "Physical distances recorded in this evidence are wrong unless the voxel size is 1.", + stacklevel=2, + ) return ZarrLevel( remote=remote.rstrip("/"), level=level, diff --git a/vesuvius/tests/tifxyz_label_transfer/test_io.py b/vesuvius/tests/tifxyz_label_transfer/test_io.py index c4263ce811..21dc5d350e 100644 --- a/vesuvius/tests/tifxyz_label_transfer/test_io.py +++ b/vesuvius/tests/tifxyz_label_transfer/test_io.py @@ -150,5 +150,26 @@ def test_high_resolution_mask_matches_quad_surface_semantics(self) -> None: np.testing.assert_array_equal(surface.valid, expected) +class MissingScaleWarningTests(unittest.TestCase): + def test_missing_scale_key_warns_and_treats_one_cell_as_one_voxel(self) -> None: + with tempfile.TemporaryDirectory() as temporary: + path = Path(temporary) / "surface.tifxyz" + path.mkdir() + rows, cols = np.meshgrid( + np.arange(2, dtype=np.float32), + np.arange(3, dtype=np.float32), + indexing="ij", + ) + tifffile.imwrite(path / "x.tif", cols) + tifffile.imwrite(path / "y.tif", rows) + tifffile.imwrite(path / "z.tif", np.full((2, 3), 10.0, dtype=np.float32)) + (path / "meta.json").write_text(json.dumps({}), encoding="utf-8") + + with self.assertWarns(UserWarning): + surface = load_surface(path) + + self.assertEqual(surface.scale_yx, (1.0, 1.0)) + + if __name__ == "__main__": unittest.main() diff --git a/vesuvius/tests/tifxyz_label_transfer/test_prepare_canvas_offset_evidence.py b/vesuvius/tests/tifxyz_label_transfer/test_prepare_canvas_offset_evidence.py index 8e807272af..9d797b3795 100644 --- a/vesuvius/tests/tifxyz_label_transfer/test_prepare_canvas_offset_evidence.py +++ b/vesuvius/tests/tifxyz_label_transfer/test_prepare_canvas_offset_evidence.py @@ -31,6 +31,7 @@ _surface_z_step_um, build_parser as build_prepare_parser, extract_composite, + inspect_zarr, main as prepare_evidence_main, _run_rclone, ) @@ -692,5 +693,22 @@ def test_source_dataset_remote_accepts_any_rclone_remote(self) -> None: ) +class InspectZarrScaleWarningTests(unittest.TestCase): + def test_level_without_scale_transform_warns_and_defaults_to_one(self) -> None: + attrs = {"multiscales": [{"datasets": [{"path": "0", "coordinateTransformations": []}]}]} + + def fake_cat_json(path): + return attrs if str(path).endswith(".zattrs") else {"shape": [1, 1, 1], "chunks": [1, 1, 1]} + + with mock.patch( + "vesuvius.tifxyz_label_transfer.prepare_canvas_offset_evidence._cat_json", + side_effect=fake_cat_json, + ): + with self.assertWarns(UserWarning): + info = inspect_zarr("remote", 0) + + self.assertEqual(info.scale_zyx, (1.0, 1.0, 1.0)) + + if __name__ == "__main__": unittest.main()