fix(py): normalize generated axes to the spec order (t, c, z, y, x)

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)

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]

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)
 

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."""

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

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
 

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,19 +166,21 @@ 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."""
         arr = dask.array.zeros((256, 256, 16), chunks=(128, 128, 16), dtype=np.uint8)
         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:

py/test/test_tiff_pyramid_rgb.py

@@ -77,30 +77,30 @@ 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}"
         )
 
     # 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}"
         )