diff --git a/.gitignore b/.gitignore index d1f43ac..3bc6b6f 100644 --- a/.gitignore +++ b/.gitignore @@ -71,3 +71,5 @@ docs/_build/ # Pyenv .python-version + +/tests/files/* diff --git a/Cargo.toml b/Cargo.toml index 9d73c95..65d26db 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "ndbioimage" -version = "2025.2.0" +version = "2025.2.1" edition = "2021" authors = ["Wim Pomp "] license = "GPL-3.0-or-later" @@ -13,7 +13,7 @@ exclude = ["/tests"] # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html [lib] name = "ndbioimage" -crate-type = ["cdylib"] +crate-type = ["cdylib", "rlib"] [dependencies] anyhow = "1.0.95" @@ -31,9 +31,9 @@ optional = true rayon = "1.10.0" [build-dependencies] -anyhow = { version = "1.0.95"} -j4rs = { version = "0.22", features = [] } -retry = { version = "2.0.0"} +anyhow = "1.0.95" +j4rs = "0.22.0" +retry = "2.0.0" [features] python = ["dep:pyo3", "dep:numpy"] diff --git a/build.rs b/build.rs index c67c341..e1495e7 100644 --- a/build.rs +++ b/build.rs @@ -1,21 +1,58 @@ -// copied from https://github.com/AzHicham/bioformats-rs - +#[cfg(not(feature = "python"))] use j4rs::{errors::J4RsError, JvmBuilder, MavenArtifact, MavenArtifactRepo, MavenSettings}; + +#[cfg(not(feature = "python"))] use retry::{delay, delay::Exponential, retry}; +#[cfg(feature = "python")] +use j4rs::Jvm; + fn main() -> anyhow::Result<()> { println!("cargo:rerun-if-changed=build.rs"); - if std::env::var("DOCS_RS").is_ok() { - Ok(()) - } else { - Ok(retry( + #[cfg(not(feature = "python"))] + if std::env::var("DOCS_RS").is_err() { + retry( Exponential::from_millis(1000).map(delay::jitter).take(4), deploy_java_artifacts, - )?) + )? } + + #[cfg(feature = "python")] + { + let py_src_path = std::env::current_dir()?.join("py").join("ndbioimage"); + let py_jassets_path = py_src_path.join("jassets"); + let py_deps_path = py_src_path.join("deps"); + if py_jassets_path.exists() { + std::fs::remove_dir_all(&py_jassets_path)?; + } + if py_deps_path.exists() { + std::fs::remove_dir_all(&py_deps_path)?; + } + + Jvm::copy_j4rs_libs_under(py_src_path.to_str().unwrap())?; + + // rename else maturin will ignore them + for file in std::fs::read_dir(&py_deps_path)? { + let f = file?.path().to_str().unwrap().to_string(); + if !f.ends_with("_") { + std::fs::rename(&f, std::format!("{f}_"))?; + } + } + + // remove so we don't include too much accidentally + for file in std::fs::read_dir(&py_jassets_path)? { + let f = file?.path(); + if !f.file_name().unwrap().to_str().unwrap().starts_with("j4rs") { + std::fs::remove_file(&f)?; + } + } + } + + Ok(()) } +#[cfg(not(feature = "python"))] fn deploy_java_artifacts() -> Result<(), J4RsError> { let jvm = JvmBuilder::new() .with_maven_settings(MavenSettings::new(vec![MavenArtifactRepo::from( @@ -23,10 +60,10 @@ fn deploy_java_artifacts() -> Result<(), J4RsError> { )])) .build()?; - jvm.deploy_artifact(&MavenArtifact::from("ome:bioformats_package:8.0.1"))?; + jvm.deploy_artifact(&MavenArtifact::from("ome:bioformats_package:8.1.0"))?; #[cfg(feature = "gpl-formats")] - jvm.deploy_artifact(&MavenArtifact::from("ome:formats-gpl:8.0.1"))?; + jvm.deploy_artifact(&MavenArtifact::from("ome:formats-gpl:8.1.0"))?; Ok(()) } diff --git a/py/ndbioimage/__init__.py b/py/ndbioimage/__init__.py new file mode 100755 index 0000000..8ec622d --- /dev/null +++ b/py/ndbioimage/__init__.py @@ -0,0 +1,1357 @@ +from __future__ import annotations + +import multiprocessing +import os +import re +import warnings +from abc import ABC, ABCMeta, abstractmethod +from argparse import ArgumentParser +from collections import OrderedDict +from datetime import datetime +from functools import cached_property, wraps +from importlib.metadata import version +from itertools import product +from operator import truediv +from pathlib import Path +from typing import Any, Callable, Generator, Iterable, Optional, Sequence, TypeVar + +import numpy as np +import yaml +from numpy.typing import ArrayLike, DTypeLike +from ome_types import OME, from_xml, model, ureg +from pint import set_application_registry +from tiffwrite import FrameInfo, IJTiffParallel +from tqdm.auto import tqdm + +from .transforms import Transform, Transforms # noqa: F401 + +try: + __version__ = version(Path(__file__).parent.name) +except Exception: # noqa + __version__ = 'unknown' + +try: + with open(Path(__file__).parent.parent / '.git' / 'HEAD') as g: + head = g.read().split(':')[1].strip() + with open(Path(__file__).parent.parent / '.git' / head) as h: + __git_commit_hash__ = h.read().rstrip('\n') +except Exception: # noqa + __git_commit_hash__ = 'unknown' + +ureg.formatter.default_format = '~P' +set_application_registry(ureg) +warnings.filterwarnings('ignore', 'Reference to unknown ID') +Number = int | float | np.integer | np.floating + + +class ReaderNotFoundError(Exception): + pass + + +class TransformTiff(IJTiffParallel): + """ transform frames in a parallel process to speed up saving """ + def __init__(self, image: Imread, *args: Any, **kwargs: Any) -> None: + self.image = image + super().__init__(*args, **kwargs) + + def parallel(self, frame: tuple[int, int, int]) -> tuple[FrameInfo]: + return (np.asarray(self.image(*frame)), 0, 0, 0), + + +class DequeDict(OrderedDict): + def __init__(self, maxlen: int = None, *args: Any, **kwargs: Any) -> None: + self.maxlen = maxlen + super().__init__(*args, **kwargs) + + def __setitem__(self, *args: Any, **kwargs: Any) -> None: + super().__setitem__(*args, **kwargs) + self.truncate() + + def truncate(self) -> None: + if self.maxlen is not None: + while len(self) > self.maxlen: + self.popitem(False) + + def update(self, *args: Any, **kwargs: Any) -> None: + super().update(*args, **kwargs) # type: ignore + self.truncate() + + +def find(obj: Sequence[Any], **kwargs: Any) -> Any: + for item in obj: + try: + if all([getattr(item, key) == value for key, value in kwargs.items()]): + return item + except AttributeError: + pass + + +R = TypeVar('R') + + +def try_default(fun: Callable[..., R], default: Any, *args: Any, **kwargs: Any) -> R: + try: + return fun(*args, **kwargs) + except Exception: # noqa + return default + + +def bioformats_ome(path: str | Path) -> OME: + from .rs import Reader + reader = Reader(str(path), 0) + return from_xml(reader.get_ome_xml(), parser='lxml') + + +class Shape(tuple): + def __new__(cls, shape: Sequence[int] | Shape, axes: str = 'yxczt') -> Shape: + if isinstance(shape, Shape): + axes = shape.axes # type: ignore + new = super().__new__(cls, shape) + new.axes = axes.lower() + return new # type: ignore + + def __getitem__(self, n: int | str) -> int | tuple[int]: + if isinstance(n, str): + if len(n) == 1: + return self[self.axes.find(n.lower())] if n.lower() in self.axes else 1 + else: + return tuple(self[i] for i in n) # type: ignore + return super().__getitem__(n) + + @cached_property + def yxczt(self) -> tuple[int, int, int, int, int]: + return tuple(self[i] for i in 'yxczt') # type: ignore + + +class OmeCache(DequeDict): + """ prevent (potentially expensive) rereading of ome data by caching """ + + instance = None + + def __new__(cls) -> OmeCache: + if cls.instance is None: + cls.instance = super().__new__(cls) + return cls.instance + + def __init__(self) -> None: + super().__init__(64) + + def __reduce__(self) -> tuple[type, tuple]: + return self.__class__, () + + def __getitem__(self, path: Path | str | tuple) -> OME: + if isinstance(path, tuple): + return super().__getitem__(path) + else: + return super().__getitem__(self.path_and_lstat(path)) + + def __setitem__(self, path: Path | str | tuple, value: OME) -> None: + if isinstance(path, tuple): + super().__setitem__(path, value) + else: + super().__setitem__(self.path_and_lstat(path), value) + + def __contains__(self, path: Path | str | tuple) -> bool: + if isinstance(path, tuple): + return super().__contains__(path) + else: + return super().__contains__(self.path_and_lstat(path)) + + @staticmethod + def path_and_lstat(path: str | Path) -> tuple[Path, Optional[os.stat_result], Optional[os.stat_result]]: + path = Path(path) + return (path, (path.lstat() if path.exists() else None), + (path.with_suffix('.ome.xml').lstat() if path.with_suffix('.ome.xml').exists() else None)) + + +def get_positions(path: str | Path) -> Optional[list[int]]: + subclass = AbstractReader.get_subclass(path) + return subclass.get_positions(AbstractReader.split_path_series(path)[0]) + + +class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): + """ class to read image files, while taking good care of important metadata, + currently optimized for .czi files, but can open anything that bioformats can handle + path: path to the image file + optional: + axes: order of axes, default: cztyx, but omitting any axes with lenght 1 + dtype: datatype to be used when returning frames + + modify images on the fly with a decorator function: + define a function which takes an instance of this object, one image frame, + and the coordinates c, z, t as arguments, and one image frame as return + >> imread.frame_decorator = fun + then use imread as usually + + Examples: + >> im = Imread('/path/to/file.image', axes='czt) + >> im + << shows summary + >> im.shape + << (15, 26, 1000, 1000) + >> im.axes + << 'ztyx' + >> plt.imshow(im[1, 0]) + << plots frame at position z=1, t=0 (python type indexing) + >> plt.imshow(im[:, 0].max('z')) + << plots max-z projection at t=0 + >> im.pxsize_um + << 0.09708737864077668 image-plane pixel size in um + >> im.laserwavelengths + << [642, 488] + >> im.laserpowers + << [0.02, 0.0005] in % + + See __init__ and other functions for more ideas. + + Subclassing: + Subclass AbstractReader to add more file types. A subclass should always have at least the following + methods: + staticmethod _can_open(path): returns True when the subclass can open the image in path + __frame__(self, c, z, t): this should return a single frame at channel c, slice z and time t + optional open(self): code to be run during initialization, e.g. to open a file handle + optional close(self): close the file in a proper way + optional class field priority: subclasses with lower priority will be tried first, default = 99 + optional get_ome(self) -> OME: return an OME structure with metadata, + if not present bioformats will be used to generate an OME + Any other method can be overridden as needed + """ + + isclosed: Optional[bool] + channel_names: Optional[list[str]] + series: Optional[int] + pxsize_um: Optional[float] + deltaz_um: Optional[float] + exposuretime_s: Optional[list[float]] + timeinterval: Optional[float] + binning: Optional[list[int]] + laserwavelengths: Optional[list[tuple[float]]] + laserpowers: Optional[list[tuple[float]]] + objective: Optional[model.Objective] + magnification: Optional[float] + tubelens: Optional[model.Objective] + filter: Optional[str] + powermode: Optional[str] + collimator: Optional[str] + tirfangle: Optional[list[float]] + gain: Optional[list[float]] + pcf: Optional[list[float]] + path: Path + __frame__: Callable[[int, int, int], np.ndarray] + + @staticmethod + def get_subclass(path: Path | str | Any): + if len(AbstractReader.__subclasses__()) == 0: + raise Exception('Restart python kernel please!') + path, _ = AbstractReader.split_path_series(path) + for subclass in sorted(AbstractReader.__subclasses__(), key=lambda subclass_: subclass_.priority): + if subclass._can_open(path): # noqa + do_not_pickle = (AbstractReader.do_not_pickle,) if isinstance(AbstractReader.do_not_pickle, str) \ + else AbstractReader.do_not_pickle + subclass_do_not_pickle = (subclass.do_not_pickle,) if isinstance(subclass.do_not_pickle, str) \ + else subclass.do_not_pickle if hasattr(subclass, 'do_not_pickle') else () + subclass.do_not_pickle = set(do_not_pickle).union(set(subclass_do_not_pickle)) + return subclass + raise ReaderNotFoundError(f'No reader found for {path}.') + + + def __new__(cls, path: Path | str | Imread | Any = None, dtype: DTypeLike = None, axes: str = None) -> Imread: + if cls is not Imread: + return super().__new__(cls) + if isinstance(path, Imread): + return path + subclass = cls.get_subclass(path) + do_not_pickle = (AbstractReader.do_not_pickle,) if isinstance(AbstractReader.do_not_pickle, str) \ + else AbstractReader.do_not_pickle + subclass_do_not_pickle = (subclass.do_not_pickle,) if isinstance(subclass.do_not_pickle, str) \ + else subclass.do_not_pickle if hasattr(subclass, 'do_not_pickle') else () + subclass.do_not_pickle = set(do_not_pickle).union(set(subclass_do_not_pickle)) + return super().__new__(subclass) + + def __init__(self, *args: Any, **kwargs: Any): + def parse(base: Imread = None, # noqa + slice: tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray] = None, # noqa + shape: tuple[int, ...] = (0, 0, 0, 0, 0), # noqa + dtype: DTypeLike = None, # noqa + frame_decorator: Callable[[Imread, np.ndarray, int, int, int], np.ndarray] = None # noqa + ) -> tuple[Any, ...]: + return base, slice, shape, dtype, frame_decorator + + base, slice, shape, dtype, frame_decorator = parse(*args, **kwargs) # noqa + self.base = base or self + self.slice = slice + self._shape = Shape(shape) + self.dtype = dtype + self.frame_decorator = frame_decorator + self.transform = Transforms() + self.flags = dict(C_CONTIGUOUS=False, F_CONTIGUOUS=False, OWNDATA=False, WRITEABLE=False, + ALIGNED=False, WRITEBACKIFCOPY=False, UPDATEIFCOPY=False) + + def __call__(self, c: int = None, z: int = None, t: int = None, x: int = None, y: int = None) -> np.ndarray: + """ same as im[] but allowing keyword axes, but slices need to made with slice() or np.s_ """ + return self[{k: slice(v) if v is None else v for k, v in dict(c=c, z=z, t=t, x=x, y=y).items()}] + + def __copy__(self) -> Imread: + return self.copy() + + def __contains__(self, item: Number) -> bool: + def unique_yield(a: Iterable[Any], b: Iterable[Any]) -> Generator[Any, None, None]: + for k in a: + yield k + for k in b: + if k not in a: + yield k + + for idx in unique_yield([key[:3] for key in self.cache.keys()], + product(range(self.shape['c']), range(self.shape['z']), range(self.shape['t']))): + yxczt = (slice(None), slice(None)) + idx + in_idx = tuple(yxczt['yxczt'.find(i)] for i in self.axes) + if item in np.asarray(self[in_idx]): + return True + return False + + def __enter__(self) -> Imread: + return self + + def __exit__(self, *args: Any, **kwargs: Any) -> None: + if not self.isclosed: + self.isclosed = True + if hasattr(self, 'close'): + self.close() + + def __getitem__(self, n: int | Sequence[int] | Sequence[slice] | slice | type(Ellipsis) | + dict[str, int | Sequence[int] | Sequence[slice] | slice | type(Ellipsis)] + ) -> Number | Imread | np.ndarray: + """ slice like a numpy array but return an Imread instance """ + if self.isclosed: + raise OSError('file is closed') + if isinstance(n, (slice, Number)): # None = : + n = (n,) + elif isinstance(n, type(Ellipsis)): + n = (None,) * len(self.axes) + elif isinstance(n, dict): # allow im[dict(z=0)] etc. + n = [n.get(i, slice(None)) for i in self.axes] + n = list(n) + + # deal with ... + ell = [i for i, e in enumerate(n) if isinstance(e, type(Ellipsis))] + if len(ell) > 1: + raise IndexError('an index can only have a single ellipsis (...)') + if len(ell): + if len(n) > self.ndim: + n.remove(Ellipsis) + else: + n[ell[0]] = None + while len(n) < self.ndim: + n.insert(ell[0], None) + while len(n) < self.ndim: + n.append(None) + + axes_idx = [self.shape.axes.find(i) for i in 'yxczt'] + n = [n[j] if 0 <= j < len(n) else None for j in axes_idx] # reorder n + + new_slice = [] + for s, e in zip(self.slice, n): + if e is None: + new_slice.append(s) + else: + new_slice.append(s[e]) + + # TODO: check output dimensionality when requested shape in some dimension is 1 + if all([isinstance(s, Number) or a < 0 and s.size == 1 for s, a in zip(new_slice, axes_idx)]): + return self.block(*new_slice).item() + else: + new = View(self) + new.slice = new_slice + new._shape = Shape([1 if isinstance(s, Number) else len(s) for s in new_slice]) + new.axes = ''.join(j for j in self.axes if j in [i for i, s in zip('yxczt', new_slice) + if not isinstance(s, Number)]) + return new + + def __getstate__(self) -> dict[str: Any]: + return ({key: value for key, value in self.__dict__.items() if key not in self.do_not_pickle} | + {'cache_size': self.cache.maxlen}) + + def __len__(self) -> int: + return self.shape[0] + + def __repr__(self) -> str: + return self.summary + + def __setstate__(self, state: dict[str, Any]) -> None: + """ What happens during unpickling """ + self.__dict__.update({key: value for key, value in state.items() if key != 'cache_size'}) + if isinstance(self, AbstractReader): + self.open() + self.cache = DequeDict(state.get('cache_size', 16)) + + def __str__(self) -> str: + return str(self.path) + + # @property + # def __array_interface__(self) -> dict[str, Any]: + # return dict(shape=tuple(self.shape), typestr=self.dtype.str, version=3, data=self.tobytes()) + + def tobytes(self) -> bytes: + return self.flatten().tobytes() + + def __array__(self, dtype: DTypeLike = None, copy: bool = None) -> np.ndarray: + if copy is False: + raise ValueError("`copy=False` isn't supported. A copy is always created.") + block = self.block(*self.slice) + axes_idx = [self.shape.axes.find(i) for i in 'yxczt'] + axes_squeeze = tuple({i for i, j in enumerate(axes_idx) if j == -1}.union( + {i for i, j in enumerate(self.slice) if isinstance(j, Number)})) + block = block.squeeze(axes_squeeze) + if dtype is not None: + block = block.astype(dtype) + if block.ndim == 0: + return block.item() + axes = ''.join(j for i, j in enumerate('yxczt') if i not in axes_squeeze) + return block.transpose([axes.find(i) for i in self.shape.axes if i in axes]) + + def __array_arg_fun__(self, fun: Callable[[ArrayLike, Optional[int]], Number | np.ndarray], + axis: int | str = None, out: np.ndarray = None) -> Number | np.ndarray: + """ frame-wise application of np.argmin and np.argmax """ + if axis is None: + value = arg = None + for idx in product(range(self.shape['c']), range(self.shape['z']), range(self.shape['t'])): + yxczt = (slice(None), slice(None)) + idx + in_idx = tuple(yxczt['yxczt'.find(i)] for i in self.axes) + new = np.asarray(self[in_idx]) + new_arg = np.unravel_index(fun(new), new.shape) # type: ignore + new_value = new[new_arg] + if value is None: + arg = new_arg + idx + value = new_value + else: + i = fun((value, new_value)) # type: ignore + arg = (arg, new_arg + idx)[i] + value = (value, new_value)[i] + axes = ''.join(i for i in self.axes if i in 'yx') + 'czt' + arg = np.ravel_multi_index([arg[axes.find(i)] for i in self.axes], self.shape) + if out is None: + return arg + else: + out.itemset(arg) + return out + else: + if isinstance(axis, str): + axis_str, axis_idx = axis, self.axes.index(axis) + else: + axis_str, axis_idx = self.axes[axis], axis + if axis_str not in self.axes: + raise IndexError(f'Axis {axis_str} not in {self.axes}.') + out_shape = list(self.shape) + out_axes = list(self.axes) + out_shape.pop(axis_idx) + out_axes.pop(axis_idx) + if out is None: + out = np.zeros(out_shape, int) + if axis_str in 'yx': + for idx in product(range(self.shape['c']), range(self.shape['z']), range(self.shape['t'])): + yxczt = (slice(None), slice(None)) + idx + out_idx = tuple(yxczt['yxczt'.find(i)] for i in out_axes) + in_idx = tuple(yxczt['yxczt'.find(i)] for i in self.axes) + new = self[in_idx] + out[out_idx] = fun(np.asarray(new), new.axes.find(axis_str)) + else: + value = np.zeros(out.shape, self.dtype) + for idx in product(range(self.shape['c']), range(self.shape['z']), range(self.shape['t'])): + yxczt = (slice(None), slice(None)) + idx + out_idx = tuple(yxczt['yxczt'.find(i)] for i in out_axes) + in_idx = tuple(yxczt['yxczt'.find(i)] for i in self.axes) + new_value = self[in_idx] + new_arg = np.full_like(new_value, idx['czt'.find(axis_str)]) + if idx['czt'.find(axis_str)] == 0: + value[out_idx] = new_value + out[out_idx] = new_arg + else: + old_value = value[out_idx] + i = fun((old_value, new_value), 0) + value[out_idx] = np.where(i, new_value, old_value) + out[out_idx] = np.where(i, new_arg, out[out_idx]) + return out + + def __array_fun__(self, funs: Sequence[Callable[[ArrayLike], Number | np.ndarray]], axis: int | str = None, + dtype: DTypeLike = None, out: np.ndarray = None, keepdims: bool = False, + initials: list[Number | np.ndarray] = None, where: bool | int | np.ndarray = True, + ffuns: Sequence[Callable[[ArrayLike], np.ndarray]] = None, + cfun: Callable[..., np.ndarray] = None) -> Number | np.ndarray: + """ frame-wise application of np.min, np.max, np.sum, np.mean and their nan equivalents """ + p = re.compile(r'\d') + dtype = self.dtype if dtype is None else np.dtype(dtype) + if initials is None: + initials = [None for _ in funs] + if ffuns is None: + ffuns = [None for _ in funs] + + def ffun_(frame: ArrayLike) -> np.ndarray: + return np.asarray(frame) + ffuns = [ffun_ if ffun is None else ffun for ffun in ffuns] + if cfun is None: + def cfun(*res): # noqa + return res[0] + + # TODO: smarter transforms + if axis is None: + for idx in product(range(self.shape['c']), range(self.shape['z']), range(self.shape['t'])): + yxczt = (slice(None), slice(None)) + idx + in_idx = tuple(yxczt['yxczt'.find(i)] for i in self.axes) + w = where if where is None or isinstance(where, bool) else where[in_idx] + initials = [fun(np.asarray(ffun(self[in_idx])), initial=initial, where=w) # type: ignore + for fun, ffun, initial in zip(funs, ffuns, initials)] + res = cfun(*initials) + res = (np.round(res) if dtype.kind in 'ui' else res).astype(p.sub('', dtype.name)) + if keepdims: + res = np.array(res, dtype, ndmin=self.ndim) + if out is None: + return res + else: + out.itemset(res) + return out + else: + if isinstance(axis, str): + axis_str, axis_idx = axis, self.axes.index(axis) + else: + axis_idx = axis % self.ndim + axis_str = self.axes[axis_idx] + if axis_str not in self.axes: + raise IndexError(f'Axis {axis_str} not in {self.axes}.') + out_shape = list(self.shape) + out_axes = list(self.axes) + if not keepdims: + out_shape.pop(axis_idx) + out_axes.pop(axis_idx) + if out is None: + out = np.zeros(out_shape, dtype) + if axis_str in 'yx': + yx = 'yx' if self.axes.find('x') > self.axes.find('y') else 'yx' + frame_ax = yx.find(axis_str) + for idx in product(range(self.shape['c']), range(self.shape['z']), range(self.shape['t'])): + yxczt = (slice(None), slice(None)) + idx + out_idx = tuple(yxczt['yxczt'.find(i)] for i in out_axes) + in_idx = tuple(yxczt['yxczt'.find(i)] for i in self.axes) + w = where if where is None or isinstance(where, bool) else where[in_idx] + res = cfun(*[fun(ffun(self[in_idx]), frame_ax, initial=initial, where=w) # type: ignore + for fun, ffun, initial in zip(funs, ffuns, initials)]) + out[out_idx] = (np.round(res) if out.dtype.kind in 'ui' else res).astype(p.sub('', dtype.name)) + else: + tmps = [np.zeros(out_shape) for _ in ffuns] + for idx in product(range(self.shape['c']), range(self.shape['z']), range(self.shape['t'])): + yxczt = (slice(None), slice(None)) + idx + out_idx = tuple(yxczt['yxczt'.find(i)] for i in out_axes) + in_idx = tuple(yxczt['yxczt'.find(i)] for i in self.axes) + + if idx['czt'.find(axis_str)] == 0: + w = where if where is None or isinstance(where, bool) else (where[in_idx],) + for tmp, fun, ffun, initial in zip(tmps, funs, ffuns, initials): + tmp[out_idx] = fun((ffun(self[in_idx]),), 0, initial=initial, where=w) # type: ignore + else: + w = where if where is None or isinstance(where, bool) else \ + (np.ones_like(where[in_idx]), where[in_idx]) + for tmp, fun, ffun in zip(tmps, funs, ffuns): + tmp[out_idx] = fun((tmp[out_idx], ffun(self[in_idx])), 0, where=w) # type: ignore + out[...] = (np.round(cfun(*tmps)) if out.dtype.kind in 'ui' else + cfun(*tmps)).astype(p.sub('', dtype.name)) + return out + + @property + def axes(self) -> str: + return self.shape.axes + + @axes.setter + def axes(self, value: str) -> None: + shape = self.shape[value] + if isinstance(shape, Number): + shape = (shape,) + self._shape = Shape(shape, value) + + @property + def dtype(self) -> np.dtype: + return self._dtype + + @dtype.setter + def dtype(self, value: DTypeLike) -> None: + self._dtype = np.dtype(value) + + @cached_property + def extrametadata(self) -> Optional[Any]: + if isinstance(self.path, Path): + if self.path.with_suffix('.pzl2').exists(): + pname = self.path.with_suffix('.pzl2') + elif self.path.with_suffix('.pzl').exists(): + pname = self.path.with_suffix('.pzl') + else: + return None + try: + return self.get_config(pname) + except Exception: # noqa + return None + return None + + @property + def ndim(self) -> int: + return len(self.shape) + + @property + def size(self) -> int: + return np.prod(self.shape) # type: ignore + + @property + def shape(self) -> Shape: + return self._shape + + @shape.setter + def shape(self, value: Shape | tuple[int, ...]) -> None: + if isinstance(value, Shape): + self._shape = value + else: + self._shape = Shape([value['yxczt'.find(i.lower())] for i in self.axes], self.axes) + + @property + def summary(self) -> str: + """ gives a helpful summary of the recorded experiment """ + s = [f'path/filename: {self.path}', + f'series/pos: {self.series}', + f"reader: {self.base.__class__.__module__.split('.')[-1]}"] + s.extend((f'dtype: {self.dtype}', + f'shape ({self.axes}):'.ljust(15) + f"{' x '.join(str(i) for i in self.shape)}")) + if self.pxsize_um: + s.append(f'pixel size: {1000 * self.pxsize_um:.2f} nm') + if self.zstack and self.deltaz_um: + s.append(f'z-interval: {1000 * self.deltaz_um:.2f} nm') + if self.exposuretime_s and not all(e is None for e in self.exposuretime_s): + s.append(f'exposuretime: {self.exposuretime_s[0]:.2f} s') + if self.timeseries and self.timeinterval: + s.append(f'time interval: {self.timeinterval:.3f} s') + if self.binning: + s.append('binning: {}x{}'.format(*self.binning)) + if self.laserwavelengths: + s.append('laser colors: ' + ' | '.join([' & '.join(len(w) * ('{:.0f}',)).format(*w) + for w in self.laserwavelengths]) + ' nm') + if self.laserpowers: + s.append('laser powers: ' + ' | '.join([' & '.join(len(p) * ('{:.3g}',)).format(*[100 * i for i in p]) + for p in self.laserpowers]) + ' %') + if self.objective and self.objective.model: + s.append(f'objective: {self.objective.model}') + if self.magnification: + s.append(f'magnification: {self.magnification}x') + if self.tubelens and self.tubelens.model: + s.append(f'tubelens: {self.tubelens.model}') + if self.filter: + s.append(f'filterset: {self.filter}') + if self.powermode: + s.append(f'powermode: {self.powermode}') + if self.collimator: + s.append('collimator: ' + (' {}' * len(self.collimator)).format(*self.collimator)) + if self.tirfangle: + s.append('TIRF angle: ' + (' {:.2f}°' * len(self.tirfangle)).format(*self.tirfangle)) + if self.gain: + s.append('gain: ' + (' {:.0f}' * len(self.gain)).format(*self.gain)) + if self.pcf: + s.append('pcf: ' + (' {:.2f}' * len(self.pcf)).format(*self.pcf)) + return '\n'.join(s) + + @property + def T(self) -> Imread: # noqa + return self.transpose() # type: ignore + + @cached_property + def timeseries(self) -> bool: + return self.shape['t'] > 1 + + @cached_property + def zstack(self) -> bool: + return self.shape['z'] > 1 + + @wraps(np.ndarray.argmax) + def argmax(self, *args, **kwargs): + return self.__array_arg_fun__(np.argmax, *args, **kwargs) + + @wraps(np.ndarray.argmin) + def argmin(self, *args, **kwargs): + return self.__array_arg_fun__(np.argmin, *args, **kwargs) + + @wraps(np.ndarray.max) + def max(self, axis=None, out=None, keepdims=False, initial=None, where=True, **_): + return self.__array_fun__([np.max], axis, None, out, keepdims, [initial], where) + + @wraps(np.ndarray.mean) + def mean(self, axis=None, dtype=None, out=None, keepdims=False, *, where=True, **_): + dtype = dtype or float + n = np.prod(self.shape) if axis is None else self.shape[axis] + + def sfun(frame: ArrayLike) -> np.ndarray: + return np.asarray(frame).astype(float) + + def cfun(res: np.ndarray) -> np.ndarray: + return res / n + + return self.__array_fun__([np.sum], axis, dtype, out, keepdims, None, where, [sfun], cfun) + + @wraps(np.ndarray.min) + def min(self, axis=None, out=None, keepdims=False, initial=None, where=True, **_): + return self.__array_fun__([np.min], axis, None, out, keepdims, [initial], where) + + @wraps(np.moveaxis) + def moveaxis(self, source, destination): + raise NotImplementedError('moveaxis is not implemented') + + @wraps(np.nanmax) + def nanmax(self, axis=None, out=None, keepdims=False, initial=None, where=True, **_): + return self.__array_fun__([np.nanmax], axis, None, out, keepdims, [initial], where) + + @wraps(np.nanmean) + def nanmean(self, axis=None, dtype=None, out=None, keepdims=False, *, where=True, **_): + dtype = dtype or float + + def sfun(frame): + return np.asarray(frame).astype(float) + + def nfun(frame): + return np.invert(np.isnan(frame)) + + return self.__array_fun__([np.nansum, np.sum], axis, dtype, out, keepdims, None, where, (sfun, nfun), truediv) + + @wraps(np.nanmin) + def nanmin(self, axis=None, out=None, keepdims=False, initial=None, where=True, **_): + return self.__array_fun__([np.nanmin], axis, None, out, keepdims, [initial], where) + + @wraps(np.nansum) + def nansum(self, axis=None, dtype=None, out=None, keepdims=False, initial=None, where=True, **_): + return self.__array_fun__([np.nansum], axis, dtype, out, keepdims, [initial], where) + + @wraps(np.nanstd) + def nanstd(self, axis=None, dtype=None, out=None, ddof=0, keepdims=None, *, where=None): + return self.nanvar(axis, dtype, out, ddof, keepdims, where=where, std=True) + + @wraps(np.nanvar) + def nanvar(self, axis=None, dtype=None, out=None, ddof=0, keepdims=None, *, where=True, std=False): + dtype = dtype or float + + def sfun(frame): + return np.asarray(frame).astype(float) + + def s2fun(frame): + return np.asarray(frame).astype(float) ** 2 + + def nfun(frame): + return np.invert(np.isnan(frame)) + + if std: + def cfun(s, s2, n): + return np.sqrt((s2 - s ** 2 / n) / (n - ddof)) + else: + def cfun(s, s2, n): + return (s2 - s ** 2 / n) / (n - ddof) + return self.__array_fun__([np.nansum, np.nansum, np.sum], axis, dtype, out, keepdims, None, where, + (sfun, s2fun, nfun), cfun) + + @wraps(np.ndarray.flatten) + def flatten(self, *args, **kwargs): + return np.asarray(self).flatten(*args, **kwargs) + + @wraps(np.ndarray.reshape) + def reshape(self, *args, **kwargs): + return np.asarray(self).reshape(*args, **kwargs) + + @wraps(np.ndarray.squeeze) + def squeeze(self, axes=None): + new = self.copy() + if axes is None: + axes = tuple(i for i, j in enumerate(new.shape) if j == 1) + elif isinstance(axes, Number): + axes = (axes,) + else: + axes = tuple(new.axes.find(ax) if isinstance(ax, str) else ax for ax in axes) + if any([new.shape[ax] != 1 for ax in axes]): + raise ValueError('cannot select an axis to squeeze out which has size not equal to one') + new.axes = ''.join(j for i, j in enumerate(new.axes) if i not in axes) + return new + + @wraps(np.ndarray.std) + def std(self, axis=None, dtype=None, out=None, ddof=0, keepdims=None, *, where=True): + return self.var(axis, dtype, out, ddof, keepdims, where=where, std=True) # type: ignore + + @wraps(np.ndarray.sum) + def sum(self, axis=None, dtype=None, out=None, keepdims=False, initial=None, where=True, **_): + return self.__array_fun__([np.sum], axis, dtype, out, keepdims, [initial], where) + + @wraps(np.ndarray.swapaxes) + def swapaxes(self, axis1, axis2): + new = self.copy() + axes = new.axes + if isinstance(axis1, str): + axis1 = axes.find(axis1) + if isinstance(axis2, str): + axis2 = axes.find(axis2) + axes[axis1], axes[axis2] = axes[axis2], axes[axis1] + new.axes = axes + return new + + @wraps(np.ndarray.transpose) + def transpose(self, *axes): + new = self.copy() + if not axes: + new.axes = new.axes[::-1] + else: + new.axes = ''.join(ax if isinstance(ax, str) else new.axes[ax] for ax in axes) + return new + + @wraps(np.ndarray.var) + def var(self, axis=None, dtype=None, out=None, ddof=0, keepdims=None, *, where=True, std=False): + dtype = dtype or float + n = np.prod(self.shape) if axis is None else self.shape[axis] + + def sfun(frame): + return np.asarray(frame).astype(float) + + def s2fun(frame): + return np.asarray(frame).astype(float) ** 2 + + if std: + def cfun(s, s2): + return np.sqrt((s2 - s ** 2 / n) / (n - ddof)) + else: + def cfun(s, s2): + return (s2 - s ** 2 / n) / (n - ddof) + return self.__array_fun__([np.sum, np.sum], axis, dtype, out, keepdims, None, where, (sfun, s2fun), cfun) + + def asarray(self) -> np.ndarray: + return self.__array__() + + @wraps(np.ndarray.astype) + def astype(self, dtype, *_, **__): + new = self.copy() + new.dtype = dtype + return new + + def block(self, y: int | Sequence[int] = None, x: int | Sequence[int] = None, + c: int | Sequence[int] = None, z: int | Sequence[int] = None, + t: int | Sequence[int] = None) -> np.ndarray: + """ returns 5D block of frames """ + y, x, c, z, t = (np.arange(self.shape[i]) if e is None else np.array(e, ndmin=1) + for i, e in zip('yxczt', (y, x, c, z, t))) + d = np.empty((len(y), len(x), len(c), len(z), len(t)), self.dtype) + for (ci, cj), (zi, zj), (ti, tj) in product(enumerate(c), enumerate(z), enumerate(t)): + d[:, :, ci, zi, ti] = self.frame(cj, zj, tj)[y][:, x] # type: ignore + return d + + def copy(self) -> View: + return View(self) + + def data(self, c: int | Sequence[int] = 0, z: int | Sequence[int] = 0, t: int | Sequence[int] = 0) -> np.ndarray: + """ returns 3D stack of frames """ + c, z, t = (np.arange(self.shape[i]) if e is None else np.array(e, ndmin=1) for i, e in zip('czt', (c, z, t))) + return np.dstack([self.frame(ci, zi, ti) for ci, zi, ti in product(c, z, t)]) + + def frame(self, c: int = 0, z: int = 0, t: int = 0) -> np.ndarray: + """ returns single 2D frame """ + c = self.get_channel(c) + c %= self.base.shape['c'] + z %= self.base.shape['z'] + t %= self.base.shape['t'] + + # cache last n (default 16) frames in memory for speed (~250x faster) + key = (c, z, t, self.transform, self.frame_decorator) + if self.cache.maxlen and key in self.cache: + self.cache.move_to_end(key) + f = self.cache[key] + else: + f = self.transform[self.channel_names[c], t].frame(self.__frame__(c, z, t)) + if self.frame_decorator is not None: + f = self.frame_decorator(self, f, c, z, t) + if self.cache.maxlen: + self.cache[key] = f + if self.dtype is not None: + return f.copy().astype(self.dtype) + else: + return f.copy() + + def get_channel(self, channel_name: str | int) -> int: + if not isinstance(channel_name, str): + return channel_name + else: + c = [i for i, c in enumerate(self.channel_names) if c.lower().startswith(channel_name.lower())] + assert len(c) > 0, f'Channel {c} not found in {self.channel_names}' + assert len(c) < 2, f'Channel {c} not unique in {self.channel_names}' + return c[0] + + @staticmethod + def get_config(file: Path | str) -> Any: + """ Open a yml config file """ + loader = yaml.SafeLoader + loader.add_implicit_resolver( + r'tag:yaml.org,2002:float', + re.compile(r'''^(?: + [-+]?([0-9][0-9_]*)\.[0-9_]*(?:[eE][-+]?[0-9]+)? + |[-+]?([0-9][0-9_]*)([eE][-+]?[0-9]+) + |\.[0-9_]+(?:[eE][-+][0-9]+)? + |[-+]?[0-9][0-9_]*(?::[0-5]?[0-9])+\.[0-9_]* + |[-+]?\\.(?:inf|Inf|INF) + |\.(?:nan|NaN|NAN))$''', re.X), + list(r'-+0123456789.')) + with open(file) as f: + return yaml.load(f, loader) + + def get_czt(self, c: int | Sequence[int], z: int | Sequence[int], + t: int | Sequence[int]) -> tuple[list[int], list[int], list[int]]: + czt = [] + for i, n in zip('czt', (c, z, t)): + if n is None: + czt.append(list(range(self.shape[i]))) + elif isinstance(n, range): + if n.stop < 0: + stop = n.stop % self.shape[i] + elif n.stop > self.shape[i]: + stop = self.shape[i] + else: + stop = n.stop + czt.append(list(range(n.start % self.shape[i], stop, n.step))) + elif isinstance(n, Number): + czt.append([n % self.shape[i]]) + else: + czt.append([k % self.shape[i] for k in n]) + return [self.get_channel(c) for c in czt[0]], *czt[1:3] # type: ignore + + @staticmethod + 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] + + @staticmethod + def fix_ome(ome: OME) -> OME: + # fix ome if necessary + for image in ome.images: + try: + if image.pixels.physical_size_z is None and len(set([plane.the_z + for plane in image.pixels.planes])) > 1: + z = np.array([(plane.position_z * ureg.Quantity(plane.position_z_unit.value).to(ureg.m).magnitude, + plane.the_z) + for plane in image.pixels.planes if plane.the_c == 0 and plane.the_t == 0]) + i = np.argsort(z[:, 1]) + image.pixels.physical_size_z = np.nanmean(np.true_divide(*np.diff(z[i], axis=0).T)) * 1e6 + image.pixels.physical_size_z_unit = 'µm' # type: ignore + except Exception: # noqa + pass + return ome + + @staticmethod + 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')) + + def get_ome(self) -> OME: + """ overload this """ + return self.bioformats_ome(self.path) + + @cached_property + def ome(self) -> OME: + cache = OmeCache() + if self.path 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] + + def is_noise(self, volume: ArrayLike = None) -> bool: + """ True if volume only has noise """ + if volume is None: + volume = self + fft = np.fft.fftn(volume) + corr = np.fft.fftshift(np.fft.ifftn(fft * fft.conj()).real / np.sum(volume ** 2)) + return 1 - corr[tuple([0] * corr.ndim)] < 0.0067 + + @staticmethod + def kill_vm() -> None: + JVM().kill_vm() + + def new(self, *args: Any, **kwargs: Any) -> View: + warnings.warn('Imread.new has been deprecated, use Imread.view instead.', DeprecationWarning, 2) + return self.view(*args, **kwargs) + + def save_as_movie(self, fname: Path | str = None, + c: int | Sequence[int] = None, z: int | Sequence[int] = None, # noqa + t: str | int | Sequence[int] = None, # noqa + colors: tuple[str] = None, brightnesses: tuple[float] = None, + scale: int = None, bar: bool = True) -> None: + """ saves the image as a mp4 or mkv file """ + from matplotlib.colors import to_rgb + from skvideo.io import FFmpegWriter + + if t is None: + t = np.arange(self.shape['t']) + elif isinstance(t, str): + t = eval(f"np.arange(self.shape['t'])[{t}]") + elif np.isscalar(t): + t = (t,) + + def get_ab(tyx: Imread, p: tuple[float, float] = (1, 99)) -> tuple[float, float]: + s = tyx.flatten() + s = s[s > 0] + a, b = np.percentile(s, p) + if a == b: + a, b = np.min(s), np.max(s) + if a == b: + a, b = 0, 1 + return a, b + + def cframe(frame: ArrayLike, color: str, a: float, b: float, scale: float = 1) -> np.ndarray: # noqa + color = to_rgb(color) + frame = (frame - a) / (b - a) + frame = np.dstack([255 * frame * i for i in color]) + return np.clip(np.round(frame), 0, 255).astype('uint8') + + ab = list(zip(*[get_ab(i) for i in self.transpose('cztyx')])) # type: ignore + colors = colors or ('r', 'g', 'b')[:self.shape['c']] + max(0, self.shape['c'] - 3) * ('w',) + brightnesses = brightnesses or (1,) * self.shape['c'] + scale = scale or 1 + shape_x = 2 * ((self.shape['x'] * scale + 1) // 2) + shape_y = 2 * ((self.shape['y'] * scale + 1) // 2) + + with FFmpegWriter( + str(fname).format(name=self.path.stem, path=str(self.path.parent)), + outputdict={'-vcodec': 'libx264', '-preset': 'veryslow', '-pix_fmt': 'yuv420p', '-r': '7', + '-vf': f'setpts={25 / 7}*PTS,scale={shape_x}:{shape_y}:flags=neighbor'} + ) as movie: + im = self.transpose('tzcyx') # type: ignore + for ti in tqdm(t, desc='Saving movie', disable=not bar): + movie.writeFrame(np.max([cframe(yx, c, a, b / s, scale) + for yx, a, b, c, s in zip(im[ti].max('z'), *ab, colors, brightnesses)], 0)) + + def save_as_tiff(self, fname: Path | str = None, c: int | Sequence[int] = None, z: int | Sequence[int] = None, + t: int | Sequence[int] = None, split: bool = False, bar: bool = True, pixel_type: str = 'uint16', + **kwargs: Any) -> None: + """ saves the image as a tif file + split: split channels into different files """ + fname = Path(str(fname).format(name=self.path.stem, path=str(self.path.parent))) + if fname is None: + fname = self.path.with_suffix('.tif') + if fname == self.path: + raise FileExistsError(f'File {fname} exists already.') + if not isinstance(fname, Path): + fname = Path(fname) + if split: + for i in range(self.shape['c']): + if self.timeseries: + self.save_as_tiff(fname.with_name(f'{fname.stem}_C{i:01d}').with_suffix('.tif'), i, 0, None, False, + bar, pixel_type) + else: + self.save_as_tiff(fname.with_name(f'{fname.stem}_C{i:01d}').with_suffix('.tif'), i, None, 0, False, + bar, pixel_type) + else: + n = [c, z, t] + for i, ax in enumerate('czt'): + if n[i] is None: + n[i] = range(self.shape[ax]) + elif not isinstance(n[i], (tuple, list)): + n[i] = (n[i],) + + shape = [len(i) for i in n] + with TransformTiff(self, fname.with_suffix('.tif'), dtype=pixel_type, + pxsize=self.pxsize_um, deltaz=self.deltaz_um, **kwargs) as tif: + for i, m in tqdm(zip(product(*[range(s) for s in shape]), product(*n)), # noqa + total=np.prod(shape), desc='Saving tiff', disable=not bar): + tif.save(m, *i) + + def with_transform(self, channels: bool = True, drift: bool = False, file: Path | str = None, + bead_files: Sequence[Path | str] = ()) -> View: + """ returns a view where channels and/or frames are registered with an affine transformation + channels: True/False register channels using bead_files + drift: True/False register frames to correct drift + file: load registration from file with name file, default: transform.yml in self.path.parent + bead_files: files used to register channels, default: files in self.path.parent, + with names starting with 'beads' + """ + view = self.view() + if file is None: + file = Path(view.path.parent) / 'transform.yml' + else: + file = Path(file) + if not bead_files: + try: + bead_files = Transforms.get_bead_files(view.path.parent) + except Exception: # noqa + if not file.exists(): + raise Exception('No transform file and no bead file found.') + bead_files = () + + if channels: + try: + view.transform = Transforms.from_file(file, T=drift) + except Exception: # noqa + view.transform = Transforms().with_beads(view.cyllens, bead_files) + if drift: + view.transform = view.transform.with_drift(view) + view.transform.save(file.with_suffix('.yml')) + view.transform.save_channel_transform_tiff(bead_files, file.with_suffix('.tif')) + elif drift: + try: + view.transform = Transforms.from_file(file, C=False) + except Exception: # noqa + view.transform = Transforms().with_drift(self) + view.transform.adapt(view.frameoffset, view.shape.yxczt, view.channel_names) + return view + + def set_cache_size(self, cache_size: int) -> None: + assert isinstance(cache_size, int) and cache_size >= 0 + self.cache.maxlen = cache_size + self.cache.truncate() + + @staticmethod + def split_path_series(path: Path | str) -> tuple[Path, int]: + if isinstance(path, str): + path = Path(path) + if isinstance(path, Path) and path.name.startswith('Pos') and path.name.lstrip('Pos').isdigit(): + return path.parent, int(path.name.lstrip('Pos')) + return path, 0 + + def view(self, *args: Any, **kwargs: Any) -> View: + return View(self, *args, **kwargs) + + +class View(Imread, ABC): + def __init__(self, base: Imread, dtype: DTypeLike = None) -> None: + super().__init__(base.base, base.slice, base.shape, dtype or base.dtype, base.frame_decorator) + self.transform = base.transform + + def __getattr__(self, item: str) -> Any: + if not hasattr(self.base, item): + raise AttributeError(f'{self.__class__} object has no attribute {item}') + return self.base.__getattribute__(item) + + +class AbstractReader(Imread, metaclass=ABCMeta): + priority = 99 + do_not_pickle = 'cache' + ureg = ureg + + @staticmethod + @abstractmethod + def _can_open(path: Path | str) -> bool: + """ Override this method, and return true when the subclass can open the file """ + return False + + @staticmethod + def get_positions(path: str | Path) -> Optional[list[int]]: + return None + + @abstractmethod + def __frame__(self, c: int, z: int, t: int) -> np.ndarray: + """ Override this, return the frame at c, z, t """ + return np.random.randint(0, 255, self.shape['yx']) + + def open(self) -> None: + """ Optionally override this, open file handles etc. + filehandles cannot be pickled and should be marked such by setting do_not_pickle = 'file_handle_name' """ + return + + def close(self) -> None: + """ Optionally override this, close file handles etc. """ + return + + def __init__(self, path: Path | str | Imread | Any = None, dtype: DTypeLike = None, axes: str = None) -> None: + if isinstance(path, Imread): + return + super().__init__() + self.isclosed = False + if isinstance(path, str): + path = Path(path) + self.path, self.series = self.split_path_series(path) + if isinstance(path, Path) and path.exists(): + self.title = self.path.name + self.acquisitiondate = datetime.fromtimestamp(self.path.stat().st_mtime).strftime('%y-%m-%dT%H:%M:%S') + else: # ndarray + self.title = self.__class__.__name__ + self.acquisitiondate = 'now' + + self.reader = None + self.pcf = None + self.powermode = None + self.collimator = None + self.tirfangle = None + self.cyllens = ['None', 'None'] + self.duolink = 'None' + self.detector = [0, 1] + self.track = [0] + self.cache = DequeDict(16) + self.frameoffset = 0, 0 # how far apart the centers of frame and sensor are + + self.open() + # extract some metadata from ome + instrument = self.ome.instruments[0] if self.ome.instruments else None + image = self.ome.images[self.series if len(self.ome.images) > 1 else 0] + pixels = image.pixels + self.shape = pixels.size_y, pixels.size_x, pixels.size_c, pixels.size_z, pixels.size_t + self.base_shape = Shape((pixels.size_y, pixels.size_x, pixels.size_c, pixels.size_z, pixels.size_t), 'yxczt') + self.dtype = pixels.type.value if dtype is None else dtype + self.pxsize = pixels.physical_size_x_quantity + try: + self.exposuretime = tuple(find(image.pixels.planes, the_c=c).exposure_time_quantity + for c in range(self.shape['c'])) + except AttributeError: + self.exposuretime = () + + if self.zstack: + self.deltaz = image.pixels.physical_size_z_quantity + self.deltaz_um = None if self.deltaz is None else self.deltaz.to(self.ureg.um).m + else: + self.deltaz = self.deltaz_um = None + if image.objective_settings: + self.objective = find(instrument.objectives, id=image.objective_settings.id) + else: + self.objective = None + try: + t0 = find(image.pixels.planes, the_c=0, the_t=0, the_z=0).delta_t + t1 = find(image.pixels.planes, the_c=0, the_t=self.shape['t'] - 1, the_z=0).delta_t + self.timeinterval = (t1 - t0) / (self.shape['t'] - 1) if self.shape['t'] > 1 and t1 > t0 else None + except AttributeError: + self.timeinterval = None + try: + self.binning = [int(i) for i in image.pixels.channels[0].detector_settings.binning.value.split('x')] + if self.pxsize is not None: + self.pxsize *= self.binning[0] + except (AttributeError, IndexError, ValueError): + self.binning = None + self.channel_names = [channel.name for channel in image.pixels.channels] + self.channel_names += [chr(97 + i) for i in range(len(self.channel_names), self.shape['c'])] + self.cnamelist = self.channel_names + try: + optovars = [objective for objective in instrument.objectives if 'tubelens' in objective.id.lower()] + except AttributeError: + optovars = [] + if len(optovars) == 0: + self.tubelens = None + else: + self.tubelens = optovars[0] + if self.objective: + if self.tubelens: + self.magnification = self.objective.nominal_magnification * self.tubelens.nominal_magnification + else: + self.magnification = self.objective.nominal_magnification + self.NA = self.objective.lens_na + else: + self.magnification = None + self.NA = None + + self.gain = [find(instrument.detectors, id=channel.detector_settings.id).amplification_gain + for channel in image.pixels.channels + if channel.detector_settings + and find(instrument.detectors, id=channel.detector_settings.id).amplification_gain] + self.laserwavelengths = [(channel.excitation_wavelength_quantity.to(self.ureg.nm).m,) + for channel in pixels.channels if channel.excitation_wavelength_quantity] + self.laserpowers = try_default(lambda: [(1 - channel.light_source_settings.attenuation,) + for channel in pixels.channels], []) + self.filter = try_default( # type: ignore + lambda: [find(instrument.filter_sets, id=channel.filter_set_ref.id).model + for channel in image.pixels.channels], None) + self.pxsize_um = None if self.pxsize is None else self.pxsize.to(self.ureg.um).m + self.exposuretime_s = [None if i is None else i.to(self.ureg.s).m for i in self.exposuretime] + + if axes is None: + self.axes = ''.join(i for i in 'cztyx' if self.shape[i] > 1) + elif axes.lower() == 'full': + self.axes = 'cztyx' + else: + self.axes = axes + self.slice = [np.arange(s, dtype=int) for s in self.shape.yxczt] + + m = self.extrametadata + if m is not None: + try: + self.cyllens = m['CylLens'] + self.duolink = m['DLFilterSet'].split(' & ')[m['DLFilterChannel']] + if 'FeedbackChannels' in m: + self.feedback = m['FeedbackChannels'] + else: + self.feedback = m['FeedbackChannel'] + except Exception: # noqa + self.cyllens = ['None', 'None'] + self.duolink = 'None' + self.feedback = [] + try: + self.cyllenschannels = np.where([self.cyllens[self.detector[c]].lower() != 'none' + for c in range(self.shape['c'])])[0].tolist() + except Exception: # noqa + pass + try: + s = int(re.findall(r'_(\d{3})_', self.duolink)[0]) * ureg.nm + except Exception: # noqa + s = 561 * ureg.nm + try: + sigma = [] + for c, d in enumerate(self.detector): + emission = (np.hstack(self.laserwavelengths[c]) + 22) * ureg.nm + sigma.append([emission[emission > s].max(initial=0), emission[emission < s].max(initial=0)][d]) + sigma = np.hstack(sigma) + sigma[sigma == 0] = 600 * ureg.nm + sigma /= 2 * self.NA * self.pxsize + self.sigma = sigma.magnitude.tolist() # type: ignore + except Exception: # noqa + self.sigma = [2] * self.shape['c'] + if not self.NA: + self.immersionN = 1 + elif 1.5 < self.NA: + self.immersionN = 1.661 + elif 1.3 < self.NA < 1.5: + self.immersionN = 1.518 + elif 1 < self.NA < 1.3: + self.immersionN = 1.33 + else: + self.immersionN = 1 + + p = re.compile(r'(\d+):(\d+)$') + try: + self.track, self.detector = zip(*[[int(i) for i in p.findall(find( + self.ome.images[self.series].pixels.channels, id=f'Channel:{c}').detector_settings.id)[0]] + for c in range(self.shape['c'])]) + except Exception: # noqa + pass + + +def main() -> None: + parser = ArgumentParser(description='Display info and save as tif') + parser.add_argument('-v', '--version', action='version', version=__version__) + parser.add_argument('file', help='image_file', type=str, nargs='*') + parser.add_argument('-w', '--write', help='path to tif/movie out, {folder}, {name} and {ext} take this from file in', + type=str, default=None) + parser.add_argument('-o', '--extract_ome', help='extract ome to xml file', action='store_true') + parser.add_argument('-r', '--register', help='register channels', action='store_true') + parser.add_argument('-c', '--channel', help='channel', type=int, default=None) + parser.add_argument('-z', '--zslice', help='z-slice', type=int, default=None) + parser.add_argument('-t', '--time', help='time (frames) in python slicing notation', type=str, default=None) + parser.add_argument('-s', '--split', help='split channels', action='store_true') + parser.add_argument('-f', '--force', help='force overwrite', action='store_true') + parser.add_argument('-C', '--movie-colors', help='colors for channels in movie', type=str, nargs='*') + parser.add_argument('-B', '--movie-brightnesses', help='scale brightness of each channel', + type=float, nargs='*') + parser.add_argument('-S', '--movie-scale', help='upscale movie xy size, int', type=float) + args = parser.parse_args() + + for file in tqdm(args.file, desc='operating on files', disable=len(args.file) == 1): + file = Path(file) + with Imread(file) as im: # noqa + if args.register: + im = im.with_transform() # noqa + if args.write: + write = Path(args.write.format(folder=str(file.parent), name=file.stem, ext=file.suffix)).absolute() # noqa + write.parent.mkdir(parents=True, exist_ok=True) + if write.exists() and not args.force: + print(f'File {args.write} exists already, add the -f flag if you want to overwrite it.') + elif write.suffix in ('.mkv', '.mp4'): + im.save_as_movie(write, args.channel, args.zslice, args.time, args.movie_colors, + args.movie_brightnesses, args.movie_scale, bar=len(args.file) == 1) + else: + im.save_as_tiff(write, args.channel, args.zslice, args.time, args.split, bar=len(args.file) == 1) + if args.extract_ome: + with open(im.path.with_suffix('.ome.xml'), 'w') as f: + f.write(im.ome.to_xml()) + if len(args.file) == 1: + print(im.summary) + + +from .readers import * # noqa diff --git a/py/ndbioimage/readers/__init__.py b/py/ndbioimage/readers/__init__.py new file mode 100644 index 0000000..d4ca3bc --- /dev/null +++ b/py/ndbioimage/readers/__init__.py @@ -0,0 +1 @@ +__all__ = 'bfread', 'cziread', 'fijiread', 'ndread', 'seqread', 'tifread', 'metaseriesread' diff --git a/py/ndbioimage/readers/bfread.py b/py/ndbioimage/readers/bfread.py new file mode 100644 index 0000000..bb36314 --- /dev/null +++ b/py/ndbioimage/readers/bfread.py @@ -0,0 +1,37 @@ +from __future__ import annotations + +from abc import ABC +from pathlib import Path + +import numpy as np + +from .. import rs +from .. import AbstractReader + +for file in (Path(__file__).parent / "deps").glob("*_"): + file.rename(str(file)[:-1]) + +if not list((Path(__file__).parent / "jassets").glob("bioformats*.jar")): + rs.download_bioformats(True) + + +class Reader(AbstractReader, ABC): + priority = 99 # panic and open with BioFormats + do_not_pickle = 'reader' + + @staticmethod + def _can_open(path: Path) -> bool: + try: + _ = rs.Reader(path) + return True + except Exception: + return False + + def open(self) -> None: + self.reader = rs.Reader(str(self.path), int(self.series)) + + def __frame__(self, c: int, z: int, t: int) -> np.ndarray: + return self.reader.get_frame(int(c), int(z), int(t)) + + def close(self) -> None: + self.reader.close() diff --git a/py/ndbioimage/readers/cziread.py b/py/ndbioimage/readers/cziread.py new file mode 100644 index 0000000..c345d7b --- /dev/null +++ b/py/ndbioimage/readers/cziread.py @@ -0,0 +1,606 @@ +import re +import warnings +from abc import ABC +from functools import cached_property +from io import BytesIO +from itertools import product +from pathlib import Path +from typing import Any, Callable, Optional, TypeVar + +import czifile +import imagecodecs +import numpy as np +from lxml import etree +from ome_types import OME, model +from tifffile import repeat_nd + +from .. import AbstractReader + +try: + # TODO: use zoom from imagecodecs implementation when available + from scipy.ndimage.interpolation import zoom +except ImportError: + try: + from ndimage.interpolation import zoom + except ImportError: + zoom = None + + +Element = TypeVar('Element') + + +def zstd_decode(data: bytes) -> bytes: # noqa + """ decode zstd bytes, copied from BioFormats ZeissCZIReader """ + def read_var_int(stream: BytesIO) -> int: # noqa + a = stream.read(1)[0] + if a & 128: + b = stream.read(1)[0] + if b & 128: + c = stream.read(1)[0] + return (c << 14) | ((b & 127) << 7) | (a & 127) + return (b << 7) | (a & 127) + return a & 255 + + try: + with BytesIO(data) as stream: + size_of_header = read_var_int(stream) + high_low_unpacking = False + while stream.tell() < size_of_header: + chunk_id = read_var_int(stream) + # only one chunk ID defined so far + if chunk_id == 1: + high_low_unpacking = (stream.read(1)[0] & 1) == 1 + else: + raise ValueError(f'Invalid chunk id: {chunk_id}') + pointer = stream.tell() + except Exception: # noqa + high_low_unpacking = False + pointer = 0 + + decoded = imagecodecs.zstd_decode(data[pointer:]) + if high_low_unpacking: + second_half = len(decoded) // 2 + return bytes([decoded[second_half + i // 2] if i % 2 else decoded[i // 2] for i in range(len(decoded))]) + else: + return decoded + + +def data(self, raw: bool = False, resize: bool = True, order: int = 0) -> np.ndarray: + """Read image data from file and return as numpy array.""" + DECOMPRESS = czifile.czifile.DECOMPRESS # noqa + DECOMPRESS[5] = imagecodecs.zstd_decode + DECOMPRESS[6] = zstd_decode + + de = self.directory_entry + fh = self._fh + if raw: + with fh.lock: + fh.seek(self.data_offset) + data = fh.read(self.data_size) # noqa + return data + if de.compression: + # if de.compression not in DECOMPRESS: + # raise ValueError('compression unknown or not supported') + with fh.lock: + fh.seek(self.data_offset) + data = fh.read(self.data_size) # noqa + data = DECOMPRESS[de.compression](data) # noqa + if de.compression == 2: + # LZW + data = np.fromstring(data, de.dtype) # noqa + elif de.compression in (5, 6): + # ZSTD + data = np.frombuffer(data, de.dtype) # noqa + else: + dtype = np.dtype(de.dtype) + with fh.lock: + fh.seek(self.data_offset) + data = fh.read_array(dtype, self.data_size // dtype.itemsize) # noqa + + data = data.reshape(de.stored_shape) # noqa + if de.compression != 4 and de.stored_shape[-1] in (3, 4): + if de.stored_shape[-1] == 3: + # BGR -> RGB + data = data[..., ::-1] # noqa + else: + # BGRA -> RGBA + tmp = data[..., 0].copy() + data[..., 0] = data[..., 2] + data[..., 2] = tmp + if de.stored_shape == de.shape or not resize: + return data + + # sub / supersampling + factors = [j / i for i, j in zip(de.stored_shape, de.shape)] + factors = [(int(round(f)) if abs(f - round(f)) < 0.0001 else f) + for f in factors] + + # use repeat if possible + if order == 0 and all(isinstance(f, int) for f in factors): + data = repeat_nd(data, factors).copy() # noqa + data.shape = de.shape + return data + + # remove leading dimensions with size 1 for speed + shape = list(de.stored_shape) + i = 0 + for s in shape: + if s != 1: + break + i += 1 + shape = shape[i:] + factors = factors[i:] + data.shape = shape + + # resize RGB components separately for speed + if zoom is None: + raise ImportError("cannot import 'zoom' from scipy or ndimage") + if shape[-1] in (3, 4) and factors[-1] == 1.0: + factors = factors[:-1] + old = data + data = np.empty(de.shape, de.dtype[-2:]) # noqa + for i in range(shape[-1]): + data[..., i] = zoom(old[..., i], zoom=factors, order=order) + else: + data = zoom(data, zoom=factors, order=order) # noqa + + data.shape = de.shape + return data + + +# monkeypatch zstd into czifile +czifile.czifile.SubBlockSegment.data = data + + +class Reader(AbstractReader, ABC): + priority = 0 + do_not_pickle = 'reader', 'filedict' + + @staticmethod + def _can_open(path: Path) -> bool: + return isinstance(path, Path) and path.suffix == '.czi' + + def open(self) -> None: + self.reader = czifile.CziFile(self.path) + filedict = {} + 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] + if len(filedict) == 0: + raise FileNotFoundError(f'Series {self.series} not found in {self.path}.') + self.filedict = filedict # noqa + + def close(self) -> None: + self.reader.close() + + def get_ome(self) -> OME: + return OmeParse.get_ome(self.reader, self.filedict) + + def __frame__(self, c: int = 0, z: int = 0, t: int = 0) -> np.ndarray: + 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} + 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 = tuple(index[self.reader.axes.index(i)] for i in 'YX') + 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)] + + +class OmeParse: + size_x: int + size_y: int + size_c: int + size_z: int + size_t: int + + nm = model.UnitsLength.NANOMETER + um = model.UnitsLength.MICROMETER + + @classmethod + def get_ome(cls, reader: czifile.CziFile, filedict: dict[tuple[int, int, int], Any]) -> OME: + new = cls(reader, filedict) + new.parse() + return new.ome + + def __init__(self, reader: czifile.CziFile, filedict: dict[tuple[int, int, int], Any]) -> None: + self.reader = reader + self.filedict = filedict + xml = reader.metadata() + self.attachments = {i.attachment_entry.name: i.attachment_entry.data_segment() + for i in reader.attachments()} + self.tree = etree.fromstring(xml) + self.metadata = self.tree.find('Metadata') + version = self.metadata.find('Version') + if version is not None: + self.version = version.text + else: + self.version = self.metadata.find('Experiment').attrib['Version'] + + self.ome = OME() + self.information = self.metadata.find('Information') + self.display_setting = self.metadata.find('DisplaySetting') + self.experiment = self.metadata.find('Experiment') + self.acquisition_block = self.experiment.find('ExperimentBlocks').find('AcquisitionBlock') + self.instrument = self.information.find('Instrument') + self.image = self.information.find('Image') + + if self.version == '1.0': + self.experiment = self.metadata.find('Experiment') + self.acquisition_block = self.experiment.find('ExperimentBlocks').find('AcquisitionBlock') + self.multi_track_setup = self.acquisition_block.find('MultiTrackSetup') + else: + self.experiment = None + self.acquisition_block = None + self.multi_track_setup = None + + def parse(self) -> None: + self.get_experimenters() + self.get_instruments() + self.get_detectors() + self.get_objectives() + self.get_tubelenses() + self.get_light_sources() + self.get_filters() + self.get_pixels() + self.get_channels() + self.get_planes() + self.get_annotations() + + @staticmethod + def text(item: Optional[Element], default: str = "") -> str: + return default if item is None else item.text + + @staticmethod + def def_list(item: Any) -> list[Any]: + return [] if item is None else item + + @staticmethod + def try_default(fun: Callable[[Any, ...], Any] | type, default: Any = None, *args: Any, **kwargs: Any) -> Any: + try: + return fun(*args, **kwargs) + except Exception: # noqa + return default + + def get_experimenters(self) -> None: + if self.version == '1.0': + self.ome.experimenters = [ + model.Experimenter(id='Experimenter:0', + user_name=self.information.find('User').find('DisplayName').text)] + elif self.version in ('1.1', '1.2'): + self.ome.experimenters = [ + model.Experimenter(id='Experimenter:0', + user_name=self.information.find('Document').find('UserName').text)] + + def get_instruments(self) -> None: + if self.version == '1.0': + self.ome.instruments.append(model.Instrument(id=self.instrument.attrib['Id'])) + elif self.version in ('1.1', '1.2'): + for _ in self.instrument.find('Microscopes'): + self.ome.instruments.append(model.Instrument(id='Instrument:0')) + + def get_detectors(self) -> None: + if self.version == '1.0': + for detector in self.instrument.find('Detectors'): + try: + detector_type = model.Detector_Type(self.text(detector.find('Type')).upper() or "") + except ValueError: + detector_type = model.Detector_Type.OTHER + + self.ome.instruments[0].detectors.append( + model.Detector( + id=detector.attrib['Id'], model=self.text(detector.find('Manufacturer').find('Model')), + amplification_gain=float(self.text(detector.find('AmplificationGain'))), + gain=float(self.text(detector.find('Gain'))), zoom=float(self.text(detector.find('Zoom'))), + type=detector_type + )) + elif self.version in ('1.1', '1.2'): + for detector in self.instrument.find('Detectors'): + try: + detector_type = model.Detector_Type(self.text(detector.find('Type')).upper() or "") + except ValueError: + detector_type = model.Detector_Type.OTHER + + self.ome.instruments[0].detectors.append( + model.Detector( + id=detector.attrib['Id'].replace(' ', ''), + model=self.text(detector.find('Manufacturer').find('Model')), + type=detector_type + )) + + def get_objectives(self) -> None: + for objective in self.instrument.find('Objectives'): + self.ome.instruments[0].objectives.append( + model.Objective( + id=objective.attrib['Id'], + model=self.text(objective.find('Manufacturer').find('Model')), + immersion=self.text(objective.find('Immersion')), # type: ignore + lens_na=float(self.text(objective.find('LensNA'))), + nominal_magnification=float(self.text(objective.find('NominalMagnification'))))) + + def get_tubelenses(self) -> None: + if self.version == '1.0': + for idx, tube_lens in enumerate({self.text(track_setup.find('TubeLensPosition')) + for track_setup in self.multi_track_setup}): + try: + nominal_magnification = float(re.findall(r'\d+[,.]\d*', tube_lens)[0].replace(',', '.')) + except Exception: # noqa + nominal_magnification = 1.0 + + self.ome.instruments[0].objectives.append( + model.Objective(id=f'Objective:Tubelens:{idx}', model=tube_lens, + nominal_magnification=nominal_magnification)) + elif self.version in ('1.1', '1.2'): + for tubelens in self.def_list(self.instrument.find('TubeLenses')): + try: + nominal_magnification = float(re.findall(r'\d+(?:[,.]\d*)?', + tubelens.attrib['Name'])[0].replace(',', '.')) + except Exception: # noqa + nominal_magnification = 1.0 + + self.ome.instruments[0].objectives.append( + model.Objective( + id=f"Objective:{tubelens.attrib['Id']}", + model=tubelens.attrib['Name'], + nominal_magnification=nominal_magnification)) + + def get_light_sources(self) -> None: + if self.version == '1.0': + for light_source in self.def_list(self.instrument.find('LightSources')): + try: + if light_source.find('LightSourceType').find('Laser') is not None: + self.ome.instruments[0].lasers.append( + model.Laser( + id=light_source.attrib['Id'], + model=self.text(light_source.find('Manufacturer').find('Model')), + power=float(self.text(light_source.find('Power'))), + wavelength=float( + self.text(light_source.find('LightSourceType').find('Laser').find('Wavelength'))))) + except AttributeError: + pass + elif self.version in ('1.1', '1.2'): + for light_source in self.def_list(self.instrument.find('LightSources')): + try: + if light_source.find('LightSourceType').find('Laser') is not None: + self.ome.instruments[0].lasers.append( + model.Laser( + id=f"LightSource:{light_source.attrib['Id']}", + power=float(self.text(light_source.find('Power'))), + wavelength=float(light_source.attrib['Id'][-3:]))) # TODO: follow Id reference + except (AttributeError, ValueError): + pass + + def get_filters(self) -> None: + if self.version == '1.0': + for idx, filter_ in enumerate({self.text(beam_splitter.find('Filter')) + for track_setup in self.multi_track_setup + for beam_splitter in track_setup.find('BeamSplitters')}): + self.ome.instruments[0].filter_sets.append( + model.FilterSet(id=f'FilterSet:{idx}', model=filter_) + ) + + def get_pixels(self) -> None: + x_min = min([f.start[f.axes.index('X')] for f in self.filedict[0, 0, 0]]) + y_min = min([f.start[f.axes.index('Y')] for f in self.filedict[0, 0, 0]]) + x_max = max([f.start[f.axes.index('X')] + f.shape[f.axes.index('X')] for f in self.filedict[0, 0, 0]]) + y_max = max([f.start[f.axes.index('Y')] + f.shape[f.axes.index('Y')] for f in self.filedict[0, 0, 0]]) + 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') + image = self.information.find('Image') + pixel_type = self.text(image.find('PixelType'), 'Gray16') + if pixel_type.startswith('Gray'): + pixel_type = 'uint' + pixel_type[4:] + objective_settings = image.find('ObjectiveSettings') + + self.ome.images.append( + model.Image( + id='Image:0', + name=f"{self.text(self.information.find('Document').find('Name'))} #1", + pixels=model.Pixels( + id='Pixels:0', size_x=self.size_x, size_y=self.size_y, + size_c=self.size_c, size_z=self.size_z, size_t=self.size_t, + dimension_order='XYCZT', type=pixel_type, # type: ignore + significant_bits=int(self.text(image.find('ComponentBitCount'))), + big_endian=False, interleaved=False, metadata_only=True), # type: ignore + experimenter_ref=model.ExperimenterRef(id='Experimenter:0'), + instrument_ref=model.InstrumentRef(id='Instrument:0'), + objective_settings=model.ObjectiveSettings( + id=objective_settings.find('ObjectiveRef').attrib['Id'], + medium=self.text(objective_settings.find('Medium')), # type: ignore + refractive_index=float(self.text(objective_settings.find('RefractiveIndex')))), + stage_label=model.StageLabel( + name=f'Scene position #0', + x=self.positions[0], x_unit=self.um, + y=self.positions[1], y_unit=self.um, + z=self.positions[2], z_unit=self.um))) + + for distance in self.metadata.find('Scaling').find('Items'): + if distance.attrib['Id'] == 'X': + self.ome.images[0].pixels.physical_size_x = float(self.text(distance.find('Value'))) * 1e6 + elif distance.attrib['Id'] == 'Y': + self.ome.images[0].pixels.physical_size_y = float(self.text(distance.find('Value'))) * 1e6 + elif self.size_z > 1 and distance.attrib['Id'] == 'Z': + self.ome.images[0].pixels.physical_size_z = float(self.text(distance.find('Value'))) * 1e6 + + @cached_property + def positions(self) -> tuple[float, float, Optional[float]]: + if self.version == '1.0': + scenes = self.image.find('Dimensions').find('S').find('Scenes') + positions = scenes[0].find('Positions')[0] + return float(positions.attrib['X']), float(positions.attrib['Y']), float(positions.attrib['Z']) + elif self.version in ('1.1', '1.2'): + try: # TODO + scenes = self.image.find('Dimensions').find('S').find('Scenes') + center_position = [float(pos) for pos in self.text(scenes[0].find('CenterPosition')).split(',')] + except AttributeError: + center_position = [0, 0] + return center_position[0], center_position[1], None + + @cached_property + def channels_im(self) -> dict: + return {channel.attrib['Id']: channel for channel in self.image.find('Dimensions').find('Channels')} + + @cached_property + def channels_ds(self) -> dict: + return {channel.attrib['Id']: channel for channel in self.display_setting.find('Channels')} + + @cached_property + def channels_ts(self) -> dict: + return {detector.attrib['Id']: track_setup + for track_setup in + self.experiment.find('ExperimentBlocks').find('AcquisitionBlock').find('MultiTrackSetup') + for detector in track_setup.find('Detectors')} + + def get_channels(self) -> None: + if self.version == '1.0': + for idx, (key, channel) in enumerate(self.channels_im.items()): + detector_settings = channel.find('DetectorSettings') + laser_scan_info = channel.find('LaserScanInfo') + detector = detector_settings.find('Detector') + try: + binning = model.Binning(self.text(detector_settings.find('Binning'))) + except ValueError: + binning = model.Binning.OTHER + + filterset = self.text(self.channels_ts[key].find('BeamSplitters')[0].find('Filter')) + filterset_idx = [filterset.model for filterset in self.ome.instruments[0].filter_sets].index(filterset) + + light_sources_settings = channel.find('LightSourcesSettings') + # no space in ome for multiple lightsources simultaneously + if len(light_sources_settings) > idx: + light_source_settings = light_sources_settings[idx] + else: + light_source_settings = light_sources_settings[0] + light_source_settings = model.LightSourceSettings( + id=light_source_settings.find('LightSource').attrib['Id'], + attenuation=float(self.text(light_source_settings.find('Attenuation'))), + wavelength=float(self.text(light_source_settings.find('Wavelength'))), + wavelength_unit=self.nm) + + self.ome.images[0].pixels.channels.append( + model.Channel( + id=f'Channel:{idx}', + name=channel.attrib['Name'], + acquisition_mode=self.text(channel.find('AcquisitionMode')), # type: ignore + color=model.Color(self.text(self.channels_ds[channel.attrib['Id']].find('Color'), 'white')), + detector_settings=model.DetectorSettings(id=detector.attrib['Id'], binning=binning), + # emission_wavelength=text(channel.find('EmissionWavelength')), # TODO: fix + excitation_wavelength=light_source_settings.wavelength, + filter_set_ref=model.FilterSetRef(id=self.ome.instruments[0].filter_sets[filterset_idx].id), + illumination_type=self.text(channel.find('IlluminationType')), # type: ignore + light_source_settings=light_source_settings, + samples_per_pixel=int(self.text(laser_scan_info.find('Averaging'))))) + elif self.version in ('1.1', '1.2'): + for idx, (key, channel) in enumerate(self.channels_im.items()): + detector_settings = channel.find('DetectorSettings') + laser_scan_info = channel.find('LaserScanInfo') + detector = detector_settings.find('Detector') + try: + color = model.Color(self.text(self.channels_ds[channel.attrib['Id']].find('Color'), 'white')) + except Exception: # noqa + color = None + try: + if (i := self.text(channel.find('EmissionWavelength'))) != '0': + emission_wavelength = float(i) + else: + emission_wavelength = None + except Exception: # noqa + emission_wavelength = None + if laser_scan_info is not None: + samples_per_pixel = int(self.text(laser_scan_info.find('Averaging'), '1')) + else: + samples_per_pixel = 1 + try: + binning = model.Binning(self.text(detector_settings.find('Binning'))) + except ValueError: + binning = model.Binning.OTHER + + light_sources_settings = channel.find('LightSourcesSettings') + # no space in ome for multiple lightsources simultaneously + if light_sources_settings is not None: + light_source_settings = light_sources_settings[0] + light_source_settings = model.LightSourceSettings( + id='LightSource:' + '_'.join([light_source_settings.find('LightSource').attrib['Id'] + for light_source_settings in light_sources_settings]), + attenuation=self.try_default(float, None, self.text(light_source_settings.find('Attenuation'))), + wavelength=self.try_default(float, None, self.text(light_source_settings.find('Wavelength'))), + wavelength_unit=self.nm) + else: + light_source_settings = None + + self.ome.images[0].pixels.channels.append( + model.Channel( + id=f'Channel:{idx}', + name=channel.attrib['Name'], + acquisition_mode=self.text(channel.find('AcquisitionMode')).replace( # type: ignore + 'SingleMoleculeLocalisation', 'SingleMoleculeImaging'), + color=color, + detector_settings=model.DetectorSettings( + id=detector.attrib['Id'].replace(' ', ""), + binning=binning), + emission_wavelength=emission_wavelength, + excitation_wavelength=self.try_default(float, None, + self.text(channel.find('ExcitationWavelength'))), + # filter_set_ref=model.FilterSetRef(id=ome.instruments[0].filter_sets[filterset_idx].id), + illumination_type=self.text(channel.find('IlluminationType')), # type: ignore + light_source_settings=light_source_settings, + samples_per_pixel=samples_per_pixel)) + + def get_planes(self) -> None: + try: + exposure_times = [float(self.text(channel.find('LaserScanInfo').find('FrameTime'))) + for channel in self.channels_im.values()] + except Exception: # noqa + exposure_times = [None] * len(self.channels_im) + delta_ts = self.attachments['TimeStamps'].data() + dt = np.diff(delta_ts) + if len(dt) and np.std(dt) / np.mean(dt) > 0.02: + dt = np.median(dt[dt > 0]) + delta_ts = dt * np.arange(len(delta_ts)) + warnings.warn(f'delta_t is inconsistent, using median value: {dt}') + + for t, z, c in product(range(self.size_t), range(self.size_z), range(self.size_c)): + self.ome.images[0].pixels.planes.append( + model.Plane(the_c=c, the_z=z, the_t=t, delta_t=delta_ts[t], + exposure_time=exposure_times[c], + position_x=self.positions[0], position_x_unit=self.um, + position_y=self.positions[1], position_y_unit=self.um, + position_z=self.positions[2], position_z_unit=self.um)) + + def get_annotations(self) -> None: + idx = 0 + for layer in [] if (ml := self.metadata.find('Layers')) is None else ml: + rectangle = layer.find('Elements').find('Rectangle') + if rectangle is not None: + geometry = rectangle.find('Geometry') + roi = model.ROI(id=f'ROI:{idx}', description=self.text(layer.find('Usage'))) + roi.union.append( + model.Rectangle( + id='Shape:0:0', + height=float(self.text(geometry.find('Height'))), + width=float(self.text(geometry.find('Width'))), + x=float(self.text(geometry.find('Left'))), + y=float(self.text(geometry.find('Top'))))) + self.ome.rois.append(roi) + self.ome.images[0].roi_refs.append(model.ROIRef(id=f'ROI:{idx}')) + idx += 1 diff --git a/py/ndbioimage/readers/fijiread.py b/py/ndbioimage/readers/fijiread.py new file mode 100644 index 0000000..f11c0e7 --- /dev/null +++ b/py/ndbioimage/readers/fijiread.py @@ -0,0 +1,59 @@ +from abc import ABC +from itertools import product +from pathlib import Path +from struct import unpack +from warnings import warn + +import numpy as np +from ome_types import model +from tifffile import TiffFile + +from .. import AbstractReader + + +class Reader(AbstractReader, ABC): + """ Can read some tif files written with Fiji which are broken because Fiji didn't finish writing. """ + priority = 90 + do_not_pickle = 'reader' + + @staticmethod + def _can_open(path): + if isinstance(path, Path) and path.suffix in ('.tif', '.tiff'): + with TiffFile(path) as tif: + return tif.is_imagej and not tif.is_bigtiff + else: + return False + + def __frame__(self, c, z, t): # Override this, return the frame at c, z, t + self.reader.filehandle.seek(self.offset + t * self.count) + return np.reshape(unpack(self.fmt, self.reader.filehandle.read(self.count)), self.base_shape['yx']) + + def open(self): + warn(f'File {self.path.name} is probably damaged, opening with fijiread.') + self.reader = TiffFile(self.path) + assert self.reader.pages[0].compression == 1, 'Can only read uncompressed tiff files.' + assert self.reader.pages[0].samplesperpixel == 1, 'Can only read 1 sample per pixel.' + self.offset = self.reader.pages[0].dataoffsets[0] # noqa + self.count = self.reader.pages[0].databytecounts[0] # noqa + self.bytes_per_sample = self.reader.pages[0].bitspersample // 8 # noqa + self.fmt = self.reader.byteorder + self.count // self.bytes_per_sample * 'BHILQ'[self.bytes_per_sample - 1] # noqa + + def close(self): + self.reader.close() + + def get_ome(self): + size_y, size_x = self.reader.pages[0].shape + size_c, size_z = 1, 1 + size_t = int(np.floor((self.reader.filehandle.size - self.reader.pages[0].dataoffsets[0]) / self.count)) + pixel_type = model.PixelType(self.reader.pages[0].dtype.name) + ome = model.OME() + ome.instruments.append(model.Instrument()) + ome.images.append( + model.Image( + pixels=model.Pixels( + size_c=size_c, size_z=size_z, size_t=size_t, size_x=size_x, size_y=size_y, + dimension_order='XYCZT', type=pixel_type), + objective_settings=model.ObjectiveSettings(id='Objective:0'))) + for c, z, t in product(range(size_c), range(size_z), range(size_t)): + ome.images[0].pixels.planes.append(model.Plane(the_c=c, the_z=z, the_t=t, delta_t=0)) + return ome diff --git a/py/ndbioimage/readers/metaseriesread.py b/py/ndbioimage/readers/metaseriesread.py new file mode 100644 index 0000000..f8087e4 --- /dev/null +++ b/py/ndbioimage/readers/metaseriesread.py @@ -0,0 +1,80 @@ +import re +from abc import ABC +from pathlib import Path +from typing import Optional + +import tifffile +from ome_types import model +from ome_types.units import _quantity_property # noqa + +from .. import AbstractReader + + +class Reader(AbstractReader, ABC): + priority = 20 + do_not_pickle = 'last_tif' + + @staticmethod + def _can_open(path): + return isinstance(path, Path) and (path.is_dir() or + (path.parent.is_dir() and path.name.lower().startswith('pos'))) + + @staticmethod + def get_positions(path: str | Path) -> Optional[list[int]]: + pat = re.compile(rf's(\d)_t\d+\.(tif|TIF)$') + return sorted({int(m.group(1)) for file in Path(path).iterdir() if (m := pat.search(file.name))}) + + def get_ome(self): + ome = model.OME() + tif = self.get_tif(0) + metadata = tif.metaseries_metadata + size_z = len(tif.pages) + page = tif.pages[0] + shape = {axis.lower(): size for axis, size in zip(page.axes, page.shape)} + size_x, size_y = shape['x'], shape['y'] + + ome.instruments.append(model.Instrument()) + + size_c = 1 + size_t = max(self.filedict.keys()) + 1 + pixel_type = f"uint{metadata['PlaneInfo']['bits-per-pixel']}" + ome.images.append( + model.Image( + pixels=model.Pixels( + size_c=size_c, size_z=size_z, size_t=size_t, + size_x=size_x, size_y=size_y, + dimension_order='XYCZT', type=pixel_type), + objective_settings=model.ObjectiveSettings(id='Objective:0'))) + return ome + + def open(self): + pat = re.compile(rf's{self.series}_t\d+\.(tif|TIF)$') + filelist = sorted([file for file in self.path.iterdir() if pat.search(file.name)]) + pattern = re.compile(r't(\d+)$') + self.filedict = {int(pattern.search(file.stem).group(1)) - 1: file for file in filelist} + if len(self.filedict) == 0: + raise FileNotFoundError + self.last_tif = 0, tifffile.TiffFile(self.filedict[0]) + + def close(self) -> None: + self.last_tif[1].close() + + def get_tif(self, t: int = None): + last_t, tif = self.last_tif + if (t is None or t == last_t) and not tif.filehandle.closed: + return tif + else: + tif.close() + tif = tifffile.TiffFile(self.filedict[t]) + self.last_tif = t, tif + return tif + + def __frame__(self, c=0, z=0, t=0): + tif = self.get_tif(t) + page = tif.pages[z] + if page.axes.upper() == 'YX': + return page.asarray() + elif page.axes.upper() == 'XY': + return page.asarray().T + else: + raise NotImplementedError(f'reading axes {page.axes} is not implemented') diff --git a/py/ndbioimage/readers/ndread.py b/py/ndbioimage/readers/ndread.py new file mode 100644 index 0000000..5a2b216 --- /dev/null +++ b/py/ndbioimage/readers/ndread.py @@ -0,0 +1,53 @@ +from abc import ABC +from itertools import product + +import numpy as np +from ome_types import model + +from .. import AbstractReader + + +class Reader(AbstractReader, ABC): + priority = 20 + + @staticmethod + def _can_open(path): + return isinstance(path, np.ndarray) and 1 <= path.ndim <= 5 + + def get_ome(self): + def shape(size_x=1, size_y=1, size_c=1, size_z=1, size_t=1): # noqa + return size_x, size_y, size_c, size_z, size_t + size_x, size_y, size_c, size_z, size_t = shape(*self.array.shape) + try: + pixel_type = model.PixelType(self.array.dtype.name) + except ValueError: + if self.array.dtype.name.startswith('int'): + pixel_type = model.PixelType('int32') + else: + pixel_type = model.PixelType('float') + + ome = model.OME() + ome.instruments.append(model.Instrument()) + ome.images.append( + model.Image( + pixels=model.Pixels( + size_c=size_c, size_z=size_z, size_t=size_t, size_x=size_x, size_y=size_y, + dimension_order='XYCZT', type=pixel_type), + objective_settings=model.ObjectiveSettings(id='Objective:0'))) + for c, z, t in product(range(size_c), range(size_z), range(size_t)): + ome.images[0].pixels.planes.append(model.Plane(the_c=c, the_z=z, the_t=t, delta_t=0)) + return ome + + def open(self): + if isinstance(self.path, np.ndarray): + self.array = np.array(self.path) + while self.array.ndim < 5: + self.array = np.expand_dims(self.array, -1) # noqa + self.path = 'numpy array' + + def __frame__(self, c, z, t): + frame = self.array[:, :, c, z, t] + if self.axes.find('y') > self.axes.find('x'): + return frame.T + else: + return frame diff --git a/py/ndbioimage/readers/seqread.py b/py/ndbioimage/readers/seqread.py new file mode 100644 index 0000000..4d74683 --- /dev/null +++ b/py/ndbioimage/readers/seqread.py @@ -0,0 +1,146 @@ +import re +from abc import ABC +from datetime import datetime +from itertools import product +from pathlib import Path + +import tifffile +import yaml +from ome_types import model +from ome_types.units import _quantity_property # noqa + +from .. import AbstractReader + + +def lazy_property(function, field, *arg_fields): + def lazy(self): + if self.__dict__.get(field) is None: + self.__dict__[field] = function(*[getattr(self, arg_field) for arg_field in arg_fields]) + try: + self.model_fields_set.add(field) + except Exception: # noqa + pass + return self.__dict__[field] + return property(lazy) + + +class Plane(model.Plane): + """ Lazily retrieve delta_t from metadata """ + def __init__(self, t0, file, **kwargs): # noqa + super().__init__(**kwargs) + # setting fields here because they would be removed by ome_types/pydantic after class definition + setattr(self.__class__, 'delta_t', lazy_property(self.get_delta_t, 'delta_t', 't0', 'file')) + setattr(self.__class__, 'delta_t_quantity', _quantity_property('delta_t')) + self.__dict__['t0'] = t0 # noqa + self.__dict__['file'] = file # noqa + + @staticmethod + def get_delta_t(t0, file): + with tifffile.TiffFile(file) as tif: + info = yaml.safe_load(tif.pages[0].tags[50839].value['Info']) + return float((datetime.strptime(info['Time'], '%Y-%m-%d %H:%M:%S %z') - t0).seconds) + + +class Reader(AbstractReader, ABC): + priority = 10 + + @staticmethod + def _can_open(path): + pat = re.compile(r'(?:\d+-)?Pos.*', re.IGNORECASE) + return (isinstance(path, Path) and path.is_dir() and + (pat.match(path.name) or any(file.is_dir() and pat.match(file.stem) for file in path.iterdir()))) + + def get_ome(self): + ome = model.OME() + with tifffile.TiffFile(self.filedict[0, 0, 0]) as tif: + metadata = {key: yaml.safe_load(value) for key, value in tif.pages[0].tags[50839].value.items()} + ome.experimenters.append( + model.Experimenter(id='Experimenter:0', user_name=metadata['Info']['Summary']['UserName'])) + objective_str = metadata['Info']['ZeissObjectiveTurret-Label'] + ome.instruments.append(model.Instrument()) + ome.instruments[0].objectives.append( + model.Objective( + id='Objective:0', manufacturer='Zeiss', model=objective_str, + nominal_magnification=float(re.findall(r'(\d+)x', objective_str)[0]), + lens_na=float(re.findall(r'/(\d\.\d+)', objective_str)[0]), + immersion=model.Objective_Immersion.OIL if 'oil' in objective_str.lower() else None)) + tubelens_str = metadata['Info']['ZeissOptovar-Label'] + ome.instruments[0].objectives.append( + model.Objective( + id='Objective:Tubelens:0', manufacturer='Zeiss', model=tubelens_str, + nominal_magnification=float(re.findall(r'\d?\d*[,.]?\d+(?=x$)', tubelens_str)[0].replace(',', '.')))) + ome.instruments[0].detectors.append( + model.Detector( + id='Detector:0', amplification_gain=100)) + ome.instruments[0].filter_sets.append( + model.FilterSet(id='FilterSet:0', model=metadata['Info']['ZeissReflectorTurret-Label'])) + + pxsize = metadata['Info']['PixelSizeUm'] + pxsize_cam = 6.5 if 'Hamamatsu' in metadata['Info']['Core-Camera'] else None + if pxsize == 0: + pxsize = pxsize_cam / ome.instruments[0].objectives[0].nominal_magnification + pixel_type = metadata['Info']['PixelType'].lower() + if pixel_type.startswith('gray'): + pixel_type = 'uint' + pixel_type[4:] + else: + pixel_type = 'uint16' # assume + + size_c, size_z, size_t = (max(i) + 1 for i in zip(*self.filedict.keys())) + t0 = datetime.strptime(metadata['Info']['Time'], '%Y-%m-%d %H:%M:%S %z') + ome.images.append( + model.Image( + pixels=model.Pixels( + size_c=size_c, size_z=size_z, size_t=size_t, + size_x=metadata['Info']['Width'], size_y=metadata['Info']['Height'], + dimension_order='XYCZT', # type: ignore + type=pixel_type, physical_size_x=pxsize, physical_size_y=pxsize, + physical_size_z=metadata['Info']['Summary']['z-step_um']), + objective_settings=model.ObjectiveSettings(id='Objective:0'))) + + for c, z, t in product(range(size_c), range(size_z), range(size_t)): + ome.images[0].pixels.planes.append( + Plane(t0, self.filedict[c, z, t], + the_c=c, the_z=z, the_t=t, exposure_time=metadata['Info']['Exposure-ms'] / 1000)) + + # compare channel names from metadata with filenames + pattern_c = re.compile(r'img_\d{3,}_(.*)_\d{3,}$', re.IGNORECASE) + for c in range(size_c): + ome.images[0].pixels.channels.append( + model.Channel( + id=f'Channel:{c}', name=pattern_c.findall(self.filedict[c, 0, 0].stem)[0], + detector_settings=model.DetectorSettings( + id='Detector:0', binning=metadata['Info']['Hamamatsu_sCMOS-Binning']), + filter_set_ref=model.FilterSetRef(id='FilterSet:0'))) + return ome + + def open(self): + # /some_path/Pos4: path = /some_path, series = 4 + # /some_path/5-Pos_001_005: path = /some_path/5-Pos_001_005, series = 0 + if re.match(r'(?:\d+-)?Pos.*', self.path.name, re.IGNORECASE) is None: + pat = re.compile(rf'^(?:\d+-)?Pos{self.series}$', re.IGNORECASE) + files = sorted(file for file in self.path.iterdir() if pat.match(file.name)) + if len(files): + path = files[0] + else: + raise FileNotFoundError(self.path / pat.pattern) + else: + path = self.path + + pat = re.compile(r'^img_\d{3,}.*\d{3,}.*\.tif$', re.IGNORECASE) + filelist = sorted([file for file in path.iterdir() if pat.search(file.name)]) + with tifffile.TiffFile(self.path / filelist[0]) as tif: + metadata = {key: yaml.safe_load(value) for key, value in tif.pages[0].tags[50839].value.items()} + + # compare channel names from metadata with filenames + cnamelist = metadata['Info']['Summary']['ChNames'] + cnamelist = [c for c in cnamelist if any([c in f.name for f in filelist])] + + pattern_c = re.compile(r'img_\d{3,}_(.*)_\d{3,}$', re.IGNORECASE) + pattern_z = re.compile(r'(\d{3,})$') + pattern_t = re.compile(r'img_(\d{3,})', re.IGNORECASE) + self.filedict = {(cnamelist.index(pattern_c.findall(file.stem)[0]), # noqa + int(pattern_z.findall(file.stem)[0]), + int(pattern_t.findall(file.stem)[0])): file for file in filelist} + + def __frame__(self, c=0, z=0, t=0): + return tifffile.imread(self.path / self.filedict[(c, z, t)]) diff --git a/py/ndbioimage/readers/tifread.py b/py/ndbioimage/readers/tifread.py new file mode 100644 index 0000000..5d3c033 --- /dev/null +++ b/py/ndbioimage/readers/tifread.py @@ -0,0 +1,85 @@ +from abc import ABC +from functools import cached_property +from itertools import product +from pathlib import Path + +import numpy as np +import tifffile +import yaml +from ome_types import model + +from .. import AbstractReader + + +class Reader(AbstractReader, ABC): + priority = 0 + do_not_pickle = 'reader' + + @staticmethod + def _can_open(path): + if isinstance(path, Path) and path.suffix in ('.tif', '.tiff'): + with tifffile.TiffFile(path) as tif: + return tif.is_imagej and tif.pages[-1]._nextifd() == 0 # noqa + else: + return False + + @cached_property + def metadata(self): + return {key: yaml.safe_load(value) if isinstance(value, str) else value + for key, value in self.reader.imagej_metadata.items()} + + def get_ome(self): + page = self.reader.pages[0] + size_y = page.imagelength + size_x = page.imagewidth + if self.p_ndim == 3: + size_c = page.samplesperpixel + size_t = self.metadata.get('frames', 1) # // C + else: + size_c = self.metadata.get('channels', 1) + size_t = self.metadata.get('frames', 1) + size_z = self.metadata.get('slices', 1) + if 282 in page.tags and 296 in page.tags and page.tags[296].value == 1: + f = page.tags[282].value + pxsize = f[1] / f[0] + else: + pxsize = None + + dtype = page.dtype.name + if dtype not in ('int8', 'int16', 'int32', 'uint8', 'uint16', 'uint32', + 'float', 'double', 'complex', 'double-complex', 'bit'): + dtype = 'float' + + interval_t = self.metadata.get('interval', 0) + + ome = model.OME() + ome.instruments.append(model.Instrument(id='Instrument:0')) + ome.instruments[0].objectives.append(model.Objective(id='Objective:0')) + ome.images.append( + model.Image( + id='Image:0', + pixels=model.Pixels( + id='Pixels:0', + size_c=size_c, size_z=size_z, size_t=size_t, size_x=size_x, size_y=size_y, + dimension_order='XYCZT', type=dtype, # type: ignore + physical_size_x=pxsize, physical_size_y=pxsize), + objective_settings=model.ObjectiveSettings(id='Objective:0'))) + for c, z, t in product(range(size_c), range(size_z), range(size_t)): + ome.images[0].pixels.planes.append(model.Plane(the_c=c, the_z=z, the_t=t, delta_t=interval_t * t)) + return ome + + def open(self): + self.reader = tifffile.TiffFile(self.path) + page = self.reader.pages[0] + self.p_ndim = page.ndim # noqa + if self.p_ndim == 3: + self.p_transpose = [i for i in [page.axes.find(j) for j in 'SYX'] if i >= 0] # noqa + + def close(self): + self.reader.close() + + def __frame__(self, c, z, t): + if self.p_ndim == 3: + return np.transpose(self.reader.asarray(z + t * self.base_shape['z']), self.p_transpose)[c] + else: + return self.reader.asarray(c + z * self.base_shape['c'] + t * self.base_shape['c'] * self.base_shape['z']) diff --git a/py/ndbioimage/rs.py b/py/ndbioimage/rs.py new file mode 100644 index 0000000..37ad889 --- /dev/null +++ b/py/ndbioimage/rs.py @@ -0,0 +1,9 @@ +from .ndbioimage_rs import Reader, download_bioformats # noqa +from pathlib import Path + + +if not list((Path(__file__).parent / "jassets").glob("bioformats*.jar")): + download_bioformats(True) + +for file in (Path(__file__).parent / "deps").glob("*_"): + file.rename(str(file)[:-1]) diff --git a/py/ndbioimage/transform.txt b/py/ndbioimage/transform.txt new file mode 100644 index 0000000..a177cc0 --- /dev/null +++ b/py/ndbioimage/transform.txt @@ -0,0 +1,7 @@ +#Insight Transform File V1.0 +#Transform 0 +Transform: CompositeTransform_double_2_2 +#Transform 1 +Transform: AffineTransform_double_2_2 +Parameters: 1 0 0 1 0 0 +FixedParameters: 255.5 255.5 diff --git a/py/ndbioimage/transforms.py b/py/ndbioimage/transforms.py new file mode 100644 index 0000000..4d73a9a --- /dev/null +++ b/py/ndbioimage/transforms.py @@ -0,0 +1,462 @@ +import warnings +from copy import deepcopy +from pathlib import Path + +import numpy as np +import yaml +from parfor import Chunks, pmap +from skimage import filters +from tiffwrite import IJTiffFile +from tqdm.auto import tqdm + +try: + # best if SimpleElastix is installed: https://simpleelastix.readthedocs.io/GettingStarted.html + import SimpleITK as sitk # noqa +except ImportError: + sitk = None + +try: + from pandas import DataFrame, Series, concat +except ImportError: + DataFrame, Series, concat = None, None, None + + +if hasattr(yaml, 'full_load'): + yamlload = yaml.full_load +else: + yamlload = yaml.load + + +class Transforms(dict): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.default = Transform() + + @classmethod + def from_file(cls, file, C=True, T=True): + with open(Path(file).with_suffix('.yml')) as f: + return cls.from_dict(yamlload(f), C, T) + + @classmethod + def from_dict(cls, d, C=True, T=True): + new = cls() + for key, value in d.items(): + if isinstance(key, str) and C: + new[key.replace(r'\:', ':').replace('\\\\', '\\')] = Transform.from_dict(value) + elif T: + new[key] = Transform.from_dict(value) + return new + + @classmethod + def from_shifts(cls, shifts): + new = cls() + for key, shift in shifts.items(): + new[key] = Transform.from_shift(shift) + return new + + def __mul__(self, other): + new = Transforms() + if isinstance(other, Transforms): + for key0, value0 in self.items(): + for key1, value1 in other.items(): + new[key0 + key1] = value0 * value1 + return new + elif other is None: + return self + else: + for key in self.keys(): + new[key] = self[key] * other + return new + + def asdict(self): + return {key.replace('\\', '\\\\').replace(':', r'\:') if isinstance(key, str) else key: value.asdict() + for key, value in self.items()} + + def __getitem__(self, item): + return np.prod([self[i] for i in item[::-1]]) if isinstance(item, tuple) else super().__getitem__(item) + + def __missing__(self, key): + return self.default + + def __getstate__(self): + return self.__dict__ + + def __setstate__(self, state): + self.__dict__.update(state) + + def __hash__(self): + return hash(frozenset((*self.__dict__.items(), *self.items()))) + + def save(self, file): + with open(Path(file).with_suffix('.yml'), 'w') as f: + yaml.safe_dump(self.asdict(), f, default_flow_style=None) + + def copy(self): + return deepcopy(self) + + def adapt(self, origin, shape, channel_names): + def key_map(a, b): + def fun(b, key_a): + for key_b in b: + if key_b in key_a or key_a in key_b: + return key_a, key_b + + return {n[0]: n[1] for key_a in a if (n := fun(b, key_a))} + + for value in self.values(): + value.adapt(origin, shape) + self.default.adapt(origin, shape) + transform_channels = {key for key in self.keys() if isinstance(key, str)} + if set(channel_names) - transform_channels: + mapping = key_map(channel_names, transform_channels) + warnings.warn(f'The image file and the transform do not have the same channels,' + f' creating a mapping: {mapping}') + for key_im, key_t in mapping.items(): + self[key_im] = self[key_t] + + @property + def inverse(self): + # TODO: check for C@T + inverse = self.copy() + for key, value in self.items(): + inverse[key] = value.inverse + return inverse + + def coords_pandas(self, array, channel_names, columns=None): + if isinstance(array, DataFrame): + return concat([self.coords_pandas(row, channel_names, columns) for _, row in array.iterrows()], axis=1).T + elif isinstance(array, Series): + key = [] + if 'C' in array: + key.append(channel_names[int(array['C'])]) + if 'T' in array: + key.append(int(array['T'])) + return self[tuple(key)].coords(array, columns) + else: + raise TypeError('Not a pandas DataFrame or Series.') + + def with_beads(self, cyllens, bead_files): + assert len(bead_files) > 0, 'At least one file is needed to calculate the registration.' + transforms = [self.calculate_channel_transforms(file, cyllens) for file in bead_files] + for key in {key for transform in transforms for key in transform.keys()}: + new_transforms = [transform[key] for transform in transforms if key in transform] + if len(new_transforms) == 1: + self[key] = new_transforms[0] + else: + self[key] = Transform() + self[key].parameters = np.mean([t.parameters for t in new_transforms], 0) + self[key].dparameters = (np.std([t.parameters for t in new_transforms], 0) / + np.sqrt(len(new_transforms))).tolist() + return self + + @staticmethod + def get_bead_files(path): + from . import Imread + files = [] + for file in path.iterdir(): + if file.name.lower().startswith('beads'): + try: + with Imread(file): + files.append(file) + except Exception: + pass + files = sorted(files) + if not files: + raise Exception('No bead file found!') + checked_files = [] + for file in files: + try: + if file.is_dir(): + file /= 'Pos0' + with Imread(file): # check for errors opening the file + checked_files.append(file) + except (Exception,): + continue + if not checked_files: + raise Exception('No bead file found!') + return checked_files + + @staticmethod + def calculate_channel_transforms(bead_file, cyllens): + """ When no channel is not transformed by a cylindrical lens, assume that the image is scaled by a factor 1.162 + in the horizontal direction """ + from . import Imread + + with Imread(bead_file, axes='zcyx') as im: # noqa + max_ims = im.max('z') + goodch = [c for c, max_im in enumerate(max_ims) if not im.is_noise(max_im)] + if not goodch: + goodch = list(range(len(max_ims))) + untransformed = [c for c in range(im.shape['c']) if cyllens[im.detector[c]].lower() == 'none'] + + good_and_untrans = sorted(set(goodch) & set(untransformed)) + if good_and_untrans: + masterch = good_and_untrans[0] + else: + masterch = goodch[0] + transform = Transform() + if not good_and_untrans: + matrix = transform.matrix + matrix[0, 0] = 0.86 + transform.matrix = matrix + transforms = Transforms() + for c in tqdm(goodch, desc='Calculating channel transforms'): # noqa + if c == masterch: + transforms[im.channel_names[c]] = transform + else: + transforms[im.channel_names[c]] = Transform.register(max_ims[masterch], max_ims[c]) * transform + return transforms + + @staticmethod + def save_channel_transform_tiff(bead_files, tiffile): + from . import Imread + n_channels = 0 + for file in bead_files: + with Imread(file) as im: + n_channels = max(n_channels, im.shape['c']) + with IJTiffFile(tiffile) as tif: + for t, file in enumerate(bead_files): + with Imread(file) as im: + with Imread(file).with_transform() as jm: + for c in range(im.shape['c']): + tif.save(np.hstack((im(c=c, t=0).max('z'), jm(c=c, t=0).max('z'))), c, 0, t) + + def with_drift(self, im): + """ Calculate shifts relative to the first frame + divide the sequence into groups, + compare each frame to the frame in the middle of the group and compare these middle frames to each other + """ + im = im.transpose('tzycx') + t_groups = [list(chunk) for chunk in Chunks(range(im.shape['t']), size=round(np.sqrt(im.shape['t'])))] + t_keys = [int(np.round(np.mean(t_group))) for t_group in t_groups] + t_pairs = [(int(np.round(np.mean(t_group))), frame) for t_group in t_groups for frame in t_group] + t_pairs.extend(zip(t_keys, t_keys[1:])) + fmaxz_keys = {t_key: filters.gaussian(im[t_key].max('z'), 5) for t_key in t_keys} + + def fun(t_key_t, im, fmaxz_keys): + t_key, t = t_key_t + if t_key == t: + return 0, 0 + else: + fmaxz = filters.gaussian(im[t].max('z'), 5) + return Transform.register(fmaxz_keys[t_key], fmaxz, 'translation').parameters[4:] + + shifts = np.array(pmap(fun, t_pairs, (im, fmaxz_keys), desc='Calculating image shifts.')) + shift_keys_cum = np.zeros(2) + for shift_keys, t_group in zip(np.vstack((-shifts[0], shifts[im.shape['t']:])), t_groups): + shift_keys_cum += shift_keys + shifts[t_group] += shift_keys_cum + + for i, shift in enumerate(shifts[:im.shape['t']]): + self[i] = Transform.from_shift(shift) + return self + + +class Transform: + def __init__(self): + if sitk is None: + self.transform = None + else: + self.transform = sitk.ReadTransform(str(Path(__file__).parent / 'transform.txt')) + self.dparameters = [0., 0., 0., 0., 0., 0.] + self.shape = [512., 512.] + self.origin = [255.5, 255.5] + self._last, self._inverse = None, None + + def __reduce__(self): + return self.from_dict, (self.asdict(),) + + def __repr__(self): + return self.asdict().__repr__() + + def __str__(self): + return self.asdict().__str__() + + @classmethod + def register(cls, fix, mov, kind=None): + """ kind: 'affine', 'translation', 'rigid' """ + if sitk is None: + raise ImportError('SimpleElastix is not installed: ' + 'https://simpleelastix.readthedocs.io/GettingStarted.html') + new = cls() + kind = kind or 'affine' + new.shape = fix.shape + fix, mov = new.cast_image(fix), new.cast_image(mov) + # TODO: implement RigidTransform + tfilter = sitk.ElastixImageFilter() + tfilter.LogToConsoleOff() + tfilter.SetFixedImage(fix) + tfilter.SetMovingImage(mov) + tfilter.SetParameterMap(sitk.GetDefaultParameterMap(kind)) + tfilter.Execute() + transform = tfilter.GetTransformParameterMap()[0] + if kind == 'affine': + new.parameters = [float(t) for t in transform['TransformParameters']] + new.shape = [float(t) for t in transform['Size']] + new.origin = [float(t) for t in transform['CenterOfRotationPoint']] + elif kind == 'translation': + new.parameters = [1.0, 0.0, 0.0, 1.0] + [float(t) for t in transform['TransformParameters']] + new.shape = [float(t) for t in transform['Size']] + new.origin = [(t - 1) / 2 for t in new.shape] + else: + raise NotImplementedError(f'{kind} tranforms not implemented (yet)') + new.dparameters = 6 * [np.nan] + return new + + @classmethod + def from_shift(cls, shift): + return cls.from_array(np.array(((1, 0, shift[0]), (0, 1, shift[1]), (0, 0, 1)))) + + @classmethod + def from_array(cls, array): + new = cls() + new.matrix = array + return new + + @classmethod + def from_file(cls, file): + with open(Path(file).with_suffix('.yml')) as f: + return cls.from_dict(yamlload(f)) + + @classmethod + def from_dict(cls, d): + new = cls() + new.origin = None if d['CenterOfRotationPoint'] is None else [float(i) for i in d['CenterOfRotationPoint']] + new.parameters = ((1., 0., 0., 1., 0., 0.) if d['TransformParameters'] is None else + [float(i) for i in d['TransformParameters']]) + new.dparameters = ([(0., 0., 0., 0., 0., 0.) if i is None else float(i) for i in d['dTransformParameters']] + if 'dTransformParameters' in d else 6 * [np.nan] and d['dTransformParameters'] is not None) + new.shape = None if d['Size'] is None else [None if i is None else float(i) for i in d['Size']] + return new + + def __mul__(self, other): # TODO: take care of dmatrix + result = self.copy() + if isinstance(other, Transform): + result.matrix = self.matrix @ other.matrix + result.dmatrix = self.dmatrix @ other.matrix + self.matrix @ other.dmatrix + else: + result.matrix = self.matrix @ other + result.dmatrix = self.dmatrix @ other + return result + + def is_unity(self): + return self.parameters == [1, 0, 0, 1, 0, 0] + + def copy(self): + return deepcopy(self) + + @staticmethod + def cast_image(im): + if not isinstance(im, sitk.Image): + im = sitk.GetImageFromArray(np.asarray(im)) + return im + + @staticmethod + def cast_array(im): + if isinstance(im, sitk.Image): + im = sitk.GetArrayFromImage(im) + return im + + @property + def matrix(self): + return np.array(((*self.parameters[:2], self.parameters[4]), + (*self.parameters[2:4], self.parameters[5]), + (0, 0, 1))) + + @matrix.setter + def matrix(self, value): + value = np.asarray(value) + self.parameters = [*value[0, :2], *value[1, :2], *value[:2, 2]] + + @property + def dmatrix(self): + return np.array(((*self.dparameters[:2], self.dparameters[4]), + (*self.dparameters[2:4], self.dparameters[5]), + (0, 0, 0))) + + @dmatrix.setter + def dmatrix(self, value): + value = np.asarray(value) + self.dparameters = [*value[0, :2], *value[1, :2], *value[:2, 2]] + + @property + def parameters(self): + if self.transform is not None: + return list(self.transform.GetParameters()) + else: + return [1., 0., 0., 1., 0., 0.] + + @parameters.setter + def parameters(self, value): + if self.transform is not None: + value = np.asarray(value) + self.transform.SetParameters(value.tolist()) + + @property + def origin(self): + if self.transform is not None: + return self.transform.GetFixedParameters() + + @origin.setter + def origin(self, value): + if self.transform is not None: + value = np.asarray(value) + self.transform.SetFixedParameters(value.tolist()) + + @property + def inverse(self): + if self.is_unity(): + return self + if self._last is None or self._last != self.asdict(): + self._last = self.asdict() + self._inverse = Transform.from_dict(self.asdict()) + self._inverse.transform = self._inverse.transform.GetInverse() + self._inverse._last = self._inverse.asdict() + self._inverse._inverse = self + return self._inverse + + def adapt(self, origin, shape): + self.origin -= np.array(origin) + (self.shape - np.array(shape)[:2]) / 2 + self.shape = shape[:2] + + def asdict(self): + return {'CenterOfRotationPoint': self.origin, 'Size': self.shape, 'TransformParameters': self.parameters, + 'dTransformParameters': np.nan_to_num(self.dparameters, nan=1e99).tolist()} + + def frame(self, im, default=0): + if self.is_unity(): + return im + else: + if sitk is None: + raise ImportError('SimpleElastix is not installed: ' + 'https://simpleelastix.readthedocs.io/GettingStarted.html') + dtype = im.dtype + im = im.astype('float') + intp = sitk.sitkBSpline if np.issubdtype(dtype, np.floating) else sitk.sitkNearestNeighbor + return self.cast_array(sitk.Resample(self.cast_image(im), self.transform, intp, default)).astype(dtype) + + def coords(self, array, columns=None): + """ Transform coordinates in 2 column numpy array, + or in pandas DataFrame or Series objects in columns ['x', 'y'] + """ + if self.is_unity(): + return array.copy() + elif DataFrame is not None and isinstance(array, (DataFrame, Series)): + columns = columns or ['x', 'y'] + array = array.copy() + if isinstance(array, DataFrame): + array[columns] = self.coords(np.atleast_2d(array[columns].to_numpy())) + elif isinstance(array, Series): + array[columns] = self.coords(np.atleast_2d(array[columns].to_numpy()))[0] + return array + else: # somehow we need to use the inverse here to get the same effect as when using self.frame + return np.array([self.inverse.transform.TransformPoint(i.tolist()) for i in np.asarray(array)]) + + def save(self, file): + """ save the parameters of the transform calculated + with affine_registration to a yaml file + """ + if not file[-3:] == 'yml': + file += '.yml' + with open(file, 'w') as f: + yaml.safe_dump(self.asdict(), f, default_flow_style=None) diff --git a/pyproject.toml b/pyproject.toml index 1c91e01..cfee4b5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,9 +1,9 @@ [build-system] -requires = ["maturin>=1.8,<2.0"] +requires = ["maturin>=1.8.2"] build-backend = "maturin" [project] -name = "ndbioimage_rs" +name = "ndbioimage" requires-python = ">=3.10" classifiers = [ "Programming Language :: Rust", @@ -11,5 +11,12 @@ classifiers = [ "Programming Language :: Python :: Implementation :: PyPy", ] dynamic = ["version"] + [tool.maturin] -features = ["pyo3/extension-module"] +python-source = "py" +features = ["pyo3/extension-module", "python", "gpl-formats"] +module-name = "ndbioimage.ndbioimage_rs" +exclude = ["py/ndbioimage/jassets/*", "py/ndbioimage/deps/*"] + +[tool.isort] +line_length = 119 diff --git a/src/bioformats.rs b/src/bioformats.rs index a9e5c17..9f705f6 100644 --- a/src/bioformats.rs +++ b/src/bioformats.rs @@ -10,11 +10,51 @@ thread_local! { /// Ensure 1 jvm per thread fn jvm() -> Rc { JVM.with(|cell| { - cell.get_or_init(move || Rc::new(JvmBuilder::new().build().expect("Failed to build JVM"))) - .clone() + cell.get_or_init(move || { + #[cfg(feature = "python")] + let path = crate::py::ndbioimage_file().unwrap(); + + #[cfg(not(feature = "python"))] + let path = std::env::current_exe() + .unwrap() + .parent() + .unwrap() + .to_path_buf(); + + let class_path = path.parent().unwrap(); + Rc::new( + JvmBuilder::new() + .with_base_path(class_path.to_str().unwrap()) + .build() + .expect("Failed to build JVM"), + ) + }) + .clone() }) } +#[cfg(feature = "python")] +pub(crate) fn download_bioformats(gpl_formats: bool) -> Result<()> { + let path = crate::py::ndbioimage_file().unwrap(); + let class_path = path.parent().unwrap(); + let jvm = JvmBuilder::new() + .with_base_path(class_path.to_str().unwrap()) + .with_maven_settings(j4rs::MavenSettings::new(vec![ + j4rs::MavenArtifactRepo::from( + "openmicroscopy::https://artifacts.openmicroscopy.org/artifactory/ome.releases", + ), + ])) + .build()?; + + jvm.deploy_artifact(&j4rs::MavenArtifact::from("ome:bioformats_package:8.1.0"))?; + + if gpl_formats { + jvm.deploy_artifact(&j4rs::MavenArtifact::from("ome:formats-gpl:8.1.0"))?; + } + + Ok(()) +} + macro_rules! method_return { ($R:ty$(|c)?) => { Result<$R> }; () => { Result<()> }; @@ -113,7 +153,7 @@ impl ImageReader { Ok(jvm() .chain(&mds)? .cast("loci.formats.ome.OMEPyramidStore")? - .invoke("dumpXML", &[])? + .invoke("dumpXML", InvocationArg::empty())? .to_rust()?) } diff --git a/src/lib.rs b/src/lib.rs index 21464b8..71f1c24 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -215,7 +215,7 @@ impl Reader { } /// Get ome metadata as xml string - pub fn ome_xml(&self) -> Result { + pub fn get_ome_xml(&self) -> Result { self.image_reader.ome_xml() } @@ -396,7 +396,7 @@ mod tests { fn ome_xml() -> Result<()> { let file = "Experiment-2029.czi"; let reader = open(file)?; - let xml = reader.ome_xml()?; + let xml = reader.get_ome_xml()?; println!("{}", xml); Ok(()) } diff --git a/src/py.rs b/src/py.rs index abaca35..db73466 100644 --- a/src/py.rs +++ b/src/py.rs @@ -1,7 +1,7 @@ +use crate::bioformats::download_bioformats; use crate::{Frame, Reader}; -use numpy::{IntoPyArray, PyArrayMethods, ToPyArray}; +use numpy::ToPyArray; use pyo3::prelude::*; -use pyo3::BoundObject; use std::path::PathBuf; #[pyclass(subclass)] @@ -41,11 +41,35 @@ impl PyReader { Frame::DOUBLE(arr) => arr.to_pyarray(py).into_any(), }) } + + fn get_ome_xml(&self) -> PyResult { + let reader = Reader::new(&self.path, self.series)?; // TODO: prevent making a new Reader each time + Ok(reader.get_ome_xml()?) + } +} + +pub(crate) fn ndbioimage_file() -> anyhow::Result { + let file = Python::with_gil(|py| { + py.import("ndbioimage") + .unwrap() + .filename() + .unwrap() + .to_string() + }); + Ok(PathBuf::from(file)) } #[pymodule] #[pyo3(name = "ndbioimage_rs")] fn ndbioimage_rs(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_class::()?; + + #[pyfn(m)] + #[pyo3(name = "download_bioformats")] + fn py_download_bioformats(gpl_formats: bool) -> PyResult<()> { + download_bioformats(gpl_formats)?; + Ok(()) + } + Ok(()) }