from_images: python-native czi support and timepoint and channel multiplexing (#822)

* add python-native czi support and timepoint and channel multiplexing across layers

* added truncate_rgba_to_rgb arg and fixed imageIO reader channel handling

* changelog link to PR

* readd largest_segment_id

* apply PR feedback part 1

* apply PR feedback part 2: add comments, better var names, minor refactoring

* add more comments

* add max_layers argument

Author: jstriebel

webknossos/Changelog.md

@@ -13,8 +13,11 @@ For upgrade instructions, please check the respective *Breaking Changes* section
 [Commits](https://github.com/scalableminds/webknossos-libs/compare/v0.10.24...HEAD)
 
 ### Breaking Changes
+- `Dataset.from_images` now adds a layer per timepoint and per channel (if the data doesn't have 1 or 3 channels). [#822](https://github.com/scalableminds/webknossos-libs/pull/822)
 
 ### Added
+- Added python-native CZI support for `Dataset.from_images` or `dataset.add_layer_from_images`, without using bioformats. [#822](https://github.com/scalableminds/webknossos-libs/pull/822)
+- `dataset.add_layer_from_images` can add a layer per timepoint and per channel when passing `allow_multiple_layers=True`. [#822](https://github.com/scalableminds/webknossos-libs/pull/822)
 
 ### Changed
 

webknossos/poetry.lock

@@ -60,13 +60,15 @@ imagecodecs = { version = ">=2021.11.20", optional = true }
 JPype1 = { version = "^1.3.0", optional = true }
 pims = { version = "^0.6.0", optional = true }
 tifffile = { version = ">=2021.11.2", optional = true }
+pylibCZIrw = { version = "^3.2", optional = true }
 
 [tool.poetry.extras]
 pims = ["pims"]
 tifffile = ["pims", "tifffile"]
 imagecodecs = ["pims", "imagecodecs"]
 bioformats = ["pims","JPype1"]
-all = ["pims","tifffile","imagecodecs","JPype1",]
+czi = ["pims","pylibCZIrw"]
+all = ["pims","tifffile","imagecodecs","JPype1","pylibCZIrw",]
 
 [tool.poetry.dev-dependencies]
 # autoflake

webknossos/pyproject.toml

@@ -1,8 +1,9 @@
 from pathlib import Path
+from shutil import copy
 from tempfile import NamedTemporaryFile, TemporaryDirectory
 from time import gmtime, strftime
 from typing import Any, Dict, List, Tuple, Union
-from zipfile import ZipFile
+from zipfile import BadZipFile, ZipFile
 
 import httpx
 import numpy as np
@@ -118,7 +119,7 @@ def test_repo_images(
     return ds
 
 
-def download_and_unpack(url: str, out_path: Path) -> None:
+def download_and_unpack(url: str, out_path: Path, filename: str) -> None:
     with NamedTemporaryFile() as download_file:
         with httpx.stream("GET", url) as response:
             total = int(response.headers["Content-Length"])
@@ -130,8 +131,12 @@ def download_and_unpack(url: str, out_path: Path) -> None:
                     progress.update(
                         download_task, completed=response.num_bytes_downloaded
                     )
-        with ZipFile(download_file, "r") as zip_file:
-            zip_file.extractall(out_path)
+        try:
+            with ZipFile(download_file, "r") as zip_file:
+                zip_file.extractall(out_path)
+        except BadZipFile:
+            out_path.mkdir(parents=True, exist_ok=True)
+            copy(download_file.name, out_path / filename)
 
 
 BIOFORMATS_ARGS = [
@@ -142,6 +147,7 @@ def download_and_unpack(url: str, out_path: Path) -> None:
         "uint8",
         3,
         (320, 240, 108),
+        1,
     ),
     (
         "https://samples.scif.io/wtembryo.zip",
@@ -150,6 +156,7 @@ def download_and_unpack(url: str, out_path: Path) -> None:
         "uint8",
         3,
         (240, 320, 1),
+        1,
     ),
     (
         "https://samples.scif.io/HEART.zip",
@@ -158,6 +165,7 @@ def download_and_unpack(url: str, out_path: Path) -> None:
         "uint8",
         1,
         (512, 512, 30),
+        1,
     ),
     (
         "https://samples.scif.io/sdub.zip",
@@ -166,6 +174,16 @@ def download_and_unpack(url: str, out_path: Path) -> None:
         "uint8",
         1,
         (192, 128, 9),
+        1,
+    ),
+    (
+        "https://samples.scif.io/sdub.zip",
+        "sdub*.pic",
+        {"allow_multiple_layers": True},
+        "uint8",
+        1,
+        (192, 128, 9),
+        12,
     ),
     (
         "https://samples.scif.io/test-avi.zip",
@@ -174,12 +192,13 @@ def download_and_unpack(url: str, out_path: Path) -> None:
         "uint8",
         3,
         (206, 218, 36),
+        1,
     ),
 ]
 
 
 @pytest.mark.parametrize(
-    "url, filename, kwargs, dtype, num_channels, size", BIOFORMATS_ARGS
+    "url, filename, kwargs, dtype, num_channels, size, num_layers", BIOFORMATS_ARGS
 )
 def test_bioformats(
     tmp_path: Path,
@@ -189,9 +208,10 @@ def test_bioformats(
     dtype: str,
     num_channels: int,
     size: Tuple[int, int, int],
+    num_layers: int,
 ) -> wk.Dataset:
     unzip_path = tmp_path / "unzip"
-    download_and_unpack(url, unzip_path)
+    download_and_unpack(url, unzip_path, filename)
     ds = wk.Dataset(tmp_path / "ds", (1, 1, 1))
     with wk.utils.get_executor_for_args(None) as executor:
         l = ds.add_layer_from_images(
@@ -205,10 +225,21 @@ def test_bioformats(
         assert l.dtype_per_channel == np.dtype(dtype)
         assert l.num_channels == num_channels
         assert l.bounding_box == wk.BoundingBox(topleft=(0, 0, 0), size=size)
+    assert len(ds.layers) == num_layers
     return ds
 
 
 TEST_IMAGES_ARGS = [
+    (
+        # published with CC0 license, taken from
+        # https://doi.org/10.6084/m9.figshare.c.3727411_D391.v1
+        "https://figshare.com/ndownloader/files/8909407",
+        "embedded_NCI_mono_matrigelcollagen_docetaxel_day10_sample10.czi",
+        {},
+        "uint16",
+        1,
+        (512, 512, 30),
+    ),
     (
         "https://samples.scif.io/test-gif.zip",
         "scifio-test.gif",
@@ -257,7 +288,7 @@ def test_test_images(
     size: Tuple[int, int, int],
 ) -> wk.Dataset:
     unzip_path = tmp_path / "unzip"
-    download_and_unpack(url, unzip_path)
+    download_and_unpack(url, unzip_path, filename)
     ds = wk.Dataset(tmp_path / "ds", (1, 1, 1))
     with wk.utils.get_executor_for_args(None) as executor:
         l_bio = ds.add_layer_from_images(

webknossos/tests/dataset/test_add_layer_from_images.py

@@ -0,0 +1,118 @@
+from contextlib import contextmanager
+from os import PathLike
+from pathlib import Path
+from typing import Iterator, List, Set
+
+import numpy as np
+
+try:
+    from pims import FramesSequenceND
+except ImportError as e:
+    raise RuntimeError(
+        "Cannot import pims, please install it e.g. using 'webknossos[all]'"
+    ) from e
+
+try:
+    from pylibCZIrw import czi as pyczi
+except ImportError as e:
+    raise RuntimeError(
+        "Cannot import pylibCZIrw, please install it e.g. using 'webknossos[czi]'"
+    ) from e
+
+PIXEL_TYPE_TO_DTYPE = {
+    "Gray8": "<u1",
+    "Gray16": "<u2",
+    "Gray32Float": "<f4",
+    "Bgr24": "<u1",
+    "Bgr48": "<u2",
+    "Bgr96Float": "<f4",
+    "Bgra32": "<4u1",
+    "Gray64ComplexFloat": "<F8",
+    "Bgr192ComplexFloat": "<F8",
+    "Gray32": "<i4",
+    "Gray64": "<i8",
+}
+
+
+class PimsCziReader(FramesSequenceND):
+    @classmethod
+    def class_exts(cls) -> Set[str]:
+        return {".czi"}
+
+    # class_priority is used in pims to pick the reader with the highest priority.
+    # Default is 10, and bioformats priority (which is the only other reader supporting czi) is 2.
+    # See http://soft-matter.github.io/pims/v0.6.1/custom_readers.html#plugging-into-pims-s-open-function
+    class_priority = 20
+
+    def __init__(self, path: PathLike, czi_channel: int = 0) -> None:
+        self.path = Path(path)
+        self.czi_channel = czi_channel
+        super().__init__()
+        with self.czi_file() as czi_file:
+            for axis, (
+                start,
+                length,
+            ) in czi_file.total_bounding_box.items():
+                axis = axis.lower()
+                if axis == "c":
+                    continue
+                assert start == 0
+                if axis not in "xy" and length == 1:
+                    # not propagating axes of length one
+                    continue
+                self._init_axis(axis, length)
+            czi_pixel_type = czi_file.get_channel_pixel_type(self.czi_channel)
+            if czi_pixel_type.startswith("Bgra"):
+                self._init_axis("c", 4)
+            elif czi_pixel_type.startswith("Bgr"):
+                self._init_axis("c", 3)
+            elif czi_pixel_type.startswith("Gray"):
+                self._init_axis("c", 1)
+            elif czi_pixel_type == "Invalid":
+                raise ValueError(
+                    f"czi_channel {self.czi_channel} does not exist in {self.path}"
+                )
+            else:
+                raise ValueError(
+                    f"Got unsupported czi pixel-type {czi_pixel_type} in {self.path}"
+                )
+
+        self._register_get_frame(self.get_frame_2D, "yxc")
+
+    @contextmanager
+    def czi_file(self) -> Iterator[pyczi.CziReader]:
+        with pyczi.open_czi(str(self.path)) as czi_file:
+            yield czi_file
+
+    def available_czi_channels(self) -> List[int]:
+        with self.czi_file() as czi_file:
+            return sorted(czi_file.pixel_types.keys())
+
+    @property  # potential @cached_property for py3.8+
+    def pixel_type(self) -> np.dtype:
+        with self.czi_file() as czi_file:
+            return np.dtype(
+                PIXEL_TYPE_TO_DTYPE[czi_file.get_channel_pixel_type(self.czi_channel)]
+            )
+
+    def get_frame_2D(self, **ind: int) -> np.ndarray:
+        plane = {k.upper(): v for k, v in ind.items()}
+
+        # safe-guard against x/y in ind argument,
+        # we always read the whole slice here:
+        plane.pop("X", None)
+        plane.pop("Y", None)
+
+        plane["C"] = self.czi_channel
+        with self.czi_file() as czi_file:
+            a = czi_file.read(plane=plane)
+            num_channels = a.shape[-1]
+            if num_channels == 3:
+                # convert from bgr to rgb
+                a = np.flip(a, axis=-1)
+            elif num_channels == 4:
+                # convert from bgra to rgba
+                a_red = a[:, :, 2].copy(order="K")
+                a[:, :, 2] = a[:, :, 0]
+                a[:, :, 0] = a_red
+            return a