diff --git a/docs/reference/cli/index.md b/docs/reference/cli/index.md index 04f83a6..5c0301a 100644 --- a/docs/reference/cli/index.md +++ b/docs/reference/cli/index.md @@ -10,7 +10,7 @@ stop. Supply tuple options as separate values, such as | Command | Purpose | | --- | --- | -| [inspect-opm](inspect-opm.md) | Print acquisition metadata. | +| [inspect-opm](inspect-opm.md) | Report raw/processed dimensions, physical sizes and camera timing. | | [process](process.md) | Calibrate, deconvolve, and deskew an acquisition. | | [fuse](fuse.md) | Register and blend processed tiles. | | [export-ome-tiff](export-ome-tiff.md) | Write fused images as compressed OME-BigTIFF. | diff --git a/docs/reference/cli/inspect-opm.md b/docs/reference/cli/inspect-opm.md index 99ab12e..ab24252 100644 --- a/docs/reference/cli/inspect-opm.md +++ b/docs/reference/cli/inspect-opm.md @@ -1,10 +1,37 @@ # inspect-opm -Print acquisition dimensions, channels, calibration, and scan geometry as JSON. -Accepts an OME-Zarr store or its containing directory. +Report available raw and processed OME-Zarr metadata as JSON, without reading +image pixels. Accepts a raw store, a processed store, or an acquisition directory. ```bash uv run inspect-opm "/path/to/acquisition" +uv run inspect-opm "/path/to/acquisition/sample_deskewed.ome.zarr" ``` +The report contains `raw` and `processed` sections. Acquisition directories and +raw paths include associated processing products when present. Processing +journals also identify products with custom names and connect a processed store +to its raw source. A processed store remains inspectable when the raw source is +unavailable; unknown calibration or timing is reported as `null`. + +Raw dimensions use `s` for **scan positions**, `y'` for distance along the tilted +camera plane, and `x` for camera columns. `sampling_um` records scan spacing and +the calibrated camera pixel pitch. `camera_field_of_view_um` gives the physical +camera-plane dimensions. `y_prime_sampling_projected_y_um` is the camera pitch +multiplied by cos(tilt angle). No raw laboratory-Z spacing is reported. + +`timing.frame_exposure_ms` lists actual camera exposures in stored channel order, +rounded to two decimal places in milliseconds. Hardware-triggered volume time is +`sum(frame_exposure_ms) × acquired scan positions`, including positions later +discarded during processing. Software arrival timestamps, requested exposures, +and processed Z depth do not determine volume time. Combined laser channels +stored in one camera image count as one exposure. + +Processed images use laboratory `z`, `y`, and `x`. Each position/series reports +its level-zero dimensions, voxel spacing, origin, channels and pyramid levels. +`volume_dimensions_um` gives **Z (µm) × Y (µm) × X (µm)** using voxel count times +recorded spacing, including any stored padding. These are the dimensions of one +volume, independent of its channel and timepoint counts. Dataset and shared +multiscale transforms are composed before reporting physical dimensions. + Append `--help` to list all options and defaults. diff --git a/src/opm_processing/dataio/acquisition.py b/src/opm_processing/dataio/acquisition.py index fbe6747..c16da6a 100644 --- a/src/opm_processing/dataio/acquisition.py +++ b/src/opm_processing/dataio/acquisition.py @@ -1061,18 +1061,25 @@ def open_acquisition_datastore( return ts.stack(arrays, axis=1) -@app.command(help="Print a metadata-only acquisition manifest as JSON.") +@app.command( + help="Print raw acquisition and processed laboratory-image metadata as JSON." +) def inspect_command( - path: Annotated[Path, typer.Argument(help="Zarr store or containing directory")], + path: Annotated[ + Path, + typer.Argument(help="Raw or processed Zarr store, or acquisition directory"), + ], ) -> None: - """Print a metadata-only acquisition manifest as JSON. + """Print raw acquisition and processed laboratory-image metadata as JSON. Parameters ---------- path : pathlib.Path - Acquisition OME-Zarr store or its containing directory. + Raw or processed OME-Zarr store, or its acquisition directory. """ - typer.echo(json.dumps(inspect_acquisition(path).to_dict(), indent=2)) + from opm_processing.dataio.inspection import inspect_dataset + + typer.echo(json.dumps(inspect_dataset(path), indent=2, allow_nan=False)) def main() -> None: diff --git a/src/opm_processing/dataio/inspection.py b/src/opm_processing/dataio/inspection.py new file mode 100644 index 0000000..148f696 --- /dev/null +++ b/src/opm_processing/dataio/inspection.py @@ -0,0 +1,295 @@ +"""Metadata-only reports for raw oblique acquisitions and laboratory images.""" + +from __future__ import annotations + +import math +from pathlib import Path +from typing import Any + +from ome_types import from_xml +from yaozarrs import open_group, v05 + +from opm_processing.dataio.acquisition import ( + AcquisitionMetadata, + _event_parts, + _initial_frame_metadata, + acquisition_stem, + inspect_acquisition, +) +from opm_processing.dataio.processing_state import ProcessingState + + +def hardware_timing( + frames: list[dict], channel_count: int, scan_positions: int +) -> dict[str, Any]: + """Calculate hardware-triggered timing from actual camera frame exposures. + + Parameters + ---------- + frames : list of dict + Camera frame records in either supported acquisition metadata format. + channel_count : int + Number of stored camera channels, including combined-laser images once. + scan_positions : int + Acquired scan positions, including planes subsequently discarded. + + Returns + ------- + dict + Channel-ordered exposures and volume time in milliseconds. Unknown + exposures remain None; software arrival times and requested camera + settings do not determine hardware timing. + """ + exposures: list[float | None] = [None] * channel_count + for frame in frames: + if frame.get("exposure_time") is not None: + exposure, factor = frame["exposure_time"], 1000 + elif frame.get("exposure_ms") is not None: + exposure, factor = frame["exposure_ms"], 1 + else: + continue + try: + value = round(float(exposure) * factor, 2) + except (TypeError, ValueError) as error: + raise ValueError("Invalid camera frame exposure") from error + if isinstance(exposure, bool) or not math.isfinite(value) or value <= 0: + raise ValueError("Camera frame exposures must be finite and positive") + index, _ = _event_parts(frame) + channel = int(index.get("c", 0)) + if not 0 <= channel < channel_count: + raise ValueError( + "Camera frame channel index is outside the stored channels" + ) + if exposures[channel] is not None and exposures[channel] != value: + raise ValueError("Camera frame exposures vary between scans or positions") + exposures[channel] = value + volume = ( + round(scan_positions * sum(exposures), 6) + if exposures and all(value is not None for value in exposures) + else None + ) + return {"frame_exposure_ms": exposures, "volume_time_ms": volume} + + +def raw_report(acquisition: AcquisitionMetadata) -> dict[str, Any]: + """Describe scan and tilted-camera coordinates without relabeling stored data. + + Parameters + ---------- + acquisition : AcquisitionMetadata + Metadata inspected from an existing raw acquisition. + + Returns + ------- + dict + Raw dimensions in TPCS Y-prime X order, camera calibration and hardware + timing. Y-prime is distance along the camera plane, not laboratory Y. + """ + report = acquisition.to_dict() + names = {"z": "s", "y": "y'"} + axes = [names.get(axis, axis) for axis in acquisition.axes] + report["axes"] = axes + report["index_sizes"] = dict(zip(axes, acquisition.shape, strict=True)) + report["acquisition_order"] = [ + names.get(axis, axis) for axis in acquisition.acquisition_order + ] + report["coordinate_frame"] = "oblique" + report["dimension_names"] = { + "s": "scan", + "y'": "tilted camera plane", + "x": "camera X", + } + pitch = acquisition.pixel_size_um + report["sampling_um"] = { + "s": acquisition.scan_axis_step_um, + "y'": pitch, + "x": pitch, + } + report["camera_field_of_view_um"] = { + names.get(axis, axis): size * pitch if pitch is not None else None + for axis, size in zip(acquisition.axes, acquisition.shape, strict=True) + if axis in ("y", "x") + } + report["y_prime_sampling_projected_y_um"] = ( + pitch * math.cos(math.radians(acquisition.angle_deg)) + if pitch is not None and acquisition.angle_deg is not None + else None + ) + count = len(acquisition.channels) * ( + 2 if acquisition.scan_position_count > 1 else 1 + ) + frames = [ + frame + for array_path in acquisition.array_paths + for frame in _initial_frame_metadata( + acquisition.path / Path(array_path).parent / "zarr.json", count + ) + ] + timing = hardware_timing( + frames, len(acquisition.channels), acquisition.scan_position_count + ) + report["timing"] = timing + for channel, exposure in zip( + report["channels"], timing["frame_exposure_ms"], strict=True + ): + channel["exposure_ms"] = exposure + return report + + +def processed_report(path: Path) -> dict[str, Any]: + """Read each processed image's laboratory sampling and physical dimensions. + + Parameters + ---------- + path : pathlib.Path + Processed root Image or Bio-Formats2Raw position collection. + + Returns + ------- + dict + Per-series dimensions, channels, ZYX spacing, origins and physical + volume dimensions. Dimensions cover N voxels, each with its recorded + spacing; they are not the center-to-center span of N minus one voxels. + """ + root = open_group(path) + layout = root.ome_metadata() + series = ( + root["OME"].ome_metadata().series if isinstance(layout, v05.Bf2Raw) else ("",) + ) + companion = path / "OME" / "METADATA.ome.xml" + ome = ( + from_xml(companion.read_text(encoding="utf-8")) if companion.is_file() else None + ) + images = [] + units = {"meter": 1e6, "millimeter": 1e3, "micrometer": 1, "nanometer": 1e-3} + for series_index, name in enumerate(series): + group = root[str(name)] if name else root + metadata = group.ome_metadata() + if not isinstance(metadata, v05.Image): + raise ValueError(f"Expected an OME-Zarr image in {path}") + multiscale = metadata.multiscales[0] + dataset = multiscale.datasets[0] + shape = group[dataset.path].metadata.shape + axes = [axis.name for axis in multiscale.axes] + scale, origin = [1.0] * len(axes), [0.0] * len(axes) + for transform in ( + *dataset.coordinateTransformations, + *(multiscale.coordinateTransformations or ()), + ): + if transform.type == "scale": + scale = [ + value * factor + for value, factor in zip(scale, transform.scale, strict=True) + ] + origin = [ + value * factor + for value, factor in zip(origin, transform.scale, strict=True) + ] + else: + origin = [ + value + offset + for value, offset in zip(origin, transform.translation, strict=True) + ] + spacing_um, origin_um, dimensions_um = {}, {}, {} + for axis, size, step, offset in zip( + multiscale.axes, shape, scale, origin, strict=True + ): + if axis.name not in "zyx": + continue + factor = units.get(axis.unit) + spacing_um[axis.name] = None if factor is None else step * factor + origin_um[axis.name] = None if factor is None else offset * factor + dimensions_um[axis.name] = None if factor is None else size * step * factor + images.append( + { + "series": str(name), + "axes": axes, + "shape": list(shape), + "index_sizes": dict(zip(axes, shape, strict=True)), + "sampling_um": spacing_um, + "origin_um": origin_um, + "volume_dimensions_um": dimensions_um, + "channels": [channel.label for channel in metadata.omero.channels] + if metadata.omero + else [ + channel.name for channel in ome.images[series_index].pixels.channels + ] + if ome is not None + else [], + "pyramid_levels": len(multiscale.datasets), + } + ) + return {"path": str(path), "coordinate_frame": "laboratory", "images": images} + + +def inspect_dataset(path: str | Path) -> dict[str, Any]: + """Inspect available raw and processed stores without opening pixel arrays. + + Parameters + ---------- + path : str or pathlib.Path + Raw store, processed store or directory containing acquisition products. + + Returns + ------- + dict + Raw acquisition report when available and associated processed images. + Processing journals associate custom output names and source paths. + """ + path = Path(path).expanduser().resolve() + direct = (path / "zarr.json").is_file() or (path / ".zattrs").is_file() + directory = path.parent if direct else path + stores = sorted(item for item in directory.glob("*.zarr") if item.is_dir()) + if direct and path not in stores: + stores.append(path) + roots = {store: open_group(store) for store in stores} + raw_paths = [store for store, root in roots.items() if "opm_v2" in root.attrs] + if direct and path in raw_paths: + raw_paths = [path] + elif direct: + raw_paths = [ + store + for store in raw_paths + if path.name.startswith(acquisition_stem(store) + "_") + ] + if len(raw_paths) > 1: + raise ValueError("Multiple raw acquisitions found; specify one raw store") + associated = set() + for sidecar in directory.glob("*.processing.json"): + state = ProcessingState.read(sidecar) + source = Path(state.document["source"]["path"]) + if not source.is_absolute(): + source = sidecar.parent / source + source = source.resolve() + products = {sidecar.parent / name for name in state.document["outputs"]} + products.update( + sidecar.parent / record[key] + for record in state.document["registration"].values() + for key in ("fused_path", "max_projection_path") + if record.get(key) + ) + if source in raw_paths or (direct and path in products): + associated.update(products) + if source.is_dir() and not raw_paths: + raw_paths = [source] + raw = raw_report(inspect_acquisition(raw_paths[0])) if raw_paths else None + if raw is not None: + prefix = acquisition_stem(raw["path"]) + "_" + associated.update( + store + for store in stores + if store.name.startswith(prefix) and "opm_v2" not in roots[store].attrs + ) + elif not direct: + associated.update(store for store in stores if store not in raw_paths) + if direct and path not in raw_paths: + associated = {path} + processed = [] + for store in sorted(associated): + if store.is_dir(): + product = processed_report(store) + product["timing"] = raw["timing"] if raw is not None else None + processed.append(product) + if raw is None and not processed: + raise ValueError(f"No raw or processed OME-Zarr data found in {path}") + return {"raw": raw, "processed": processed} diff --git a/tests/cpu/test_cli_workflows.py b/tests/cpu/test_cli_workflows.py index 27b9add..660ef3a 100644 --- a/tests/cpu/test_cli_workflows.py +++ b/tests/cpu/test_cli_workflows.py @@ -29,11 +29,12 @@ def test_inspection_and_exports_process_simulated_acquisition( options = cli_export_options metadata = case.dataset.metadata source = case.dataset.path.parent - manifest = json.loads( + report = json.loads( invoke_cli( "inspect-opm", source if options["depth_color"] else case.dataset.path ).output ) + manifest = report["raw"] assert manifest["shape"] == list(case.dataset.raw_data.shape) assert manifest["scan_axis_step_um"] == metadata.scan_axis_step_um assert manifest["camera_offset"] == metadata.camera_offset diff --git a/tests/cpu/test_inspection.py b/tests/cpu/test_inspection.py new file mode 100644 index 0000000..d024cc3 --- /dev/null +++ b/tests/cpu/test_inspection.py @@ -0,0 +1,205 @@ +"""Analytical metadata units and CLI inspection of simulated OPM acquisitions.""" + +import json + +import numpy as np +import pytest +import zarr +from yaozarrs import open_group + +from opm_processing.dataio.acquisition import inspect_acquisition +from opm_processing.dataio.inspection import ( + hardware_timing, + inspect_dataset, + processed_report, + raw_report, +) + + +@pytest.mark.unit +@pytest.mark.parametrize("field,factor", [("exposure_ms", 1), ("exposure_time", 1000)]) +@pytest.mark.parametrize("planes", [1, 17]) +def test_hardware_timing_uses_camera_exposures(field, factor, planes): + """Sum unequal camera exposures once per stored channel and scan position.""" + frames = [ + { + "event_index": {"c": channel}, + field: exposure / factor, + "runner_time_ms": 123456, + "event_metadata": {"Camera": {"exposure_ms": 999}}, + } + for channel, exposure in enumerate((5.002102941, 15.0)) + for _ in range(2) + ] + assert hardware_timing(frames, 2, planes) == { + "frame_exposure_ms": [5, 15], + "volume_time_ms": 20 * planes, + } + + +@pytest.mark.unit +@pytest.mark.parametrize("exposure", [None, -1, 0, float("nan"), True, "invalid"]) +def test_missing_or_invalid_hardware_exposure_is_not_inferred(exposure): + """Keep unknown camera timing unknown and reject nonphysical exposures.""" + frames = [ + {"mda_event": {"index": {"c": 0}, "metadata": {"Camera": {"exposure_ms": 10}}}} + ] + if exposure is None: + assert hardware_timing(frames, 1, 17) == { + "frame_exposure_ms": [None], + "volume_time_ms": None, + } + else: + frames[0]["exposure_ms"] = exposure + with pytest.raises(ValueError, match="exposure"): + hardware_timing(frames, 1, 17) + + +@pytest.mark.unit +@pytest.mark.parametrize( + "unit,factor", [("micrometer", 1), ("nanometer", 1000), (None, None)] +) +@pytest.mark.parametrize("global_transform", [False, True]) +def test_processed_volume_dimensions_use_recorded_lab_grid( + fused_store_factory, unit, factor, global_transform +): + """Report known voxel extents and composed transforms without reading pixels.""" + path = fused_store_factory( + np.arange(60, dtype=np.uint16).reshape(1, 1, 3, 4, 5), + spacing=(0.23, 0.115, 0.115), + origin=(1, 2, 3), + factors=(1, 2), + ) + ome = dict(open_group(path).attrs["ome"]) + multiscale = ome["multiscales"][0] + for axis in multiscale["axes"][2:]: + if unit is None: + axis.pop("unit", None) + else: + axis["unit"] = unit + for dataset in multiscale["datasets"]: + for transform in dataset["coordinateTransformations"]: + key = transform["type"] + transform[key][2:] = [value * (factor or 1) for value in transform[key][2:]] + if global_transform: + multiscale["coordinateTransformations"] = [ + {"type": "scale", "scale": [1, 1, 2, 2, 2]}, + {"type": "translation", "translation": [0, 0, *([5 * (factor or 1)] * 3)]}, + ] + zarr.open_group(path, mode="r+").attrs.update({"ome": ome}) + image = processed_report(path)["images"][0] + gain = 2 if global_transform else 1 + for axis, size, spacing, origin in zip( + "zyx", (3, 4, 5), (0.23, 0.115, 0.115), (1, 2, 3), strict=True + ): + if unit is None: + assert image["sampling_um"][axis] is None + assert image["volume_dimensions_um"][axis] is None + else: + assert image["sampling_um"][axis] == pytest.approx(gain * spacing) + assert image["volume_dimensions_um"][axis] == pytest.approx( + size * gain * spacing + ) + assert image["origin_um"][axis] == pytest.approx( + gain * origin + (5 if global_transform else 0) + ) + assert image["pyramid_levels"] == 2 + + +@pytest.mark.integration +@pytest.mark.parametrize( + "selection", + ["raw-only", "directory", "raw", "processed", "custom-processed", "processed-only"], +) +def test_inspect_cli_reports_simulated_raw_and_processed_data( + inspection_case, + selection, + invoke_cli, +): + """Inspect real processing products of an optically simulated oblique object.""" + case, collection = inspection_case.case, inspection_case.collection + if collection is not None: + assert any(np.any(array.read().result() > 0) for array in collection.arrays) + report = json.loads(invoke_cli("inspect-opm", inspection_case.selected).output) + if selection == "processed-only": + assert report["raw"] is None + else: + raw = report["raw"] + assert raw["axes"] == ["t", "p", "c", "s", "y'", "x"] + assert raw["index_sizes"]["s"] == case.raw_data.shape[3] + assert "z" not in raw["index_sizes"] + assert raw["sampling_um"] == { + "s": case.metadata.scan_axis_step_um, + "y'": case.metadata.pixel_size_um, + "x": case.metadata.pixel_size_um, + } + assert raw["y_prime_sampling_projected_y_um"] == pytest.approx( + case.metadata.pixel_size_um * np.sqrt(3) / 2 + ) + assert raw["timing"] == { + "frame_exposure_ms": [10], + "volume_time_ms": 10 * case.raw_data.shape[3], + } + assert len(report["processed"]) == int(selection != "raw-only") + if collection is not None: + assert report["processed"][0]["timing"] == ( + None if selection == "processed-only" else report["raw"]["timing"] + ) + assert len(report["processed"][0]["images"]) == case.raw_data.shape[1] + for image, array in zip( + report["processed"][0]["images"], collection.arrays, strict=True + ): + assert image["axes"] == list("tczyx") + assert image["shape"] == list(array.shape) + assert image["channels"] == ["488nm"] + for axis, size in zip("zyx", array.shape[-3:], strict=True): + assert image["volume_dimensions_um"][axis] == pytest.approx( + size * case.metadata.pixel_size_um + ) + + +@pytest.mark.unit +def test_missing_dataset_is_reported(tmp_path): + """Report absent input instead of inventing acquisition metadata.""" + with pytest.raises(ValueError, match="No raw or processed"): + inspect_dataset(tmp_path) + + +@pytest.mark.unit +def test_raw_report_prefers_actual_camera_exposures(current_opm_v2_stage_scan): + """Ignore stale requested exposure and file-writing timestamps in raw timing.""" + path = current_opm_v2_stage_scan + for position in ("0", "1"): + group = zarr.open_group(path / position, mode="r+") + writer = dict(group.attrs["ome_writers"]) + for frame in writer["frame_metadata"]: + frame["event_metadata"]["Camera"]["exposure_ms"] = 999 + frame["delta_t"] = 123456 + group.attrs.update({"ome_writers": writer}) + report = raw_report(inspect_acquisition(path)) + assert [channel["exposure_ms"] for channel in report["channels"]] == [10, 15] + assert report["timing"] == {"frame_exposure_ms": [10, 15], "volume_time_ms": 75} + assert report["camera_field_of_view_um"] == {"y'": 4 * 0.115, "x": 5 * 0.115} + + +@pytest.mark.unit +@pytest.mark.parametrize( + "problem", ["missing-channel", "varying-exposure", "out-of-range-channel"] +) +def test_hardware_timing_requires_consistent_channel_measurements(problem): + """Reject contradictory records and leave incomplete channel timing unknown.""" + frames = [{"event_index": {"c": 0}, "exposure_ms": 5}] + if problem == "missing-channel": + assert hardware_timing(frames, 2, 17) == { + "frame_exposure_ms": [5, None], + "volume_time_ms": None, + } + else: + frames.append( + { + "event_index": {"c": 2 if problem == "out-of-range-channel" else 0}, + "exposure_ms": 10, + } + ) + with pytest.raises(ValueError, match=r"channel|exposures vary"): + hardware_timing(frames, 2, 17) diff --git a/tests/fixtures/acquisition.py b/tests/fixtures/acquisition.py index dfe3812..8001424 100644 --- a/tests/fixtures/acquisition.py +++ b/tests/fixtures/acquisition.py @@ -604,6 +604,74 @@ def processing_options(): } +@pytest.fixture(params=("mirror", "stage")) +def inspection_case( + request, + simulated_camera, + acquisition_factory, + processing_options, + selection, + small_cpu_kernels, +): + """Prepare CLI inputs from a simulated object and real deskew processing. + + Parameters + ---------- + request : pytest.FixtureRequest + Mirror or stage scan mode. + simulated_camera : callable + Shared physical forward model for a known fluorescent object. + acquisition_factory : callable + Native writer for calibrated opm-v2 acquisition metadata and pixels. + processing_options : dict + Shared settings preserving the laboratory output grid. + selection : str + Raw, directory, processed, custom-name or unavailable-source scenario. + small_cpu_kernels : None + Production CPU kernels configured for these small simulated volumes. + + Returns + ------- + types.SimpleNamespace + Simulated acquisition, actual processed collection and CLI input path. + """ + from opm_processing.dataio.acquisition import acquisition_stem + from opm_processing.dataio.processing_state import ProcessingState + from opm_processing.process import process + + observed = simulated_camera((16, 32, 32)) + pixels = np.broadcast_to(observed.raw, (1, 2, 1, *observed.raw.shape)).copy() + if request.param == "stage": + pixels = pixels[:, :, :, ::-1].copy() + case = acquisition_factory( + pixels, + name="inspection", + mode=request.param, + acquired_plane_count=pixels.shape[3], + plane_period_s=0.01, + ) + path = case.path + output = path.parent / f"{acquisition_stem(path)}_deskewed.ome.zarr" + collection = None + if selection != "raw-only": + process(path, **processing_options) + if selection == "custom-processed": + output = output.rename(path.parent / "custom-output.ome.zarr") + state = ProcessingState.create(path.parent / "custom.processing.json", path) + state.initialize_run(output, configuration={}, overwrite=True) + collection = open_position_collection(output) + if selection == "processed-only": + path.rename(path.parent / "offline-raw") + selected = ( + path.parent + if selection == "directory" + else output + if "processed" in selection + else path + ) + return SimpleNamespace(case=case, collection=collection, selected=selected) + + @pytest.fixture def calibrated_mirror_timelapse( current_single_position_mirror_timelapse, simulated_camera