From bba24f2156fffd9e79df4c0e48b72f54276eaf85 Mon Sep 17 00:00:00 2001 From: "w.pomp" Date: Thu, 26 Mar 2026 16:39:16 +0100 Subject: [PATCH] - add support for tiled czi's --- ndbioimage/__init__.py | 33 +++++---- ndbioimage/readers/cziread.py | 123 +++++++++++++++++++++++++++------- ndbioimage/readers/seqread.py | 4 +- pyproject.toml | 2 +- 4 files changed, 122 insertions(+), 40 deletions(-) mode change 100755 => 100644 ndbioimage/readers/cziread.py diff --git a/ndbioimage/__init__.py b/ndbioimage/__init__.py index c0f31de..d1f5a37 100755 --- a/ndbioimage/__init__.py +++ b/ndbioimage/__init__.py @@ -87,6 +87,7 @@ def find(obj: Sequence[Any], **kwargs: Any) -> Any: return item except AttributeError: pass + return None R = TypeVar("R") @@ -154,23 +155,26 @@ class OmeCache(DequeDict): def __reduce__(self) -> tuple[type, tuple]: return self.__class__, () - def __getitem__(self, path: Path | str | tuple) -> OME: + def __getitem__(self, path_and_series: tuple[Path | str | tuple, int]) -> OME: + path, series = path_and_series if isinstance(path, tuple): - return super().__getitem__(path) + return super().__getitem__((path, series)) else: - return super().__getitem__(self.path_and_lstat(path)) + return super().__getitem__((self.path_and_lstat(path), series)) - def __setitem__(self, path: Path | str | tuple, value: OME) -> None: + def __setitem__(self, path_and_series: tuple[Path | str | tuple, int], value: OME) -> None: + path, series = path_and_series if isinstance(path, tuple): - super().__setitem__(path, value) + super().__setitem__((path, series), value) else: - super().__setitem__(self.path_and_lstat(path), value) + super().__setitem__((self.path_and_lstat(path), series), value) - def __contains__(self, path: Path | str | tuple) -> bool: + def __contains__(self, path_and_series: tuple[Path | str | tuple, int]) -> bool: + path, series = path_and_series if isinstance(path, tuple): - return super().__contains__(path) + return super().__contains__((path, series)) else: - return super().__contains__(self.path_and_lstat(path)) + return super().__contains__((self.path_and_lstat(path), series)) @staticmethod def path_and_lstat(path: str | Path) -> tuple[Path, Optional[os.stat_result], Optional[os.stat_result]]: @@ -1025,7 +1029,7 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): return [self.get_channel(c) for c in czt[0]], *czt[1:3] # type: ignore @staticmethod - def bioformats_ome(path: [str, Path]) -> OME: + def bioformats_ome(path: str | Path) -> OME: """Use java BioFormats to make an ome metadata structure.""" with multiprocessing.get_context("spawn").Pool(1) as pool: return pool.map(bioformats_ome, (path,))[0] @@ -1057,10 +1061,11 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): return ome @staticmethod - def read_ome(path: [str, Path]) -> Optional[OME]: + def read_ome(path: str | Path) -> Optional[OME]: path = Path(path) if path.with_suffix(".ome.xml").exists(): return OME.from_xml(path.with_suffix(".ome.xml")) + return None def get_ome(self) -> OME: """overload this""" @@ -1069,12 +1074,12 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): @cached_property def ome(self) -> OME: cache = OmeCache() - if self.path not in cache: + if (self.path, self.series) not in cache: ome = self.read_ome(self.path) if ome is None: ome = self.get_ome() - cache[self.path] = self.fix_ome(ome) - return cache[self.path] + cache[self.path, self.series] = self.fix_ome(ome) + return cache[self.path, self.series] def is_noise(self, volume: ArrayLike = None) -> bool: """True if volume only has noise""" diff --git a/ndbioimage/readers/cziread.py b/ndbioimage/readers/cziread.py old mode 100755 new mode 100644 index cb4dbdd..962a735 --- a/ndbioimage/readers/cziread.py +++ b/ndbioimage/readers/cziread.py @@ -14,7 +14,7 @@ from lxml import etree from ome_types import OME, model from tifffile import repeat_nd -from .. import AbstractReader +from .. import AbstractReader, ureg try: # TODO: use zoom from imagecodecs implementation when available @@ -152,6 +152,14 @@ def data(self, raw: bool = False, resize: bool = True, order: int = 0) -> np.nda czifile.czifile.SubBlockSegment.data = data +def xml_walk(tree, elements): + element, *elements = elements + if elements: + return [j for i in tree.findall(element) for j in xml_walk(i, elements)] + else: + return tree.findall(element) + + class Reader(AbstractReader, ABC): priority = 0 do_not_pickle = "reader", "filedict" @@ -163,16 +171,49 @@ class Reader(AbstractReader, ABC): def open(self) -> None: self.reader = czifile.CziFile(self.path) filedict = {} + syx = set() + si = self.reader.axes.index("S") if "S" in self.reader.axes else None + ci = self.reader.axes.index("C") if "C" in self.reader.axes else None + zi = self.reader.axes.index("Z") if "Z" in self.reader.axes else None + ti = self.reader.axes.index("T") if "T" in self.reader.axes else None + yi = self.reader.axes.index("Y") if "Y" in self.reader.axes else None + xi = self.reader.axes.index("X") if "X" in self.reader.axes else None + for directory_entry in self.reader.filtered_subblock_directory: idx = self.get_index(directory_entry, self.reader.start) - if "S" not in self.reader.axes or self.series in range(*idx[self.reader.axes.index("S")]): - for c in range(*idx[self.reader.axes.index("C")]): - for z in range(*idx[self.reader.axes.index("Z")]): - for t in range(*idx[self.reader.axes.index("T")]): - if (c, z, t) in filedict: - filedict[c, z, t].append(directory_entry) - else: - filedict[c, z, t] = [directory_entry] + syx.add((0 if si is None else idx[si], idx[yi], idx[xi])) + + if self.tiles != (1, 1): + assert len({s for s, *_ in list(syx)}) == 1, "multiple tiled series not supported" + x, y = np.array(list(syx))[:, 1:, 0].T + x = np.unique(x) + y = np.unique(y) + _, bx = np.histogram(x, self.tiles[0]) + _, by = np.histogram(y, self.tiles[1]) + by, bx = list(product([(i, j) for i, j in zip(by, by[1:])], [(i, j) for i, j in zip(bx, bx[1:])]))[ + self.series + ] + for directory_entry in self.reader.filtered_subblock_directory: + idx = self.get_index(directory_entry, self.reader.start) + if bx[0] <= idx[xi][0] <= bx[1] and by[0] <= idx[yi][0] <= by[1]: + for cj in (0,) if ci is None else range(*idx[ci]): + for zj in (0,) if zi is None else range(*idx[zi]): + for tj in (0,) if ti is None else range(*idx[ti]): + if (cj, zj, tj) in filedict: + filedict[cj, zj, tj].append(directory_entry) + else: + filedict[cj, zj, tj] = [directory_entry] + else: + for directory_entry in self.reader.filtered_subblock_directory: + idx = self.get_index(directory_entry, self.reader.start) + if si is None or self.series == idx[si][0]: + for cj in (0,) if ci is None else range(*idx[ci]): + for zj in (0,) if zi is None else range(*idx[zi]): + for tj in (0,) if ti is None else range(*idx[ti]): + if (cj, zj, tj) in filedict: + filedict[cj, zj, tj].append(directory_entry) + else: + filedict[cj, zj, tj] = [directory_entry] if len(filedict) == 0: raise FileNotFoundError(f"Series {self.series} not found in {self.path}.") self.filedict = filedict # noqa @@ -187,28 +228,51 @@ class Reader(AbstractReader, ABC): f = np.zeros(self.base_shape["yx"], self.dtype) if (c, z, t) in self.filedict: directory_entries = self.filedict[c, z, t] - x_min = min([f.start[f.axes.index("X")] for f in directory_entries]) - y_min = min([f.start[f.axes.index("Y")] for f in directory_entries]) - xy_min = {"X": x_min, "Y": y_min} + start = np.min([directory_entry.start for directory_entry in directory_entries], 0) for directory_entry in directory_entries: subblock = directory_entry.data_segment() tile = subblock.data(resize=True, order=0) - axes_min = [xy_min.get(ax, 0) for ax in directory_entry.axes] - index = [ - slice(i - j - m, i - j + k) - for i, j, k, m in zip(directory_entry.start, self.reader.start, tile.shape, axes_min) - ] + index = [slice(i - j, i - j + k) for i, j, k in zip(directory_entry.start, start, tile.shape)] index = tuple(index[self.reader.axes.index(i)] for i in "YX") - try: - f[index] = tile.squeeze() - except ValueError: - f = tile.squeeze() + f[index] = tile.squeeze() return f @staticmethod def get_index(directory_entry: czifile.DirectoryEntryDV, start: tuple[int]) -> list[tuple[int, int]]: return [(i - j, i - j + k) for i, j, k in zip(directory_entry.start, start, directory_entry.shape)] + @cached_property + def tiles(self): + columns = 1 + rows = 1 + xml = self.reader.metadata() + tree = etree.fromstring(xml) + tile_regions = xml_walk( + tree, + ( + "Metadata", + "Experiment", + "ExperimentBlocks", + "AcquisitionBlock", + "SubDimensionSetups", + "RegionsSetup", + "SampleHolder", + "TileRegions", + "TileRegion", + ), + ) + for tile_region in tile_regions: + used = tile_region.find("IsUsedForAcquisition") + if used is not None and used.text.lower() == "true": + c = tile_region.find("Columns") + if c is not None: + columns = int(c.text) + r = tile_region.find("Rows") + if r is not None: + rows = int(r.text) + break + return columns, rows + class OmeParse: size_x: int @@ -432,7 +496,7 @@ class OmeParse: self.size_x = x_max - x_min self.size_y = y_max - y_min self.size_c, self.size_z, self.size_t = ( - self.reader.shape[self.reader.axes.index(directory_entry)] for directory_entry in "CZT" + self.reader.shape[self.reader.axes.index(axis)] if axis in self.reader.axes else 1 for axis in "CZT" ) image = self.information.find("Image") pixel_type = self.text(image.find("PixelType"), "Gray16") @@ -498,6 +562,8 @@ class OmeParse: except AttributeError: center_position = [0, 0] return center_position[0], center_position[1], None + else: + raise NotImplementedError(f"unknown czi version: {self.version}") @cached_property def channels_im(self) -> dict: @@ -642,7 +708,16 @@ class OmeParse: delta_ts = dt * np.arange(len(delta_ts)) warnings.warn(f"delta_t is inconsistent, using median value: {dt}") + pxsize_x = self.ome.images[0].pixels.physical_size_x_quantity.to(ureg.um).magnitude + pxsize_y = self.ome.images[0].pixels.physical_size_y_quantity.to(ureg.um).magnitude + for t, z, c in product(range(self.size_t), range(self.size_z), range(self.size_c)): + x_min = ( + min([f.start[f.axes.index("X")] for f in self.filedict[c, z, t]]) if (c, z, t) in self.filedict else 0 + ) + y_min = ( + min([f.start[f.axes.index("Y")] for f in self.filedict[c, z, t]]) if (c, z, t) in self.filedict else 0 + ) self.ome.images[0].pixels.planes.append( model.Plane( the_c=c, @@ -650,9 +725,9 @@ class OmeParse: the_t=t, delta_t=delta_ts[t], exposure_time=exposure_times[min(c, len(exposure_times) - 1)] if len(exposure_times) > 0 else None, - position_x=self.positions[0], + position_x=self.positions[0] + x_min * pxsize_x, position_x_unit=self.um, - position_y=self.positions[1], + position_y=self.positions[1] + y_min * pxsize_y, position_y_unit=self.um, position_z=self.positions[2], position_z_unit=self.um, diff --git a/ndbioimage/readers/seqread.py b/ndbioimage/readers/seqread.py index 5fd83ed..da05430 100644 --- a/ndbioimage/readers/seqread.py +++ b/ndbioimage/readers/seqread.py @@ -114,7 +114,9 @@ class Reader(AbstractReader, ABC): type=pixel_type, physical_size_x=pxsize, physical_size_y=pxsize, - physical_size_z=metadata["Info"]["Summary"]["z-step_um"] if "Summary" in metadata["Info"] else None, + physical_size_z=metadata["Info"]["Summary"]["z-step_um"] + if "Summary" in metadata["Info"] + else None, ), objective_settings=model.ObjectiveSettings(id="Objective:0"), ) diff --git a/pyproject.toml b/pyproject.toml index 32d2d48..9c6a005 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "ndbioimage" -version = "2026.3.2" +version = "2026.3.3" description = "Bio image reading, metadata and some affine registration." authors = [ { name = "W. Pomp", email = "w.pomp@nki.nl" }