Add support for OME-NGFF ZARR inputs

README.rst

@@ -52,7 +52,8 @@ Sample data is available and described on Zenodo:
    :target: https://doi.org/10.5281/zenodo.4304786
 
 Any image and transform file format supported by `SimpleITK's IO <https://simpleitk.readthedocs.io/en/master/IO.html>`_
-can be use by sitk-ibex. The 3D `nrrd` format, and `txt` transform file extension are recommended.
+can be used by sitk-ibex. The NRRD or`NGFF <https://ngff.openmicroscopy.org/latest/>`_ image formats, and `txt` transform file
+extension are recommended.
 
 
 Example
@@ -82,6 +83,32 @@ Then apply the registration transform by resampling all channels of the the inpu
  python -m sitkibex resample "spleen_panel2.nrrd@CD4 AF594" spleen_panel3.nrrd tx_p2_to_p3.txt \
         -o spleen_onto_p2_panel3.nrrd
 
+Additional Example
+------------------
+
+Additional sample data:
+
+.. image:: https://zenodo.org/badge/DOI/10.5281/zenodo.4632320.svg
+   :target: https://doi.org/10.5281/zenodo.4632320
+
+The sample Imaris files can be converted to OME-NGFF ZARR with
+`bioformats2raw <https://github.com/glencoesoftware/bioformats2raw/releases>`_. The ims files contains one series, and
+for simplicity, the structure is generated without a series index in the hierarchy with the following commands::
+
+ bioformats2raw  --series 0 --scale-format-string '%2$d/'  Human_Spleen_Panel1.ims Human_Spleen_Panel1.zarr
+ bioformats2raw  --series 0 --scale-format-string '%2$d/'  Human_Spleen_Panel2.ims Human_Spleen_Panel2.zarr
+ bioformats2raw  --series 0 --scale-format-string '%2$d/'  Human_Spleen_Panel2.ims Human_Spleen_Panel3.zarr
+
+These images will be registered based on the common "Hoechst". The name of the channels are lost in this conversion from
+Imaris to OME-NGFF ZARR. Some conversions produce "omera" metadata in the ZARR file which contains channel labels, which
+can be used. When the channel labels are unavailable, the channel index can be used such as the following commands::
+
+ python -m sitkibex registration --affine "Human_Spleen_Panel2.zarr@3" "Human_Spleen_Panel1.zarr@2" tx_p2_to_p1.txt
+ python -m sitkibex registration --affine "Human_Spleen_Panel2.zarr@3" "Human_Spleen_Panel3.zarr@4" tx_p2_to_p3.txt
+
+The quick 2D visualization can be run similarly to the NRRD example. The OME-NGFF ZARR files are not supported for
+writing, so the resample command can produce NRRD files as well.
+
 
 How to Cite
 -----------

requirements.txt

@@ -1,3 +1,4 @@
 SimpleITK >= 2.0.0
 click
-numpy
+numpy
+zarr

setup.py

@@ -46,5 +46,8 @@
     python_requires=">=3.5",
     install_requires=requirements,
     tests_require=dev_requirements,
+    extras_require={
+        "zarr": ["zarr"],
+    },
     setup_requires=["setuptools_scm"],
 )

sitkibex/cli.py

@@ -150,8 +150,8 @@ def cli(**kwargs):
     help="Use wall-clock instead of a fixed seed for random initialization.",
 )
 @click.option("--samples-per-parameter", default=5000, type=int, show_default=True)
-@click.argument("fixed_image", type=ImagePathChannel(exists=True, dir_okay=False, resolve_path=True))
-@click.argument("moving_image", type=ImagePathChannel(exists=True, dir_okay=False, resolve_path=True))
+@click.argument("fixed_image", type=ImagePathChannel(exists=True, dir_okay=True, resolve_path=True))
+@click.argument("moving_image", type=ImagePathChannel(exists=True, dir_okay=True, resolve_path=True))
 @click.argument("output_transform", type=click.Path(exists=False, resolve_path=True))
 def reg_cli(fixed_image, moving_image, output_transform, **kwargs):
     """Perform registration to solve for an OUTPUT_TRANSFORM mapping points from the FIXED_IMAGE to the MOVING_IMAGE."""
@@ -218,8 +218,8 @@ def r(img, channel_name, bin_xy):
     help="filename for output image, if not provided the SimpleITK Show method is called",
     type=click.Path(exists=False, resolve_path=True),
 )
-@click.argument("fixed_image", type=ImagePathChannel(exists=True, dir_okay=False, resolve_path=True))
-@click.argument("moving_image", type=ImagePathChannel(exists=True, dir_okay=False, resolve_path=True))
+@click.argument("fixed_image", type=ImagePathChannel(exists=True, dir_okay=True, resolve_path=True))
+@click.argument("moving_image", type=ImagePathChannel(exists=True, dir_okay=True, resolve_path=True))
 @click.argument("transform", required=False, type=click.Path(exists=True, dir_okay=False, resolve_path=True))
 def resample_cli(fixed_image, moving_image, transform, **kwargs):
     """Create new image by transforming the MOVING_IMAGE onto the FIXED_IMAGE.
@@ -255,7 +255,7 @@ def binner(image):
 
     moving_img = binner(im_read_channel(moving_image, moving_channel_name))
 
-    if fixed_channel_name is None:
+    if fixed_channel_name is None and os.path.isfile(fixed_image):
         reader = sitk.ImageFileReader()
         reader.SetFileName(fixed_image)
         reader.ReadImageInformation()

sitkibex/io.py

@@ -18,6 +18,16 @@
 import os.path
 import logging
 from .xml_info import XMLInfo, OMEInfo
+from pathlib import Path
+import numpy as np
+
+try:
+    import zarr
+
+    _has_zarr = True
+except ImportError:
+    _has_zarr = False
+
 
 _logger = logging.getLogger(__name__)
 
@@ -33,6 +43,58 @@
 }
 
 
+def _zarr_read_channel(filename: Path, channel=None) -> sitk.Image:
+    """
+    Read a channel from a zarr file.
+    """
+
+    store = zarr.DirectoryStore(filename)
+    zarr_group = zarr.open_group(store=store, mode="r")
+
+    axes = zarr_group.attrs["multiscales"][0]["axes"]
+    ds_attr = zarr_group.attrs["multiscales"][0]["datasets"][0]
+
+    arr = zarr_group[ds_attr["path"]]
+    _logger.debug(arr)
+
+    spacing = ds_attr["coordinateTransformations"][0]["scale"]
+
+    axes_names = [ax["name"].upper() for ax in axes]
+
+    if axes_names != list("TCZYX"):
+        raise ValueError(f"Only TCZYX axes are supported, not {axes_names}.")
+
+    if arr.shape[0] > 1:
+        raise ValueError("Only single time point is supported.")
+
+    if channel is None:
+        arr = np.moveaxis(arr[0, ...], 0, -1)
+        img = sitk.GetImageFromArray(arr.astype(arr.dtype.newbyteorder("=")), isVector=True)
+    else:
+
+        if isinstance(channel, int):
+            channel_number = channel
+        else:
+            if "omero" not in zarr_group.attrs or "channels" not in zarr_group.attrs["omero"]:
+                raise ValueError("No OMERO metadata. Only integer channel numbers are supported.")
+            omero_channel_labels = [c["label"] for c in zarr_group.attrs["omero"]["channels"]]
+
+            if channel not in omero_channel_labels:
+                raise Exception(f'Channel name "{channel}" is not in OMERO marker list: {omero_channel_labels}.')
+
+            channel_number = omero_channel_labels.index(channel)
+
+        if channel_number >= arr.shape[1] or channel_number < 0:
+            raise ValueError("Channel number is out of range.")
+
+        arr = arr[0, channel_number, ...]
+        img = sitk.GetImageFromArray(arr.astype(arr.dtype.newbyteorder("=")), isVector=False)
+
+    img.SetSpacing(spacing[2:])
+
+    return img
+
+
 def im_read_channel(filename, channel=None):  # noqa: C901
     """
     Read a channel from an image file.
@@ -48,8 +110,15 @@ def im_read_channel(filename, channel=None):  # noqa: C901
 
     """
 
+    filename = Path(filename)
+    if filename.is_dir() and (filename / ".zattrs").exists():
+        if _has_zarr:
+            return _zarr_read_channel(filename, channel)
+        else:
+            raise ImportError("zarr is not installed.")
+
     reader = sitk.ImageFileReader()
-    reader.SetFileName(filename)
+    reader.SetFileName(str(filename))
 
     if channel is None:
         _logger.info('Reading whole "{}" image file.'.format(filename))
@@ -69,7 +138,6 @@ def im_read_channel(filename, channel=None):  # noqa: C901
         channel_number = channel
 
     else:
-
         IMAGE_DESCRIPTION = "ImageDescription"
         IMARIS_CHANNEL_INFORMATION = "imaris_channels_information"