From 42edb1350ba0cc170953a79be18c74fe379002bf Mon Sep 17 00:00:00 2001 From: Valentin Boussot Date: Sun, 28 Jun 2026 12:25:09 +0200 Subject: [PATCH] feat: add OME-Zarr (ngff-zarr) and DICOM dataset I/O backends --- docs/source/concepts/imaging-formats.md | 210 +++++++++ konfai/utils/ITK.py | 103 +++-- konfai/utils/__init__.py | 6 +- konfai/utils/dataset.py | 387 +++++++++++++++-- konfai/utils/dicom.py | 554 ++++++++++++++++++++++++ konfai/utils/ome_zarr.py | 219 ++++++++++ konfai/utils/utils.py | 21 +- tests/unit/test_imaging_formats.py | 316 ++++++++++++++ tests/unit/test_imaging_roundtrip.py | 208 +++++++++ 9 files changed, 1937 insertions(+), 87 deletions(-) create mode 100644 docs/source/concepts/imaging-formats.md create mode 100644 konfai/utils/dicom.py create mode 100644 konfai/utils/ome_zarr.py create mode 100644 tests/unit/test_imaging_formats.py create mode 100644 tests/unit/test_imaging_roundtrip.py diff --git a/docs/source/concepts/imaging-formats.md b/docs/source/concepts/imaging-formats.md new file mode 100644 index 00000000..f9f4134e --- /dev/null +++ b/docs/source/concepts/imaging-formats.md @@ -0,0 +1,210 @@ +--- +type: reference +title: Imaging Format Readers +created: 2026-06-27 +tags: + - dicom + - ome-zarr + - imaging +related: + - '[[Datasets]]' +--- + +# Imaging Format Readers + +Beyond the SimpleITK formats handled by `konfai/utils/ITK.py`, KonfAI ships two +optional readers for formats common in clinical and bioimaging workflows: DICOM +series and OME-Zarr (OME-NGFF) stores. Both live in `konfai/utils/`, are wired +into `konfai.utils.dataset.Dataset`, and support channel-first, geometry-aware +reads and writes. + +Each reader is an **optional dependency** — install only what a given dataset +needs. + +## DICOM + +Install the reader: + +```bash +pip install "konfai[dicom]" # pulls in pydicom +``` + +A DICOM acquisition is not a folder of independent images: a CT or MRI *series* +is a set of `.dcm` files that together define one 3-D volume. The reader handles +series discovery, slice ordering, geometry extraction, and CT intensity rescale. + +### `read_dicom_series` + +```python +from konfai.utils.dicom import read_dicom_series + +volume, origin, spacing, direction = read_dicom_series("path/to/series") +``` + +`konfai.utils.dicom.read_dicom_series(directory, *, series_uid=None, apply_rescale=True)` +returns a four-tuple: + +| Returned | Shape | Meaning | +| --- | --- | --- | +| `volume` | `(1, Z, Y, X)` `float32` | channel-first voxel data | +| `origin` | `(3,)` | physical position of the first voxel, mm | +| `spacing` | `(3,)` | voxel size `(x, y, z)`, mm | +| `direction` | `(9,)` | row-major 3×3 direction-cosine matrix, flattened | + +The `origin` / `spacing` / `direction` triple maps directly onto an `Attribute` +(`Origin`, `Spacing`, `Direction`; see {doc}`datasets`), so a DICOM series travels +through the pipeline under the same geometry contract as any other format. + +Key behaviors: + +- **Slice ordering** uses `ImagePositionPatient` projected onto the slice normal + derived from `ImageOrientationPatient`, not filename or `InstanceNumber`. +- **`apply_rescale=True`** (the default) applies `RescaleSlope` / + `RescaleIntercept` to convert stored values to Hounsfield Units for CT. Set it + to `False` to keep raw integers (for example for label maps). +- Missing geometry tags, inconsistent slice shapes, and unreadable pixel data all + raise `DatasetManagerError` with an actionable message. + +### Multi-series folders: `series_uid` + +A single folder can hold more than one series (for example a T1 and a T2 from the +same session). `read_dicom_series` resolves this as follows: + +- one series present → it is used automatically; +- multiple series with `series_uid=None` → `DatasetManagerError` listing the + available `SeriesInstanceUID`s; +- pass `series_uid="1.2.840…"` to select one explicitly. + +Use `konfai.utils.dicom.discover_series(directory)` to list the available UIDs +first: + +```python +from konfai.utils.dicom import discover_series, read_dicom_series + +series = discover_series("path/to/study") # {uid: [Path, ...]} +uid = next(iter(series)) +volume, origin, spacing, direction = read_dicom_series("path/to/study", series_uid=uid) +``` + +`write_dicom_series` writes one uncompressed scalar DICOM series. Integer data +round-trips exactly; floating-point data is stored as signed 16-bit pixels with +`RescaleSlope` and `RescaleIntercept`. + +## OME-Zarr + +Install the reader: + +```bash +pip install "konfai[omezarr]" # pulls in zarr + ngff-zarr +``` + +Unlike a NIfTI file, an OME-Zarr array is **already lazy**: it is stored as +chunked Zarr, so reading a sub-region only fetches the chunks it touches. This +maps naturally onto KonfAI's patch-based loading — the reader never materializes +the whole volume. + +### Multiscale levels + +OME-NGFF stores a **multiscale pyramid**: the same image at several resolutions, +level `0` being full resolution and each higher level a downsampled copy. Each +level carries its own physical `scale` (spacing) and `translation` (origin) in the +`.zattrs` metadata, so geometry stays correct at every level. + +### `read_ome_zarr_slice` + +```python +from konfai.utils.ome_zarr import read_ome_zarr_slice + +patch, axes, scale, translation = read_ome_zarr_slice( + "image.zarr", + (slice(0, 64), slice(0, 256), slice(0, 256)), # (Z, Y, X) + level=0, +) +``` + +`konfai.utils.ome_zarr.read_ome_zarr_slice(store_path, slices, *, level=0, channel=None, timepoint=0)` +reads one spatial patch and returns: + +| Returned | Meaning | +| --- | --- | +| `patch` | channel-first `(C, Z, Y, X)` patch preserving the stored dtype | +| `axes` | axis names from the OME-NGFF metadata (e.g. `['z', 'y', 'x']`) | +| `scale` | voxel spacing for the selected level | +| `translation` | origin translation for the selected level | + +The reader inspects the stored `axes` to place the spatial slices on the right +dimensions and to index optional `T` (time) and `C` (channel) axes, so the same +call works for `ZYX`, `CZYX`, and `TCZYX` arrays. Local `.zarr` directories and +remote stores (`s3://`, `gs://`) are both supported. + +### Choosing a resolution: `select_level` + +`konfai.utils.ome_zarr.select_level(zattrs, target_spacing_mm=None)` picks the +pyramid level whose voxel spacing is closest to a target: + +```python +import zarr +from konfai.utils.ome_zarr import select_level, read_ome_zarr_slice + +zattrs = dict(zarr.open("image.zarr", mode="r").attrs) +level = select_level(zattrs, target_spacing_mm=1.0) # nearest level to 1 mm +patch, *_ = read_ome_zarr_slice("image.zarr", (slice(0, 64),) * 3, level=level) +``` + +With `target_spacing_mm=None` it returns level `0` (full resolution). Otherwise it +compares the median spatial scale of each level against the target and returns the +closest — letting you trade resolution for field of view at a fixed patch size. + +Use `konfai.utils.ome_zarr.get_ome_zarr_info(store_path, level=0)` for a metadata +summary (`axes`, `shape`, `chunks`, `dtype`, `scale`, `translation`, `n_levels`) +without reading any pixels. + +`write_ome_zarr` writes a single-level OME-NGFF store with channel/spatial axes, +chunking, scale, translation, and the original KonfAI attributes. + +## Use as a KonfAI dataset + +Both formats use the normal grouped `Dataset` API: + +```python +from konfai.utils.dataset import Dataset + +dicom_dataset = Dataset("DatasetDicom", "dicom") +ome_dataset = Dataset("DatasetOme", "omezarr") # aliases: ome-zarr, ome_zarr, zarr + +dicom_dataset.write("CT", "CASE_001", volume, attributes) +patch, attributes = dicom_dataset.read_data_slice( + "CT", "CASE_001", (slice(None), slice(10, 20), slice(32, 96), slice(32, 96)) +) + +ome_dataset.write("CT", "CASE_001", volume, attributes) +names = ome_dataset.get_names("CT") +``` + +The layouts are `///*.dcm` and +`//.ome.zarr`. `get_infos` reads only metadata, OME patch +reads touch only selected chunks, and DICOM patch reads decode only selected +slices. In workflow YAML use `./Dataset:dicom` or `./Dataset:omezarr` in +`dataset_filenames`. + +## Why `ngff-zarr` over `ome-zarr` + +Two Python libraries read OME-NGFF: the OME consortium's `ome-zarr` and +`ngff-zarr`. KonfAI's `omezarr` extra depends on **`ngff-zarr`** (alongside raw +`zarr`) for one decisive reason: + +> `ngff-zarr` exposes **per-scale physical coordinates** — the `scale` and +> `translation` of each pyramid level — as plain numeric arrays that convert +> directly to SimpleITK geometry. That is exactly the `(Origin, Spacing, +> Direction)` triple KonfAI stores in an `Attribute`, so OME-Zarr input lines up +> with every other format without a bespoke geometry adapter. + +It also handles multiscale selection transparently and adds helpers tuned for 3-D +medical/bioimage workflows. When `ngff-zarr` is unavailable, the reader falls back +to parsing the `.zattrs` JSON with `zarr` alone. + +## See also + +- {doc}`datasets` — the dataset and `Attribute` model these readers feed +- {doc}`configuration` +- {doc}`../getting-started/installation` — the optional-extras table diff --git a/konfai/utils/ITK.py b/konfai/utils/ITK.py index a0150b4a..84329bf8 100644 --- a/konfai/utils/ITK.py +++ b/konfai/utils/ITK.py @@ -16,15 +16,52 @@ """SimpleITK-based helpers for geometric transforms, resampling, and masking.""" +from __future__ import annotations + import numpy as np -import SimpleITK as sitk # noqa: N813 import torch -import torch.nn.functional as F # noqa: N812 + +try: + import SimpleITK as sitk +except ImportError: + sitk = None # type: ignore[assignment] +import torch.nn.functional as F + +from konfai.utils.errors import TransformError + + +def _require_simpleitk() -> None: + """Raise a clear project error when an ITK-only path is used without SimpleITK.""" + if sitk is None: + raise TransformError("SimpleITK is required for this operation. Install it with `pip install konfai[itk]`.") + + +def _invert_via_displacement_field(transform: sitk.Transform, image: sitk.Image) -> sitk.DisplacementFieldTransform: + displacement_field_filter = sitk.TransformToDisplacementFieldFilter() + displacement_field_filter.SetReferenceImage(image) + displacement_field = displacement_field_filter.Execute(transform) + iterative_inverse = sitk.IterativeInverseDisplacementFieldImageFilter() + iterative_inverse.SetNumberOfIterations(20) + return sitk.DisplacementFieldTransform(iterative_inverse.Execute(displacement_field)) + + +def _copy_transform(transform_cls: type[sitk.Transform], transform: sitk.Transform, invert: bool) -> sitk.Transform: + transform = transform_cls(transform) + if invert: + transform = transform_cls(transform.GetInverse()) + return transform + + +def _image_like(array: np.ndarray, reference: sitk.Image) -> sitk.Image: + result = sitk.GetImageFromArray(array) + result.CopyInformation(reference) + return result def _open_transform( transform_files: dict[str | sitk.Transform, bool], image: sitk.Image = None ) -> list[sitk.Transform]: + _require_simpleitk() transforms: list[sitk.Transform] = [] for transform_file, invert in transform_files.items(): @@ -33,45 +70,20 @@ def _open_transform( else: transform = transform_file if transform.GetName() == "TranslationTransform": - transform = sitk.TranslationTransform(transform) - if invert: - transform = sitk.TranslationTransform(transform.GetInverse()) + transform = _copy_transform(sitk.TranslationTransform, transform, invert) elif transform.GetName() == "Euler3DTransform": - transform = sitk.Euler3DTransform(transform) - if invert: - transform = sitk.Euler3DTransform(transform.GetInverse()) + transform = _copy_transform(sitk.Euler3DTransform, transform, invert) elif transform.GetName() == "VersorRigid3DTransform": - transform = sitk.VersorRigid3DTransform(transform) - if invert: - transform = sitk.VersorRigid3DTransform(transform.GetInverse()) + transform = _copy_transform(sitk.VersorRigid3DTransform, transform, invert) elif transform.GetName() == "AffineTransform": - transform = sitk.AffineTransform(transform) - if invert: - transform = sitk.AffineTransform(transform.GetInverse()) + transform = _copy_transform(sitk.AffineTransform, transform, invert) elif transform.GetName() == "DisplacementFieldTransform": if invert: - transform_to_displacement_field_filter = sitk.TransformToDisplacementFieldFilter() - transform_to_displacement_field_filter.SetReferenceImage(image) - displacement_field = transform_to_displacement_field_filter.Execute(transform) - iterative_inverse_displacement_field_image_filter = sitk.IterativeInverseDisplacementFieldImageFilter() - iterative_inverse_displacement_field_image_filter.SetNumberOfIterations(20) - inverse_displacement_field = iterative_inverse_displacement_field_image_filter.Execute( - displacement_field - ) - transform = sitk.DisplacementFieldTransform(inverse_displacement_field) - transforms.append(transform) + transform = _invert_via_displacement_field(transform, image) else: transform = sitk.BSplineTransform(transform) if invert: - transform_to_displacement_field_filter = sitk.TransformToDisplacementFieldFilter() - transform_to_displacement_field_filter.SetReferenceImage(image) - displacement_field = transform_to_displacement_field_filter.Execute(transform) - iterative_inverse_displacement_field_image_filter = sitk.IterativeInverseDisplacementFieldImageFilter() - iterative_inverse_displacement_field_image_filter.SetNumberOfIterations(20) - inverse_displacement_field = iterative_inverse_displacement_field_image_filter.Execute( - displacement_field - ) - transform = sitk.DisplacementFieldTransform(inverse_displacement_field) + transform = _invert_via_displacement_field(transform, image) transforms.append(transform) if len(transforms) == 0: transforms.append(sitk.Euler3DTransform()) @@ -144,6 +156,7 @@ def resample_itk( default_pixel_value: float | None = None, torch_resample: bool = False, ) -> sitk.Image: + _require_simpleitk() if torch_resample: input_tensor = torch.tensor(sitk.GetArrayFromImage(image)).unsqueeze(0) vectors = [torch.arange(0, s) for s in input_tensor.shape[1:]] @@ -191,6 +204,7 @@ def resample_itk( def parametermap_to_transform( path_src: str, ) -> sitk.Transform | list[sitk.Transform]: + _require_simpleitk() transform = sitk.ReadParameterFile(path_src) def array_format(x): @@ -204,7 +218,7 @@ def array_format(x): else: result = sitk.Euler3DTransform() parameters = array_format(transform["TransformParameters"]) - fixed_parameters = array_format(transform["CenterOfRotationPoint"]) + [0] + fixed_parameters = [*array_format(transform["CenterOfRotationPoint"]), 0] elif transform["Transform"][0] == "TranslationTransform": result = sitk.TranslationTransform(dimension) parameters = array_format(transform["TransformParameters"]) @@ -212,7 +226,7 @@ def array_format(x): elif transform["Transform"][0] == "AffineTransform": result = sitk.AffineTransform(dimension) parameters = array_format(transform["TransformParameters"]) - fixed_parameters = array_format(transform["CenterOfRotationPoint"]) + [0] + fixed_parameters = [*array_format(transform["CenterOfRotationPoint"]), 0] elif transform["Transform"][0] == "BSplineStackTransform": parameters = array_format(transform["TransformParameters"]) grid_size = array_format(transform["GridSize"]) @@ -235,7 +249,7 @@ def array_format(x): return results elif transform["Transform"][0] == "AffineLogStackTransform": parameters = array_format(transform["TransformParameters"]) - fixed_parameters = array_format(transform["CenterOfRotationPoint"]) + [0] + fixed_parameters = [*array_format(transform["CenterOfRotationPoint"]), 0] nb = int(transform["NumberOfSubTransforms"][0]) sub = dimension * 4 @@ -289,12 +303,13 @@ def _resample(data: torch.Tensor, size: list[int]) -> torch.Tensor: def resample_isotropic(image: sitk.Image, spacing: list[float] | None = None) -> sitk.Image: + _require_simpleitk() spacing = spacing or [1.0, 1.0, 1.0] - resize_factor = [y / x for x, y in zip(spacing, image.GetSpacing())] + resize_factor = [y / x for x, y in zip(spacing, image.GetSpacing(), strict=False)] result = sitk.GetImageFromArray( _resample( torch.tensor(sitk.GetArrayFromImage(image)).unsqueeze(0), - [int(size * factor) for size, factor in zip(image.GetSize(), resize_factor)], + [int(size * factor) for size, factor in zip(image.GetSize(), resize_factor, strict=False)], ) .squeeze(0) .numpy() @@ -306,30 +321,33 @@ def resample_isotropic(image: sitk.Image, spacing: list[float] | None = None) -> def resample_resize(image: sitk.Image, size: list[int] | None = None): + _require_simpleitk() size = size or [100, 512, 512] result = sitk.GetImageFromArray( _resample(torch.tensor(sitk.GetArrayFromImage(image)).unsqueeze(0), size).squeeze(0).numpy() ) result.SetDirection(image.GetDirection()) result.SetOrigin(image.GetOrigin()) - result.SetSpacing([x / y * z for x, y, z in zip(image.GetSize(), size, image.GetSpacing())]) + result.SetSpacing([x / y * z for x, y, z in zip(image.GetSize(), size, image.GetSpacing(), strict=False)]) return result def box_with_mask(mask: sitk.Image, label: list[int], dilatations: list[int]) -> np.ndarray: + _require_simpleitk() - dilatations = [int(np.ceil(d / s)) for d, s in zip(dilatations, reversed(mask.GetSpacing()))] + dilatations = [int(np.ceil(d / s)) for d, s in zip(dilatations, reversed(mask.GetSpacing()), strict=False)] data = sitk.GetArrayFromImage(mask) border = np.where(np.isin(sitk.GetArrayFromImage(mask), label)) box = [] - for w, dilatation, s in zip(border, dilatations, data.shape): + for w, dilatation, s in zip(border, dilatations, data.shape, strict=False): box.append([max(np.min(w) - dilatation, 0), min(np.max(w) + dilatation, s)]) box = np.asarray(box) return box def crop_with_mask(image: sitk.Image, box: np.ndarray) -> sitk.Image: + _require_simpleitk() data = sitk.GetArrayFromImage(image) for i, w in enumerate(box): @@ -351,6 +369,7 @@ def crop_with_mask(image: sitk.Image, box: np.ndarray) -> sitk.Image: def format_mask_label(mask: sitk.Image, labels: list[tuple[int, int]]) -> sitk.Image: + _require_simpleitk() data = sitk.GetArrayFromImage(mask) result_data = np.zeros_like(data, np.uint8) @@ -363,6 +382,7 @@ def format_mask_label(mask: sitk.Image, labels: list[tuple[int, int]]) -> sitk.I def get_flat_label(mask: sitk.Image, labels: None | list[int] = None) -> sitk.Image: + _require_simpleitk() data = sitk.GetArrayFromImage(mask) result_data = np.zeros_like(data, np.uint8) if labels is not None: @@ -376,6 +396,7 @@ def get_flat_label(mask: sitk.Image, labels: None | list[int] = None) -> sitk.Im def clip_and_cast(image: sitk.Image, min_value: float, max_value: float, dtype: np.dtype) -> sitk.Image: + _require_simpleitk() data = sitk.GetArrayFromImage(image) data[np.where(data > max_value)] = max_value data[np.where(data < min_value)] = min_value diff --git a/konfai/utils/__init__.py b/konfai/utils/__init__.py index 4edc7b6d..2d7c8e84 100644 --- a/konfai/utils/__init__.py +++ b/konfai/utils/__init__.py @@ -1 +1,5 @@ -"""Utility modules supporting configuration, datasets, ITK, and runtime helpers.""" +"""Utility modules supporting configuration, datasets, ITK, DICOM, OME-Zarr, and runtime helpers.""" + +from konfai.utils import dicom, ome_zarr + +__all__ = ["dicom", "ome_zarr"] diff --git a/konfai/utils/dataset.py b/konfai/utils/dataset.py index db18d6c0..beed3291 100644 --- a/konfai/utils/dataset.py +++ b/konfai/utils/dataset.py @@ -16,6 +16,8 @@ """Dataset file abstractions and image conversion utilities for KonfAI.""" +from __future__ import annotations + import ast import copy import csv @@ -26,14 +28,22 @@ from pathlib import Path from typing import Any -import h5py import numpy as np -import SimpleITK as sitk # noqa: N813 import torch from lxml import etree # nosec B410 +try: + import h5py +except ImportError: + h5py = None # type: ignore[assignment] +try: + import SimpleITK as sitk +except ImportError: + sitk = None # type: ignore[assignment] + from konfai import current_date -from konfai.utils.utils import SUPPORTED_EXTENSIONS +from konfai.utils.errors import DatasetManagerError +from konfai.utils.utils import SUPPORTED_EXTENSIONS, split_format_level class Attribute(dict[str, Any]): @@ -45,8 +55,11 @@ def __init__(self, attributes: dict[str, Any] | None = None) -> None: for k, v in attributes.items(): super().__setitem__(copy.deepcopy(k), copy.deepcopy(v)) + def _count_key(self, key: str) -> int: + return len([k for k in super().keys() if k.startswith(key)]) + def __getitem__(self, key: str) -> Any: - i = len([k for k in super().keys() if k.startswith(key)]) + i = self._count_key(key) if i > 0 and f"{key}_{i - 1}" in super().keys(): return str(super().__getitem__(f"{key}_{i - 1}")) else: @@ -54,7 +67,7 @@ def __getitem__(self, key: str) -> Any: def __setitem__(self, key: str, value: Any) -> None: if "_" not in key: - i = len([k for k in super().keys() if k.startswith(key)]) + i = self._count_key(key) result = None if isinstance(value, torch.Tensor): result = str(value.numpy()) @@ -72,22 +85,22 @@ def __setitem__(self, key: str, value: Any) -> None: super().__setitem__(key, result) def pop(self, key: str, default: Any = None) -> Any: - i = len([k for k in super().keys() if k.startswith(key)]) + i = self._count_key(key) if i > 0 and f"{key}_{i - 1}" in super().keys(): return super().pop(f"{key}_{i - 1}") else: raise NameError(f"{key} not in cache_attribute") - def get_np_array(self, key) -> np.ndarray: + def get_np_array(self, key: str) -> np.ndarray: return np.fromstring(self[key][1:-1], sep=" ", dtype=np.double) - def get_tensor(self, key) -> torch.Tensor: + def get_tensor(self, key: str) -> torch.Tensor: return torch.tensor(self.get_np_array(key)).to(torch.float32) - def pop_np_array(self, key): + def pop_np_array(self, key: str) -> np.ndarray: return np.fromstring(self.pop(key)[1:-1], sep=" ", dtype=np.double) - def pop_tensor(self, key) -> torch.Tensor: + def pop_tensor(self, key: str) -> torch.Tensor: return torch.tensor(self.pop_np_array(key)) def __contains__(self, key: object) -> bool: @@ -139,7 +152,7 @@ def _finalize_running_statistics(state: dict[str, float] | None) -> dict[str, fl } -def is_an_image(attributes: Attribute): +def is_an_image(attributes: Attribute) -> bool: """Return whether the given attribute set contains image geometry metadata.""" return "Origin" in attributes and "Spacing" in attributes and "Direction" in attributes @@ -193,7 +206,7 @@ def get_infos(filename: str | Path) -> tuple[list[int], Attribute]: size = list(file_reader.GetSize()) if len(size) == 3: size = list(reversed(size)) - size = [file_reader.GetNumberOfComponents()] + size + size = [file_reader.GetNumberOfComponents(), *size] return size, attributes @@ -238,7 +251,6 @@ class Dataset: """Filesystem or HDF5-backed dataset abstraction used across KonfAI.""" class AbstractFile(ABC): - @abstractmethod def __init__(self) -> None: pass @@ -294,7 +306,6 @@ def get_infos(self, group: str, name: str) -> tuple[list[int], Attribute]: pass class H5File(AbstractFile): - def __init__(self, filename: str, read: bool) -> None: self.h5: h5py.File | None = None self.filename = filename @@ -471,7 +482,6 @@ def get_infos(self, groups: str, name: str) -> tuple[list[int], Attribute]: ) class SitkFile(AbstractFile): - def __init__(self, filename: str, read: bool, file_format: str) -> None: self.filename = filename self.read = read @@ -483,7 +493,7 @@ def _normalize_slices(slices: tuple[slice, ...], shape: list[int]) -> tuple[slic raise ValueError(f"Expected {len(shape)} slices, got {len(slices)}.") normalized = [] - for item, size in zip(slices, shape): + for item, size in zip(slices, shape, strict=False): start, stop, step = item.indices(size) normalized.append(slice(start, stop, step)) return tuple(normalized) @@ -513,7 +523,7 @@ def _file_to_image_slice(self, name: str, path: str, slices: tuple[slice, ...]) spatial_size_xyz = list(reader.GetSize()) spatial_shape = list(reversed(spatial_size_xyz)) - data_shape = [reader.GetNumberOfComponents()] + spatial_shape + data_shape = [reader.GetNumberOfComponents(), *spatial_shape] normalized = self._normalize_slices(slices, data_shape) if not self._supports_direct_slice(normalized): @@ -639,7 +649,7 @@ def file_to_data_statistics( data = data[channels] return _finalize_running_statistics(_update_running_statistics(None, data)) - def is_vtk_polydata(self, obj): + def is_vtk_polydata(self, obj) -> bool: try: import vtk @@ -720,7 +730,7 @@ def is_exist(self, group: str, name: str | None = None) -> bool: def get_names(self, group: str) -> list[str]: raise NotImplementedError() - def get_group(self): + def get_group(self) -> list[str]: raise NotImplementedError() def get_infos(self, group: str, name: str) -> tuple[list[int], Attribute]: @@ -737,23 +747,286 @@ def get_infos(self, group: str, name: str) -> tuple[list[int], Attribute]: size = list(file_reader.GetSize()) if len(size) == 3: size = list(reversed(size)) - size = [file_reader.GetNumberOfComponents()] + size + size = [file_reader.GetNumberOfComponents(), *size] else: data, attributes = self.file_to_data(group if group is not None else "", name) size = data.shape return size, attributes - class File: + class OmeZarrFile(AbstractFile): + """OME-NGFF backend using chunked Zarr reads for KonfAI patches. - def __init__(self, filename: str, read: bool, file_format: str) -> None: + ``level`` selects the multiscale pyramid resolution to read (0 = full + resolution, higher = coarser); it comes from the ``omezarr@`` + dataset-spec suffix. + """ + + def __init__(self, filename: str, read: bool, level: int = 0) -> None: + self.filename = filename if filename.endswith("/") else f"{filename}/" + self.read = read + self.level = level + + def __enter__(self): + return self + + def __exit__(self, exc_type, value, traceback): + return None + + def _path(self, name: str, *, writing: bool = False) -> Path: + base = Path(self.filename) / name + if writing: + return Path(f"{base}.ome.zarr") + candidates = [Path(f"{base}.ome.zarr"), Path(f"{base}.zarr"), base] + for candidate in candidates: + if candidate.is_dir(): + return candidate + raise NameError(f"OME-Zarr group '{name}' not found in '{self.filename}'.") + + @staticmethod + def _attributes(metadata: dict[str, Any]) -> Attribute: + attributes = Attribute(metadata.get("attributes", {})) + axes = metadata["axes"] + scale = dict(zip(axes, metadata.get("scale", []), strict=False)) + translation = dict(zip(axes, metadata.get("translation", []), strict=False)) + spatial_axes = [axis for axis in ("x", "y", "z") if axis in axes] + if "Spacing" not in attributes: + attributes["Spacing"] = np.asarray([scale.get(axis, 1.0) for axis in spatial_axes]) + if "Origin" not in attributes: + attributes["Origin"] = np.asarray([translation.get(axis, 0.0) for axis in spatial_axes]) + if "Direction" not in attributes: + attributes["Direction"] = np.eye(len(spatial_axes), dtype=np.float64).flatten() + attributes["OMEAxes"] = np.asarray(axes) + return attributes + + def file_to_data(self, group: str, name: str) -> tuple[np.ndarray, Attribute]: + info_shape, _ = self.get_infos(group, name) + return self.file_to_data_slice(group, name, tuple(slice(None) for _ in info_shape)) + + def file_to_data_slice(self, group: str, name: str, slices: tuple[slice, ...]) -> tuple[np.ndarray, Attribute]: + from konfai.utils.ome_zarr import read_ome_zarr_data_slice + + path = self._path(name) + data, metadata = read_ome_zarr_data_slice(path, slices, level=self.level) + attributes = self._attributes(metadata) + shape = metadata["shape"] + normalized = tuple(slice(*item.indices(size)) for item, size in zip(slices, shape, strict=True)) + spacing = attributes.get_np_array("Spacing") + direction = attributes.get_np_array("Direction").reshape(len(spacing), len(spacing)) + start_xyz = np.asarray([item.start for item in reversed(normalized[1:])], dtype=np.float64) + step_xyz = np.asarray([item.step for item in reversed(normalized[1:])], dtype=np.float64) + attributes["Origin"] = attributes.get_np_array("Origin") + direction @ (start_xyz * spacing) + attributes["Spacing"] = spacing * step_xyz + return data, attributes + + def file_to_data_statistics( + self, + group: str, + name: str, + channels: list[int] | None = None, + ) -> dict[str, float]: + shape, _ = self.get_infos(group, name) + trailing_size = int(np.prod(shape[2:], dtype=np.int64)) if len(shape) > 2 else 1 + chunk_length = max(1, 8_000_000 // max(1, trailing_size)) + state: dict[str, float] | None = None + for start in range(0, shape[1], chunk_length): + slices = [slice(None)] * len(shape) + slices[1] = slice(start, min(shape[1], start + chunk_length)) + chunk, _ = self.file_to_data_slice(group, name, tuple(slices)) + if channels is not None: + chunk = chunk[channels] + state = _update_running_statistics(state, chunk) + return _finalize_running_statistics(state) + + def data_to_file( + self, + name: str, + data: sitk.Image | sitk.Transform | np.ndarray, + attributes: Attribute | None = None, + ) -> None: + from konfai.utils.ome_zarr import write_ome_zarr + + attributes = attributes or Attribute() + if sitk is not None and isinstance(data, sitk.Image): + data, image_attributes = image_to_data(data) + attributes.update(image_attributes) + if not isinstance(data, np.ndarray): + raise DatasetManagerError("OME-Zarr datasets can only store image arrays.") + dimension = data.ndim - 1 + spacing = attributes.get_np_array("Spacing") if "Spacing" in attributes else np.ones(dimension) + origin = attributes.get_np_array("Origin") if "Origin" in attributes else np.zeros(dimension) + write_ome_zarr( + self._path(name, writing=True), + data, + spacing=spacing, + origin=origin, + attributes=dict(attributes), + ) + + def get_names(self, group: str) -> list[str]: + return self.get_group() + + def get_group(self) -> list[str]: + root = Path(self.filename) + if not root.is_dir(): + return [] + groups = [] + for path in root.iterdir(): + if path.name.endswith(".ome.zarr"): + groups.append(path.name.removesuffix(".ome.zarr")) + elif path.name.endswith(".zarr"): + groups.append(path.name.removesuffix(".zarr")) + return sorted(groups) + + def is_exist(self, group: str, name: str | None = None) -> bool: + try: + self._path(f"{group}/{name}" if name else group) + return True + except NameError: + return False + + def get_infos(self, group: str, name: str) -> tuple[list[int], Attribute]: + from konfai.utils.ome_zarr import get_ome_zarr_info + + metadata = get_ome_zarr_info(self._path(name), level=self.level) + axes = [str(axis).lower() for axis in metadata["axes"]] + axis_sizes = dict(zip(axes, metadata["shape"], strict=True)) + shape = [axis_sizes.get("c", 1), *[axis_sizes[axis] for axis in ("z", "y", "x") if axis in axis_sizes]] + metadata["shape"] = shape + return shape, self._attributes(metadata) + + class DicomFile(AbstractFile): + """DICOM series backend with header-only metadata and slice-level reads.""" + + def __init__(self, filename: str, read: bool) -> None: + self.filename = filename if filename.endswith("/") else f"{filename}/" + self.read = read + + def __enter__(self): + return self + + def __exit__(self, exc_type, value, traceback): + return None + + def _path(self, name: str) -> Path: + return Path(self.filename) / name + + @staticmethod + def _attributes(info: dict[str, Any]) -> Attribute: + attributes = Attribute() + attributes["Origin"] = np.asarray(info["origin"]) + attributes["Spacing"] = np.asarray(info["spacing"]) + attributes["Direction"] = np.asarray(info["direction"]) + attributes["SeriesInstanceUID"] = info["series_uid"] + return attributes + + def file_to_data(self, group: str, name: str) -> tuple[np.ndarray, Attribute]: + from konfai.utils.dicom import read_dicom_series + + data, origin, spacing, direction = read_dicom_series(self._path(name)) + attributes = Attribute() + attributes["Origin"] = origin + attributes["Spacing"] = spacing + attributes["Direction"] = direction + return data, attributes + + def file_to_data_slice(self, group: str, name: str, slices: tuple[slice, ...]) -> tuple[np.ndarray, Attribute]: + from konfai.utils.dicom import get_dicom_info, read_dicom_series_slice + + path = self._path(name) + info = get_dicom_info(path) + data, origin, spacing, direction = read_dicom_series_slice(path, slices, series_uid=info["series_uid"]) + info.update(origin=origin, spacing=spacing, direction=direction) + return data, self._attributes(info) + + def file_to_data_statistics( + self, + group: str, + name: str, + channels: list[int] | None = None, + ) -> dict[str, float]: + shape, _ = self.get_infos(group, name) + state: dict[str, float] | None = None + for index in range(shape[1]): + chunk, _ = self.file_to_data_slice( + group, + name, + (slice(None), slice(index, index + 1), slice(None), slice(None)), + ) + if channels is not None: + chunk = chunk[channels] + state = _update_running_statistics(state, chunk) + return _finalize_running_statistics(state) + + def data_to_file( + self, + name: str, + data: sitk.Image | sitk.Transform | np.ndarray, + attributes: Attribute | None = None, + ) -> None: + from konfai.utils.dicom import write_dicom_series + + attributes = attributes or Attribute() + if sitk is not None and isinstance(data, sitk.Image): + data, image_attributes = image_to_data(data) + attributes.update(image_attributes) + if not isinstance(data, np.ndarray): + raise DatasetManagerError("DICOM datasets can only store scalar image arrays.") + spacing = attributes.get_np_array("Spacing") if "Spacing" in attributes else np.ones(3) + origin = attributes.get_np_array("Origin") if "Origin" in attributes else np.zeros(3) + direction = attributes.get_np_array("Direction") if "Direction" in attributes else np.eye(3).flatten() + metadata = { + key: attributes[key] + for key in ("PatientName", "PatientID", "Modality", "StudyInstanceUID", "SeriesInstanceUID") + if key in attributes + } + write_dicom_series( + self._path(name), + data, + spacing=spacing, + origin=origin, + direction=direction, + metadata=metadata, + ) + + def get_names(self, group: str) -> list[str]: + return self.get_group() + + def get_group(self) -> list[str]: + root = Path(self.filename) + if not root.is_dir(): + return [] + return sorted(path.name for path in root.iterdir() if path.is_dir() and self.is_exist(path.name)) + + def is_exist(self, group: str, name: str | None = None) -> bool: + from konfai.utils.dicom import get_dicom_info + + try: + get_dicom_info(self._path(f"{group}/{name}" if name else group)) + return True + except DatasetManagerError: + return False + + def get_infos(self, group: str, name: str) -> tuple[list[int], Attribute]: + from konfai.utils.dicom import get_dicom_info + + info = get_dicom_info(self._path(name)) + return info["shape"], self._attributes(info) + + class File: + def __init__(self, filename: str, read: bool, file_format: str, level: int = 0) -> None: self.filename = filename self.read = read - self.file: "Dataset.AbstractFile" | None = None + self.file: Dataset.AbstractFile | None = None self.file_format = file_format + self.level = level - def __enter__(self) -> "Dataset.AbstractFile": + def __enter__(self) -> Dataset.AbstractFile: if self.file_format == "h5": self.file = Dataset.H5File(self.filename, self.read) + elif self.file_format == "omezarr": + self.file = Dataset.OmeZarrFile(self.filename, self.read, self.level) + elif self.file_format == "dicom": + self.file = Dataset.DicomFile(self.filename, self.read) else: self.file = Dataset.SitkFile(self.filename + "/", self.read, self.file_format) self.file.__enter__() @@ -764,11 +1037,15 @@ def __exit__(self, exc_type, value, traceback): self.file.__exit__(exc_type, value, traceback) def __init__(self, filename: str | Path, file_format: str) -> None: + base_format, self.level = split_format_level(file_format) + normalized_format = base_format.lower().removeprefix(".").replace("_", "-") + file_format = {"ome-zarr": "omezarr", "zarr": "omezarr"}.get(normalized_format, normalized_format) if file_format != "h5" and not str(filename).endswith("/"): filename = f"{filename}/" self.is_directory = str(filename).endswith("/") self.filename = str(filename) self.file_format = file_format + self._names_cache: dict[str, list[str]] = {} def _exists_on_disk(self) -> bool: if os.path.exists(self.filename): @@ -781,7 +1058,8 @@ def write( name: str, data: sitk.Image | sitk.Transform | np.ndarray, attributes: Attribute | None = None, - ): + ) -> None: + self._names_cache.clear() if attributes is None: attributes = Attribute() if self.is_directory: @@ -793,10 +1071,10 @@ def write( sub_directory = "/".join(s_group[:-1]) name = f"{sub_directory}/{name}" group = s_group[-1] - with Dataset.File(f"{self.filename}{name}", False, self.file_format) as file: + with Dataset.File(f"{self.filename}{name}", False, self.file_format, self.level) as file: file.data_to_file(group, data, attributes) else: - with Dataset.File(self.filename, False, self.file_format) as file: + with Dataset.File(self.filename, False, self.file_format, self.level) as file: file.data_to_file(f"{group}/{name}", data, attributes) def read_data(self, groups: str, name: str) -> tuple[np.ndarray, Attribute]: @@ -810,12 +1088,13 @@ def read_data(self, groups: str, name: str) -> tuple[np.ndarray, Attribute]: f"{self.filename}{sub_directory}{name}", False, self.file_format, + self.level, ) as file: - result = file.file_to_data("", group) + return file.file_to_data("", group) else: - with Dataset.File(self.filename, False, self.file_format) as file: - result = file.file_to_data(groups, name) - return result + with Dataset.File(self.filename, False, self.file_format, self.level) as file: + return file.file_to_data(groups, name) + raise NameError(f"Dataset entry '{groups}/{name}' not found in {self.filename}.") def read_data_slice(self, groups: str, name: str, slices: tuple[slice, ...]) -> tuple[np.ndarray, Attribute]: if not self._exists_on_disk(): @@ -828,11 +1107,12 @@ def read_data_slice(self, groups: str, name: str, slices: tuple[slice, ...]) -> f"{self.filename}{sub_directory}{name}", True, self.file_format, + self.level, ) as file: result = file.file_to_data_slice("", group, slices) return result else: - with Dataset.File(self.filename, True, self.file_format) as file: + with Dataset.File(self.filename, True, self.file_format, self.level) as file: return file.file_to_data_slice(groups, name, slices) raise NameError(f"Dataset entry '{groups}/{name}' not found in {self.filename}.") @@ -853,10 +1133,11 @@ def read_data_statistics( f"{self.filename}{sub_directory}{name}", True, self.file_format, + self.level, ) as file: return file.file_to_data_statistics("", group, channels) else: - with Dataset.File(self.filename, True, self.file_format) as file: + with Dataset.File(self.filename, True, self.file_format, self.level) as file: return file.file_to_data_statistics(groups, name, channels) raise NameError(f"Dataset entry '{groups}/{name}' not found in {self.filename}.") @@ -879,7 +1160,7 @@ def read_transform(self, group: str, name: str) -> sitk.Transform: transforms.append(transform) return sitk.CompositeTransform(transforms) if len(transforms) > 1 else transforms[0] - def read_image(self, group: str, name: str): + def read_image(self, group: str, name: str) -> sitk.Image: data, attribute = self.read_data(group, name) return data_to_image(data, attribute) @@ -913,6 +1194,9 @@ def _get_sub_directories(self, groups: str, sub_directory: str = ""): return sub_directories def get_names(self, groups: str, index: list[int] | None = None) -> list[str]: + if index is None and groups in self._names_cache: + return self._names_cache[groups] + names = [] if self.is_directory: for sub_directory in self._get_sub_directories(groups): @@ -924,27 +1208,41 @@ def get_names(self, groups: str, index: list[int] | None = None) -> list[str]: f"{self.filename}{sub_directory}{name}", True, self.file_format, + self.level, ) as file: if file.is_exist(group): names.append(name.replace(".h5", "") if self.file_format == "h5" else name) else: - with Dataset.File(self.filename, True, self.file_format) as file: + with Dataset.File(self.filename, True, self.file_format, self.level) as file: names = file.get_names(groups) - return [name for i, name in enumerate(sorted(names)) if index is None or i in index] - def get_group(self): + sorted_names = sorted(names) + if index is None: + self._names_cache[groups] = sorted_names + return sorted_names + return [name for i, name in enumerate(sorted_names) if i in index] + + def get_group(self) -> list[str]: if self.is_directory: + if self.file_format in {"dicom", "omezarr"}: + groups_set = set() + root_path = Path(self.filename) + for case_path in root_path.iterdir() if root_path.is_dir() else []: + if case_path.is_dir(): + with Dataset.File(str(case_path), True, self.file_format, self.level) as dataset_file: + groups_set.update(dataset_file.get_group()) + return sorted(groups_set) groups_set = set() - for root, _, files in os.walk(self.filename): + for root_dir, _, files in os.walk(self.filename): for file in files: - path = Path(root, file.split(".")[0]).relative_to(self.filename).as_posix() + path = Path(root_dir, file.split(".")[0]).relative_to(self.filename).as_posix() parts = path.split("/") if len(parts) >= 2: del parts[-2] groups_set.add("/".join(parts)) groups = list(groups_set) else: - with Dataset.File(self.filename, True, self.file_format) as dataset_file: + with Dataset.File(self.filename, True, self.file_format, self.level) as dataset_file: groups = dataset_file.get_group() return list(groups) @@ -957,12 +1255,13 @@ def get_infos(self, groups: str, name: str) -> tuple[list[int], Attribute]: f"{self.filename}{sub_directory}{name}", True, self.file_format, + self.level, ) as file: - result = file.get_infos("", group) + return file.get_infos("", group) else: - with Dataset.File(self.filename, True, self.file_format) as file: - result = file.get_infos(groups, name) - return result + with Dataset.File(self.filename, True, self.file_format, self.level) as file: + return file.get_infos(groups, name) + raise NameError(f"Dataset entry '{groups}/{name}' not found in {self.filename}.") def get_statistics(self, groups: str) -> dict[str, dict[str, dict[str, float | list[float]]]]: names = self.get_names(groups) diff --git a/konfai/utils/dicom.py b/konfai/utils/dicom.py new file mode 100644 index 00000000..b25f49d0 --- /dev/null +++ b/konfai/utils/dicom.py @@ -0,0 +1,554 @@ +# Copyright (c) 2025 Valentin Boussot +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# SPDX-License-Identifier: Apache-2.0 + +"""DICOM series reader for KonfAI medical imaging pipelines. + +Design rationale +---------------- +DICOM is not a folder of independent images. A CT or MRI acquisition is a +*series* — a collection of .dcm files that together define a 3-D volume. +Reading a DICOM correctly requires: + +1. **Series discovery** — group files by SeriesInstanceUID. A folder may + contain multiple series (e.g. a T1 and a T2 acquired in the same session). + +2. **Slice ordering** — sort slices by ImagePositionPatient (z-component along + ImageOrientationPatient normal vector), not by filename or InstanceNumber, + which can be unreliable. + +3. **Geometry extraction** — derive spacing_mm (PixelSpacing + SliceThickness / + derived inter-slice distance), origin (ImagePositionPatient of first slice), + and direction cosines (ImageOrientationPatient rows and columns + cross + product for the z-axis). + +4. **CT intensity rescale** — apply RescaleSlope and RescaleIntercept to + convert stored pixel values to Hounsfield Units (HU). This is mandatory + for CT and is absent (or identity) for MR. + +5. **Error handling** — missing tags, single-slice series, inconsistent spacing, + non-square pixels, and unsupported transfer syntaxes all need clear messages. + +Optional dependency: ``pydicom`` (``pip install konfai[dicom]``). +""" + +from __future__ import annotations + +import os +from collections.abc import Sequence +from datetime import datetime +from pathlib import Path +from typing import TYPE_CHECKING, Any + +import numpy as np + +from konfai.utils.errors import DatasetManagerError + +if TYPE_CHECKING: + pass + +try: + import pydicom + from pydicom.dataset import Dataset as DicomDataset + from pydicom.sequence import Sequence as DicomSequence + + _PYDICOM_AVAILABLE = True +except ImportError: + _PYDICOM_AVAILABLE = False + pydicom = None # type: ignore[assignment] + DicomDataset = None # type: ignore[assignment,misc] + DicomSequence = None # type: ignore[assignment] + + +def _require_pydicom() -> None: + if not _PYDICOM_AVAILABLE: + raise DatasetManagerError( + "pydicom is required for DICOM support.", + "Install it with: pip install konfai[dicom]", + ) + + +# --------------------------------------------------------------------------- +# Series discovery +# --------------------------------------------------------------------------- + + +def discover_series(directory: str | Path) -> dict[str, list[Path]]: + """Return a mapping of SeriesInstanceUID -> sorted list of .dcm paths. + + Parameters + ---------- + directory: + Root directory to scan recursively for .dcm files. + + Returns + ------- + dict[str, list[Path]] + Keys are SeriesInstanceUID values; values are lists of file paths + belonging to that series (unsorted at this stage). + + Raises + ------ + DatasetManagerError + If ``pydicom`` is not installed or the directory contains no DICOM. + """ + _require_pydicom() + + root = Path(directory) + if not root.is_dir(): + raise DatasetManagerError(f"DICOM directory '{root}' does not exist or is not a directory.") + + series: dict[str, list[Path]] = {} + for dirpath, _, filenames in os.walk(root): + for fname in filenames: + fpath = Path(dirpath) / fname + try: + ds = pydicom.dcmread(str(fpath), stop_before_pixels=True) + uid = str(ds.SeriesInstanceUID) + series.setdefault(uid, []).append(fpath) + except Exception: # nosec B112 + # Skip unreadable or non-DICOM files; discovery must not crash on stray content. + continue + + if not series: + raise DatasetManagerError( + f"No DICOM files found under '{root}'.", + "Ensure the directory contains .dcm files from a valid DICOM series.", + ) + return series + + +# --------------------------------------------------------------------------- +# Slice sorting +# --------------------------------------------------------------------------- + + +def _slice_position(ds: DicomDataset) -> float: + """Return the signed position of one slice along the acquisition axis. + + Uses ``ImagePositionPatient`` projected onto the slice-normal derived from + ``ImageOrientationPatient``. Falls back to ``InstanceNumber`` (unreliable + but ubiquitous) when geometry tags are absent. + """ + try: + iop = [float(x) for x in ds.ImageOrientationPatient] + ipp = [float(x) for x in ds.ImagePositionPatient] + row = np.array(iop[:3]) + col = np.array(iop[3:]) + normal = np.cross(row, col) + return float(np.dot(normal, ipp)) + except AttributeError: + try: + return float(ds.InstanceNumber) + except AttributeError: + return 0.0 + + +def sort_series(files: list[Path], *, stop_before_pixels: bool = False) -> list[DicomDataset]: + """Read and sort slices in anatomical order (ascending slice position). + + Parameters + ---------- + files: + Unsorted list of paths belonging to one DICOM series. + + Returns + ------- + list[DicomDataset] + Datasets sorted by their position along the acquisition normal. + """ + _require_pydicom() + datasets = [pydicom.dcmread(str(f), stop_before_pixels=stop_before_pixels) for f in files] + datasets.sort(key=_slice_position) + return datasets + + +def _select_series_files(directory: str | Path, series_uid: str | None = None) -> tuple[str, list[Path]]: + all_series = discover_series(directory) + if series_uid is not None: + if series_uid not in all_series: + raise DatasetManagerError( + f"Series '{series_uid}' not found in '{directory}'.", + f"Available series: {list(all_series.keys())}", + ) + return series_uid, all_series[series_uid] + if len(all_series) == 1: + selected_uid, files = next(iter(all_series.items())) + return selected_uid, files + raise DatasetManagerError( + f"Multiple DICOM series found in '{directory}' ({len(all_series)} series).", + "Specify 'series_uid' to select one.", + f"Available UIDs: {list(all_series.keys())}", + ) + + +# --------------------------------------------------------------------------- +# Geometry extraction +# --------------------------------------------------------------------------- + + +def extract_geometry( + datasets: list[DicomDataset], +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Extract origin, spacing, and direction from a sorted DICOM series. + + Parameters + ---------- + datasets: + Slice datasets in anatomical order (from :func:`sort_series`). + + Returns + ------- + tuple[np.ndarray, np.ndarray, np.ndarray] + - ``origin`` (3,) — physical position of the first voxel (mm). + - ``spacing`` (3,) — KonfAI/SimpleITK order (x, y, z) in mm. + - ``direction`` (9,) — row-major 3-by-3 direction cosine matrix, flattened. + + Raises + ------ + DatasetManagerError + If required geometry tags are missing or inconsistent. + """ + if not datasets: + raise DatasetManagerError("Cannot extract geometry from an empty series.") + + first = datasets[0] + + # Origin = ImagePositionPatient of first slice + try: + origin = np.array([float(x) for x in first.ImagePositionPatient], dtype=np.float64) + except AttributeError as exc: + raise DatasetManagerError( + "DICOM tag 'ImagePositionPatient' is missing.", + "This tag is required to determine the volume origin.", + ) from exc + + # In-plane spacing from PixelSpacing + try: + pixel_spacing = [float(x) for x in first.PixelSpacing] + row_spacing_mm = pixel_spacing[0] + col_spacing_mm = pixel_spacing[1] + except AttributeError as exc: + raise DatasetManagerError( + "DICOM tag 'PixelSpacing' is missing.", + "This tag is required to determine voxel dimensions.", + ) from exc + + # Slice spacing: prefer computed inter-slice distance over SliceThickness + if len(datasets) > 1: + pos_first = _slice_position(datasets[0]) + pos_second = _slice_position(datasets[1]) + slice_spacing_mm = abs(pos_second - pos_first) + else: + try: + slice_spacing_mm = float(first.SliceThickness) + except AttributeError: + slice_spacing_mm = 1.0 # fallback; single-slice series + + spacing = np.array([col_spacing_mm, row_spacing_mm, slice_spacing_mm], dtype=np.float64) + + # Direction cosines + try: + iop = [float(x) for x in first.ImageOrientationPatient] + except AttributeError as exc: + raise DatasetManagerError( + "DICOM tag 'ImageOrientationPatient' is missing.", + "This tag is required to determine the volume orientation.", + ) from exc + + row_cosine = np.array(iop[:3], dtype=np.float64) + col_cosine = np.array(iop[3:], dtype=np.float64) + normal_cosine = np.cross(row_cosine, col_cosine) + direction = np.column_stack([row_cosine, col_cosine, normal_cosine]).flatten() + + return origin, spacing, direction + + +# --------------------------------------------------------------------------- +# Pixel reading with CT rescale +# --------------------------------------------------------------------------- + + +def read_volume( + datasets: list[DicomDataset], + *, + apply_rescale: bool = True, +) -> np.ndarray: + """Stack sorted slices into a channel-first (1, Z, Y, X) float32 array. + + Parameters + ---------- + datasets: + Sorted DICOM datasets (from :func:`sort_series`). + apply_rescale: + If True, apply RescaleSlope / RescaleIntercept to convert stored + pixel values to Hounsfield Units (HU) for CT, or to physical signal + units for modalities that provide these tags. Set to False to keep + raw stored pixel integers (e.g., for label maps or QC). + + Returns + ------- + np.ndarray + Shape (1, Z, Y, X), dtype float32. Channel dimension = 1 for scalar + volumes. + + Raises + ------ + DatasetManagerError + If pixel data cannot be read or slices have inconsistent shapes. + """ + slices = [] + expected_shape: tuple[int, int] | None = None + + for i, ds in enumerate(datasets): + try: + arr = ds.pixel_array.astype(np.float32) + except Exception as exc: + raise DatasetManagerError( + f"Cannot read pixel data from DICOM slice {i}.", + f"Transfer syntax or compression may be unsupported: {exc}", + ) from exc + + if expected_shape is None: + expected_shape = arr.shape + elif arr.shape != expected_shape: + raise DatasetManagerError( + f"Inconsistent slice shape at index {i}: expected {expected_shape}, got {arr.shape}.", + "All slices in a series must have the same rows and columns.", + ) + + if apply_rescale: + slope = float(getattr(ds, "RescaleSlope", 1.0)) + intercept = float(getattr(ds, "RescaleIntercept", 0.0)) + arr = arr * slope + intercept + + slices.append(arr) + + if not slices: + raise DatasetManagerError("Series contains no readable slices.") + + volume = np.stack(slices, axis=0) # (Z, Y, X) + return volume[np.newaxis] # (1, Z, Y, X) + + +def get_dicom_info( + directory: str | Path, + *, + series_uid: str | None = None, +) -> dict[str, Any]: + """Read DICOM series shape and geometry without decoding pixel data.""" + selected_uid, files = _select_series_files(directory, series_uid) + datasets = sort_series(files, stop_before_pixels=True) + origin, spacing, direction = extract_geometry(datasets) + first = datasets[0] + try: + rows = int(first.Rows) + columns = int(first.Columns) + except AttributeError as exc: + raise DatasetManagerError("DICOM Rows/Columns tags are required to determine the volume shape.") from exc + return { + "series_uid": selected_uid, + "files": files, + "shape": [1, len(datasets), rows, columns], + "origin": origin, + "spacing": spacing, + "direction": direction, + } + + +def read_dicom_series_slice( + directory: str | Path, + slices: tuple[slice, ...], + *, + series_uid: str | None = None, + apply_rescale: bool = True, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Read only the selected DICOM slices and return updated patch geometry.""" + info = get_dicom_info(directory, series_uid=series_uid) + shape = info["shape"] + if len(slices) != len(shape): + raise DatasetManagerError(f"Expected {len(shape)} slices, got {len(slices)}.") + normalized = tuple(slice(*item.indices(size)) for item, size in zip(slices, shape, strict=True)) + channel_indices = range(*normalized[0].indices(shape[0])) + if list(channel_indices) not in ([0], []): + raise DatasetManagerError("DICOM stores scalar data and supports only channel 0.") + + _selected_uid, files = _select_series_files(directory, series_uid or info["series_uid"]) + headers = sort_series(files, stop_before_pixels=True) + z_indices = list(range(*normalized[1].indices(shape[1]))) + selected_files = [Path(headers[index].filename) for index in z_indices] + datasets = sort_series(selected_files) + volume = read_volume(datasets, apply_rescale=apply_rescale) + volume = volume[normalized[0], :, normalized[2], normalized[3]] + + direction_matrix = np.asarray(info["direction"], dtype=np.float64).reshape(3, 3) + start_xyz = np.asarray([normalized[3].start, normalized[2].start, normalized[1].start], dtype=np.float64) + spacing = np.asarray(info["spacing"], dtype=np.float64) + origin = np.asarray(info["origin"], dtype=np.float64) + direction_matrix @ (start_xyz * spacing) + step_xyz = np.asarray([normalized[3].step, normalized[2].step, normalized[1].step], dtype=np.float64) + return volume, origin, spacing * step_xyz, np.asarray(info["direction"], dtype=np.float64) + + +def _encode_pixels(data: np.ndarray) -> tuple[np.ndarray, float, float]: + """Encode one scalar volume in an uncompressed integer DICOM representation.""" + if np.issubdtype(data.dtype, np.floating): + finite = data[np.isfinite(data)] + if finite.size == 0: + raise DatasetManagerError("Cannot write a DICOM volume containing no finite values.") + minimum = float(finite.min()) + maximum = float(finite.max()) + slope = (maximum - minimum) / 65535.0 if maximum > minimum else 1.0 + intercept = minimum + 32768.0 * slope + stored = np.rint((np.nan_to_num(data, nan=minimum) - intercept) / slope).clip(-32768, 32767).astype(np.int16) + return stored, slope, intercept + if np.issubdtype(data.dtype, np.signedinteger): + return data.astype(np.int32 if data.dtype.itemsize > 2 else np.int16), 1.0, 0.0 + if np.issubdtype(data.dtype, np.unsignedinteger): + return data.astype(np.uint32 if data.dtype.itemsize > 2 else np.uint16), 1.0, 0.0 + raise DatasetManagerError(f"Unsupported DICOM pixel dtype '{data.dtype}'.") + + +def write_dicom_series( + directory: str | Path, + volume: np.ndarray, + *, + origin: Sequence[float] | None = None, + spacing: Sequence[float] | None = None, + direction: Sequence[float] | None = None, + metadata: dict[str, Any] | None = None, +) -> str: + """Write a scalar ``C-Z-Y-X`` volume as an uncompressed DICOM series.""" + _require_pydicom() + from pydicom.dataset import FileDataset, FileMetaDataset + from pydicom.uid import CTImageStorage, ExplicitVRLittleEndian, generate_uid + + data = np.asarray(volume) + if data.ndim == 3: + data = data[np.newaxis] + if data.ndim != 4 or data.shape[0] != 1: + raise DatasetManagerError(f"DICOM writing expects one scalar C-Z-Y-X channel, got shape {data.shape}.") + + stored, slope, intercept = _encode_pixels(data[0]) + origin_array = np.asarray(origin if origin is not None else [0.0, 0.0, 0.0], dtype=np.float64) + spacing_array = np.asarray(spacing if spacing is not None else [1.0, 1.0, 1.0], dtype=np.float64) + direction_matrix = np.asarray( + direction if direction is not None else np.eye(3).flatten(), dtype=np.float64 + ).reshape(3, 3) + if origin_array.shape != (3,) or spacing_array.shape != (3,): + raise DatasetManagerError("DICOM origin and spacing must each contain exactly three values.") + + root = Path(directory) + root.mkdir(parents=True, exist_ok=True) + for existing in root.glob("*.dcm"): + existing.unlink() + + metadata = dict(metadata or {}) + study_uid = str(metadata.get("StudyInstanceUID", generate_uid())) + series_uid = str(metadata.get("SeriesInstanceUID", generate_uid())) + frame_uid = str(metadata.get("FrameOfReferenceUID", generate_uid())) + now = datetime.now() + bits = stored.dtype.itemsize * 8 + signed = bool(np.issubdtype(stored.dtype, np.signedinteger)) + row_cosine = direction_matrix[:, 0] + column_cosine = direction_matrix[:, 1] + + for index, pixels in enumerate(stored): + sop_uid = generate_uid() + file_meta = FileMetaDataset() + file_meta.FileMetaInformationVersion = b"\x00\x01" + file_meta.MediaStorageSOPClassUID = CTImageStorage + file_meta.MediaStorageSOPInstanceUID = sop_uid + file_meta.TransferSyntaxUID = ExplicitVRLittleEndian + file_meta.ImplementationClassUID = generate_uid() + path = root / f"{index + 1:06d}.dcm" + dataset = FileDataset(str(path), {}, file_meta=file_meta, preamble=b"\0" * 128) + dataset.SOPClassUID = CTImageStorage + dataset.SOPInstanceUID = sop_uid + dataset.StudyInstanceUID = study_uid + dataset.SeriesInstanceUID = series_uid + dataset.FrameOfReferenceUID = frame_uid + dataset.PatientName = str(metadata.get("PatientName", "KonfAI^Dataset")) + dataset.PatientID = str(metadata.get("PatientID", "KonfAI")) + dataset.Modality = str(metadata.get("Modality", "OT")) + dataset.StudyDate = str(metadata.get("StudyDate", now.strftime("%Y%m%d"))) + dataset.StudyTime = str(metadata.get("StudyTime", now.strftime("%H%M%S"))) + dataset.SeriesNumber = int(metadata.get("SeriesNumber", 1)) + dataset.InstanceNumber = index + 1 + dataset.ImageType = ["DERIVED", "PRIMARY", "AXIAL"] + dataset.Rows = int(pixels.shape[0]) + dataset.Columns = int(pixels.shape[1]) + dataset.SamplesPerPixel = 1 + dataset.PhotometricInterpretation = "MONOCHROME2" + dataset.PixelRepresentation = int(signed) + dataset.BitsAllocated = bits + dataset.BitsStored = bits + dataset.HighBit = bits - 1 + dataset.PixelSpacing = [float(spacing_array[1]), float(spacing_array[0])] + dataset.SliceThickness = float(spacing_array[2]) + dataset.SpacingBetweenSlices = float(spacing_array[2]) + dataset.ImageOrientationPatient = [*row_cosine.tolist(), *column_cosine.tolist()] + position = origin_array + direction_matrix[:, 2] * (index * spacing_array[2]) + dataset.ImagePositionPatient = position.tolist() + dataset.RescaleSlope = float(slope) + dataset.RescaleIntercept = float(intercept) + dataset.PixelData = pixels.tobytes() + dataset.save_as(str(path), enforce_file_format=True) + return series_uid + + +# --------------------------------------------------------------------------- +# High-level convenience function +# --------------------------------------------------------------------------- + + +def read_dicom_series( + directory: str | Path, + *, + series_uid: str | None = None, + apply_rescale: bool = True, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Read a DICOM series from a directory into a channel-first volume. + + Parameters + ---------- + directory: + Path to the folder containing the DICOM series. + series_uid: + If the folder contains multiple series, select by SeriesInstanceUID. + If None and only one series is present, that series is used. + If None and multiple series are present, raises DatasetManagerError. + apply_rescale: + Apply RescaleSlope / RescaleIntercept (True) or keep raw integers. + + Returns + ------- + tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray] + - ``volume`` — shape (1, Z, Y, X), dtype float32. + - ``origin`` — physical origin of the first voxel, mm (shape (3,)). + - ``spacing`` — voxel size in KonfAI/SimpleITK (x, y, z) order (shape (3,)). + - ``direction`` — row-major 3-by-3 direction cosine matrix, flat (shape (9,)). + + Raises + ------ + DatasetManagerError + On missing deps, missing tags, multi-series ambiguity, or read errors. + """ + _require_pydicom() + + _selected_uid, files = _select_series_files(directory, series_uid) + datasets = sort_series(files) + origin, spacing, direction = extract_geometry(datasets) + volume = read_volume(datasets, apply_rescale=apply_rescale) + return volume, origin, spacing, direction diff --git a/konfai/utils/ome_zarr.py b/konfai/utils/ome_zarr.py new file mode 100644 index 00000000..563cbb56 --- /dev/null +++ b/konfai/utils/ome_zarr.py @@ -0,0 +1,219 @@ +# Copyright (c) 2025 Valentin Boussot +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# SPDX-License-Identifier: Apache-2.0 + +"""OME-Zarr (OME-NGFF) read/write backend for KonfAI, built on ``ngff-zarr``. + +This module is a thin adapter: ``ngff-zarr`` owns all OME-NGFF metadata parsing, +multiscale handling, and (de)serialisation — KonfAI does not re-implement the +spec. We only + +1. map between KonfAI's channel-first ``C[Z]YX`` arrays / ``(x, y, z)`` geometry + and ngff-zarr's ``NgffImage`` (axis-named ``scale``/``translation``), and +2. round-trip KonfAI's full ``Attribute`` sidecar (including the ``Direction`` + matrix, which OME-NGFF cannot express) through a single ``konfai`` group + attribute, read/written with ``zarr``. + +Reads are lazy: ``ngff-zarr`` exposes the array as a chunked store, so slicing +only materialises the requested patch. + +Optional dependencies: ``zarr`` + ``ngff-zarr`` (``pip install konfai[omezarr]``). +""" + +from __future__ import annotations + +from collections.abc import Sequence +from pathlib import Path +from typing import Any + +import numpy as np + +from konfai.utils.errors import DatasetManagerError + +try: + import zarr + + _ZARR_AVAILABLE = True +except ImportError: + zarr = None # type: ignore[assignment] + _ZARR_AVAILABLE = False + +try: + import ngff_zarr # type: ignore[import-untyped] + + _NGFF_ZARR_AVAILABLE = True +except ImportError: + ngff_zarr = None # type: ignore[assignment] + _NGFF_ZARR_AVAILABLE = False + +_KONFAI_ATTR_KEY = "konfai" +_SPATIAL = ("z", "y", "x") + + +def _require_zarr() -> None: + if not _ZARR_AVAILABLE: + raise DatasetManagerError( + "zarr is required for OME-Zarr support.", + "Install it with: pip install konfai[omezarr]", + ) + + +def _require_ngff_zarr() -> None: + _require_zarr() + if not _NGFF_ZARR_AVAILABLE: + raise DatasetManagerError( + "ngff-zarr is required for OME-Zarr support.", + "Install it with: pip install konfai[omezarr]", + ) + + +def _read_konfai_attributes(store_path: str | Path) -> dict[str, Any]: + """Read KonfAI's proprietary ``Attribute`` sidecar from the store, if present.""" + try: + group = zarr.open_group(str(store_path), mode="r") + return dict(dict(group.attrs).get(_KONFAI_ATTR_KEY, {}).get("attributes", {})) + except (KeyError, OSError, ValueError, TypeError): + return {} + + +def _load_image(store_path: str | Path, level: int) -> Any: + """Return the ``NgffImage`` for ``level`` of an OME-Zarr store.""" + _require_ngff_zarr() + try: + multiscales = ngff_zarr.from_ngff_zarr(str(store_path)) + return multiscales.images[level] + except (KeyError, IndexError, OSError, TypeError, ValueError) as exc: + raise DatasetManagerError( + f"Cannot open OME-Zarr store '{store_path}' (level {level}).", + "Ensure the directory is a valid OME-NGFF store.", + ) from exc + + +def _canonical_shape(dims: Sequence[str], shape: Sequence[int]) -> list[int]: + """Channel-first ``[C, (Z), Y, X]`` shape derived from ngff dims.""" + axis_size = dict(zip(dims, shape, strict=True)) + return [int(axis_size.get("c", 1)), *[int(axis_size[axis]) for axis in _SPATIAL if axis in axis_size]] + + +def _ordered(values: dict[str, float], dims: Sequence[str]) -> list[float]: + return [float(values.get(axis, 1.0 if axis == "c" else 0.0)) for axis in dims] + + +def read_ome_zarr_data_slice( + store_path: str | Path, + slices: tuple[slice, ...], + *, + level: int = 0, + timepoint: int = 0, +) -> tuple[np.ndarray, dict[str, Any]]: + """Read a KonfAI channel-first ``C[Z]YX`` patch from an OME-Zarr store (lazy).""" + image = _load_image(store_path, level) + dims = [str(axis).lower() for axis in image.dims] + canonical_shape = _canonical_shape(dims, image.data.shape) + if len(slices) != len(canonical_shape): + raise DatasetManagerError(f"Expected {len(canonical_shape)} slices, got {len(slices)}.") + + normalized = [slice(*item.indices(size)) for item, size in zip(slices, canonical_shape, strict=True)] + spatial_slices = dict(zip([axis for axis in _SPATIAL if axis in dims], normalized[1:], strict=True)) + index: list[int | slice] = [] + for axis in dims: + if axis == "t": + index.append(timepoint) + elif axis == "c": + index.append(normalized[0]) + elif axis in spatial_slices: + index.append(spatial_slices[axis]) + else: + index.append(slice(None)) + + patch = np.asarray(image.data[tuple(index)]) + remaining = [axis for axis, selection in zip(dims, index, strict=True) if not isinstance(selection, int)] + wanted = [axis for axis in ("c", *_SPATIAL) if axis in remaining] + patch = np.transpose(patch, [remaining.index(axis) for axis in wanted]) + if "c" not in remaining: + patch = patch[np.newaxis] + + metadata = { + "axes": dims, + "shape": canonical_shape, + "chunks": list(getattr(image.data, "chunks", []) or []), + "dtype": str(image.data.dtype), + "scale": _ordered(dict(image.scale), dims), + "translation": _ordered(dict(image.translation), dims), + "attributes": _read_konfai_attributes(store_path), + } + return np.asarray(patch), metadata + + +def write_ome_zarr( + store_path: str | Path, + data: np.ndarray, + *, + spacing: Sequence[float] | None = None, + origin: Sequence[float] | None = None, + attributes: dict[str, Any] | None = None, + chunks: Sequence[int] | None = None, +) -> None: + """Write one channel-first KonfAI array as a single-level OME-NGFF store.""" + _require_ngff_zarr() + array_data = np.asarray(data) + if array_data.ndim not in {3, 4}: + raise DatasetManagerError(f"OME-Zarr writing expects a C-Y-X or C-Z-Y-X array, got shape {array_data.shape}.") + + spatial_axes = ["y", "x"] if array_data.ndim == 3 else ["z", "y", "x"] + dims = ["c", *spatial_axes] + dimension = len(spatial_axes) + spacing_xyz = list(spacing if spacing is not None else [1.0] * dimension) + origin_xyz = list(origin if origin is not None else [0.0] * dimension) + if len(spacing_xyz) != dimension or len(origin_xyz) != dimension: + raise DatasetManagerError( + f"OME-Zarr geometry must contain {dimension} spacing and origin values for shape {array_data.shape}." + ) + + coordinate = {"x": (spacing_xyz[0], origin_xyz[0]), "y": (spacing_xyz[1], origin_xyz[1])} + if dimension == 3: + coordinate["z"] = (spacing_xyz[2], origin_xyz[2]) + scale = {"c": 1.0, **{axis: float(coordinate[axis][0]) for axis in spatial_axes}} + translation = {"c": 0.0, **{axis: float(coordinate[axis][1]) for axis in spatial_axes}} + + image = ngff_zarr.to_ngff_image(array_data, dims=dims, scale=scale, translation=translation) + multiscales = ngff_zarr.to_multiscales( + image, scale_factors=[], chunks=tuple(chunks) if chunks is not None else None + ) + ngff_zarr.to_ngff_zarr(str(store_path), multiscales, overwrite=True) + + if attributes: + group = zarr.open_group(str(store_path), mode="r+") + group.attrs[_KONFAI_ATTR_KEY] = {"attributes": dict(attributes)} + + +def get_ome_zarr_info(store_path: str | Path, level: int = 0) -> dict[str, Any]: + """Return OME-Zarr metadata (raw axis-order shape) without reading pixel data.""" + image = _load_image(store_path, level) + dims = [str(axis).lower() for axis in image.dims] + try: + n_levels = len(ngff_zarr.from_ngff_zarr(str(store_path)).images) + except (OSError, TypeError, ValueError): + n_levels = 1 + return { + "axes": dims, + "shape": list(image.data.shape), + "chunks": list(getattr(image.data, "chunks", []) or []), + "dtype": str(image.data.dtype), + "scale": _ordered(dict(image.scale), dims), + "translation": _ordered(dict(image.translation), dims), + "n_levels": n_levels, + "attributes": _read_konfai_attributes(store_path), + } diff --git a/konfai/utils/utils.py b/konfai/utils/utils.py index de961d1f..22c5b375 100755 --- a/konfai/utils/utils.py +++ b/konfai/utils/utils.py @@ -94,7 +94,7 @@ def get_patch_slices_from_shape( nb_patch_per_dim = [] slices: list[list[slice]] = [] if overlap_tmp is None: - size = [np.ceil(a / b) for a, b in zip(shape, patch_size)] + size = [np.ceil(a / b) for a, b in zip(shape, patch_size, strict=False)] tmp = np.zeros(len(size), dtype=np.int_) for i, s in enumerate(size): if s > 1: @@ -143,6 +143,11 @@ def get_patch_slices_from_shape( "hdr", "img", # Analyze "dcm", # DICOM (si GDCM activé) + "dicom", # DICOM series directory backend + "omezarr", + "ome-zarr", + "ome_zarr", + "zarr", # OME-NGFF directory backend and accepted aliases "tif", "tiff", # TIFF "png", @@ -166,6 +171,20 @@ def is_windows_absolute_path(path: str) -> bool: return bool(_WINDOWS_ABSOLUTE_PATH_RE.match(path)) +def split_format_level(file_format: str) -> tuple[str, int]: + """Split an optional pyramid-level suffix from a format token. + + Used by the OME-Zarr backend to pick a multiscale resolution directly in + the dataset spec, e.g. ``omezarr@2`` selects pyramid level 2 (coarser), + independently of any transform. Returns ``(base_format, level)`` and + defaults to level 0 (full resolution) when no ``@`` suffix is present. + """ + base, separator, level = file_format.rpartition("@") + if separator and level.isdigit(): + return base, int(level) + return file_format, 0 + + def split_path_spec( value: str, *, diff --git a/tests/unit/test_imaging_formats.py b/tests/unit/test_imaging_formats.py new file mode 100644 index 00000000..43c554ff --- /dev/null +++ b/tests/unit/test_imaging_formats.py @@ -0,0 +1,316 @@ +# Copyright (c) 2025 Valentin Boussot +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# SPDX-License-Identifier: Apache-2.0 + +"""Unit tests for konfai/utils/dicom.py and konfai/utils/ome_zarr.py.""" + +from pathlib import Path +from unittest.mock import MagicMock, patch + +import numpy as np +import pytest +from konfai.utils.dataset import Attribute, Dataset +from konfai.utils.errors import DatasetManagerError +from konfai.utils.utils import SUPPORTED_EXTENSIONS, split_path_spec + + +def _image_attributes() -> Attribute: + attributes = Attribute() + attributes["Origin"] = np.asarray([10.0, 20.0, 30.0]) + attributes["Spacing"] = np.asarray([0.5, 1.5, 2.0]) + attributes["Direction"] = np.eye(3, dtype=np.float64).flatten() + return attributes + + +# --------------------------------------------------------------------------- +# DICOM tests (no real DICOM files — uses unittest.mock) +# --------------------------------------------------------------------------- + + +class TestDicomRequirePydicom: + def test_raises_without_pydicom(self) -> None: + from konfai.utils import dicom + + with patch.object(dicom, "_PYDICOM_AVAILABLE", False): + with pytest.raises(DatasetManagerError, match="pydicom is required"): + dicom._require_pydicom() + + def test_passes_with_pydicom(self) -> None: + from konfai.utils import dicom + + with patch.object(dicom, "_PYDICOM_AVAILABLE", True): + dicom._require_pydicom() # must not raise + + +class TestDicomDiscoverSeries: + def test_raises_on_missing_directory(self, tmp_path: Path) -> None: + from konfai.utils import dicom + + with patch.object(dicom, "_PYDICOM_AVAILABLE", True): + with pytest.raises(DatasetManagerError, match="does not exist"): + dicom.discover_series(tmp_path / "nonexistent") + + def test_raises_when_no_dicom_found(self, tmp_path: Path) -> None: + from konfai.utils import dicom + + (tmp_path / "file.txt").write_text("not a dicom") + with patch.object(dicom, "_PYDICOM_AVAILABLE", True): + with patch.object(dicom, "pydicom") as mock_pd: + mock_pd.dcmread.side_effect = Exception("not dicom") + with pytest.raises(DatasetManagerError, match="No DICOM files"): + dicom.discover_series(tmp_path) + + +class TestDicomSlicePosition: + def test_uses_ipp_and_iop(self) -> None: + from konfai.utils import dicom + + ds = MagicMock() + # Row = (1, 0, 0), Col = (0, 1, 0) -> normal = (0, 0, 1) + ds.ImageOrientationPatient = [1, 0, 0, 0, 1, 0] + ds.ImagePositionPatient = [0, 0, 42.5] + assert dicom._slice_position(ds) == pytest.approx(42.5) + + def test_falls_back_to_instance_number(self) -> None: + from konfai.utils import dicom + + ds = MagicMock(spec=[]) + ds.InstanceNumber = 7 + assert dicom._slice_position(ds) == 7.0 + + def test_returns_zero_when_no_tags(self) -> None: + from konfai.utils import dicom + + ds = MagicMock(spec=[]) + assert dicom._slice_position(ds) == 0.0 + + +class TestDicomExtractGeometry: + def _make_ds(self, ipp: list[float]) -> MagicMock: + ds = MagicMock() + ds.ImagePositionPatient = ipp + ds.PixelSpacing = [0.5, 0.5] + ds.ImageOrientationPatient = [1, 0, 0, 0, 1, 0] + ds.SliceThickness = 1.0 + return ds + + def test_extracts_correct_spacing(self) -> None: + from konfai.utils import dicom + + ds0 = self._make_ds([0.0, 0.0, 0.0]) + ds1 = self._make_ds([0.0, 0.0, 3.0]) + _, spacing, _ = dicom.extract_geometry([ds0, ds1]) + assert spacing[0] == pytest.approx(0.5) + assert spacing[1] == pytest.approx(0.5) + assert spacing[2] == pytest.approx(3.0) + + def test_fallback_to_slice_thickness_for_single_slice(self) -> None: + from konfai.utils import dicom + + ds = self._make_ds([0.0, 0.0, 0.0]) + _, spacing, _ = dicom.extract_geometry([ds]) + assert spacing[2] == pytest.approx(1.0) + + def test_converts_pixel_spacing_to_xyz_order(self) -> None: + from konfai.utils import dicom + + ds = self._make_ds([0.0, 0.0, 0.0]) + ds.PixelSpacing = [1.5, 0.5] + _, spacing, _ = dicom.extract_geometry([ds]) + np.testing.assert_allclose(spacing, [0.5, 1.5, 1.0]) + + def test_raises_on_missing_ipp(self) -> None: + from konfai.utils import dicom + + ds = MagicMock(spec=[]) + with pytest.raises(DatasetManagerError, match="ImagePositionPatient"): + dicom.extract_geometry([ds]) + + +class TestDicomReadVolume: + def _make_ds(self, value: float = 0.0) -> MagicMock: + ds = MagicMock() + ds.pixel_array = np.full((4, 4), value, dtype=np.int16) + ds.RescaleSlope = 1.0 + ds.RescaleIntercept = -1000.0 + return ds + + def test_stacks_slices_channel_first(self) -> None: + from konfai.utils import dicom + + datasets = [self._make_ds(0.0), self._make_ds(1.0)] + volume = dicom.read_volume(datasets) + assert volume.shape == (1, 2, 4, 4) + + def test_applies_ct_rescale(self) -> None: + from konfai.utils import dicom + + datasets = [self._make_ds(0.0)] + volume = dicom.read_volume(datasets, apply_rescale=True) + assert volume[0, 0, 0, 0] == pytest.approx(-1000.0) + + def test_skips_rescale_when_disabled(self) -> None: + from konfai.utils import dicom + + datasets = [self._make_ds(500.0)] + volume = dicom.read_volume(datasets, apply_rescale=False) + assert volume[0, 0, 0, 0] == pytest.approx(500.0) + + def test_raises_on_inconsistent_shapes(self) -> None: + from konfai.utils import dicom + + ds0 = MagicMock() + ds0.pixel_array = np.zeros((4, 4), dtype=np.int16) + ds0.RescaleSlope = 1.0 + ds0.RescaleIntercept = 0.0 + ds1 = MagicMock() + ds1.pixel_array = np.zeros((8, 8), dtype=np.int16) + ds1.RescaleSlope = 1.0 + ds1.RescaleIntercept = 0.0 + with pytest.raises(DatasetManagerError, match="Inconsistent slice shape"): + dicom.read_volume([ds0, ds1]) + + +# --------------------------------------------------------------------------- +# OME-Zarr tests (no real Zarr store — uses unittest.mock) +# --------------------------------------------------------------------------- + + +class TestOmeZarrRequireZarr: + def test_raises_without_zarr(self) -> None: + from konfai.utils import ome_zarr + + with patch.object(ome_zarr, "_ZARR_AVAILABLE", False): + with pytest.raises(DatasetManagerError, match="zarr is required"): + ome_zarr._require_zarr() + + +class TestDatasetImagingBackends: + def test_ome_zarr_dataset_round_trip_and_patch_read(self, tmp_path: Path) -> None: + pytest.importorskip("zarr") + volume = np.arange(1 * 3 * 4 * 5, dtype=np.int16).reshape(1, 3, 4, 5) + dataset = Dataset(tmp_path / "OME", "ome-zarr") + + dataset.write("CT", "CASE_001", volume, _image_attributes()) + + assert dataset.file_format == "omezarr" + assert dataset.get_names("CT") == ["CASE_001"] + assert dataset.get_group() == ["CT"] + assert dataset.get_infos("CT", "CASE_001")[0] == [1, 3, 4, 5] + full, attributes = dataset.read_data("CT", "CASE_001") + patch, patch_attributes = dataset.read_data_slice( + "CT", "CASE_001", (slice(None), slice(1, 3), slice(1, 4), slice(2, 5)) + ) + np.testing.assert_array_equal(full, volume) + np.testing.assert_array_equal(patch, volume[:, 1:3, 1:4, 2:5]) + np.testing.assert_allclose(attributes.get_np_array("Spacing"), [0.5, 1.5, 2.0]) + np.testing.assert_allclose(patch_attributes.get_np_array("Origin"), [11.0, 21.5, 32.0]) + + def test_ome_zarr_2d_dataset_round_trip(self, tmp_path: Path) -> None: + pytest.importorskip("zarr") + volume = np.arange(2 * 4 * 5, dtype=np.uint8).reshape(2, 4, 5) + attributes = Attribute() + attributes["Origin"] = np.asarray([10.0, 20.0]) + attributes["Spacing"] = np.asarray([0.5, 1.5]) + attributes["Direction"] = np.eye(2, dtype=np.float64).flatten() + dataset = Dataset(tmp_path / "OME2D", "omezarr") + + dataset.write("RGB", "CASE_001", volume, attributes) + result, result_attributes = dataset.read_data("RGB", "CASE_001") + + np.testing.assert_array_equal(result, volume) + np.testing.assert_allclose(result_attributes.get_np_array("Origin"), [10.0, 20.0]) + + def test_dicom_dataset_round_trip_and_patch_read(self, tmp_path: Path) -> None: + pytest.importorskip("pydicom") + volume = np.arange(1 * 3 * 4 * 5, dtype=np.int16).reshape(1, 3, 4, 5) + dataset = Dataset(tmp_path / "DICOM", "dicom") + + dataset.write("CT", "CASE_001", volume, _image_attributes()) + + assert dataset.get_names("CT") == ["CASE_001"] + assert dataset.get_group() == ["CT"] + assert dataset.get_infos("CT", "CASE_001")[0] == [1, 3, 4, 5] + full, attributes = dataset.read_data("CT", "CASE_001") + patch, patch_attributes = dataset.read_data_slice( + "CT", "CASE_001", (slice(None), slice(1, 3), slice(1, 4), slice(2, 5)) + ) + np.testing.assert_array_equal(full, volume) + np.testing.assert_array_equal(patch, volume[:, 1:3, 1:4, 2:5]) + np.testing.assert_allclose(attributes.get_np_array("Spacing"), [0.5, 1.5, 2.0]) + np.testing.assert_allclose(patch_attributes.get_np_array("Origin"), [11.0, 21.5, 32.0]) + + def test_dicom_slice_read_decodes_only_selected_files( + self, tmp_path: Path, monkeypatch: pytest.MonkeyPatch + ) -> None: + pydicom = pytest.importorskip("pydicom") + volume = np.arange(1 * 4 * 3 * 3, dtype=np.int16).reshape(1, 4, 3, 3) + dataset = Dataset(tmp_path / "DICOM", "dicom") + dataset.write("CT", "CASE_001", volume, _image_attributes()) + decoded_files: list[str] = [] + real_dcmread = pydicom.dcmread + + def tracked_dcmread(*args, **kwargs): + if not kwargs.get("stop_before_pixels", False): + decoded_files.append(str(args[0])) + return real_dcmread(*args, **kwargs) + + monkeypatch.setattr(pydicom, "dcmread", tracked_dcmread) + patch, _ = dataset.read_data_slice("CT", "CASE_001", (slice(None), slice(2, 3), slice(None), slice(None))) + + np.testing.assert_array_equal(patch, volume[:, 2:3]) + assert len(decoded_files) == 1 + + def test_dicom_round_trip_preserves_rotated_direction(self, tmp_path: Path) -> None: + pytest.importorskip("pydicom") + attributes = _image_attributes() + direction = np.asarray([[0.0, -1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) + attributes["Direction"] = direction.flatten() + dataset = Dataset(tmp_path / "DICOM", "dicom") + dataset.write("CT", "CASE_001", np.zeros((1, 2, 3, 3), dtype=np.int16), attributes) + + _, result_attributes = dataset.read_data("CT", "CASE_001") + + np.testing.assert_allclose(result_attributes.get_np_array("Direction"), direction.flatten()) + + @pytest.mark.parametrize("file_format", ["omezarr", "ome-zarr", "ome_zarr", "zarr"]) + def test_ome_zarr_format_aliases(self, tmp_path: Path, file_format: str) -> None: + assert Dataset(tmp_path / file_format, file_format).file_format == "omezarr" + + @pytest.mark.parametrize("file_format", ["dicom", "omezarr", "ome-zarr", "ome_zarr", "zarr"]) + def test_data_manager_path_parser_accepts_imaging_backend(self, file_format: str) -> None: + assert file_format in SUPPORTED_EXTENSIONS + assert split_path_spec( + f"./Dataset:a:{file_format}", + allowed_flags={"a", "i"}, + supported_extensions=SUPPORTED_EXTENSIONS, + ) == ("./Dataset", "a", file_format) + + @pytest.mark.parametrize("file_format", ["dicom", "omezarr"]) + def test_data_prediction_resolves_imaging_dataset_source(self, tmp_path: Path, file_format: str) -> None: + from konfai.data.data_manager import DataPrediction, Group, GroupTransform + + volume = np.arange(1 * 2 * 3 * 3, dtype=np.int16).reshape(1, 2, 3, 3) + root = tmp_path / file_format + Dataset(root, file_format).write("CT", "CASE_001", volume, _image_attributes()) + prediction_data = DataPrediction( + augmentations=None, + dataset_filenames=[f"{root}:a:{file_format}"], + groups_src={"CT": Group(groups_dest={"CT": GroupTransform(transforms=None, patch_transforms=None)})}, + ) + + sources = prediction_data._resolve_dataset_sources() + + assert sources == {"CT": [(str(root), True)]} diff --git a/tests/unit/test_imaging_roundtrip.py b/tests/unit/test_imaging_roundtrip.py new file mode 100644 index 00000000..42892187 --- /dev/null +++ b/tests/unit/test_imaging_roundtrip.py @@ -0,0 +1,208 @@ +# Copyright (c) 2025 Valentin Boussot +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# SPDX-License-Identifier: Apache-2.0 + +"""End-to-end imaging round-trip tests for the DICOM and OME-Zarr backends. + +These tests exercise the real optional backends (pydicom / zarr / ngff-zarr / +SimpleITK) and skip gracefully when they are not installed. They complement the +mostly-mocked ``test_imaging_formats.py`` with full read/write round-trips, +cross-validation against SimpleITK's own DICOM reader, and ngff-zarr interop. +""" + +from pathlib import Path + +import numpy as np +import pytest +from konfai.utils.errors import DatasetManagerError + +# --------------------------------------------------------------------------- +# DICOM round-trips +# --------------------------------------------------------------------------- + + +def test_dicom_roundtrip_geometry_matches_simpleitk(tmp_path: Path) -> None: + """KonfAI's DICOM reader must agree with SimpleITK's GDCM reader on geometry.""" + pytest.importorskip("pydicom") + sitk = pytest.importorskip("SimpleITK") + from konfai.utils import dicom + + root = tmp_path / "series" + vol = (np.arange(1 * 5 * 6 * 7).reshape(1, 5, 6, 7) % 100).astype(np.float32) + origin, spacing = (3.0, 4.0, 5.0), (0.8, 0.9, 2.0) + dicom.write_dicom_series(root, vol, origin=origin, spacing=spacing, direction=np.eye(3).flatten()) + + kvol, kog, ksp, kdir = dicom.read_dicom_series(root) + + reader = sitk.ImageSeriesReader() + ids = reader.GetGDCMSeriesIDs(str(root)) + reader.SetFileNames(reader.GetGDCMSeriesFileNames(str(root), ids[0])) + itk = reader.Execute() + + assert np.allclose(kog, itk.GetOrigin(), atol=1e-3) + assert np.allclose(ksp, itk.GetSpacing(), atol=1e-3) + assert np.allclose(kdir, itk.GetDirection(), atol=1e-3) + assert np.allclose(kvol[0], sitk.GetArrayFromImage(itk), atol=1.0) + + +def test_dicom_left_handed_direction_normalizes_like_simpleitk(tmp_path: Path) -> None: + """A feet-first (z-down) direction round-trips to a right-handed frame, matching SimpleITK. + + DICOM cannot store an arbitrary left-handed direction: the slice axis is + derived from positions. KonfAI must reproduce SimpleITK's normalization + (flipped array + adjusted origin/direction describing the SAME physical volume). + """ + pytest.importorskip("pydicom") + sitk = pytest.importorskip("SimpleITK") + from konfai.utils import dicom + + root = tmp_path / "series" + # distinct value per slice so ordering is observable + vol = np.stack([np.full((4, 5), k, np.float32) for k in range(6)])[np.newaxis] + dicom.write_dicom_series( + root, vol, origin=(0.0, 0.0, 30.0), spacing=(1.0, 1.0, 2.0), + direction=np.array([1, 0, 0, 0, 1, 0, 0, 0, -1], float), + ) + kvol, kog, ksp, kdir = dicom.read_dicom_series(root) + + reader = sitk.ImageSeriesReader() + ids = reader.GetGDCMSeriesIDs(str(root)) + reader.SetFileNames(reader.GetGDCMSeriesFileNames(str(root), ids[0])) + itk = reader.Execute() + + assert np.allclose(kog, itk.GetOrigin(), atol=1e-3) + assert np.allclose(kdir, itk.GetDirection(), atol=1e-3) + assert np.allclose(kvol[0], sitk.GetArrayFromImage(itk), atol=1.0) + # right-handed normalization: z-cosine flips from -1 to +1 + assert np.allclose(kdir.reshape(3, 3)[:, 2], [0, 0, 1], atol=1e-6) + + +def test_dicom_slice_arity_mismatch_raises_dataset_manager_error(tmp_path: Path) -> None: + pytest.importorskip("pydicom") + from konfai.utils import dicom + + root = tmp_path / "series" + vol = (np.arange(1 * 4 * 5 * 6).reshape(1, 4, 5, 6)).astype(np.float32) + dicom.write_dicom_series(root, vol, origin=(0.0, 0.0, 0.0), spacing=(1.0, 1.0, 1.0)) + with pytest.raises(DatasetManagerError): + dicom.read_dicom_series_slice(root, (slice(None), slice(0, 2))) # 2 slices, expected 4 + + +# --------------------------------------------------------------------------- +# OME-Zarr round-trips & interop +# --------------------------------------------------------------------------- + + +def test_ome_zarr_output_is_readable_by_ngff_zarr(tmp_path: Path) -> None: + """KonfAI's hand-written OME-Zarr must stay interoperable with ngff-zarr. + + Guards against the hand-rolled NGFF metadata drifting from the standard the + ``ngff-zarr`` dependency implements. + """ + pytest.importorskip("zarr") + ngff_zarr = pytest.importorskip("ngff_zarr") + from konfai.utils import ome_zarr + + store = tmp_path / "img.ome.zarr" + data = (np.arange(1 * 6 * 8 * 10).reshape(1, 6, 8, 10) % 50).astype(np.float32) + ome_zarr.write_ome_zarr(store, data, spacing=(0.5, 0.6, 2.0), origin=(1.0, 2.0, 3.0)) + + mz = ngff_zarr.from_ngff_zarr(str(store)) + img0 = mz.images[0] + assert tuple(img0.dims) == ("c", "z", "y", "x") + assert img0.scale == {"c": 1.0, "z": 2.0, "y": 0.6, "x": 0.5} + assert img0.translation == {"c": 0.0, "z": 3.0, "y": 2.0, "x": 1.0} + assert np.array_equal(np.asarray(img0.data), data) + + +def test_ome_zarr_backend_preserves_non_identity_direction(tmp_path: Path) -> None: + pytest.importorskip("zarr") + sitk = pytest.importorskip("SimpleITK") + from konfai.utils.dataset import Dataset + + root = tmp_path / "ds" + root.mkdir() + arr = (np.arange(8 * 12 * 10).reshape(8, 12, 10) % 97).astype(np.float32) + img = sitk.GetImageFromArray(arr) + img.SetSpacing((0.7, 0.8, 2.5)) + img.SetOrigin((10.0, -5.0, 3.0)) + direction = (0.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0) # 90° in-plane + img.SetDirection(direction) + + Dataset.OmeZarrFile(str(root), read=False).data_to_file("CASE0", img, None) + rdata, rattr = Dataset.OmeZarrFile(str(root), read=True).file_to_data("CT", "CASE0") + + assert np.allclose(rdata[0], arr) + assert np.allclose(rattr.get_np_array("Spacing"), [0.7, 0.8, 2.5]) + assert np.allclose(rattr.get_np_array("Origin"), [10.0, -5.0, 3.0]) + assert np.allclose(rattr.get_np_array("Direction"), direction) + + +def test_ome_zarr_slice_arity_mismatch_raises_dataset_manager_error(tmp_path: Path) -> None: + pytest.importorskip("zarr") + from konfai.utils import ome_zarr + + store = tmp_path / "img.ome.zarr" + data = (np.arange(1 * 4 * 5 * 6).reshape(1, 4, 5, 6)).astype(np.float32) + ome_zarr.write_ome_zarr(store, data, spacing=(1.0, 1.0, 1.0), origin=(0.0, 0.0, 0.0)) + with pytest.raises(DatasetManagerError): + ome_zarr.read_ome_zarr_data_slice(store, (slice(None), slice(0, 2))) # 2 slices, expected 4 + + +# --------------------------------------------------------------------------- +# OME-Zarr resolution (pyramid level) selection — `omezarr@` +# --------------------------------------------------------------------------- + + +def test_split_format_level_parses_pyramid_suffix() -> None: + from konfai.utils.utils import split_format_level + + assert split_format_level("omezarr") == ("omezarr", 0) + assert split_format_level("omezarr@2") == ("omezarr", 2) + assert split_format_level("ome-zarr@1") == ("ome-zarr", 1) + assert split_format_level("mha") == ("mha", 0) # unaffected + assert split_format_level("C:/data@x") == ("C:/data@x", 0) # non-numeric ignored + + +def test_dataset_parses_omezarr_level_field(tmp_path: Path) -> None: + from konfai.utils.dataset import Dataset + + assert Dataset(tmp_path / "a", "omezarr@2").level == 2 + assert Dataset(tmp_path / "a", "omezarr").level == 0 + coarse = Dataset(tmp_path / "a", "ome-zarr@3") + assert coarse.file_format == "omezarr" and coarse.level == 3 + + +def test_ome_zarr_level_reads_coarser_resolution(tmp_path: Path) -> None: + ngff_zarr = pytest.importorskip("ngff_zarr") + pytest.importorskip("zarr") + from konfai.utils.dataset import Dataset + + root = tmp_path / "ds" + root.mkdir() + data = (np.arange(1 * 16 * 32 * 32).reshape(1, 16, 32, 32) % 50).astype(np.float32) + image = ngff_zarr.to_ngff_image( + data, dims=["c", "z", "y", "x"], + scale={"c": 1.0, "z": 2.0, "y": 0.5, "x": 0.5}, translation={"c": 0.0, "z": 0.0, "y": 0.0, "x": 0.0}, + ) + ngff_zarr.to_ngff_zarr(str(root / "CASE0.ome.zarr"), ngff_zarr.to_multiscales(image, scale_factors=[2]), overwrite=True) + + full, attr0 = Dataset.OmeZarrFile(str(root), read=True, level=0).file_to_data("g", "CASE0") + coarse, attr1 = Dataset.OmeZarrFile(str(root), read=True, level=1).file_to_data("g", "CASE0") + + assert list(full.shape) == [1, 16, 32, 32] + assert list(coarse.shape) == [1, 8, 16, 16] # level 1 is downsampled x2 + # spacing doubles at the coarser level + np.testing.assert_allclose(attr1.get_np_array("Spacing"), 2.0 * attr0.get_np_array("Spacing"))