Tiffreader for nd data (#1043)

* WIP change axis recognition of tifffiles.

* Add basic custom tiff reader.

* Revert changes on _assume_color_channel method.

* Adapt test_repo_images.

* Test out _register_get_frame changes.

* Work on fixing image conversion for images without z axis.

* Update expected_bbox function to deliver a mag1 bbox.

* Working on fixing tiff_reader.

* reverting some changes at the pims image reader class.

* Working on data extraction from tiff reader to desired Pims Sequence.

* Use memmap for tiffreader to avoid OOM issue.

* Implement np memmap for pims tiff reader because nd ome tiffs are not memmappable with tifffile.memmap.

* Adapt calculation of num_channels and set default axes for get_frame.

* Update test parameters for test multiple multitiffs.

* Update changelog and update error message for unsupported conversion to WKW.

* Update webknossos/webknossos/dataset/_utils/pims_images.py

Co-authored-by: Norman Rzepka <[email protected]>

* Started to implement requested changes.

* Fix issue with zarr array creation and add zarr to supported dataformats for nd data.

* Add an entry for merge-fallback in cli docs.

* Add 2 additional testdatasets and adapt logging behaviour.

* Remove user warning from conversion test.

* Run linter.

* Add pytest.warn() to cli_test.py.

* Revert some changes that are not relevant for this PR.

* Add comment for chunk shape calculation.

* Update changelog.

---------

Co-authored-by: Norman Rzepka <[email protected]>

Author: markbader

webknossos/Changelog.md

@@ -13,6 +13,7 @@ For upgrade instructions, please check the respective _Breaking Changes_ section
 [Commits](https://github.com/scalableminds/webknossos-libs/compare/v0.14.22...HEAD)
 
 ### Breaking Changes
+- Preferring a custom tiff reader over the default PIMS reader to convert tiff files. This change enables the recognition of axis information and the support of tifffiles with more than 3 dimensions. However, it also leads to changed behavior when converting tiff files. Tiffs with axes other than c, x, y, and z, with a shape bigger than 1, are no longer supported for conversion to WKW. Please convert these files to Zarr or Zarr3 Datasets instead. [#1043](https://github.com/scalableminds/webknossos-libs/pull/1043)
 
 ### Added
 - Added a pixel level heuristic for distinguishing color and segmentation layers when importing image data with the `from_images` or `add_layer_from_images` method. [#1007](https://github.com/scalableminds/webknossos-libs/pull/1007)

webknossos/README.md

@@ -61,4 +61,5 @@ Excerpts for testing purposes have been sampled from:
 * Dow Jacobo Hossain Siletti Hudspeth (2018). **Connectomics of the zebrafish's lateral-line neuromast reveals wiring and miswiring in a simple microcircuit.** eLife. [DOI:10.7554/eLife.33988](https://elifesciences.org/articles/33988)
 * Zheng Lauritzen Perlman Robinson Nichols Milkie Torrens Price Fisher Sharifi Calle-Schuler Kmecova Ali Karsh Trautman Bogovic Hanslovsky Jefferis Kazhdan Khairy Saalfeld Fetter Bock (2018). **A Complete Electron Microscopy Volume of the Brain of Adult Drosophila melanogaster.** Cell. [DOI:10.1016/j.cell.2018.06.019](https://www.cell.com/cell/fulltext/S0092-8674(18)30787-6). License: [CC BY-NC 4.0](https://creativecommons.org/licenses/by-nc/4.0/)
 * Bosch Ackels Pacureanu et al (2022). **Functional and multiscale 3D structural investigation of brain tissue through correlative in vivo physiology, synchrotron microtomography and volume electron microscopy.** Nature Communications. [DOI:10.1038/s41467-022-30199-6](https://www.nature.com/articles/s41467-022-30199-6)
-* Hanke, M., Baumgartner, F. J., Ibe, P., Kaule, F. R., Pollmann, S., Speck, O., Zinke, W. & Stadler, J. (2014). **A high-resolution 7-Tesla fMRI dataset from complex natural stimulation with an audio movie. Scientific Data, 1:140003.** [DOI:10.1038/sdata.2014.3](http://www.nature.com/articles/sdata20143)
+* Hanke, M., Baumgartner, F. J., Ibe, P., Kaule, F. R., Pollmann, S., Speck, O., Zinke, W. & Stadler, J. (2014). **A high-resolution 7-Tesla fMRI dataset from complex natural stimulation with an audio movie. Scientific Data, 1:140003.** [DOI:10.1038/sdata.2014.3](http://www.nature.com/articles/sdata20143)
+* Sample OME-TIFF files (c) by the OME Consortium [https://downloads.openmicroscopy.org/images/OME-TIFF/2016-06/bioformats-artificial/](https://downloads.openmicroscopy.org/images/OME-TIFF/2016-06/bioformats-artificial/)

webknossos/examples/convert_4d_tiff.py

@@ -10,7 +10,6 @@ def main() -> None:
         "testoutput/4D_series",
         voxel_size=(10, 10, 10),
         data_format="zarr3",
-        use_bioformats=True,
     )
 
     # Access the first color layer and the Mag 1 view of this layer

webknossos/testdata/4D/multi_channel_z_series/multi-channel-z-series.ome.tif

@@ -50,7 +50,6 @@ def test_compare_nd_tifffile(tmp_path: Path) -> None:
         layer_name="color",
         category="color",
         topleft=(2, 55, 100, 100),
-        use_bioformats=True,
         data_format="zarr3",
         chunk_shape=(8, 8, 8),
         chunks_per_shard=(8, 8, 8),
@@ -69,13 +68,14 @@ def test_compare_nd_tifffile(tmp_path: Path) -> None:
 
 
 REPO_IMAGES_ARGS: List[
-    Tuple[Union[str, List[Path]], Dict[str, Any], str, int, Tuple[int, int, int]]
+    Tuple[Union[str, List[Path]], Dict[str, Any], str, int, int, Tuple[int, ...]]
 ] = [
     (
         "testdata/tiff/test.*.tiff",
         {"category": "segmentation"},
         "uint8",
         1,
+        1,
         (265, 265, 257),
     ),
     (
@@ -87,34 +87,39 @@ def test_compare_nd_tifffile(tmp_path: Path) -> None:
         {},
         "uint8",
         1,
+        1,
         (265, 265, 3),
     ),
     (
         "testdata/rgb_tiff/test_rgb.tif",
         {"mag": 2},
         "uint8",
-        3,
-        (64, 64, 2),
+        1,
+        1,
+        (64, 64, 6),
     ),
     (
         "testdata/rgb_tiff",
-        {"mag": 2, "channel": 1, "dtype": "uint32"},
+        {"mag": 2, "channel": 0, "dtype": "uint32"},
         "uint32",
         1,
-        (64, 64, 2),
+        1,
+        (64, 64, 6),
     ),
     (
         "testdata/temca2/*/*/*.jpg",
         {"flip_x": True, "batch_size": 2048},
         "uint8",
         1,
+        1,
         (1024, 1024, 12),
     ),
     (
         "testdata/temca2",
         {"flip_z": True, "batch_size": 2048},
         "uint8",
         1,
+        1,
         # The topmost folder contains an extra image,
         # which is included here as well, but not in
         # the glob pattern above. Therefore z is +1.
@@ -125,25 +130,73 @@ def test_compare_nd_tifffile(tmp_path: Path) -> None:
         {"flip_y": True},
         "uint8",
         1,
+        1,
         (2970, 2521, 4),
     ),
-    ("testdata/various_tiff_formats/test_CS.tif", {}, "uint8", 3, (128, 128, 320)),
-    ("testdata/various_tiff_formats/test_C.tif", {}, "uint8", 1, (128, 128, 320)),
+    (
+        "testdata/various_tiff_formats/test_CS.tif",
+        {"data_format": "zarr3", "allow_multiple_layers": True},
+        "uint8",
+        1,
+        5,
+        (3, 64, 128, 128),
+    ),
+    (
+        "testdata/various_tiff_formats/test_C.tif",
+        {"allow_multiple_layers": True},
+        "uint8",
+        1,
+        5,
+        (128, 128, 64),
+    ),
     # same as test_C.tif above, but as a single file in a folder:
-    ("testdata/single_multipage_tiff_folder", {}, "uint8", 1, (128, 128, 320)),
-    ("testdata/various_tiff_formats/test_I.tif", {}, "uint32", 1, (64, 128, 64)),
-    ("testdata/various_tiff_formats/test_S.tif", {}, "uint16", 3, (128, 128, 64)),
+    (
+        "testdata/single_multipage_tiff_folder",
+        {"allow_multiple_layers": True},
+        "uint8",
+        1,
+        5,
+        (128, 128, 64),
+    ),
+    ("testdata/various_tiff_formats/test_I.tif", {}, "uint32", 1, 1, (64, 128, 64)),
+    (
+        "testdata/various_tiff_formats/test_S.tif",
+        {"data_format": "zarr3"},
+        "uint16",
+        1,
+        1,
+        (3, 64, 128, 128),
+    ),
+    (
+        "testdata/4D/single_channel/single-channel.ome.tiff",
+        {},
+        "int8",
+        1,
+        1,
+        (439, 167, 1),
+    ),
+    (
+        "testdata/4D/multi_channel_z_series/multi-channel-z-series.ome.tif",
+        {"allow_multiple_layers": True},
+        "int8",
+        1,
+        3,
+        (439, 167, 5),
+    ),
 ]
 
 
-@pytest.mark.parametrize("path, kwargs, dtype, num_channels, size", REPO_IMAGES_ARGS)
+@pytest.mark.parametrize(
+    "path, kwargs, dtype, num_channels, num_layers, size", REPO_IMAGES_ARGS
+)
 def test_repo_images(
     tmp_path: Path,
     path: str,
     kwargs: Dict,
     dtype: str,
     num_channels: int,
-    size: Tuple[int, int, int],
+    num_layers: int,
+    size: Tuple[int, ...],
 ) -> wk.Dataset:
     with wk.utils.get_executor_for_args(None) as executor:
         ds = wk.Dataset(tmp_path, (1, 1, 1))
@@ -157,7 +210,8 @@ def test_repo_images(
         )
         assert layer.dtype_per_channel == np.dtype(dtype)
         assert layer.num_channels == num_channels
-        assert layer.bounding_box == wk.BoundingBox(topleft=(0, 0, 0), size=size)
+        assert len(ds.layers) == num_layers
+        assert layer.bounding_box.size.to_tuple() == size
         if isinstance(layer, wk.SegmentationLayer):
             assert layer.largest_segment_id is not None
             assert layer.largest_segment_id > 0
@@ -396,7 +450,7 @@ def test_test_images(
         else:
             assert l_bio.dtype_per_channel == np.dtype(dtype)
             assert l_bio.num_channels == num_channels
-            assert l_bio.bounding_box == wk.BoundingBox(topleft=(0, 0, 0), size=size)
+            assert l_bio.bounding_box.size.to_tuple() == size
         l_normal = ds.add_layer_from_images(
             path,
             layer_name="normal_" + layer_name,
@@ -407,7 +461,7 @@ def test_test_images(
         )
         assert l_normal.dtype_per_channel == np.dtype(dtype)
         assert l_normal.num_channels == num_channels
-        assert l_normal.bounding_box == wk.BoundingBox(topleft=(0, 0, 0), size=size)
+        assert l_normal.bounding_box.size.to_tuple() == size
         if l_bio is not None:
             assert np.array_equal(
                 l_bio.get_finest_mag().read(), l_normal.get_finest_mag().read()

webknossos/testdata/4D/single_channel/single-channel.ome.tiff

@@ -44,15 +44,24 @@ def test_multiple_multitiffs(tmp_path: Path) -> None:
         TESTDATA_DIR / "various_tiff_formats",
         tmp_path,
         (1, 1, 1),
+        data_format="zarr3",
         layer_name="tiffs",
     )
-    assert len(ds.layers) == 4
+    assert len(ds.layers) == 12
 
     expected_dtype_channels_size_per_layer = {
-        "tiffs_test_CS.tif": ("uint8", 3, (128, 128, 320)),
-        "tiffs_test_C.tif": ("uint8", 1, (128, 128, 320)),
+        "tiffs_test_CS.tif__channel0": ("uint8", 1, (3, 64, 128, 128)),
+        "tiffs_test_CS.tif__channel1": ("uint8", 1, (3, 64, 128, 128)),
+        "tiffs_test_CS.tif__channel2": ("uint8", 1, (3, 64, 128, 128)),
+        "tiffs_test_CS.tif__channel3": ("uint8", 1, (3, 64, 128, 128)),
+        "tiffs_test_CS.tif__channel4": ("uint8", 1, (3, 64, 128, 128)),
+        "tiffs_test_C.tif__channel0": ("uint8", 1, (128, 128, 64)),
+        "tiffs_test_C.tif__channel1": ("uint8", 1, (128, 128, 64)),
+        "tiffs_test_C.tif__channel2": ("uint8", 1, (128, 128, 64)),
+        "tiffs_test_C.tif__channel3": ("uint8", 1, (128, 128, 64)),
+        "tiffs_test_C.tif__channel4": ("uint8", 1, (128, 128, 64)),
         "tiffs_test_I.tif": ("uint32", 1, (64, 128, 64)),
-        "tiffs_test_S.tif": ("uint16", 3, (128, 128, 64)),
+        "tiffs_test_S.tif": ("uint16", 1, (3, 64, 128, 128)),
     }
 
     for layer_name, layer in ds.layers.items():

webknossos/tests/dataset/test_add_layer_from_images.py

@@ -593,7 +593,11 @@ def create(cls, path: Path, array_info: ArrayInfo) -> "ZarritaArray":
             ArrayV2.create(
                 store=path,
                 shape=(array_info.shape),
-                chunks=(array_info.num_channels,) + array_info.chunk_shape.to_tuple(),
+                chunks=(array_info.num_channels,)
+                + tuple(
+                    getattr(array_info.chunk_shape, axis, 1)
+                    for axis in array_info.dimension_names[1:]
+                ),  # The chunk shape consists of the number of channels, the x, y, and z dimensions of a chunk, and 1 for all other dimensions.
                 dtype=array_info.voxel_type,
                 compressor=(
                     {"id": "blosc", "cname": "zstd", "clevel": 5}

webknossos/tests/dataset/test_from_images.py

@@ -27,7 +27,6 @@
 from webknossos.geometry.bounding_box import BoundingBox
 from webknossos.geometry.nd_bounding_box import NDBoundingBox
 
-# pylint: disable=unused-import
 try:
     from .pims_czi_reader import PimsCziReader
 except ImportError:
@@ -48,6 +47,11 @@
 except ImportError:
     pass
 
+try:
+    from .pims_tiff_reader import PimsTiffReader  # noqa: F401 unused-import
+except ImportError:
+    pass
+
 
 from ...geometry.vec_int import VecInt
 from ..mag_view import MagView
@@ -120,7 +124,7 @@ def __init__(
 
         ## attributes that will be set in __init__()
         # _bundle_axes
-        self._iter_axes = None
+        self._iter_axes: List[str] = []
         self._iter_loop_size = None
         self._possible_layers = {}
 
@@ -216,18 +220,15 @@ def __init__(
                 if len(images.shape) == 2:
                     # Assume yx
                     self._bundle_axes = ["y", "x"]
-                    self._iter_axes = []
                 elif len(images.shape) == 3:
                     # Assume yxc, cyx or zyx
                     if _assume_color_channel(images.shape[2], images.dtype):
                         self._bundle_axes = ["y", "x", "c"]
-                        self._iter_axes = []
                     elif images.shape[0] == 1 or (
                         _allow_channels_first
                         and _assume_color_channel(images.shape[0], images.dtype)
                     ):
                         self._bundle_axes = ["c", "y", "x"]
-                        self._iter_axes = []
                     else:
                         self._bundle_axes = ["y", "x"]
                         self._iter_axes = ["z"]
@@ -277,18 +278,16 @@ def __init__(
         #########################
 
         with self._open_images() as images:
-            try:
-                c_index = self._bundle_axes.index("c")
-                if isinstance(images, list):
-                    images_shape = (len(images),) + cast(
-                        pims.FramesSequence, images[0]
-                    ).shape
+            if "c" in self._bundle_axes:
+                if isinstance(images, pims.FramesSequenceND):
+                    self.num_channels = images.sizes.get("c", 1)
+                elif isinstance(images, list):
+                    self.num_channels = cast(pims.FramesSequence, images[0]).shape[
+                        self._bundle_axes.index("c")
+                    ]
                 else:
-                    images_shape = images.shape  # pylint: disable=no-member
-
-                self.num_channels = images_shape[c_index + 1]
-
-            except ValueError:
+                    self.num_channels = images.shape[self._bundle_axes.index("c") + 1]
+            else:
                 self.num_channels = 1
 
         self._first_n_channels = None
@@ -302,7 +301,7 @@ def __init__(
                 self._possible_layers["channel"] = [0, 1]
                 self.num_channels = 1
                 self._channel = 0
-            elif self.num_channels > 3:
+            elif self.num_channels >= 3:
                 self._possible_layers["channel"] = list(range(0, self.num_channels))
                 self.num_channels = 3
                 self._first_n_channels = 3
@@ -483,13 +482,15 @@ def _open_images(
                         images.bundle_axes = self._bundle_axes
                         images.iter_axes = self._iter_axes
                 else:
-                    if self._timepoint is not None:
-                        images = images[self._timepoint]
-                        if self._iter_axes and "t" in self._iter_axes:
-                            self._iter_axes.remove("t")
-                    if self._iter_axes == []:
-                        # add outer list to wrap 2D images as 3D-like structure
-                        images = [images]
+                    if hasattr(self, "_bundle_axes"):
+                        # first part of __init__() has happened
+                        if self._timepoint is not None:
+                            images = images[self._timepoint]
+                            if "t" in self._iter_axes:
+                                self._iter_axes.remove("t")
+                        if not self._iter_axes:
+                            # add outer list to wrap 2D images as 3D-like structure
+                            images = [images]
                 yield images
 
     def copy_to_view(
@@ -519,7 +520,7 @@ def copy_to_view(
         max_value = 0
 
         with self._open_images() as images:
-            if self._iter_axes is not None and self._iter_loop_size is not None:
+            if self._iter_axes and self._iter_loop_size is not None:
                 # select the range of images that represents one xyz combination in the mag_view
                 lower_bounds = sum(
                     self._iter_loop_size[axis_name]
@@ -529,7 +530,7 @@ def copy_to_view(
                 upper_bounds = lower_bounds + mag_view.bounding_box.get_shape("z")
                 images = images[lower_bounds:upper_bounds]
             if self._flip_z:
-                images = images[::-1]  # pylint: disable=unsubscriptable-object
+                images = images[::-1]
 
             with mag_view.get_buffered_slice_writer(
                 # Previously only z_start and its end were important, now the slice writer needs to know
@@ -600,7 +601,7 @@ def expected_bbox(self) -> NDBoundingBox:
                     ).shape
 
                 else:
-                    images_shape = images.shape  # pylint: disable=no-member
+                    images_shape = images.shape
                 if len(images_shape) == 3:
                     axes = ("z", "y", "x")
                 else:
@@ -633,10 +634,7 @@ def expected_bbox(self) -> NDBoundingBox:
                     axes_names = (self._iter_axes or []) + [
                         axis for axis in self._bundle_axes if axis != "c"
                     ]
-                    axes_sizes = [
-                        images.sizes[axis]  # pylint: disable=no-member
-                        for axis in axes_names
-                    ]
+                    axes_sizes = [images.sizes[axis] for axis in axes_names]
                     axes_index = list(range(1, len(axes_names) + 1))
                     topleft = VecInt.zeros(tuple(axes_names))