diff --git a/py/ngff_zarr/methods/_support.py b/py/ngff_zarr/methods/_support.py index 0e369246..d6832e17 100644 --- a/py/ngff_zarr/methods/_support.py +++ b/py/ngff_zarr/methods/_support.py @@ -9,6 +9,36 @@ _spatial_dims = {"x", "y", "z"} +#: The OME-Zarr specification axis order: time, then channel, then space. +_CANONICAL_AXIS_ORDER = ("t", "c", "z", "y", "x") + + +def _canonical_axis_order(ngff_image: NgffImage) -> NgffImage: + """Return the image with dims in the spec axis order (t, c, z, y, x). + + OME-Zarr requires axes ordered by type: time, then channel, then space. + Conversion sources produce channel-last layouts -- the TIFF ``S`` (sample) + axis, ITK/ITKWasm component images, and the default dims inference for 4-D + and 5-D arrays all yield ``(..., c)``. The multiscale pipeline normalizes + them through this helper so generated metadata is spec-ordered; the data + transpose is lazy. + + An image with dims outside the canonical set is returned unchanged: axis + models not expressible before RFC-3 carry no spec ordering to normalize to. + """ + dims = list(ngff_image.dims) + new_dims = [dim for dim in _CANONICAL_AXIS_ORDER if dim in dims] + if len(new_dims) != len(dims) or tuple(new_dims) == tuple(dims): + return ngff_image + + new_order = [dims.index(dim) for dim in new_dims] + + result = copy.copy(ngff_image) + result.data = ngff_image.data.transpose(new_order) + result.dims = tuple(new_dims) + + return result + def _spatial_dims_last(ngff_image: NgffImage) -> NgffImage: dims = list(ngff_image.dims) diff --git a/py/ngff_zarr/tiff_to_ngff_image.py b/py/ngff_zarr/tiff_to_ngff_image.py index 0f4b3a7d..1a4ec8cf 100644 --- a/py/ngff_zarr/tiff_to_ngff_image.py +++ b/py/ngff_zarr/tiff_to_ngff_image.py @@ -24,6 +24,7 @@ import dask.array as da import zarr +from .methods._support import _canonical_axis_order from .ngff_image import NgffImage from .to_ngff_image import to_ngff_image @@ -812,8 +813,10 @@ def _build_multiscales_from_pyramid( base_scale = ngff_image_0.scale base_translation = ngff_image_0.translation - # Build NgffImages for all levels - images = [ngff_image_0] + # Build NgffImages for all levels. Levels are computed in the TIFF's own + # axis order and normalized to the spec order (t, c, z, y, x) on append, + # so the emitted metadata orders time, then channel, then space. + images = [_canonical_axis_order(ngff_image_0)] for path in paths[1:]: arr = root[path] level_scale: dict[str, float] = {} @@ -875,11 +878,11 @@ def _build_multiscales_from_pyramid( channel_names=ngff_image_0.channel_names, channel_colors=ome_channel_colors, ) - images.append(level_image) + images.append(_canonical_axis_order(level_image)) # Build Metadata (axes, datasets, coordinate transforms) axes = [] - for dim in ngff_image_0.dims: + for dim in images[0].dims: unit = None if ngff_image_0.axes_units and dim in ngff_image_0.axes_units: unit = ngff_image_0.axes_units[dim] diff --git a/py/ngff_zarr/to_multiscales.py b/py/ngff_zarr/to_multiscales.py index 3941812e..7af29947 100644 --- a/py/ngff_zarr/to_multiscales.py +++ b/py/ngff_zarr/to_multiscales.py @@ -35,7 +35,7 @@ _downsample_itkwasm, ) from .methods._metadata import get_method_metadata -from .methods._support import _spatial_dims +from .methods._support import _canonical_axis_order, _spatial_dims from .multiscales import NgffMultiscales from .ngff_image import NgffImage from .rfc4 import AnatomicalOrientation, orientation_from_name @@ -496,7 +496,10 @@ def to_multiscales( output automatically whenever it is present; no separate opt-in is required. :type orientation: str or mapping of str to AnatomicalOrientation, optional - :return: NgffImage for each resolution and NGFF multiscales metadata + :return: NgffImage for each resolution and NGFF multiscales metadata. + Axes are normalized to the OME-Zarr specification order -- time, then + channel, then space (t, c, z, y, x) -- with a lazy transpose when the + input image orders them differently. :rtype : NgffMultiscales """ ngff_image = data if isinstance(data, NgffImage) else to_ngff_image(data) @@ -566,13 +569,24 @@ def to_multiscales( if "t" in ngff_image.dims: default_chunks["t"] = 1 - da_out_chunks = tuple(out_chunks[d] for d in ngff_image.dims) if not isinstance(ngff_image.data, DaskArray): if isinstance(ngff_image.data, (ZarrArray, str, MutableMapping)): ngff_image.data = dask.array.from_zarr(ngff_image.data) else: ngff_image.data = dask.array.from_array(ngff_image.data) + # OME-Zarr orders axes time, then channel, then space. Channel-last input + # (the TIFF S axis, ITK component images, the 4-D/5-D default dims) is + # normalized with a lazy transpose so the generated metadata and every + # scale are spec-ordered. + ngff_image = _canonical_axis_order(ngff_image) + # Re-key the dim-keyed chunk mappings to follow the (possibly reordered) + # dims; _ngff_image_scale_factors asserts this ordering. + out_chunks = {dim: out_chunks[dim] for dim in ngff_image.dims} + default_chunks = {dim: default_chunks[dim] for dim in ngff_image.dims} + + da_out_chunks = tuple(out_chunks[d] for d in ngff_image.dims) + if isinstance(scale_factors, int): scale_factors = _ngff_image_scale_factors(ngff_image, scale_factors, out_chunks) diff --git a/py/test/test_bin_shrink_map_blocks_fast_path.py b/py/test/test_bin_shrink_map_blocks_fast_path.py index d236a284..263e07e8 100644 --- a/py/test/test_bin_shrink_map_blocks_fast_path.py +++ b/py/test/test_bin_shrink_map_blocks_fast_path.py @@ -142,7 +142,9 @@ def test_zyxc_vector_mode(self): image, scale_factors=[2], method=Methods.ITKWASM_BIN_SHRINK ) assert len(multiscales.images) == 2 - assert multiscales.images[1].data.shape == (8, 16, 16, 3) + # Channels lead in the spec axis order (c, z, y, x) + assert multiscales.images[1].dims == ("c", "z", "y", "x") + assert multiscales.images[1].data.shape == (3, 8, 16, 16) store = MemoryStore() to_ngff_zarr(store, multiscales) @@ -157,12 +159,13 @@ def test_zyxc_many_channels(self): image, scale_factors=[2], method=Methods.ITKWASM_BIN_SHRINK ) assert len(multiscales.images) == 2 - # Verify channels are preserved - assert multiscales.images[1].data.shape[-1] == 12 + # Channels lead in the spec axis order (c, z, y, x) and are preserved + assert multiscales.images[1].dims == ("c", "z", "y", "x") + assert multiscales.images[1].data.shape[0] == 12 # c # Verify spatial dimensions are downsampled - assert multiscales.images[1].data.shape[0] == 4 # z - assert multiscales.images[1].data.shape[1] == 16 # y - assert multiscales.images[1].data.shape[2] == 16 # x + assert multiscales.images[1].data.shape[1] == 4 # z + assert multiscales.images[1].data.shape[2] == 16 # y + assert multiscales.images[1].data.shape[3] == 16 # x def test_dict_scale_factors(self): """Dict-based scale factors on the fast path.""" diff --git a/py/test/test_canonical_axis_order.py b/py/test/test_canonical_axis_order.py new file mode 100644 index 00000000..22e4ee2d --- /dev/null +++ b/py/test/test_canonical_axis_order.py @@ -0,0 +1,77 @@ +# SPDX-FileCopyrightText: Copyright (c) Fideus Labs LLC +# SPDX-License-Identifier: MIT +"""The multiscale pipeline normalizes axes to the spec order (t, c, z, y, x).""" + +import dask.array +import numpy as np +from ngff_zarr import to_multiscales, to_ngff_image +from ngff_zarr.methods._support import _canonical_axis_order +from ngff_zarr.ngff_image import NgffImage + + +def test_channel_last_input_is_normalized(): + """A (y, x, c) image comes out as (c, y, x) with matching metadata.""" + data = np.random.randint(0, 255, (32, 64, 3), dtype=np.uint8) + image = to_ngff_image(data, dims=("y", "x", "c")) + + multiscales = to_multiscales(image, scale_factors=[]) + + assert multiscales.images[0].dims == ("c", "y", "x") + assert multiscales.images[0].data.shape == (3, 32, 64) + axes = multiscales.metadata.coordinateSystems[0].axes + assert [ax.name for ax in axes] == ["c", "y", "x"] + + # The transpose is a relabeling, not a data change. + original = np.moveaxis(data, -1, 0) + np.testing.assert_array_equal(np.asarray(multiscales.images[0].data), original) + + +def test_default_dims_inference_is_normalized(): + """The 4-D default dims (z, y, x, c) are normalized to (c, z, y, x).""" + data = np.zeros((8, 16, 16, 2), dtype=np.uint8) + + multiscales = to_multiscales(data, scale_factors=[]) + + assert multiscales.images[0].dims == ("c", "z", "y", "x") + assert multiscales.images[0].data.shape == (2, 8, 16, 16) + + +def test_canonical_input_is_unchanged(): + """An image already in spec order passes through untouched.""" + data = dask.array.zeros((2, 8, 16, 16), dtype=np.uint8) + image = to_ngff_image(data, dims=("c", "z", "y", "x")) + + multiscales = to_multiscales(image, scale_factors=[]) + + assert multiscales.images[0].dims == ("c", "z", "y", "x") + axes = multiscales.metadata.coordinateSystems[0].axes + assert [ax.name for ax in axes] == ["c", "z", "y", "x"] + + +def test_scale_and_translation_follow_the_axes(): + """Dim-keyed metadata keeps its values across the reorder.""" + data = np.zeros((32, 64, 3), dtype=np.uint8) + image = to_ngff_image( + data, dims=("y", "x", "c"), scale={"y": 2.0, "x": 4.0}, translation={"y": 1.0} + ) + + multiscales = to_multiscales(image, scale_factors=[]) + + out = multiscales.images[0] + assert out.scale["y"] == 2.0 + assert out.scale["x"] == 4.0 + assert out.translation["y"] == 1.0 + sequence = multiscales.metadata.datasets[0].coordinateTransformations[0] + assert sequence.transformations[0].scale == [1.0, 2.0, 4.0] # (c, y, x) order + + +def test_unknown_dims_are_left_untouched(): + """Axis models outside {t, c, z, y, x} have no spec order to normalize to.""" + data = dask.array.zeros((4, 5, 6), dtype=np.uint8) + image = NgffImage( + data=data, dims=("b", "a", "x"), scale={"x": 1.0}, translation={"x": 0.0} + ) + + result = _canonical_axis_order(image) + + assert result is image diff --git a/py/test/test_issue_436.py b/py/test/test_issue_436.py index f39387e4..b2e9c0fe 100644 --- a/py/test/test_issue_436.py +++ b/py/test/test_issue_436.py @@ -143,13 +143,13 @@ def test_multiseries_tiff_to_ozx_creates_requested_file(): # The auxiliary series is written alongside, preserving the .ozx format. assert (tmp / "out_Thumbnail.ozx").exists() - # The primary .ozx must hold the full-resolution Baseline pyramid. + # The primary .ozx must hold the full-resolution Baseline pyramid, + # normalized from the TIFF's YXS order to the spec axis order (c, y, x). store = zarr.storage.ZipStore(str(out), mode="r") root = zarr.open_group(store, mode="r") scale0 = root["scale0"] arr = scale0[next(iter(scale0.keys()))] - assert arr.shape[0] == base.shape[0] - assert arr.shape[1] == base.shape[1] + assert arr.shape == (base.shape[2], base.shape[0], base.shape[1]) # Real pixel data, not an all-zero / blank array. assert int(np.asarray(arr[:]).max()) > 0 diff --git a/py/test/test_large_image_chunking.py b/py/test/test_large_image_chunking.py index b486e428..728784d4 100644 --- a/py/test/test_large_image_chunking.py +++ b/py/test/test_large_image_chunking.py @@ -71,10 +71,12 @@ def test_3d_rgb_input_chunks_preserved(self): image = to_ngff_image(arr, dims=("z", "y", "x", "c")) multiscales = to_multiscales(image, scale_factors=[]) - out_data = multiscales.images[0].data - c_idx = list(image.dims).index("c") - y_idx = list(image.dims).index("y") - x_idx = list(image.dims).index("x") + # Output dims are normalized to the spec axis order (c, z, y, x) + out_image = multiscales.images[0] + out_data = out_image.data + c_idx = list(out_image.dims).index("c") + y_idx = list(out_image.dims).index("y") + x_idx = list(out_image.dims).index("x") # Channels should stay together assert out_data.chunksize[c_idx] == 3 # Spatial should stay >= 512 @@ -164,9 +166,10 @@ def test_rgb_channels_kept_together(self): image = to_ngff_image(arr, dims=("y", "x", "c")) multiscales = to_multiscales(image, scale_factors=[]) - out_data = multiscales.images[0].data - c_idx = list(image.dims).index("c") - assert out_data.chunksize[c_idx] == 3 + # Output dims are normalized to the spec axis order (c, y, x) + out_image = multiscales.images[0] + c_idx = list(out_image.dims).index("c") + assert out_image.data.chunksize[c_idx] == 3 def test_multichannel_kept_together(self): """Multi-channel data (e.g. 16 channels) should keep channel chunks.""" @@ -174,9 +177,10 @@ def test_multichannel_kept_together(self): image = to_ngff_image(arr, dims=("y", "x", "c")) multiscales = to_multiscales(image, scale_factors=[]) - out_data = multiscales.images[0].data - c_idx = list(image.dims).index("c") - assert out_data.chunksize[c_idx] == 16 + # Output dims are normalized to the spec axis order (c, y, x) + out_image = multiscales.images[0] + c_idx = list(out_image.dims).index("c") + assert out_image.data.chunksize[c_idx] == 16 class TestCaching2DStrips: diff --git a/py/test/test_tiff_pyramid_rgb.py b/py/test/test_tiff_pyramid_rgb.py index 11952973..fb0221f1 100644 --- a/py/test/test_tiff_pyramid_rgb.py +++ b/py/test/test_tiff_pyramid_rgb.py @@ -77,17 +77,17 @@ def test_pyramidal_rgb_tiff_channel_consistency(): f"Expected 3 pyramid levels, got {len(result.images)}" ) - # Verify all levels have consistent dimensions + # Verify all levels are normalized to the spec axis order (c, y, x) for i, img in enumerate(result.images): assert img.dims == ( + "c", "y", "x", - "c", - ), f"Level {i} should have dims ('y', 'x', 'c'), got {img.dims}" + ), f"Level {i} should have dims ('c', 'y', 'x'), got {img.dims}" # Check shape consistency expected_size = 256 // (2**i) # 256, 128, 64 - expected_shape = (expected_size, expected_size, 3) + expected_shape = (3, expected_size, expected_size) assert img.data.shape == expected_shape, ( f"Level {i} should have shape {expected_shape}, got {img.data.shape}" ) @@ -95,12 +95,12 @@ def test_pyramidal_rgb_tiff_channel_consistency(): # Verify data can be computed without errors for i, img in enumerate(result.images): # Compute a small slice to ensure lazy loading works - slice_data = img.data[:4, :4, :].compute() + slice_data = img.data[:, :4, :4].compute() assert slice_data.shape == ( + 3, 4, 4, - 3, - ), f"Level {i} slice should have shape (4, 4, 3), got {slice_data.shape}" + ), f"Level {i} slice should have shape (3, 4, 4), got {slice_data.shape}" assert slice_data.dtype == np.uint8, ( f"Level {i} should preserve uint8 dtype, got {slice_data.dtype}" ) diff --git a/py/test/test_to_ngff_zarr_itkwasm.py b/py/test/test_to_ngff_zarr_itkwasm.py index 59d81ed6..f9889390 100644 --- a/py/test/test_to_ngff_zarr_itkwasm.py +++ b/py/test/test_to_ngff_zarr_itkwasm.py @@ -46,10 +46,10 @@ def test_downsample_zycx(): multiscales = to_multiscales(image, scale_factors=[2, 4], chunks=32) store = MemoryStore() to_ngff_zarr(store, multiscales) - assert multiscales.images[0].dims[0] == "z" - assert multiscales.images[0].dims[2] == "c" - assert multiscales.images[1].data.shape[0] == 16 - assert multiscales.images[1].data.shape[2] == 2 + assert multiscales.images[0].dims[0] == "c" + assert multiscales.images[0].dims[1] == "z" + assert multiscales.images[1].data.shape[0] == 2 + assert multiscales.images[1].data.shape[1] == 16 def test_downsample_cxyz(): @@ -85,11 +85,11 @@ def test_downsample_tzycx(): store = MemoryStore() to_ngff_zarr(store, multiscales) assert multiscales.images[0].dims[0] == "t" - assert multiscales.images[0].dims[1] == "z" - assert multiscales.images[0].dims[3] == "c" + assert multiscales.images[0].dims[1] == "c" + assert multiscales.images[0].dims[2] == "z" assert multiscales.images[1].data.shape[0] == 2 - assert multiscales.images[1].data.shape[1] == 16 - assert multiscales.images[1].data.shape[3] == 2 + assert multiscales.images[1].data.shape[1] == 2 + assert multiscales.images[1].data.shape[2] == 16 def test_downsample_tcxyz(): @@ -280,12 +280,13 @@ def test_itkwasm_gaussian_many_channels(): ) assert len(multiscales.images) == 2 - # Verify channels are preserved - assert multiscales.images[1].data.shape[-1] == 16 + # Channels lead in the spec axis order (c, z, y, x) and are preserved + assert multiscales.images[1].dims == ("c", "z", "y", "x") + assert multiscales.images[1].data.shape[0] == 16 # c # Verify spatial dimensions are downsampled - assert multiscales.images[1].data.shape[0] == 4 # z - assert multiscales.images[1].data.shape[1] == 16 # y - assert multiscales.images[1].data.shape[2] == 16 # x + assert multiscales.images[1].data.shape[1] == 4 # z + assert multiscales.images[1].data.shape[2] == 16 # y + assert multiscales.images[1].data.shape[3] == 16 # x def test_itkwasm_gaussian_few_channels(): @@ -299,12 +300,13 @@ def test_itkwasm_gaussian_few_channels(): ) assert len(multiscales.images) == 2 - # Verify channels are preserved - assert multiscales.images[1].data.shape[-1] == 3 + # Channels lead in the spec axis order (c, z, y, x) and are preserved + assert multiscales.images[1].dims == ("c", "z", "y", "x") + assert multiscales.images[1].data.shape[0] == 3 # c # Verify spatial dimensions are downsampled - assert multiscales.images[1].data.shape[0] == 4 # z - assert multiscales.images[1].data.shape[1] == 16 # y - assert multiscales.images[1].data.shape[2] == 16 # x + assert multiscales.images[1].data.shape[1] == 4 # z + assert multiscales.images[1].data.shape[2] == 16 # y + assert multiscales.images[1].data.shape[3] == 16 # x def test_itkwasm_bin_shrink_many_channels(): @@ -322,9 +324,10 @@ def test_itkwasm_bin_shrink_many_channels(): ) assert len(multiscales.images) == 2 - # Verify channels are preserved - assert multiscales.images[1].data.shape[-1] == 12 + # Channels lead in the spec axis order (c, z, y, x) and are preserved + assert multiscales.images[1].dims == ("c", "z", "y", "x") + assert multiscales.images[1].data.shape[0] == 12 # c # Verify spatial dimensions are downsampled - assert multiscales.images[1].data.shape[0] == 4 # z - assert multiscales.images[1].data.shape[1] == 16 # y - assert multiscales.images[1].data.shape[2] == 16 # x + assert multiscales.images[1].data.shape[1] == 4 # z + assert multiscales.images[1].data.shape[2] == 16 # y + assert multiscales.images[1].data.shape[3] == 16 # x