Merge pull request #516 from jo-mueller/clean-scalers

deprecate scaler class

Author: will-moore

docs/requirements.txt

@@ -7,3 +7,4 @@ numpy
 rangehttpserver
 scipy
 scikit-image
+Deprecated

docs/source/cli.rst

@@ -56,6 +56,36 @@ Specify a different output directory::
 
     ome_zarr download https://uk1s3.embassy.ebi.ac.uk/idr/zarr/v0.5/idr0062A/6001240_labels.zarr --output image_dir
 
+scale
+=====
+
+Use the `ome_zarr scale` command to generate a multiscale pyramid from a Zarr array.
+This creates downsampled resolutions of the input image using the specified downsampling method.
+
+Basic usage::
+
+    ome_zarr scale input.zarr output.zarr zyx
+
+This reads the input Zarr array with dimensions 'zyx' and writes a multiscale pyramid to the output directory.
+
+Options:
+
+- `--downscale`: The downsampling factor (default: 2). For example, `--downscale 2` creates scale factors of 2, 4, 8, etc.
+- `--max_layer`: Maximum number of pyramid levels to generate (default: 4)
+- `--method`: Downsampling method to use (default: "resize"). Options are:
+    - `nearest`: Nearest neighbor
+    - `resize`: Resize with anti-aliasing (default)
+    - `laplacian`: Laplacian pyramid
+    - `local_mean`: Local mean
+    - `zoom`: Scipy zoom
+- `--copy-metadata`: If specified, copies the input array metadata to the output group
+
+Example with custom options::
+
+    ome_zarr scale input.zarr output.zarr tczyx --downscale 3 --max_layer 5 --method nearest --copy-metadata
+
+
+
 create
 ======
 

docs/source/python.rst

@@ -38,6 +38,47 @@ This image can be viewed in `napari` using the
 
     $ napari test_ngff_image.zarr
 
+Building a pyramid
+------------------
+
+Multi-resolution pyramids are an integral part of ome-zarr image data
+and enable fast rendering of large images.
+The entrypoints to writing ome-zarr images in ome-zarr-py (`write_image` and `write_labels`)
+build these pyramids under the hood as delayed dask arrays based on the settings for the scaling functions and scale factors.
+
+The scale factors can be passed as a list of integers or a list of dicts::
+
+    from ome_zarr.writer import write_image
+
+    scale_factors = [2, 4, 8]
+    write_image(
+        your_data
+        path,
+        axes="zyx",
+        scale_factors=scale_factors,
+        )
+
+In this example, the downsampling will be applied in all spatial dimensions *except the z dimension*, which will be left at a scale factor of 1.
+To apply equal or custom downsampling factors along all spatial dimensions, pass the scale factors as a list of dicts, e.g.::
+
+    from ome_zarr.writer import write_image
+
+    scale_factors = [
+        {"z": 2, "y": 2, "x": 2},
+        {"z": 4, "y": 4, "x": 4},
+        {"z": 8, "y": 8, "x": 8}
+        ]
+    write_image(
+        your_data
+        path,
+        axes="zyx",
+        scale_factors=scale_factors,
+        )
+
+If you have already built a pyramid representation by other means,
+you can pass it directly to the :py:func:`ome_zarr.writer.write_multiscale` or use :py:func:`ome_zarr.writer.write_multiscale_labels`,
+which do not perform any down-sampling but just write the passed pyramid to disk with the correct metadata.
+
 Rendering settings
 ------------------
 Rendering settings can be added to an existing zarr group::

mypy.ini

@@ -2,5 +2,8 @@
 disallow_untyped_defs = False
 ignore_missing_imports = True
 
+[mypy-deprecated.*]
+ignore_missing_imports = True
+
 [mypy-ome_zarr]
 disallow_untyped_defs = True

ome_zarr/cli.py

@@ -7,7 +7,6 @@
 from .csv import csv_to_zarr
 from .data import astronaut, coins, create_zarr
 from .format import CurrentFormat, Format, format_from_version
-from .scale import Scaler
 from .utils import download as zarr_download
 from .utils import finder as bff_finder
 from .utils import info as zarr_info
@@ -72,16 +71,30 @@ def create(args: argparse.Namespace) -> None:
 
 
 def scale(args: argparse.Namespace) -> None:
-    """Wrap the :func:`~ome_zarr.scale.Scaler.scale` method."""
-    scaler = Scaler(
-        copy_metadata=args.copy_metadata,
-        downscale=args.downscale,
-        in_place=args.in_place,
-        labeled=args.labeled,
-        max_layer=args.max_layer,
+    import dask.array as da
+    import zarr
+
+    from .writer import write_image
+
+    """Wrap the :func:`~ome_zarr.scale._build_pyramid` method."""
+    base = zarr.open_array(args.input_array, mode="r")
+    scale_factors = tuple(args.downscale**i for i in range(1, args.max_layer + 1))
+
+    data = da.from_zarr(args.input_array)
+
+    write_image(
+        data,
+        args.output_directory,
+        axes=str(args.axes),
         method=args.method,
+        scale_factors=scale_factors,
     )
-    scaler.scale(args.input_array, args.output_directory)
+
+    grp = zarr.open_group(args.output_directory, mode="a")
+
+    if args.copy_metadata:
+        print(f"copying attribute keys: {list(base.attrs.keys())}")
+        grp.attrs.update(base.attrs)
 
 
 def csv_to_labels(args: argparse.Namespace) -> None:
@@ -167,23 +180,24 @@ def main(args: list[str] | None = None) -> None:
     parser_scale.add_argument("input_array")
     parser_scale.add_argument("output_directory")
     parser_scale.add_argument(
-        "--labeled",
-        action="store_true",
-        help="assert that the list of unique pixel values doesn't change",
+        "axes", type=str, help="Dimensions of input data, i.e. 'zyx' or 'tczyx'."
     )
     parser_scale.add_argument(
         "--copy-metadata",
         action="store_true",
         help="copies the array metadata to the new group",
     )
     parser_scale.add_argument(
-        "--method", choices=list(Scaler.methods()), default="nearest"
+        "--method",
+        choices=["nearest", "resize", "laplacian", "local_mean", "zoom"],
+        default="resize",
     )
     parser_scale.add_argument(
         "--in-place", action="store_true", help="if true, don't write the base array"
     )
     parser_scale.add_argument("--downscale", type=int, default=2)
     parser_scale.add_argument("--max_layer", type=int, default=4)
+
     parser_scale.set_defaults(func=scale)
 
     # csv to label properties

ome_zarr/dask_utils.py

@@ -1,3 +1,4 @@
+from collections.abc import Sequence
 from typing import Any
 
 import numpy as np
@@ -8,15 +9,37 @@
 # See https://github.com/toloudis/ome-zarr-py/pull/1
 
 
+def _better_chunksize(
+    image: da.Array, factors: np.ndarray
+) -> tuple[Sequence[int], Sequence[int]]:
+    better_chunksize = tuple(
+        np.maximum(1, np.round(np.array(image.chunksize) * factors) / factors).astype(
+            int
+        )
+    )
+
+    # If E.g. we resize image from 6675 by 0.5 to 3337, factor is 0.49992509 so each
+    # chunk of size e.g. 1000 will resize to 499. When assumbled into a new array, the
+    # array will now be of size 3331 instead of 3337 because each of 6 chunks was
+    # smaller by 1. When we compute() this, dask will read 6 chunks of 1000 and expect
+    # last chunk to be 337 but instead it will only be 331.
+    # So we use ceil() here (and in resize_block) to round 499.925 up to chunk of 500
+    block_output_shape = tuple(
+        np.ceil(np.array(better_chunksize) * factors).astype(int)
+    )
+
+    return better_chunksize, block_output_shape
+
+
 def resize(
-    image: da.Array, output_shape: tuple[int, ...], *args: Any, **kwargs: Any
+    image: da.Array, output_shape: Sequence[int], *args: Any, **kwargs: Any
 ) -> da.Array:
     r"""
     Wrapped copy of "skimage.transform.resize"
     Resize image to match a certain size.
     :type image: :class:`dask.array`
     :param image: The dask array to resize
-    :type output_shape: tuple
+    :type output_shape: Sequence[int]
     :param output_shape: The shape of the resize array
     :type \*args: list
     :param \*args: Arguments of skimage.transform.resize
@@ -27,23 +50,9 @@ def resize(
     factors = np.array(output_shape) / np.array(image.shape).astype(float)
     # Rechunk the input blocks so that the factors achieve an output
     # blocks size of full numbers.
-    better_chunksize = tuple(
-        np.maximum(1, np.round(np.array(image.chunksize) * factors) / factors).astype(
-            int
-        )
-    )
+    better_chunksize, block_output_shape = _better_chunksize(image, factors)
     image_prepared = image.rechunk(better_chunksize)
 
-    # If E.g. we resize image from 6675 by 0.5 to 3337, factor is 0.49992509 so each
-    # chunk of size e.g. 1000 will resize to 499. When assumbled into a new array, the
-    # array will now be of size 3331 instead of 3337 because each of 6 chunks was
-    # smaller by 1. When we compute() this, dask will read 6 chunks of 1000 and expect
-    # last chunk to be 337 but instead it will only be 331.
-    # So we use ceil() here (and in resize_block) to round 499.925 up to chunk of 500
-    block_output_shape = tuple(
-        np.ceil(np.array(better_chunksize) * factors).astype(int)
-    )
-
     # Map overlap
     def resize_block(image_block: da.Array, block_info: dict) -> da.Array:
         # if the input block is smaller than a 'regular' chunk (e.g. edge of image)
@@ -62,6 +71,65 @@ def resize_block(image_block: da.Array, block_info: dict) -> da.Array:
     return output.rechunk(image.chunksize).astype(image.dtype)
 
 
+def local_mean(
+    image: da.Array, output_shape: Sequence[int], *args, **kwargs
+) -> da.Array:
+    """
+    Local mean downscaling.
+
+    :type image: :class:`dask.array.Array`
+    :param image: The dask array to resize
+    :type output_shape: Sequence[int]
+    :param output_shape: The shape of the resize array
+    :return: Resized image.
+    """
+    from skimage.transform import downscale_local_mean
+
+    factors = np.array(image.shape).astype(float) / np.array(output_shape)
+    better_chunksize, block_output_shape = _better_chunksize(image, 1 / factors)
+    image_prepared = image.rechunk(better_chunksize)
+
+    def local_mean_block(image_block: da.Array) -> da.Array:
+        local_mean = downscale_local_mean(
+            image_block, tuple(factors.astype(int)), *args, **kwargs
+        )
+        return local_mean.astype(image_block.dtype)
+
+    output_slices = tuple(slice(0, d) for d in output_shape)
+    output = da.map_blocks(
+        local_mean_block, image_prepared, dtype=image.dtype, chunks=block_output_shape
+    )[output_slices]
+    return output.rechunk(image.chunksize).astype(image.dtype)
+
+
+def zoom(image: da.Array, output_shape: Sequence[int], *args, **kwargs) -> da.Array:
+    """
+    Scipy zoom downscaling by integer factors.
+
+    :type image: :class:`dask.array.Array`
+    :param image: The dask array to resize
+    :type output_shape: Sequence[int]
+    :param output_shape: The shape of the resize array
+    :return: Resized image.
+    """
+    from scipy.ndimage import zoom
+
+    factors = np.array(image.shape).astype(float) / np.array(output_shape)
+    better_chunksize, block_output_shape = _better_chunksize(image, 1 / factors)
+    image_prepared = image.rechunk(better_chunksize)
+
+    def zoom_block(image_block: da.Array) -> da.Array:
+        zoomed = zoom(image_block, 1 / factors, order=1)
+        return zoomed.astype(image_block.dtype)
+
+    output_shape = tuple(d // f for d, f in zip(image.shape, factors))
+    output_slices = tuple(slice(0, d) for d in output_shape)
+    output = da.map_blocks(
+        zoom_block, image_prepared, dtype=image.dtype, chunks=block_output_shape
+    )[output_slices]
+    return output.rechunk(image.chunksize).astype(image.dtype)
+
+
 def downscale_nearest(image: da.Array, factors: tuple[int, ...]) -> da.Array:
     """
     Primitive downscaling by integer factors using stepped slicing.