From 3dc8e6af045a07548039a6ab34d23845f9835ac8 Mon Sep 17 00:00:00 2001 From: Wim Pomp Date: Thu, 21 Aug 2025 19:45:02 +0200 Subject: [PATCH] - bump bioformats to 8.3.0 - rust: command line binary, save as mp4, save as tiff, ome metadata, more methods for View, bugfixes, less unsafe code - python: ome as dict --- Cargo.toml | 39 +- build.rs | 69 +-- py/ndbioimage/__init__.py | 1064 +++---------------------------------- py/ndbioimage/ome.py | 39 ++ pyproject.toml | 9 +- src/axes.rs | 24 +- src/bioformats.rs | 38 +- src/colors.rs | 207 ++++++++ src/lib.rs | 126 ++++- src/main.rs | 110 ++++ src/metadata.rs | 372 +++++++++++++ src/movie.rs | 258 +++++++++ src/py.rs | 375 +++++++++++-- src/reader.rs | 51 +- src/stats.rs | 10 +- src/tiff.rs | 193 +++++++ src/view.rs | 594 ++++++++++++++++----- 17 files changed, 2317 insertions(+), 1261 deletions(-) create mode 100644 py/ndbioimage/ome.py create mode 100644 src/colors.rs create mode 100644 src/main.rs create mode 100644 src/metadata.rs create mode 100644 src/movie.rs create mode 100644 src/tiff.rs diff --git a/Cargo.toml b/Cargo.toml index bdd842e..7332247 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,8 +1,8 @@ [package] name = "ndbioimage" -version = "2025.4.1" -edition = "2021" -rust-version = "1.78.0" +version = "2025.8.0" +edition = "2024" +rust-version = "1.85.1" authors = ["Wim Pomp "] license = "MIT" description = "Read bio image formats using the bio-formats java package." @@ -19,20 +19,28 @@ name = "ndbioimage" crate-type = ["cdylib", "rlib"] [dependencies] -anyhow = "1.0.98" +anyhow = { version = "1.0.99", features = ["backtrace"] } +clap = { version = "4.5.45", features = ["derive"] } +ffmpeg-sidecar = { version = "2.1.0", optional = true } itertools = "0.14.0" -indexmap = { version = "2.9.0", features = ["serde"] } +indexmap = { version = "2.0.0", features = ["serde"] } +indicatif = { version = "0.18.0", features = ["rayon"], optional = true } j4rs = "0.22.0" -ndarray = { version = "0.16.1", features = ["serde"] } +ndarray = { version = "0.16.1", features = ["serde"] } num = "0.4.3" -numpy = { version = "0.24.0", optional = true } +numpy = { version = "0.25.0", optional = true } +ordered-float = "5.0.0" +rayon = { version = "1.11.0", optional = true } serde = { version = "1.0.219", features = ["rc"] } -serde_json = { version = "1.0.140", optional = true } +serde_json = { version = "1.0.143", optional = true } serde_with = "3.12.0" -thread_local = "1.1.8" +tiffwrite = { version = "2025.5.0", optional = true} +thread_local = "1.1.9" +ome-metadata = "0.2.2" +lazy_static = "1.5.0" [dependencies.pyo3] -version = "0.24.2" +version = "0.25.1" features = ["extension-module", "abi3-py310", "generate-import-lib", "anyhow"] optional = true @@ -40,10 +48,17 @@ optional = true rayon = "1.10.0" [build-dependencies] -anyhow = "1.0.98" +anyhow = "1.0.99" j4rs = "0.22.0" +ffmpeg-sidecar = "2.1.0" retry = "2.1.0" [features] +# Enables formats for which code in bioformats with a GPL license is needed +gpl-formats = [] +# Enables python ffi using pyO3 python = ["dep:pyo3", "dep:numpy", "dep:serde_json"] -gpl-formats = [] \ No newline at end of file +# Enables writing as tiff +tiff = ["dep:tiffwrite", "dep:indicatif", "dep:rayon"] +# Enables writing as mp4 using ffmpeg +movie = ["dep:ffmpeg-sidecar"] \ No newline at end of file diff --git a/build.rs b/build.rs index 99f9ae9..364ed0f 100644 --- a/build.rs +++ b/build.rs @@ -1,5 +1,5 @@ #[cfg(not(feature = "python"))] -use j4rs::{errors::J4RsError, JvmBuilder, MavenArtifact, MavenArtifactRepo, MavenSettings}; +use j4rs::{JvmBuilder, MavenArtifact, MavenArtifactRepo, MavenSettings, errors::J4RsError}; #[cfg(not(feature = "python"))] use retry::{delay, delay::Exponential, retry}; @@ -7,44 +7,50 @@ use retry::{delay, delay::Exponential, retry}; #[cfg(feature = "python")] use j4rs::Jvm; +#[cfg(feature = "movie")] +use ffmpeg_sidecar::download::auto_download; + fn main() -> anyhow::Result<()> { println!("cargo::rerun-if-changed=build.rs"); - #[cfg(not(feature = "python"))] if std::env::var("DOCS_RS").is_err() { + #[cfg(feature = "movie")] + auto_download()?; + + #[cfg(not(feature = "python"))] 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}_"))?; + #[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)?; } - } - // 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)?; + 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)?; + } } } } @@ -55,15 +61,16 @@ fn main() -> anyhow::Result<()> { #[cfg(not(feature = "python"))] fn deploy_java_artifacts() -> Result<(), J4RsError> { let jvm = JvmBuilder::new() + .skip_setting_native_lib() .with_maven_settings(MavenSettings::new(vec![MavenArtifactRepo::from( "openmicroscopy::https://artifacts.openmicroscopy.org/artifactory/ome.releases", )])) .build()?; - jvm.deploy_artifact(&MavenArtifact::from("ome:bioformats_package:8.1.0"))?; + jvm.deploy_artifact(&MavenArtifact::from("ome:bioformats_package:8.3.0"))?; #[cfg(feature = "gpl-formats")] - jvm.deploy_artifact(&MavenArtifact::from("ome:formats-gpl:8.1.0"))?; + jvm.deploy_artifact(&MavenArtifact::from("ome:formats-gpl:8.3.0"))?; Ok(()) } diff --git a/py/ndbioimage/__init__.py b/py/ndbioimage/__init__.py index 0f20518..37d9fc4 100755 --- a/py/ndbioimage/__init__.py +++ b/py/ndbioimage/__init__.py @@ -1,27 +1,23 @@ from __future__ import annotations import os -import re import warnings from abc import ABC 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, Optional, Sequence, TypeVar import numpy as np -import yaml from numpy.typing import ArrayLike, DTypeLike -from ome_types import OME, from_xml, ureg -from pint import set_application_registry from tiffwrite import FrameInfo, IJTiffParallel from tqdm.auto import tqdm +from ome_metadata import Ome +from ome_metadata.ome_metadata_rs import Length # noqa from . import ndbioimage_rs as rs # noqa from .transforms import Transform, Transforms # noqa: F401 @@ -38,11 +34,6 @@ try: except Exception: # noqa __git_commit_hash__ = "unknown" -if hasattr(ureg, "formatter"): - ureg.formatter.default_format = "~P" -else: - ureg.default_format = "~P" -set_application_registry(ureg) warnings.filterwarnings("ignore", "Reference to unknown ID") Number = int | float | np.integer | np.floating @@ -141,13 +132,13 @@ class OmeCache(DequeDict): def __reduce__(self) -> tuple[type, tuple]: return self.__class__, () - def __getitem__(self, path: Path | str | tuple) -> OME: + 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: + def __setitem__(self, path: Path | str | tuple, value: Ome) -> None: if isinstance(path, tuple): super().__setitem__(path, value) else: @@ -180,7 +171,7 @@ def get_positions(path: str | Path) -> Optional[list[int]]: # noqa return None -class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): +class Imread(rs.View, 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 @@ -200,7 +191,7 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): << 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 + >> im.pxsize << 0.09708737864077668 image-plane pixel size in um >> im.laserwavelengths << [642, 488] @@ -208,292 +199,60 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): << [0.02, 0.0005] in % See __init__ and other functions for more ideas. + + # TODO: argmax, argmin, nanmax, nanmin, nanmean, nansum, nanstd, nanvar, std, var, squeeze """ - ureg = ureg + def __getitem__(self, item): + new = super().__getitem__(item) + return Imread(new) if isinstance(new, rs.View) else new - def close(self) -> None: - self.reader.close() + def __copy__(self): + Imread(super().__copy__()) + + def copy(self): + Imread(super().copy()) + + def astype(self): + Imread(super().astype()) + + def squeeze(self): + new = super().squeeze() + return Imread(new) if isinstance(new, rs.View) else new + + def min(self, *args, **kwargs) -> Imread | float: + new = super().min(*args, **kwargs) + return Imread(new) if isinstance(new, rs.View) else new + + def max(self, *args, **kwargs) -> Imread | float: + new = super().max(*args, **kwargs) + return Imread(new) if isinstance(new, rs.View) else new + + def mean(self, *args, **kwargs) -> Imread | float: + new = super().mean(*args, **kwargs) + return Imread(new) if isinstance(new, rs.View) else new + + def sum(self, *args, **kwargs) -> Imread | float: + new = super().sum(*args, **kwargs) + return Imread(new) if isinstance(new, rs.View) else new + + def transpose(self, *args, **kwargs) -> Imread | float: + new = super().transpose(*args, **kwargs) + return Imread(new) if isinstance(new, rs.View) else new + + def swap_axes(self, *args, **kwargs) -> Imread | float: + new = super().swap_axes(*args, **kwargs) + return Imread(new) if isinstance(new, rs.View) else new + + @property + def T(self) -> Imread | float: + return Imread(super().T) @staticmethod def get_positions(path: str | Path) -> Optional[list[int]]: # noqa # TODO return None - def __init__( - self, path: Path | str = None, dtype: DTypeLike = None, axes: str = None - ) -> None: - super().__init__() - if path is not None: - path = Path(path) - path, series = self.split_path_series(path) - self.reader = rs.View(str(path), int(series)) - 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.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.frameoffset = 0, 0 # how far apart the centers of frame and sensor are - - # 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.reader.dtype = str(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 - - 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 __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 be made with slice() or np.s_""" - return np.asarray( - 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: - raise NotImplementedError - - def __enter__(self) -> Imread: - return self - - def __exit__(self, *args: Any, **kwargs: Any) -> None: - pass - - def __getitem__( - self, - n: int - | Sequence[int] - | Sequence[slice] - | slice - | type(Ellipsis) - | dict[str, int | Sequence[int] | Sequence[slice] | slice | type(Ellipsis)], - ) -> Number | Imread: - """slice like a numpy array but return an Imread instance""" - item = self.reader.__getitem__(n) - if isinstance(item, rs.View): - new = self.new() - new.reader = item - return new - else: - return item - - def __len__(self) -> int: - return self.shape[0] - - def __repr__(self) -> str: - return self.summary - - def __str__(self) -> str: - return str(self.path) - @staticmethod def as_axis(axis): if axis is None: @@ -503,532 +262,10 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): else: return str(axis) - # @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, *_, **__) -> np.ndarray: - if dtype is not None and self.dtype != np.dtype(dtype): - reader = self.reader.as_type(str(np.dtype(dtype))) - else: - reader = self.reader - return reader.as_array() - - 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 - axes = "".join(i for i in self.axes if i in "yx") + "czt" - 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 - elif value == new_value: - a = np.ravel_multi_index( - [arg[axes.find(i)] for i in self.axes], self.shape - ) - b = np.ravel_multi_index( - [(new_arg + idx)[axes.find(i)] for i in self.axes], self.shape - ) - if b < a: - arg = new_arg + idx - else: - i = fun((value, new_value)) # type: ignore - arg = (arg, new_arg + idx)[i] - value = (value, new_value)[i] - arg = np.ravel_multi_index( - [arg[axes.find(i)] for i in self.axes], self.shape - ) - if out is None: - return arg - else: - out[:] = 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[:] = 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 - - def get_frame(self, c: int = 0, z: int = 0, t: int = 0) -> np.ndarray: - """retrieve a single frame at czt, sliced accordingly""" - return self.reader.get_frame(int(c), int(z), int(t)) - - @property - def path(self) -> Path: - return Path(self.reader.path) - - @property - def series(self) -> int: - return self.reader.series - - @property - def axes(self) -> str: - return "".join(self.reader.axes).lower() - - @axes.setter - def axes(self, value: str) -> None: - value = value.lower() - self.reader = self.reader.__getitem__( - [slice(None) if ax in value else 0 for ax in self.axes] - ).transpose([ax for ax in value.upper()]) - - @property - def dtype(self) -> np.dtype: - return np.dtype(self.reader.dtype.lower()) - - @dtype.setter - def dtype(self, value: DTypeLike) -> None: - self.reader.dtype = str(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 self.reader.ndim - - @property - def size(self) -> int: - return self.reader.size - - @property - def shape(self) -> Shape: - return Shape(self.reader.shape, "".join(self.reader.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}"] - 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, *_, **__) -> Number | Imread: - reader = self.reader.max(self.as_axis(axis)) - if isinstance(reader, rs.View): - new = self.new() - new.reader = reader - return new - else: - return reader - - @wraps(np.ndarray.mean) - def mean(self, axis=None, *_, **__) -> Number | Imread: - reader = self.reader.mean(self.as_axis(axis)) - if isinstance(reader, rs.View): - new = self.new() - new.reader = reader - return new - else: - return reader - - @wraps(np.ndarray.min) - def min(self, axis=None, *_, **__) -> Number | Imread: - reader = self.reader.min(self.as_axis(axis)) - if isinstance(reader, rs.View): - new = self.new() - new.reader = reader - return new - else: - return reader - - @wraps(np.ndarray.sum) - def sum(self, axis=None, *_, **__) -> Number | Imread: - reader = self.reader.sum(self.as_axis(axis)) - if isinstance(reader, rs.View): - new = self.new() - new.reader = reader - return new - else: - return reader - @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) -> np.ndarray: return np.asarray(self).flatten(*args, **kwargs) @@ -1037,167 +274,17 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): def reshape(self, *args, **kwargs) -> np.ndarray: return np.asarray(self).reshape(*args, **kwargs) # noqa - @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.swapaxes) - def swapaxes(self, axis1, axis2): - axis1 = self.as_axis(axis1) - axis2 = self.as_axis(axis2) - new = self.new() - new.reader = self.reader.swap_axes(axis1, axis2) - return new - - @wraps(np.ndarray.transpose) - def transpose(self, *axes): - axes = [self.as_axis(axis) for axis in axes] - new = self.new() - new.reader = self.reader.transpose(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: + def as_array(self) -> np.ndarray: return self.__array__() @wraps(np.ndarray.astype) - def astype(self, dtype, *_, **__): - new = self.copy() - new.dtype = dtype - return new - - def copy(self) -> Imread: - new = self.new() - new.reader = self.reader.copy() - return new - - def new(self) -> Imread: - new = Imread.__new__(Imread) - new.__dict__.update( - {key: value for key, value in self.__dict__.items() if key != "reader"} - ) - return new - - 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] + def astype(self, dtype: DTypeLike, *_, **__) -> Imread: + return Imread(super().astype(str(np.dtype(dtype)))) @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 fix_ome(ome: OME) -> OME: + def fix_ome(ome: Ome) -> Ome: # fix ome if necessary - for image in ome.images: + for image in ome.image: try: if ( image.pixels.physical_size_z is None @@ -1206,10 +293,7 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): z = np.array( [ ( - plane.position_z - * ureg.Quantity(plane.position_z_unit.value) - .to(ureg.m) - .magnitude, + plane.position_z_unit.convert("um", plane.position_z), plane.the_z, ) for plane in image.pixels.planes @@ -1220,24 +304,24 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): 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 + image.pixels.physical_size_z_unit = Length("um") # type: ignore except Exception: # noqa pass return ome @staticmethod - def read_ome(path: str | Path) -> Optional[OME]: + def read_ome(path: str | Path) -> Optional[Ome]: path = Path(path) # type: ignore if path.with_suffix(".ome.xml").exists(): - return OME.from_xml(path.with_suffix(".ome.xml")) + return Ome.from_xml(path.with_suffix(".ome.xml").read_text()) return None - def get_ome(self) -> OME: + def get_ome(self) -> Ome: """OME metadata structure""" - return from_xml(self.reader.get_ome_xml(), validate=False) + return Ome.from_xml(self.get_ome_xml()) @cached_property - def ome(self) -> OME: + def ome(self) -> Ome: cache = OmeCache() if self.path not in cache: ome = self.read_ome(self.path) @@ -1293,7 +377,11 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): return a, b def cframe( - frame: ArrayLike, color: str, a: float, b: float, scale: float = 1 # noqa + frame: ArrayLike, + color: str, + a: float, + b: float, + scale: float = 1, # noqa ) -> np.ndarray: color = to_rgb(color) frame = (frame - a) / (b - a) @@ -1388,8 +476,8 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): self, fname.with_suffix(".tif"), dtype=pixel_type, - pxsize=self.pxsize_um, - deltaz=self.deltaz_um, + pxsize=self.pxsize, + deltaz=self.deltaz, **kwargs, ) as tif: for i, m in tqdm( # noqa @@ -1445,18 +533,6 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): # view.transform.adapt(view.frameoffset, view.shape.yxczt, view.channel_names) # return view - @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 main() -> None: parser = ArgumentParser(description="Display info and save as tif") diff --git a/py/ndbioimage/ome.py b/py/ndbioimage/ome.py new file mode 100644 index 0000000..f21b7d4 --- /dev/null +++ b/py/ndbioimage/ome.py @@ -0,0 +1,39 @@ +from __future__ import annotations + +from ome_metadata import ome_metadata_rs as rs # noqa +from collections import UserDict, UserList + + +class Ome(UserDict): + @staticmethod + def from_xml(xml: str) -> Ome: + """Create the OME structure from an XML string""" + new = Ome() + new.update(rs.ome(str(xml))) + return new + + def __dir__(self) -> list[str]: + return list(self.keys()) + list(super().__dir__()) + + def __getattr__(self, key: str) -> Ome | OmeList | int | float | str: + try: + new = self.__getitem__(key) + except KeyError: + raise AttributeError(f"'Ome' object has no attribute '{key}'") + if isinstance(new, dict): + return Ome(**new) + elif isinstance(new, list): + return OmeList(new) + else: + return new + + +class OmeList(UserList): + def __getitem__(self, item: int) -> Ome | OmeList | int | float | str: + new = super().__getitem__(item) + if isinstance(new, dict): + return Ome(**new) + elif isinstance(new, list): + return OmeList(new) + else: + return new diff --git a/pyproject.toml b/pyproject.toml index 28f84aa..694914a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,16 +27,9 @@ classifiers = [ dynamic = ["version"] dependencies = [ "numpy", - "pandas", "tiffwrite", - "ome-types", - "pint", + "ome-metadata >= 0.2.1", "tqdm", - "lxml", - "pyyaml", - "parfor", - "SimpleITK-SimpleElastix", - "scikit-image", ] [project.optional-dependencies] diff --git a/src/axes.rs b/src/axes.rs index 03e2575..ab48b97 100644 --- a/src/axes.rs +++ b/src/axes.rs @@ -1,8 +1,9 @@ use crate::stats::MinMax; -use anyhow::{anyhow, Error, Result}; +use anyhow::{Error, Result, anyhow}; use ndarray::{Array, Dimension, Ix2, SliceInfo, SliceInfoElem}; use serde::{Deserialize, Deserializer, Serialize, Serializer}; use serde_with::{DeserializeAs, SerializeAs}; +use std::fmt::{Display, Formatter}; use std::hash::{Hash, Hasher}; use std::str::FromStr; @@ -51,6 +52,20 @@ impl FromStr for Axis { } } +impl Display for Axis { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + let s = match self { + Axis::C => "C", + Axis::Z => "Z", + Axis::T => "T", + Axis::Y => "Y", + Axis::X => "X", + Axis::New => "N", + }; + write!(f, "{}", s) + } +} + impl Ax for Axis { fn n(&self) -> usize { *self as usize @@ -60,7 +75,10 @@ impl Ax for Axis { if let Some(pos) = axes.iter().position(|a| a == self) { Ok(pos) } else { - Err(Error::msg("Axis not found in axes")) + Err(Error::msg(format!( + "Axis {:?} not found in axes {:?}", + self, axes + ))) } } @@ -96,7 +114,7 @@ impl Ax for usize { } }) .collect(); - assert!(*self < idx.len(), "self: {}, idx: {:?}", self, idx); + debug_assert!(*self < idx.len(), "self: {}, idx: {:?}", self, idx); Ok(idx[*self]) } } diff --git a/src/bioformats.rs b/src/bioformats.rs index b83e096..44442c6 100644 --- a/src/bioformats.rs +++ b/src/bioformats.rs @@ -21,9 +21,21 @@ fn jvm() -> Rc { .unwrap() .to_path_buf(); - let class_path = path.parent().unwrap(); + let class_path = if path.join("jassets").exists() { + path.as_path() + } else { + path.parent().unwrap() + }; + if !class_path.join("jassets").exists() { + panic!( + "jassets directory does not exist in {}", + class_path.display() + ); + } + Rc::new( JvmBuilder::new() + .skip_setting_native_lib() .with_base_path(class_path.to_str().unwrap()) .build() .expect("Failed to build JVM"), @@ -33,11 +45,25 @@ fn jvm() -> Rc { }) } -#[cfg(feature = "python")] -pub(crate) fn download_bioformats(gpl_formats: bool) -> Result<()> { - let path = crate::py::ndbioimage_file().unwrap(); +pub fn download_bioformats(gpl_formats: bool) -> Result<()> { + #[cfg(feature = "python")] + let path = crate::py::ndbioimage_file()?; + + #[cfg(not(feature = "python"))] + let path = std::env::current_exe() + .unwrap() + .parent() + .unwrap() + .to_path_buf(); + let class_path = path.parent().unwrap(); + let jassets = class_path.join("jassets"); + if !jassets.exists() { + std::fs::create_dir_all(jassets)?; + } + println!("installing jassets in {}", class_path.display()); let jvm = JvmBuilder::new() + .skip_setting_native_lib() .with_base_path(class_path.to_str().unwrap()) .with_maven_settings(j4rs::MavenSettings::new(vec![ j4rs::MavenArtifactRepo::from( @@ -46,10 +72,10 @@ pub(crate) fn download_bioformats(gpl_formats: bool) -> Result<()> { ])) .build()?; - jvm.deploy_artifact(&j4rs::MavenArtifact::from("ome:bioformats_package:8.1.0"))?; + jvm.deploy_artifact(&j4rs::MavenArtifact::from("ome:bioformats_package:8.3.0"))?; if gpl_formats { - jvm.deploy_artifact(&j4rs::MavenArtifact::from("ome:formats-gpl:8.1.0"))?; + jvm.deploy_artifact(&j4rs::MavenArtifact::from("ome:formats-gpl:8.3.0"))?; } Ok(()) diff --git a/src/colors.rs b/src/colors.rs new file mode 100644 index 0000000..0cfbe2e --- /dev/null +++ b/src/colors.rs @@ -0,0 +1,207 @@ +use anyhow::{Error, Result}; +use lazy_static::lazy_static; +use std::collections::HashMap; +use std::fmt::Display; +use std::str::FromStr; + +lazy_static! { + static ref COLORS: HashMap = { + HashMap::from([ + ("b".to_string(), "#0000FF".to_string()), + ("g".to_string(), "#008000".to_string()), + ("r".to_string(), "#FF0000".to_string()), + ("c".to_string(), "#00BFBF".to_string()), + ("m".to_string(), "#BF00BF".to_string()), + ("y".to_string(), "#BFBF00".to_string()), + ("k".to_string(), "#000000".to_string()), + ("w".to_string(), "#FFFFFF".to_string()), + ("aliceblue".to_string(), "#F0F8FF".to_string()), + ("antiquewhite".to_string(), "#FAEBD7".to_string()), + ("aqua".to_string(), "#00FFFF".to_string()), + ("aquamarine".to_string(), "#7FFFD4".to_string()), + ("azure".to_string(), "#F0FFFF".to_string()), + ("beige".to_string(), "#F5F5DC".to_string()), + ("bisque".to_string(), "#FFE4C4".to_string()), + ("black".to_string(), "#000000".to_string()), + ("blanchedalmond".to_string(), "#FFEBCD".to_string()), + ("blue".to_string(), "#0000FF".to_string()), + ("blueviolet".to_string(), "#8A2BE2".to_string()), + ("brown".to_string(), "#A52A2A".to_string()), + ("burlywood".to_string(), "#DEB887".to_string()), + ("cadetblue".to_string(), "#5F9EA0".to_string()), + ("chartreuse".to_string(), "#7FFF00".to_string()), + ("chocolate".to_string(), "#D2691E".to_string()), + ("coral".to_string(), "#FF7F50".to_string()), + ("cornflowerblue".to_string(), "#6495ED".to_string()), + ("cornsilk".to_string(), "#FFF8DC".to_string()), + ("crimson".to_string(), "#DC143C".to_string()), + ("cyan".to_string(), "#00FFFF".to_string()), + ("darkblue".to_string(), "#00008B".to_string()), + ("darkcyan".to_string(), "#008B8B".to_string()), + ("darkgoldenrod".to_string(), "#B8860B".to_string()), + ("darkgray".to_string(), "#A9A9A9".to_string()), + ("darkgreen".to_string(), "#006400".to_string()), + ("darkgrey".to_string(), "#A9A9A9".to_string()), + ("darkkhaki".to_string(), "#BDB76B".to_string()), + ("darkmagenta".to_string(), "#8B008B".to_string()), + ("darkolivegreen".to_string(), "#556B2F".to_string()), + ("darkorange".to_string(), "#FF8C00".to_string()), + ("darkorchid".to_string(), "#9932CC".to_string()), + ("darkred".to_string(), "#8B0000".to_string()), + ("darksalmon".to_string(), "#E9967A".to_string()), + ("darkseagreen".to_string(), "#8FBC8F".to_string()), + ("darkslateblue".to_string(), "#483D8B".to_string()), + ("darkslategray".to_string(), "#2F4F4F".to_string()), + ("darkslategrey".to_string(), "#2F4F4F".to_string()), + ("darkturquoise".to_string(), "#00CED1".to_string()), + ("darkviolet".to_string(), "#9400D3".to_string()), + ("deeppink".to_string(), "#FF1493".to_string()), + ("deepskyblue".to_string(), "#00BFFF".to_string()), + ("dimgray".to_string(), "#696969".to_string()), + ("dimgrey".to_string(), "#696969".to_string()), + ("dodgerblue".to_string(), "#1E90FF".to_string()), + ("firebrick".to_string(), "#B22222".to_string()), + ("floralwhite".to_string(), "#FFFAF0".to_string()), + ("forestgreen".to_string(), "#228B22".to_string()), + ("fuchsia".to_string(), "#FF00FF".to_string()), + ("gainsboro".to_string(), "#DCDCDC".to_string()), + ("ghostwhite".to_string(), "#F8F8FF".to_string()), + ("gold".to_string(), "#FFD700".to_string()), + ("goldenrod".to_string(), "#DAA520".to_string()), + ("gray".to_string(), "#808080".to_string()), + ("green".to_string(), "#008000".to_string()), + ("greenyellow".to_string(), "#ADFF2F".to_string()), + ("grey".to_string(), "#808080".to_string()), + ("honeydew".to_string(), "#F0FFF0".to_string()), + ("hotpink".to_string(), "#FF69B4".to_string()), + ("indianred".to_string(), "#CD5C5C".to_string()), + ("indigo".to_string(), "#4B0082".to_string()), + ("ivory".to_string(), "#FFFFF0".to_string()), + ("khaki".to_string(), "#F0E68C".to_string()), + ("lavender".to_string(), "#E6E6FA".to_string()), + ("lavenderblush".to_string(), "#FFF0F5".to_string()), + ("lawngreen".to_string(), "#7CFC00".to_string()), + ("lemonchiffon".to_string(), "#FFFACD".to_string()), + ("lightblue".to_string(), "#ADD8E6".to_string()), + ("lightcoral".to_string(), "#F08080".to_string()), + ("lightcyan".to_string(), "#E0FFFF".to_string()), + ("lightgoldenrodyellow".to_string(), "#FAFAD2".to_string()), + ("lightgray".to_string(), "#D3D3D3".to_string()), + ("lightgreen".to_string(), "#90EE90".to_string()), + ("lightgrey".to_string(), "#D3D3D3".to_string()), + ("lightpink".to_string(), "#FFB6C1".to_string()), + ("lightsalmon".to_string(), "#FFA07A".to_string()), + ("lightseagreen".to_string(), "#20B2AA".to_string()), + ("lightskyblue".to_string(), "#87CEFA".to_string()), + ("lightslategray".to_string(), "#778899".to_string()), + ("lightslategrey".to_string(), "#778899".to_string()), + ("lightsteelblue".to_string(), "#B0C4DE".to_string()), + ("lightyellow".to_string(), "#FFFFE0".to_string()), + ("lime".to_string(), "#00FF00".to_string()), + ("limegreen".to_string(), "#32CD32".to_string()), + ("linen".to_string(), "#FAF0E6".to_string()), + ("magenta".to_string(), "#FF00FF".to_string()), + ("maroon".to_string(), "#800000".to_string()), + ("mediumaquamarine".to_string(), "#66CDAA".to_string()), + ("mediumblue".to_string(), "#0000CD".to_string()), + ("mediumorchid".to_string(), "#BA55D3".to_string()), + ("mediumpurple".to_string(), "#9370DB".to_string()), + ("mediumseagreen".to_string(), "#3CB371".to_string()), + ("mediumslateblue".to_string(), "#7B68EE".to_string()), + ("mediumspringgreen".to_string(), "#00FA9A".to_string()), + ("mediumturquoise".to_string(), "#48D1CC".to_string()), + ("mediumvioletred".to_string(), "#C71585".to_string()), + ("midnightblue".to_string(), "#191970".to_string()), + ("mintcream".to_string(), "#F5FFFA".to_string()), + ("mistyrose".to_string(), "#FFE4E1".to_string()), + ("moccasin".to_string(), "#FFE4B5".to_string()), + ("navajowhite".to_string(), "#FFDEAD".to_string()), + ("navy".to_string(), "#000080".to_string()), + ("oldlace".to_string(), "#FDF5E6".to_string()), + ("olive".to_string(), "#808000".to_string()), + ("olivedrab".to_string(), "#6B8E23".to_string()), + ("orange".to_string(), "#FFA500".to_string()), + ("orangered".to_string(), "#FF4500".to_string()), + ("orchid".to_string(), "#DA70D6".to_string()), + ("palegoldenrod".to_string(), "#EEE8AA".to_string()), + ("palegreen".to_string(), "#98FB98".to_string()), + ("paleturquoise".to_string(), "#AFEEEE".to_string()), + ("palevioletred".to_string(), "#DB7093".to_string()), + ("papayawhip".to_string(), "#FFEFD5".to_string()), + ("peachpuff".to_string(), "#FFDAB9".to_string()), + ("peru".to_string(), "#CD853F".to_string()), + ("pink".to_string(), "#FFC0CB".to_string()), + ("plum".to_string(), "#DDA0DD".to_string()), + ("powderblue".to_string(), "#B0E0E6".to_string()), + ("purple".to_string(), "#800080".to_string()), + ("rebeccapurple".to_string(), "#663399".to_string()), + ("red".to_string(), "#FF0000".to_string()), + ("rosybrown".to_string(), "#BC8F8F".to_string()), + ("royalblue".to_string(), "#4169E1".to_string()), + ("saddlebrown".to_string(), "#8B4513".to_string()), + ("salmon".to_string(), "#FA8072".to_string()), + ("sandybrown".to_string(), "#F4A460".to_string()), + ("seagreen".to_string(), "#2E8B57".to_string()), + ("seashell".to_string(), "#FFF5EE".to_string()), + ("sienna".to_string(), "#A0522D".to_string()), + ("silver".to_string(), "#C0C0C0".to_string()), + ("skyblue".to_string(), "#87CEEB".to_string()), + ("slateblue".to_string(), "#6A5ACD".to_string()), + ("slategray".to_string(), "#708090".to_string()), + ("slategrey".to_string(), "#708090".to_string()), + ("snow".to_string(), "#FFFAFA".to_string()), + ("springgreen".to_string(), "#00FF7F".to_string()), + ("steelblue".to_string(), "#4682B4".to_string()), + ("tan".to_string(), "#D2B48C".to_string()), + ("teal".to_string(), "#008080".to_string()), + ("thistle".to_string(), "#D8BFD8".to_string()), + ("tomato".to_string(), "#FF6347".to_string()), + ("turquoise".to_string(), "#40E0D0".to_string()), + ("violet".to_string(), "#EE82EE".to_string()), + ("wheat".to_string(), "#F5DEB3".to_string()), + ("white".to_string(), "#FFFFFF".to_string()), + ("whitesmoke".to_string(), "#F5F5F5".to_string()), + ("yellow".to_string(), "#FFFF00".to_string()), + ("yellowgreen".to_string(), "#9ACD32".to_string()), + ]) + }; +} + +#[derive(Clone, Debug)] +pub struct Color { + r: u8, + g: u8, + b: u8, +} + +impl FromStr for Color { + type Err = Error; + + fn from_str(s: &str) -> Result { + let s = if !s.starts_with("#") { + if let Some(s) = COLORS.get(s) { + s + } else { + return Err(Error::msg(format!("invalid color: {}", s))); + } + } else { + s + }; + let r = u8::from_str_radix(&s[1..3], 16)?; + let g = u8::from_str_radix(&s[3..5], 16)?; + let b = u8::from_str_radix(&s[5..], 16)?; + Ok(Self { r, g, b }) + } +} + +impl Display for Color { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "#{:02X}{:02X}{:02X}", self.r, self.g, self.b) + } +} + +impl Color { + pub fn to_rgb(&self) -> Vec { + vec![self.r, self.g, self.b] + } +} diff --git a/src/lib.rs b/src/lib.rs index cd693ff..2ce9f52 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,22 +1,31 @@ mod bioformats; pub mod axes; +pub mod metadata; #[cfg(feature = "python")] mod py; pub mod reader; pub mod stats; pub mod view; +pub mod colors; +#[cfg(feature = "movie")] +pub mod movie; +#[cfg(feature = "tiff")] +pub mod tiff; + +pub use bioformats::download_bioformats; + #[cfg(test)] mod tests { - use crate::stats::MinMax; - use ndarray::{Array, Array4, Array5, NewAxis}; - use rayon::prelude::*; - use crate::axes::Axis; use crate::reader::{Frame, Reader}; + use crate::stats::MinMax; + use crate::view::Item; use anyhow::Result; - use ndarray::{s, Array2}; + use ndarray::{Array, Array4, Array5, NewAxis}; + use ndarray::{Array2, s}; + use rayon::prelude::*; fn open(file: &str) -> Result { let path = std::env::current_dir()? @@ -137,7 +146,6 @@ mod tests { let a = a.slice(s![0, ..5, 0, NewAxis, 100..200, ..]); let v = view.slice(s![0, ..5, 0, NewAxis, 100..200, ..])?; assert_eq!(v.shape(), a.shape()); - println!("\nshape: {:?}", v.shape()); let a = a.slice(s![NewAxis, .., .., NewAxis, .., .., NewAxis]); let v = v.slice(s![NewAxis, .., .., NewAxis, .., .., NewAxis])?; assert_eq!(v.shape(), a.shape()); @@ -160,10 +168,7 @@ mod tests { assert_eq!(view.shape(), a.shape()); let b: Array5 = view.clone().try_into()?; assert_eq!(b.shape(), a.shape()); - - println!("{:?}", view.axes()); let view = view.permute_axes(&[Axis::X, Axis::Z, Axis::Y])?; - println!("{:?}", view.axes()); let a = a.permuted_axes([4, 1, 2, 0, 3]); assert_eq!(view.shape(), a.shape()); let b: Array5 = view.clone().try_into()?; @@ -259,4 +264,107 @@ mod tests { assert_eq!(d.shape(), [1024, 1024]); Ok(()) } + + #[test] + fn item() -> Result<()> { + let file = "1xp53-01-AP1.czi"; + let reader = open(file)?; + let view = reader.view(); + let a = view.slice(s![.., 0, 0, 0, 0])?; + let b = a.slice(s![0])?; + let item = b.item::()?; + assert_eq!(item, 2); + Ok(()) + } + + #[test] + fn slice_cztyx() -> Result<()> { + let file = "1xp53-01-AP1.czi"; + let reader = open(file)?; + let view = reader.view().max_proj(Axis::Z)?.into_dyn(); + println!("view.axes: {:?}", view.get_axes()); + println!("view.slice: {:?}", view.get_slice()); + let r = view.reset_axes()?; + println!("r.axes: {:?}", r.get_axes()); + println!("r.slice: {:?}", r.get_slice()); + let a = view.slice_cztyx(s![0, 0, 0, .., ..])?; + println!("a.axes: {:?}", a.get_axes()); + println!("a.slice: {:?}", a.get_slice()); + assert_eq!(a.axes(), [Axis::Y, Axis::X]); + Ok(()) + } + + #[test] + fn reset_axes() -> Result<()> { + let file = "1xp53-01-AP1.czi"; + let reader = open(file)?; + let view = reader.view().max_proj(Axis::Z)?; + let view = view.reset_axes()?; + assert_eq!(view.axes(), [Axis::C, Axis::New, Axis::T, Axis::Y, Axis::X]); + let a = view.as_array::()?; + assert_eq!(a.ndim(), 5); + Ok(()) + } + + #[test] + fn reset_axes2() -> Result<()> { + let file = "Experiment-2029.czi"; + let reader = open(file)?; + let view = reader.view().squeeze()?; + let a = view.reset_axes()?; + assert_eq!(a.axes(), [Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X]); + Ok(()) + } + + #[test] + fn reset_axes3() -> Result<()> { + let file = "Experiment-2029.czi"; + let reader = open(file)?; + let view4 = reader.view().squeeze()?; + let view = view4.max_proj(Axis::Z)?.into_dyn(); + let slice = view.slice_cztyx(s![0, .., .., .., ..])?.into_dyn(); + let a = slice.as_array::()?; + assert_eq!(slice.shape(), [1, 10, 1280, 1280]); + assert_eq!(a.shape(), [1, 10, 1280, 1280]); + let r = slice.reset_axes()?; + let b = r.as_array::()?; + assert_eq!(r.shape(), [1, 1, 10, 1280, 1280]); + assert_eq!(b.shape(), [1, 1, 10, 1280, 1280]); + let q = slice.max_proj(Axis::C)?.max_proj(Axis::T)?; + let c = q.as_array::()?; + assert_eq!(q.shape(), [1, 1280, 1280]); + assert_eq!(c.shape(), [1, 1280, 1280]); + let p = q.reset_axes()?; + let d = p.as_array::()?; + println!("axes: {:?}", p.get_axes()); + println!("operations: {:?}", p.get_operations()); + println!("slice: {:?}", p.get_slice()); + assert_eq!(p.shape(), [1, 1, 1, 1280, 1280]); + assert_eq!(d.shape(), [1, 1, 1, 1280, 1280]); + Ok(()) + } + + #[test] + fn max() -> Result<()> { + let file = "Experiment-2029.czi"; + let reader = open(file)?; + let view = reader.view(); + let m = view.max_proj(Axis::T)?; + let a = m.as_array::()?; + assert_eq!(m.shape(), [2, 1, 1280, 1280]); + assert_eq!(a.shape(), [2, 1, 1280, 1280]); + let mc = view.max_proj(Axis::C)?; + let a = mc.as_array::()?; + assert_eq!(mc.shape(), [1, 10, 1280, 1280]); + assert_eq!(a.shape(), [1, 10, 1280, 1280]); + let mz = mc.max_proj(Axis::Z)?; + let a = mz.as_array::()?; + assert_eq!(mz.shape(), [10, 1280, 1280]); + assert_eq!(a.shape(), [10, 1280, 1280]); + let mt = mz.max_proj(Axis::T)?; + let a = mt.as_array::()?; + assert_eq!(mt.shape(), [1280, 1280]); + assert_eq!(a.shape(), [1280, 1280]); + Ok(()) + } } diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..50d0bdb --- /dev/null +++ b/src/main.rs @@ -0,0 +1,110 @@ +use anyhow::Result; +use clap::{Parser, Subcommand}; +#[cfg(feature = "movie")] +use ndbioimage::movie::MovieOptions; +use ndbioimage::reader::split_path_and_series; +#[cfg(feature = "tiff")] +use ndbioimage::tiff::TiffOptions; +use ndbioimage::view::View; +use std::path::PathBuf; + +#[derive(Parser)] +#[command(arg_required_else_help = true, version, about, long_about = None, propagate_version = true)] +struct Cli { + #[command(subcommand)] + command: Commands, +} + +#[derive(Subcommand)] +enum Commands { + /// Print some metadata + Info { + #[arg(value_name = "FILE", num_args(1..))] + file: Vec, + }, + /// Save the image as tiff file + #[cfg(feature = "tiff")] + Tiff { + #[arg(value_name = "FILE", num_args(1..))] + file: Vec, + #[arg(short, long, value_name = "COLOR", num_args(1..))] + colors: Vec, + #[arg(short, long, value_name = "OVERWRITE")] + overwrite: bool, + }, + /// Save the image as mp4 file + #[cfg(feature = "movie")] + Movie { + #[arg(value_name = "FILE", num_args(1..))] + file: Vec, + #[arg(short, long, value_name = "Velocity", default_value = "3.6")] + velocity: f64, + #[arg(short, long, value_name = "BRIGHTNESS")] + brightness: Vec, + #[arg(short, long, value_name = "SCALE", default_value = "1.0")] + scale: f64, + #[arg(short, long, value_name = "COLOR", num_args(1..))] + colors: Vec, + #[arg(short, long, value_name = "OVERWRITE")] + overwrite: bool, + }, + /// Download the BioFormats jar into the correct folder + DownloadBioFormats { + #[arg(short, long, value_name = "GPL_FORMATS")] + gpl_formats: bool, + }, +} + +pub(crate) fn main() -> Result<()> { + let cli = Cli::parse(); + match &cli.command { + Commands::Info { file } => { + for f in file { + let (path, series) = split_path_and_series(f)?; + let view = View::from_path(path, series.unwrap_or(0))?.squeeze()?; + println!("{}", view.summary()?); + } + } + #[cfg(feature = "tiff")] + Commands::Tiff { + file, + colors, + overwrite, + } => { + let mut options = TiffOptions::new(true, None, colors.clone(), *overwrite)?; + options.enable_bar()?; + for f in file { + let (path, series) = split_path_and_series(f)?; + let view = View::from_path(path, series.unwrap_or(0))?; + view.save_as_tiff(f.with_extension("tiff"), &options)?; + } + } + #[cfg(feature = "movie")] + Commands::Movie { + file, + velocity: speed, + brightness, + scale, + colors, + overwrite, + } => { + let options = MovieOptions::new( + *speed, + brightness.to_vec(), + *scale, + colors.to_vec(), + *overwrite, + )?; + for f in file { + let (path, series) = split_path_and_series(f)?; + let view = View::from_path(path, series.unwrap_or(0))?; + view.save_as_movie(f.with_extension("mp4"), &options)?; + } + } + Commands::DownloadBioFormats { gpl_formats } => { + ndbioimage::download_bioformats(*gpl_formats)? + } + } + + Ok(()) +} diff --git a/src/metadata.rs b/src/metadata.rs new file mode 100644 index 0000000..494e13d --- /dev/null +++ b/src/metadata.rs @@ -0,0 +1,372 @@ +use anyhow::{Result, anyhow}; +use itertools::Itertools; +use ome_metadata::Ome; +use ome_metadata::ome::{ + BinningType, Convert, Image, Instrument, Objective, Pixels, UnitsLength, UnitsTime, +}; + +impl Metadata for Ome { + fn get_instrument(&self) -> Option<&Instrument> { + let instrument_id = self.get_image()?.instrument_ref.as_ref()?.id.clone(); + self.instrument + .as_ref()? + .iter() + .find(|i| i.id == instrument_id) + } + + fn get_image(&self) -> Option<&Image> { + if let Some(image) = &self.image { + if !image.is_empty() { + Some(&image[0]) + } else { + None + } + } else { + None + } + } +} + +pub trait Metadata { + fn get_instrument(&self) -> Option<&Instrument>; + fn get_image(&self) -> Option<&Image>; + + fn get_pixels(&self) -> Option<&Pixels> { + if let Some(image) = self.get_image() { + Some(&image.pixels) + } else { + None + } + } + + fn get_objective(&self) -> Option<&Objective> { + let objective_id = self.get_image()?.objective_settings.as_ref()?.id.clone(); + self.get_instrument()? + .objective + .iter() + .find(|o| o.id == objective_id) + } + + fn get_tube_lens(&self) -> Option { + Some(Objective { + manufacturer: None, + model: Some("Unknown".to_string()), + serial_number: None, + lot_number: None, + id: "TubeLens:1".to_string(), + correction: None, + immersion: None, + lens_na: None, + nominal_magnification: Some(1.0), + calibrated_magnification: None, + working_distance: None, + working_distance_unit: UnitsLength::um, + iris: None, + annotation_ref: vec![], + }) + } + + /// shape of the data along cztyx axes + fn shape(&self) -> Result<(usize, usize, usize, usize, usize)> { + if let Some(pixels) = self.get_pixels() { + Ok(( + pixels.size_c as usize, + pixels.size_z as usize, + pixels.size_t as usize, + pixels.size_y as usize, + pixels.size_x as usize, + )) + } else { + Err(anyhow!("No image or pixels found")) + } + } + + /// pixel size in nm + fn pixel_size(&self) -> Result> { + if let Some(pixels) = self.get_pixels() { + match (pixels.physical_size_x, pixels.physical_size_y) { + (Some(x), Some(y)) => Ok(Some( + (pixels + .physical_size_x_unit + .convert(&UnitsLength::nm, x as f64)? + + pixels + .physical_size_y_unit + .convert(&UnitsLength::nm, y as f64)?) + / 2f64, + )), + (Some(x), None) => Ok(Some( + pixels + .physical_size_x_unit + .convert(&UnitsLength::nm, x as f64)? + .powi(2), + )), + (None, Some(y)) => Ok(Some( + pixels + .physical_size_y_unit + .convert(&UnitsLength::nm, y as f64)? + .powi(2), + )), + _ => Ok(None), + } + } else { + Ok(None) + } + } + + /// distance between planes in z-stack in nm + fn delta_z(&self) -> Result> { + if let Some(pixels) = self.get_pixels() { + if let Some(z) = pixels.physical_size_z { + return Ok(Some( + pixels + .physical_size_z_unit + .convert(&UnitsLength::nm, z as f64)?, + )); + } + } + Ok(None) + } + + /// time interval in seconds for time-lapse images + fn time_interval(&self) -> Result> { + if let Some(pixels) = self.get_pixels() { + if let Some(plane) = &pixels.plane { + if let Some(t) = plane.iter().map(|p| p.the_t).max() { + if t > 0 { + let plane_a = plane + .iter() + .find(|p| (p.the_c == 0) && (p.the_z == 0) && (p.the_t == 0)); + let plane_b = plane + .iter() + .find(|p| (p.the_c == 0) && (p.the_z == 0) && (p.the_t == t)); + if let (Some(a), Some(b)) = (plane_a, plane_b) { + if let (Some(a_t), Some(b_t)) = (a.delta_t, b.delta_t) { + return Ok(Some( + (b.delta_t_unit.convert(&UnitsTime::s, b_t as f64)? + - a.delta_t_unit.convert(&UnitsTime::s, a_t as f64)?) + .abs() + / (t as f64), + )); + } + } + } + } + } + } + Ok(None) + } + + /// exposure time for channel, z=0 and t=0 + fn exposure_time(&self, channel: usize) -> Result> { + let c = channel as i32; + if let Some(pixels) = self.get_pixels() { + if let Some(plane) = &pixels.plane { + if let Some(p) = plane + .iter() + .find(|p| (p.the_c == c) && (p.the_z == 0) && (p.the_t == 0)) + { + if let Some(t) = p.exposure_time { + return Ok(Some(p.exposure_time_unit.convert(&UnitsTime::s, t as f64)?)); + } + } + } + } + Ok(None) + } + + fn binning(&self, channel: usize) -> Option { + match self + .get_pixels()? + .channel + .get(channel)? + .detector_settings + .as_ref()? + .binning + .as_ref()? + { + BinningType::_1X1 => Some(1), + BinningType::_2X2 => Some(2), + BinningType::_4X4 => Some(4), + BinningType::_8X8 => Some(8), + BinningType::Other => None, + } + } + + fn laser_wavelengths(&self, channel: usize) -> Result> { + if let Some(pixels) = self.get_pixels() { + if let Some(channel) = pixels.channel.get(channel) { + if let Some(w) = channel.excitation_wavelength { + return Ok(Some( + channel + .excitation_wavelength_unit + .convert(&UnitsLength::nm, w as f64)?, + )); + } + } + } + Ok(None) + } + + fn laser_powers(&self, channel: usize) -> Result> { + if let Some(pixels) = self.get_pixels() { + if let Some(channel) = pixels.channel.get(channel) { + if let Some(ls) = &channel.light_source_settings { + if let Some(a) = ls.attenuation { + return if (0. ..=1.).contains(&a) { + Ok(Some(1f64 - (a as f64))) + } else { + Err(anyhow!("Invalid attenuation value")) + }; + } + } + } + } + Ok(None) + } + + fn objective_name(&self) -> Option { + Some(self.get_objective()?.model.as_ref()?.clone()) + } + + fn magnification(&self) -> Option { + Some( + (self.get_objective()?.nominal_magnification? as f64) + * (self.get_tube_lens()?.nominal_magnification? as f64), + ) + } + + fn tube_lens_name(&self) -> Option { + self.get_tube_lens()?.model.clone() + } + + fn filter_set_name(&self, channel: usize) -> Option { + let filter_set_id = self + .get_pixels()? + .channel + .get(channel)? + .filter_set_ref + .as_ref()? + .id + .clone(); + self.get_instrument() + .as_ref()? + .filter_set + .iter() + .find(|f| f.id == filter_set_id)? + .model + .clone() + } + + fn gain(&self, channel: usize) -> Option { + if let Some(pixels) = self.get_pixels() { + Some( + *pixels + .channel + .get(channel)? + .detector_settings + .as_ref()? + .gain + .as_ref()? as f64, + ) + } else { + None + } + } + + fn summary(&self) -> Result { + let size_c = if let Some(pixels) = self.get_pixels() { + pixels.channel.len() + } else { + 0 + }; + let mut s = "".to_string(); + if let Ok(Some(pixel_size)) = self.pixel_size() { + s.push_str(&format!("pixel size: {pixel_size:.2} nm\n")); + } + if let Ok(Some(delta_z)) = self.delta_z() { + s.push_str(&format!("z-interval: {delta_z:.2} nm\n")) + } + if let Ok(Some(time_interval)) = self.time_interval() { + s.push_str(&format!("time interval: {time_interval:.2} s\n")) + } + let exposure_time = (0..size_c) + .map(|c| self.exposure_time(c)) + .collect::>>()? + .into_iter() + .flatten() + .collect::>(); + if !exposure_time.is_empty() { + s.push_str(&format!( + "exposure time: {}\n", + exposure_time.into_iter().join(" | ") + )); + } + if let Some(magnification) = self.magnification() { + s.push_str(&format!("magnification: {magnification}x\n")) + } + if let Some(objective_name) = self.objective_name() { + s.push_str(&format!("objective: {objective_name}\n")) + } + if let Some(tube_lens_name) = self.tube_lens_name() { + s.push_str(&format!("tube lens: {tube_lens_name}\n")) + } + let filter_set_name = (0..size_c) + .map(|c| self.filter_set_name(c)) + .collect::>() + .into_iter() + .flatten() + .collect::>(); + if !filter_set_name.is_empty() { + s.push_str(&format!( + "filter set: {}\n", + filter_set_name.into_iter().join(" | ") + )); + } + let gain = (0..size_c) + .map(|c| self.gain(c)) + .collect::>() + .into_iter() + .flatten() + .collect::>(); + if !gain.is_empty() { + s.push_str(&format!("gain: {}\n", gain.into_iter().join(" "))); + } + let laser_wavelengths = (0..size_c) + .map(|c| self.laser_wavelengths(c)) + .collect::>>()? + .into_iter() + .flatten() + .collect::>(); + if !laser_wavelengths.is_empty() { + s.push_str(&format!( + "laser colors: {} nm\n", + laser_wavelengths.into_iter().join(" | ") + )); + } + let laser_powers = (0..size_c) + .map(|c| self.laser_powers(c)) + .collect::>>()? + .into_iter() + .flatten() + .collect::>(); + if !laser_powers.is_empty() { + s.push_str(&format!( + "laser powers: {}\n", + laser_powers.into_iter().join(" | ") + )); + } + let binning = (0..size_c) + .map(|c| self.binning(c)) + .collect::>() + .into_iter() + .flatten() + .collect::>(); + if !binning.is_empty() { + s.push_str(&format!( + "binning: {}\n", + binning.into_iter().join(" | ") + )); + } + Ok(s.to_string()) + } +} diff --git a/src/movie.rs b/src/movie.rs new file mode 100644 index 0000000..b85fbc0 --- /dev/null +++ b/src/movie.rs @@ -0,0 +1,258 @@ +use crate::axes::Axis; +use crate::colors::Color; +use crate::view::View; +use anyhow::{Result, anyhow}; +use ffmpeg_sidecar::command::FfmpegCommand; +use ffmpeg_sidecar::download::auto_download; +use ffmpeg_sidecar::event::{FfmpegEvent, LogLevel}; +use itertools::Itertools; +use ndarray::{Array2, Array3, Dimension, IxDyn, s, stack}; +use ordered_float::OrderedFloat; +use std::io::Write; +use std::path::Path; +use std::thread; + +pub struct MovieOptions { + velocity: f64, + brightness: Vec, + scale: f64, + colors: Option>>, + overwrite: bool, +} + +impl Default for MovieOptions { + fn default() -> Self { + Self { + velocity: 3.6, + brightness: Vec::new(), + scale: 1.0, + colors: None, + overwrite: false, + } + } +} + +impl MovieOptions { + pub fn new( + velocity: f64, + brightness: Vec, + scale: f64, + colors: Vec, + overwrite: bool, + ) -> Result { + let colors = if colors.is_empty() { + None + } else { + let colors = colors + .iter() + .map(|c| c.parse::()) + .collect::>>()?; + Some(colors.into_iter().map(|c| c.to_rgb()).collect()) + }; + Ok(Self { + velocity, + brightness, + scale, + colors, + overwrite, + }) + } + + pub fn set_velocity(&mut self, velocity: f64) { + self.velocity = velocity; + } + + pub fn set_brightness(&mut self, brightness: Vec) { + self.brightness = brightness; + } + + pub fn set_scale(&mut self, scale: f64) { + self.scale = scale; + } + + pub fn set_colors(&mut self, colors: &[String]) -> Result<()> { + let colors = colors + .iter() + .map(|c| c.parse::()) + .collect::>>()?; + self.colors = Some(colors.into_iter().map(|c| c.to_rgb()).collect()); + Ok(()) + } + + pub fn set_overwrite(&mut self, overwrite: bool) { + self.overwrite = overwrite; + } +} + +fn get_ab(tyx: View) -> Result<(f64, f64)> { + let s = tyx + .as_array::()? + .iter() + .filter_map(|&i| { + if i == 0.0 || !i.is_finite() { + None + } else { + Some(OrderedFloat::from(i)) + } + }) + .sorted_unstable() + .map(f64::from) + .collect::>(); + + let n = s.len(); + let mut a = s[n / 100]; + let mut b = s[n - n / 100 - 1]; + if a == b { + a = s[0]; + b = s[n - 1]; + } + if a == b { + a = 1.0; + b = 1.0; + } + Ok((a, b)) +} + +fn cframe(frame: Array2, color: &[u8], a: f64, b: f64) -> Array3 { + let frame = (frame - a) / (b - a); + let color = color + .iter() + .map(|&c| (c as f64) / 255.0) + .collect::>(); + let frame = color + .iter() + .map(|&c| (c * &frame).to_owned()) + .collect::>>(); + let view = frame.iter().map(|c| c.view()).collect::>(); + stack(ndarray::Axis(0), &view).unwrap() +} + +impl View +where + D: Dimension, +{ + pub fn save_as_movie

(&self, path: P, options: &MovieOptions) -> Result<()> + where + P: AsRef, + { + let path = path.as_ref().to_path_buf(); + if path.exists() { + if options.overwrite { + std::fs::remove_file(&path)?; + } else { + return Err(anyhow!("File {} already exists", path.display())); + } + } + let view = self.max_proj(Axis::Z)?.reset_axes()?; + let velocity = options.velocity; + let mut brightness = options.brightness.clone(); + let scale = options.scale; + let shape = view.shape(); + let size_c = shape[0]; + let size_t = shape[2]; + let size_x = shape[3]; + let size_y = shape[4]; + let shape_x = 2 * (((size_x as f64 * scale) / 2.).round() as usize); + let shape_y = 2 * (((size_y as f64 * scale) / 2.).round() as usize); + + while brightness.len() < size_c { + brightness.push(1.0); + } + let mut colors = if let Some(colors) = options.colors.as_ref() { + colors.to_vec() + } else { + Vec::new() + }; + while colors.len() < size_c { + colors.push(vec![255, 255, 255]); + } + + auto_download()?; + let mut movie = FfmpegCommand::new() + .args([ + "-f", + "rawvideo", + "-pix_fmt", + "gray", + "-s", + &format!("{}x{}", size_x, size_y), + ]) + .input("-") + .args([ + "-vcodec", + "libx264", + "-preset", + "veryslow", + "-pix_fmt", + "yuv420p", + "-r", + "7", + "-vf", + &format!("setpts={velocity}*PTS,scale={shape_x}:{shape_y}:flags=neighbor"), + ]) + .output(path.to_str().expect("path cannot be converted to string")) + .spawn()?; + let mut stdin = movie.take_stdin().unwrap(); + + let ab = (0..size_c) + .map(|c| match view.slice(s![c, .., .., .., ..]) { + Ok(slice) => get_ab(slice.into_dyn()), + Err(e) => Err(e), + }) + .collect::>>()?; + + thread::spawn(move || { + for t in 0..size_t { + let mut frame = Array3::::zeros((3, size_y, size_y)); + for c in 0..size_c { + frame = frame + + cframe( + view.get_frame(c, 0, t).unwrap(), + &colors[c], + ab[c].0, + ab[c].1 / brightness[c], + ); + } + let frame = frame.mapv(|i| { + if i < 0.0 { + 0 + } else if i > 1.0 { + 1 + } else { + (255.0 * i).round() as u8 + } + }); + let bytes: Vec<_> = frame.flatten().into_iter().collect(); + stdin.write_all(&bytes).unwrap(); + } + }); + + movie.iter()?.for_each(|e| match e { + FfmpegEvent::Log(LogLevel::Error, e) => println!("Error: {}", e), + FfmpegEvent::Progress(p) => println!("Progress: {} / 00:00:15", p.time), + _ => {} + }); + Ok(()) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::reader::Reader; + + #[test] + fn movie() -> Result<()> { + let file = "1xp53-01-AP1.czi"; + let path = std::env::current_dir()? + .join("tests") + .join("files") + .join(file); + let reader = Reader::new(&path, 0)?; + let view = reader.view(); + let mut options = MovieOptions::default(); + options.set_overwrite(true); + view.save_as_movie("/home/wim/tmp/movie.mp4", &options)?; + Ok(()) + } +} diff --git a/src/py.rs b/src/py.rs index cb64a6b..31c8a75 100644 --- a/src/py.rs +++ b/src/py.rs @@ -1,16 +1,21 @@ use crate::axes::Axis; use crate::bioformats::download_bioformats; +use crate::metadata::Metadata; use crate::reader::{PixelType, Reader}; use crate::view::{Item, View}; -use anyhow::{anyhow, Result}; +use anyhow::{Result, anyhow}; +use itertools::Itertools; use ndarray::{Ix0, IxDyn, SliceInfoElem}; use numpy::IntoPyArray; +use ome_metadata::Ome; +use pyo3::IntoPyObjectExt; +use pyo3::exceptions::PyNotImplementedError; use pyo3::prelude::*; use pyo3::types::{PyEllipsis, PyInt, PyList, PySlice, PySliceMethods, PyString, PyTuple}; -use pyo3::IntoPyObjectExt; use serde::{Deserialize, Serialize}; use serde_json::{from_str, to_string}; use std::path::PathBuf; +use std::sync::Arc; #[pyclass(module = "ndbioimage.ndbioimage_rs")] struct ViewConstructor; @@ -22,12 +27,10 @@ impl ViewConstructor { Self } - fn __getstate__(&self) -> (u8,) { - (0,) + fn __getnewargs__<'py>(&self, py: Python<'py>) -> Bound<'py, PyTuple> { + PyTuple::empty(py) } - fn __setstate__(&self, _state: (u8,)) {} - #[staticmethod] fn __call__(state: String) -> PyResult { if let Ok(new) = from_str(&state) { @@ -40,32 +43,122 @@ impl ViewConstructor { #[pyclass(subclass, module = "ndbioimage.ndbioimage_rs")] #[pyo3(name = "View")] -#[derive(Debug, Serialize, Deserialize)] +#[derive(Clone, Debug, Serialize, Deserialize)] struct PyView { view: View, dtype: PixelType, + ome: Arc, } #[pymethods] impl PyView { - #[new] - #[pyo3(signature = (path, series = 0, dtype = "uint16"))] /// new view on a file at path, open series #, open as dtype: (u)int(8/16/32) or float(32/64) - fn new(path: &str, series: usize, dtype: &str) -> PyResult { - let mut path = PathBuf::from(path); - if path.is_dir() { - for file in path.read_dir()?.flatten() { - let p = file.path(); - if file.path().is_file() & (p.extension() == Some("tif".as_ref())) { - path = p; - break; + #[new] + #[pyo3(signature = (path, series = 0, dtype = "uint16", axes = "cztyx"))] + fn new(path: Bound<'_, PyAny>, series: usize, dtype: &str, axes: &str) -> PyResult { + if path.is_instance_of::() { + Ok(path.downcast_into::()?.extract::()?) + } else { + let mut path = PathBuf::from(path.downcast_into::()?.extract::()?); + if path.is_dir() { + for file in path.read_dir()?.flatten() { + let p = file.path(); + if file.path().is_file() && (p.extension() == Some("tif".as_ref())) { + path = p; + break; + } } } + let axes = axes + .chars() + .map(|a| a.to_string().parse()) + .collect::>>()?; + let reader = Reader::new(&path, series)?; + let view = View::new_with_axes(Arc::new(reader), axes)?; + let dtype = dtype.parse()?; + let ome = Arc::new(view.get_ome()?); + Ok(Self { view, dtype, ome }) + } + } + + fn squeeze<'py>(&self, py: Python<'py>) -> PyResult> { + let view = self.view.squeeze()?; + if view.ndim() == 0 { + Ok(match self.dtype { + PixelType::I8 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::U8 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::I16 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::U16 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::I32 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::U32 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::I64 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::U64 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::I128 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::U128 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::F32 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::F64 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::F128 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + }) + } else { + PyView { + view, + dtype: self.dtype.clone(), + ome: self.ome.clone(), + } + .into_bound_py_any(py) } - Ok(Self { - view: Reader::new(&path, series as i32)?.view().into_dyn(), - dtype: dtype.parse()?, - }) } /// close the file: does nothing as this is handled automatically @@ -73,11 +166,18 @@ impl PyView { Ok(()) } - fn copy(&self) -> PyView { - PyView { + /// change the data type of the view: (u)int(8/16/32) or float(32/64) + fn as_type(&self, dtype: &str) -> PyResult { + Ok(PyView { view: self.view.clone(), - dtype: self.dtype.clone(), - } + dtype: dtype.parse()?, + ome: self.ome.clone(), + }) + } + + /// change the data type of the view: (u)int(8/16/32) or float(32/64) + fn astype(&self, dtype: &str) -> PyResult { + self.as_type(dtype) } /// slice the view and return a new view or a single number @@ -216,11 +316,33 @@ impl PyView { PyView { view, dtype: self.dtype.clone(), + ome: self.ome.clone(), } .into_bound_py_any(py) } } + #[pyo3(signature = (dtype = None))] + fn __array__<'py>(&self, py: Python<'py>, dtype: Option<&str>) -> PyResult> { + if let Some(dtype) = dtype { + self.as_type(dtype)?.as_array(py) + } else { + self.as_array(py) + } + } + + fn __contains__(&self, _item: Bound<'_, PyAny>) -> PyResult { + Err(PyNotImplementedError::new_err("contains not implemented")) + } + + fn __enter__<'p>(slf: PyRef<'p, Self>, _py: Python<'p>) -> PyResult> { + Ok(slf) + } + + fn __exit__(&self) -> PyResult<()> { + self.close() + } + fn __reduce__(&self) -> PyResult<(ViewConstructor, (String,))> { if let Ok(s) = to_string(self) { Ok((ViewConstructor, (s,))) @@ -229,6 +351,34 @@ impl PyView { } } + fn __copy__(&self) -> Self { + Self { + view: self.view.clone(), + dtype: self.dtype.clone(), + ome: self.ome.clone(), + } + } + + fn copy(&self) -> Self { + Self { + view: self.view.clone(), + dtype: self.dtype.clone(), + ome: self.ome.clone(), + } + } + + fn __len__(&self) -> PyResult { + Ok(self.view.len()) + } + + fn __repr__(&self) -> PyResult { + Ok(self.view.summary()?) + } + + fn __str__(&self) -> PyResult { + Ok(self.view.path.display().to_string()) + } + /// retrieve a single frame at czt, sliced accordingly fn get_frame<'py>( &self, @@ -240,73 +390,113 @@ impl PyView { Ok(match self.dtype { PixelType::I8 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), PixelType::U8 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), PixelType::I16 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), PixelType::U16 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), PixelType::I32 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), PixelType::U32 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), PixelType::F32 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), PixelType::F64 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), PixelType::I64 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), PixelType::U64 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), PixelType::I128 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), PixelType::U128 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), PixelType::F128 => self .view - .get_frame::(c, z, t)? + .get_frame::(c, z, t)? .into_pyarray(py) .into_any(), }) } - /// retrieve the ome metadata as an xml string + fn flatten<'py>(&self, py: Python<'py>) -> PyResult> { + Ok(match self.dtype { + PixelType::I8 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::U8 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::I16 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::U16 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::I32 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::U32 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::F32 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::F64 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::I64 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::U64 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::I128 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::U128 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::F128 => self.view.flatten::()?.into_pyarray(py).into_any(), + }) + } + + fn to_bytes(&self) -> PyResult> { + Ok(match self.dtype { + PixelType::I8 => self.view.to_bytes::()?, + PixelType::U8 => self.view.to_bytes::()?, + PixelType::I16 => self.view.to_bytes::()?, + PixelType::U16 => self.view.to_bytes::()?, + PixelType::I32 => self.view.to_bytes::()?, + PixelType::U32 => self.view.to_bytes::()?, + PixelType::F32 => self.view.to_bytes::()?, + PixelType::F64 => self.view.to_bytes::()?, + PixelType::I64 => self.view.to_bytes::()?, + PixelType::U64 => self.view.to_bytes::()?, + PixelType::I128 => self.view.to_bytes::()?, + PixelType::U128 => self.view.to_bytes::()?, + PixelType::F128 => self.view.to_bytes::()?, + }) + } + + fn tobytes(&self) -> PyResult> { + self.to_bytes() + } + + /// retrieve the ome metadata as an XML string fn get_ome_xml(&self) -> PyResult { Ok(self.view.get_ome_xml()?) } @@ -319,18 +509,14 @@ impl PyView { /// the series in the file #[getter] - fn series(&self) -> PyResult { + fn series(&self) -> PyResult { Ok(self.view.series) } /// the axes in the view #[getter] - fn axes(&self) -> Vec { - self.view - .axes() - .iter() - .map(|a| format!("{:?}", a)) - .collect() + fn axes(&self) -> String { + self.view.axes().iter().map(|a| format!("{:?}", a)).join("") } /// the shape of the view @@ -391,6 +577,7 @@ impl PyView { Ok(PyView { view, dtype: self.dtype.clone(), + ome: self.ome.clone(), }) } @@ -409,9 +596,16 @@ impl PyView { Ok(PyView { view, dtype: self.dtype.clone(), + ome: self.ome.clone(), }) } + #[allow(non_snake_case)] + #[getter] + fn T(&self) -> PyResult { + self.transpose(None) + } + /// collect data into a numpy array fn as_array<'py>(&self, py: Python<'py>) -> PyResult> { Ok(match self.dtype { @@ -431,14 +625,6 @@ impl PyView { }) } - /// change the data type of the view: (u)int(8/16/32) or float(32/64) - fn as_type(&self, dtype: String) -> PyResult { - Ok(PyView { - view: self.view.clone(), - dtype: dtype.parse()?, - }) - } - #[getter] fn get_dtype(&self) -> PyResult<&str> { Ok(match self.dtype { @@ -475,6 +661,7 @@ impl PyView { PyView { dtype: self.dtype.clone(), view: self.view.max_proj(self.get_ax(axis)?)?, + ome: self.ome.clone(), } .into_bound_py_any(py) } else { @@ -507,6 +694,7 @@ impl PyView { PyView { dtype: self.dtype.clone(), view: self.view.min_proj(self.get_ax(axis)?)?, + ome: self.ome.clone(), } .into_bound_py_any(py) } else { @@ -544,6 +732,7 @@ impl PyView { PyView { dtype, view: self.view.mean_proj(self.get_ax(axis)?)?, + ome: self.ome.clone(), } .into_bound_py_any(py) } else { @@ -580,6 +769,7 @@ impl PyView { PyView { dtype, view: self.view.sum_proj(self.get_ax(axis)?)?, + ome: self.ome.clone(), } .into_bound_py_any(py) } else { @@ -595,6 +785,83 @@ impl PyView { }) } } + + #[getter] + fn z_stack(&self) -> PyResult { + if let Some(s) = self.view.size_ax(Axis::Z) { + Ok(s > 1) + } else { + Ok(false) + } + } + + #[getter] + fn time_series(&self) -> PyResult { + if let Some(s) = self.view.size_ax(Axis::T) { + Ok(s > 1) + } else { + Ok(false) + } + } + + #[getter] + fn pixel_size(&self) -> PyResult> { + Ok(self.ome.pixel_size()?) + } + + #[getter] + fn delta_z(&self) -> PyResult> { + Ok(self.ome.delta_z()?) + } + + #[getter] + fn time_interval(&self) -> PyResult> { + Ok(self.ome.time_interval()?) + } + + fn exposure_time(&self, channel: usize) -> PyResult> { + Ok(self.ome.exposure_time(channel)?) + } + + fn binning(&self, channel: usize) -> Option { + self.ome.binning(channel) + } + + fn laser_wavelengths(&self, channel: usize) -> PyResult> { + Ok(self.ome.laser_wavelengths(channel)?) + } + + fn laser_power(&self, channel: usize) -> PyResult> { + Ok(self.ome.laser_powers(channel)?) + } + + #[getter] + fn objective_name(&self) -> Option { + self.ome.objective_name() + } + + #[getter] + fn magnification(&self) -> Option { + self.ome.magnification() + } + + #[getter] + fn tube_lens_name(&self) -> Option { + self.ome.tube_lens_name() + } + + fn filter_set_name(&self, channel: usize) -> Option { + self.ome.filter_set_name(channel) + } + + fn gain(&self, channel: usize) -> Option { + self.ome.gain(channel) + } + + /// gives a helpful summary of the recorded experiment + fn summary(&self) -> PyResult { + Ok(self.view.summary()?) + } } pub(crate) fn ndbioimage_file() -> anyhow::Result { diff --git a/src/reader.rs b/src/reader.rs index c7a6e2b..89d40e4 100644 --- a/src/reader.rs +++ b/src/reader.rs @@ -2,9 +2,10 @@ use crate::axes::Axis; use crate::bioformats; use crate::bioformats::{DebugTools, ImageReader, MetadataTools}; use crate::view::View; -use anyhow::{anyhow, Error, Result}; -use ndarray::{s, Array2, Ix5}; +use anyhow::{Error, Result, anyhow}; +use ndarray::{Array2, Ix5, s}; use num::{FromPrimitive, Zero}; +use ome_metadata::Ome; use serde::{Deserialize, Serialize}; use std::any::type_name; use std::fmt::Debug; @@ -14,6 +15,26 @@ use std::str::FromStr; use std::sync::Arc; use thread_local::ThreadLocal; +pub fn split_path_and_series

(path: P) -> Result<(PathBuf, Option)> +where + P: Into, +{ + let path = path.into(); + let file_name = path + .file_name() + .ok_or_else(|| anyhow!("No file name"))? + .to_str() + .ok_or_else(|| anyhow!("No file name"))?; + if file_name.to_lowercase().starts_with("pos") { + if let Some(series) = file_name.get(3..) { + if let Ok(series) = series.parse::() { + return Ok((path, Some(series))); + } + } + } + Ok((path, None)) +} + /// Pixel types (u)int(8/16/32) or float(32/64), (u/i)(64/128) are not included in bioformats #[allow(clippy::upper_case_acronyms)] #[derive(Clone, Debug, Serialize, Deserialize)] @@ -228,7 +249,7 @@ pub struct Reader { /// path to file pub path: PathBuf, /// which (if more than 1) of the series in the file to open - pub series: i32, + pub series: usize, /// size x (horizontal) pub size_x: usize, /// size y (vertical) @@ -254,7 +275,7 @@ impl Deref for Reader { let ome_meta = meta_data_tools.create_ome_xml_metadata().unwrap(); reader.set_metadata_store(ome_meta).unwrap(); reader.set_id(self.path.to_str().unwrap()).unwrap(); - reader.set_series(self.series).unwrap(); + reader.set_series(self.series as i32).unwrap(); reader }) } @@ -284,11 +305,14 @@ impl Debug for Reader { impl Reader { /// Create a new reader for the image file at a path, and open series #. - pub fn new(path: &Path, series: i32) -> Result { + pub fn new

(path: P, series: usize) -> Result + where + P: AsRef, + { DebugTools::set_root_level("ERROR")?; let mut reader = Reader { image_reader: ThreadLocal::default(), - path: PathBuf::from(path), + path: path.as_ref().to_path_buf(), series, size_x: 0, size_y: 0, @@ -308,6 +332,17 @@ impl Reader { Ok(reader) } + /// Get ome metadata as ome structure + pub fn get_ome(&self) -> Result { + let mut ome = self.ome_xml()?.parse::()?; + if let Some(image) = ome.image.as_ref() { + if image.len() > 1 { + ome.image = Some(vec![image[self.series].clone()]); + } + } + Ok(ome) + } + /// Get ome metadata as xml string pub fn get_ome_xml(&self) -> Result { self.ome_xml() @@ -339,7 +374,7 @@ impl Reader { /// Retrieve fame at channel c, slize z and time t. pub fn get_frame(&self, c: usize, z: usize, t: usize) -> Result { - let bytes = if self.is_rgb()? & self.is_interleaved()? { + let bytes = if self.is_rgb()? && self.is_interleaved()? { let index = self.get_index(z as i32, 0, t as i32)?; self.deinterleave(self.open_bytes(index)?, c)? } else if self.get_rgb_channel_count()? > 1 { @@ -418,7 +453,7 @@ impl Reader { } } - /// get a slicable view on the image file + /// get a sliceable view on the image file pub fn view(&self) -> View { let slice = s![ 0..self.size_c, diff --git a/src/stats.rs b/src/stats.rs index 24d128f..1ea6ca8 100644 --- a/src/stats.rs +++ b/src/stats.rs @@ -1,4 +1,4 @@ -use anyhow::{anyhow, Result}; +use anyhow::{Result, anyhow}; use ndarray::{Array, ArrayD, ArrayView, Axis, Dimension, RemoveAxis}; /// a trait to define the min, max, sum and mean operations along an axis @@ -196,6 +196,10 @@ macro_rules! impl_frame_stats_int { } impl_frame_stats_float_view!(f32, f64); -impl_frame_stats_int_view!(u8, i8, u16, i16, u32, i32, u64, i64, u128, i128, usize, isize); +impl_frame_stats_int_view!( + u8, i8, u16, i16, u32, i32, u64, i64, u128, i128, usize, isize +); impl_frame_stats_float!(f32, f64); -impl_frame_stats_int!(u8, i8, u16, i16, u32, i32, u64, i64, u128, i128, usize, isize); +impl_frame_stats_int!( + u8, i8, u16, i16, u32, i32, u64, i64, u128, i128, usize, isize +); diff --git a/src/tiff.rs b/src/tiff.rs new file mode 100644 index 0000000..2eed311 --- /dev/null +++ b/src/tiff.rs @@ -0,0 +1,193 @@ +use crate::colors::Color; +use crate::metadata::Metadata; +use crate::reader::PixelType; +use crate::stats::MinMax; +use crate::view::{Number, View}; +use anyhow::{Result, anyhow}; +use indicatif::{ProgressBar, ProgressStyle}; +use itertools::iproduct; +use ndarray::{Array0, Array1, Array2, ArrayD, Dimension}; +use rayon::prelude::*; +use std::path::Path; +use std::sync::{Arc, Mutex}; +use tiffwrite::{Bytes, Colors, Compression, IJTiffFile}; + +#[derive(Clone)] +pub struct TiffOptions { + bar: Option, + compression: Compression, + colors: Option>>, + overwrite: bool, +} + +impl Default for TiffOptions { + fn default() -> Self { + Self { + bar: None, + compression: Compression::Zstd(10), + colors: None, + overwrite: false, + } + } +} + +impl TiffOptions { + pub fn new( + bar: bool, + compression: Option, + colors: Vec, + overwrite: bool, + ) -> Result { + let mut options = Self { + bar: None, + compression: compression.unwrap_or(Compression::Zstd(10)), + colors: None, + overwrite, + }; + if bar { + options.enable_bar()?; + } + if !colors.is_empty() { + options.set_colors(&colors)?; + } + Ok(options) + } + + /// show a progress bar while saving tiff + pub fn enable_bar(&mut self) -> Result<()> { + self.bar = Some(ProgressStyle::with_template( + "{spinner:.green} [{elapsed_precise}, {percent}%] [{wide_bar:.green/lime}] {pos:>7}/{len:7} ({eta_precise}, {per_sec:<5})", + )?.progress_chars("▰▱▱")); + Ok(()) + } + + /// do not show a progress bar while saving tiff + pub fn disable_bar(&mut self) { + self.bar = None; + } + + /// save tiff with zstd compression (default) + pub fn set_zstd_compression(&mut self) { + self.compression = Compression::Zstd(10) + } + + /// save tiff with zstd compression, choose a level between 7..=22 + pub fn set_zstd_compression_level(&mut self, level: i32) { + self.compression = Compression::Zstd(level) + } + + /// save tiff with deflate compression + pub fn set_deflate_compression(&mut self) { + self.compression = Compression::Deflate + } + + pub fn set_colors(&mut self, colors: &[String]) -> Result<()> { + let colors = colors + .iter() + .map(|c| c.parse::()) + .collect::>>()?; + self.colors = Some(colors.into_iter().map(|c| c.to_rgb()).collect()); + Ok(()) + } + + pub fn set_overwrite(&mut self, overwrite: bool) { + self.overwrite = overwrite; + } +} + +impl View +where + D: Dimension, +{ + /// save as tiff with a certain type + pub fn save_as_tiff_with_type(&self, path: P, options: &TiffOptions) -> Result<()> + where + P: AsRef, + T: Bytes + Number + Send + Sync, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, + { + let path = path.as_ref().to_path_buf(); + if path.exists() { + if options.overwrite { + std::fs::remove_file(&path)?; + } else { + return Err(anyhow!("File {} already exists", path.display())); + } + } + let size_c = self.size_c(); + let size_z = self.size_z(); + let size_t = self.size_t(); + let mut tiff = IJTiffFile::new(path)?; + tiff.set_compression(options.compression.clone()); + let ome = self.get_ome()?; + tiff.px_size = ome.pixel_size()?.map(|i| i / 1e3); + tiff.time_interval = ome.time_interval()?.map(|i| i / 1e3); + tiff.delta_z = ome.delta_z()?.map(|i| i / 1e3); + tiff.comment = Some(self.ome_xml()?); + if let Some(mut colors) = options.colors.clone() { + while colors.len() < self.size_c { + colors.push(vec![255, 255, 255]); + } + tiff.colors = Colors::Colors(colors); + } + let tiff = Arc::new(Mutex::new(tiff)); + if let Some(style) = &options.bar { + let bar = ProgressBar::new((size_c as u64) * (size_z as u64) * (size_t as u64)) + .with_style(style.clone()); + iproduct!(0..size_c, 0..size_z, 0..size_t) + .collect::>() + .into_par_iter() + .map(|(c, z, t)| { + if let Ok(mut tiff) = tiff.lock() { + tiff.save(&self.get_frame::(c, z, t)?, c, z, t)?; + bar.inc(1); + Ok(()) + } else { + Err(anyhow::anyhow!("tiff is locked")) + } + }) + .collect::>()?; + bar.finish(); + } else { + iproduct!(0..size_c, 0..size_z, 0..size_t) + .collect::>() + .into_par_iter() + .map(|(c, z, t)| { + if let Ok(mut tiff) = tiff.lock() { + tiff.save(&self.get_frame::(c, z, t)?, c, z, t)?; + Ok(()) + } else { + Err(anyhow::anyhow!("tiff is locked")) + } + }) + .collect::>()?; + }; + Ok(()) + } + + /// save as tiff with whatever pixel type the view has + pub fn save_as_tiff

(&self, path: P, options: &TiffOptions) -> Result<()> + where + P: AsRef, + { + match self.pixel_type { + PixelType::I8 => self.save_as_tiff_with_type::(path, options)?, + PixelType::U8 => self.save_as_tiff_with_type::(path, options)?, + PixelType::I16 => self.save_as_tiff_with_type::(path, options)?, + PixelType::U16 => self.save_as_tiff_with_type::(path, options)?, + PixelType::I32 => self.save_as_tiff_with_type::(path, options)?, + PixelType::U32 => self.save_as_tiff_with_type::(path, options)?, + PixelType::F32 => self.save_as_tiff_with_type::(path, options)?, + PixelType::F64 => self.save_as_tiff_with_type::(path, options)?, + PixelType::I64 => self.save_as_tiff_with_type::(path, options)?, + PixelType::U64 => self.save_as_tiff_with_type::(path, options)?, + PixelType::I128 => self.save_as_tiff_with_type::(path, options)?, + PixelType::U128 => self.save_as_tiff_with_type::(path, options)?, + PixelType::F128 => self.save_as_tiff_with_type::(path, options)?, + } + + Ok(()) + } +} diff --git a/src/view.rs b/src/view.rs index 87a77e4..111ea52 100644 --- a/src/view.rs +++ b/src/view.rs @@ -1,23 +1,25 @@ -use crate::axes::{slice_info, Ax, Axis, Operation, Slice, SliceInfoElemDef}; +use crate::axes::{Ax, Axis, Operation, Slice, SliceInfoElemDef, slice_info}; +use crate::metadata::Metadata; use crate::reader::Reader; use crate::stats::MinMax; -use anyhow::{anyhow, Error, Result}; +use anyhow::{Error, Result, anyhow}; use indexmap::IndexMap; -use itertools::iproduct; +use itertools::{Itertools, iproduct}; use ndarray::{ - s, Array, Array0, Array1, Array2, ArrayD, Dimension, IntoDimension, Ix0, Ix1, Ix2, Ix5, IxDyn, - SliceArg, SliceInfoElem, + Array, Array0, Array1, Array2, ArrayD, Dimension, IntoDimension, Ix0, Ix1, Ix2, Ix5, IxDyn, + SliceArg, SliceInfoElem, s, }; -use num::{Bounded, FromPrimitive, Zero}; +use num::traits::ToBytes; +use num::{Bounded, FromPrimitive, ToPrimitive, Zero}; use serde::{Deserialize, Serialize}; use serde_with::serde_as; use std::any::type_name; use std::collections::HashMap; +use std::fmt::{Display, Formatter}; use std::iter::Sum; use std::marker::PhantomData; -use std::mem::{transmute, transmute_copy}; use std::ops::{AddAssign, Deref, Div}; -use std::path::PathBuf; +use std::path::{Path, PathBuf}; use std::sync::Arc; fn idx_bnd(idx: isize, bnd: isize) -> Result { @@ -65,8 +67,10 @@ impl Number for T where #[derive(Clone, Debug, Serialize, Deserialize)] pub struct View { reader: Arc, + /// same order as axes #[serde_as(as = "Vec")] slice: Vec, + /// always has all of cztyx with possibly some new axes added axes: Vec, operations: IndexMap, dimensionality: PhantomData, @@ -83,13 +87,79 @@ impl View { } } + #[allow(dead_code)] + pub(crate) fn new_with_axes(reader: Arc, axes: Vec) -> Result { + let mut slice = Vec::new(); + for axis in axes.iter() { + match axis { + Axis::C => slice.push(SliceInfoElem::Slice { + start: 0, + end: Some(reader.size_c as isize), + step: 1, + }), + Axis::Z => slice.push(SliceInfoElem::Slice { + start: 0, + end: Some(reader.size_z as isize), + step: 1, + }), + Axis::T => slice.push(SliceInfoElem::Slice { + start: 0, + end: Some(reader.size_t as isize), + step: 1, + }), + Axis::Y => slice.push(SliceInfoElem::Slice { + start: 0, + end: Some(reader.size_y as isize), + step: 1, + }), + Axis::X => slice.push(SliceInfoElem::Slice { + start: 0, + end: Some(reader.size_x as isize), + step: 1, + }), + Axis::New => { + slice.push(SliceInfoElem::NewAxis); + } + } + } + let mut axes = axes.clone(); + for axis in [Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X] { + if !axes.contains(&axis) { + let size = match axis { + Axis::C => reader.size_c, + Axis::Z => reader.size_z, + Axis::T => reader.size_t, + Axis::Y => reader.size_y, + Axis::X => reader.size_x, + Axis::New => 1, + }; + if size > 1 { + return Err(anyhow!( + "Axis {:?} has length {}, but was not included", + axis, + size + )); + } + slice.push(SliceInfoElem::Index(0)); + axes.push(axis); + } + } + Ok(Self { + reader, + slice, + axes, + operations: IndexMap::new(), + dimensionality: PhantomData, + }) + } + /// the file path pub fn path(&self) -> &PathBuf { &self.reader.path } /// the series in the file - pub fn series(&self) -> i32 { + pub fn series(&self) -> usize { self.reader.series } @@ -139,6 +209,11 @@ impl View { &self.axes } + #[allow(dead_code)] + pub(crate) fn get_operations(&self) -> &IndexMap { + &self.operations + } + /// the slice defining the view pub fn get_slice(&self) -> &[SliceInfoElem] { &self.slice @@ -159,6 +234,27 @@ impl View { .collect() } + /// remove axes of size 1 + pub fn squeeze(&self) -> Result> { + let view = self.clone().into_dyn(); + let slice: Vec<_> = self + .shape() + .into_iter() + .map(|s| { + if s == 1 { + SliceInfoElem::Index(0) + } else { + SliceInfoElem::Slice { + start: 0, + end: None, + step: 1, + } + } + }) + .collect(); + view.slice(slice.as_slice()) + } + pub(crate) fn op_axes(&self) -> Vec { self.operations.keys().cloned().collect() } @@ -172,11 +268,27 @@ impl View { } } + /// the number of pixels in the first dimension + pub fn len(&self) -> usize { + self.shape()[0] + } + + pub fn is_empty(&self) -> bool { + self.shape()[0] == 0 + } + /// the number of pixels in the view pub fn size(&self) -> usize { self.shape().into_iter().product() } + pub fn size_ax(&self, ax: Axis) -> Option { + self.axes() + .iter() + .position(|a| *a == ax) + .map(|i| self.shape()[i]) + } + /// the shape of the view pub fn shape(&self) -> Vec { let mut shape = Vec::::new(); @@ -202,6 +314,43 @@ impl View { shape } + pub fn size_of(&self, axis: Axis) -> usize { + if let Some(axis_position) = self.axes.iter().position(|a| *a == axis) { + match self.slice[axis_position] { + SliceInfoElem::Slice { start, end, step } => { + if let Some(e) = end { + ((e - start).max(0) / step) as usize + } else { + panic!("slice has no end") + } + } + _ => 1, + } + } else { + 1 + } + } + + pub fn size_c(&self) -> usize { + self.size_of(Axis::C) + } + + pub fn size_z(&self) -> usize { + self.size_of(Axis::Z) + } + + pub fn size_t(&self) -> usize { + self.size_of(Axis::T) + } + + pub fn size_y(&self) -> usize { + self.size_of(Axis::Y) + } + + pub fn size_x(&self) -> usize { + self.size_of(Axis::X) + } + fn shape_all(&self) -> Vec { let mut shape = Vec::::new(); for s in self.slice.iter() { @@ -227,7 +376,7 @@ impl View { slice.swap(idx0, idx1); let mut axes = self.axes.clone(); axes.swap(idx0, idx1); - Ok(View::new(self.reader.clone(), slice, axes)) + Ok(View::new(self.reader.clone(), slice, axes).with_operations(self.operations.clone())) } /// subset of gives axes will be reordered in given order @@ -244,7 +393,7 @@ impl View { slice[j] = self.slice[i]; axes[j] = self.axes[i]; } - Ok(View::new(self.reader.clone(), slice, axes)) + Ok(View::new(self.reader.clone(), slice, axes).with_operations(self.operations.clone())) } /// reverse the order of the axes @@ -253,17 +402,35 @@ impl View { self.reader.clone(), self.slice.iter().rev().cloned().collect(), self.axes.iter().rev().cloned().collect(), - )) + ) + .with_operations(self.operations.clone())) } fn operate(&self, axis: A, operation: Operation) -> Result> { let pos = axis.pos_op(&self.axes, &self.slice, &self.op_axes())?; - let mut operations = self.operations.clone(); - operations.insert(self.axes[pos], operation); - Ok( - View::new(self.reader.clone(), self.slice.clone(), self.axes.clone()) - .with_operations(operations), - ) + let ax = self.axes[pos]; + let (axes, slice, operations) = if Axis::New == ax { + let mut axes = self.axes.clone(); + let mut slice = self.slice.clone(); + axes.remove(pos); + slice.remove(pos); + (axes, slice, self.operations.clone()) + } else if self.operations.contains_key(&ax) { + if D::NDIM.is_none() { + ( + self.axes.clone(), + self.slice.clone(), + self.operations.clone(), + ) + } else { + return Err(anyhow!("axis {}: {} is already operated on!", pos, ax)); + } + } else { + let mut operations = self.operations.clone(); + operations.insert(ax, operation); + (self.axes.clone(), self.slice.clone(), operations) + }; + Ok(View::new(self.reader.clone(), slice, axes).with_operations(operations)) } /// maximum along axis @@ -300,51 +467,88 @@ impl View { let mut new_slice = Vec::new(); let mut new_axes = Vec::new(); let reader_slice = self.slice.as_slice(); - while (r_idx < reader_slice.len()) & (n_idx < info.len()) { - let n = &info[n_idx]; - let r = &reader_slice[r_idx]; - let a = &self.axes[r_idx]; - if self.operations.contains_key(a) { - new_slice.push(*r); - new_axes.push(*a); - r_idx += 1; - } else { - match (n, r) { + while (r_idx < reader_slice.len()) | (n_idx < info.len()) { + let n = info.get(n_idx); + let r = reader_slice.get(r_idx); + let a = self.axes.get(r_idx); + match a { + Some(i) if self.operations.contains_key(i) => { + new_slice.push(*r.expect("slice should exist for axes under operation")); + new_axes.push(*i); + r_idx += 1; + } + _ => match (n, r) { ( - SliceInfoElem::Slice { + Some(SliceInfoElem::Slice { start: info_start, end: info_end, step: info_step, - }, - SliceInfoElem::Slice { start, end, step }, + }), + Some(SliceInfoElem::Slice { start, end, step }), ) => { let new_start = start + info_start; - let new_end = match (info_end, end) { - (Some(m), Some(n)) => *n.min(&(start + info_step * m)), - (None, Some(n)) => *n, - _ => panic!("slice has no end"), + let end = end.expect("slice has no end"); + let new_end = if let Some(m) = info_end { + end.min(start + info_step * m) + } else { + end }; let new_step = (step * info_step).abs(); + if new_start > end { + return Err(anyhow!( + "Index {} out of bounds {}", + info_start, + (end - start) / step + )); + } new_slice.push(SliceInfoElem::Slice { start: new_start, end: Some(new_end), step: new_step, }); - new_axes.push(*a); + new_axes.push(*a.expect("axis should exist when slice exists")); n_idx += 1; r_idx += 1; } - (SliceInfoElem::Index(k), SliceInfoElem::Slice { start, end, step }) => { - if *k < 0 { - new_slice.push(SliceInfoElem::Index(end.unwrap_or(0) + step.abs() * k)) + ( + Some(SliceInfoElem::Index(k)), + Some(SliceInfoElem::Slice { start, end, step }), + ) => { + let i = if *k < 0 { + end.unwrap_or(0) + step.abs() * k } else { - new_slice.push(SliceInfoElem::Index(start + step.abs() * k)); + start + step.abs() * k + }; + let end = end.expect("slice has no end"); + if i >= end { + return Err(anyhow!( + "Index {} out of bounds {}", + i, + (end - start) / step + )); } - new_axes.push(*a); + new_slice.push(SliceInfoElem::Index(i)); + new_axes.push(*a.expect("axis should exist when slice exists")); n_idx += 1; r_idx += 1; } - (SliceInfoElem::NewAxis, SliceInfoElem::NewAxis) => { + (Some(SliceInfoElem::Slice { start, .. }), Some(SliceInfoElem::NewAxis)) => { + if *start != 0 { + return Err(anyhow!("Index {} out of bounds 1", start)); + } + new_slice.push(SliceInfoElem::NewAxis); + new_axes.push(Axis::New); + n_idx += 1; + r_idx += 1; + } + (Some(SliceInfoElem::Index(k)), Some(SliceInfoElem::NewAxis)) => { + if *k != 0 { + return Err(anyhow!("Index {} out of bounds 1", k)); + } + n_idx += 1; + r_idx += 1; + } + (Some(SliceInfoElem::NewAxis), Some(SliceInfoElem::NewAxis)) => { new_slice.push(SliceInfoElem::NewAxis); new_slice.push(SliceInfoElem::NewAxis); new_axes.push(Axis::New); @@ -352,23 +556,22 @@ impl View { n_idx += 1; r_idx += 1; } - (_, SliceInfoElem::NewAxis) => { - new_slice.push(SliceInfoElem::NewAxis); - new_axes.push(Axis::New); - n_idx += 1; - r_idx += 1; - } - (SliceInfoElem::NewAxis, _) => { + (Some(SliceInfoElem::NewAxis), _) => { new_slice.push(SliceInfoElem::NewAxis); new_axes.push(Axis::New); n_idx += 1; } - (_, SliceInfoElem::Index(k)) => { + (_, Some(SliceInfoElem::Index(k))) => { new_slice.push(SliceInfoElem::Index(*k)); - new_axes.push(*a); + new_axes.push(*a.expect("axis should exist when slice exists")); r_idx += 1; } - } + _ => { + panic!("unreachable"); + // n_idx += 1; + // r_idx += 1; + } + }, } } debug_assert_eq!(r_idx, reader_slice.len()); @@ -382,37 +585,50 @@ impl View { .with_operations(self.operations.clone())) } + /// resets axes to cztyx order, with all 5 axes present, + /// inserts new axes in place of axes under operation (max_proj etc.) + pub fn reset_axes(&self) -> Result> { + let mut axes = Vec::new(); + let mut slice = Vec::new(); + + for ax in [Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X] { + axes.push(ax); + let s = self.slice[ax.pos(&self.axes, &self.slice)?]; + match s { + SliceInfoElem::Slice { .. } => slice.push(s), + SliceInfoElem::Index(i) => slice.push(SliceInfoElem::Slice { + start: i, + end: Some(i + 1), + step: 1, + }), + SliceInfoElem::NewAxis => { + panic!("slice should not be NewAxis when axis is one of cztyx") + } + } + if self.operations.contains_key(&ax) { + axes.push(Axis::New); + slice.push(SliceInfoElem::NewAxis) + } + } + Ok(View::new(self.reader.clone(), slice, axes).with_operations(self.operations.clone())) + } + /// slice, but slice elements are in cztyx order, all cztyx must be given, /// but axes not present in view will be ignored, view axes are reordered in cztyx order pub fn slice_cztyx(&self, info: I) -> Result> where I: SliceArg, { - let axes = self.axes(); - let slice: Vec<_> = info - .as_ref() - .iter() - .zip(self.axes.iter()) - .filter_map(|(&s, ax)| if axes.contains(ax) { Some(s) } else { None }) - .collect(); - let new_axes: Vec<_> = [Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X] - .into_iter() - .filter(|ax| axes.contains(ax)) - .collect(); - self.clone() - .into_dyn() - .permute_axes(&new_axes)? - .slice(slice.as_slice())? - .into_dimensionality() + self.reset_axes()?.slice(info)?.into_dimensionality() } /// the pixel intensity at a given index pub fn item_at(&self, index: &[isize]) -> Result where T: Number, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, { let slice: Vec<_> = index.iter().map(|s| SliceInfoElem::Index(*s)).collect(); let view = self.clone().into_dyn().slice(slice.as_slice())?; @@ -424,9 +640,9 @@ impl View { pub fn as_array(&self) -> Result> where T: Number, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, { Ok(self.as_array_dyn()?.into_dimensionality()?) } @@ -435,9 +651,9 @@ impl View { pub fn as_array_dyn(&self) -> Result> where T: Number, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, { let mut op_xy = IndexMap::new(); if let Some((&ax, op)) = self.operations.first() { @@ -461,18 +677,15 @@ impl View { for (s, a) in self.slice.iter().zip(&self.axes) { match s { SliceInfoElem::Slice { start, end, step } => { - if let Some(e) = end { - if !op_xy.contains_key(a) && !op_czt.contains_key(a) { - shape_out.push(((e - start).max(0) / step) as usize); - slice.push(SliceInfoElem::Slice { - start: 0, - end: None, - step: 1, - }); - ax_out.push(*a); - } - } else { - panic!("slice has no end") + let end = end.expect("slice has no end"); + if !op_xy.contains_key(a) && !op_czt.contains_key(a) { + shape_out.push(((end - start).max(0) / step) as usize); + slice.push(SliceInfoElem::Slice { + start: 0, + end: None, + step: 1, + }); + ax_out.push(*a); } } SliceInfoElem::Index(_) => {} @@ -542,7 +755,9 @@ impl View { let mut axes_out_idx = [None; 5]; for (i, ax) in ax_out.iter().enumerate() { - axes_out_idx[*ax as usize] = Some(i); + if *ax < Axis::New { + axes_out_idx[*ax as usize] = Some(i); + } } for (c, z, t) in iproduct!(&slice_reader[0], &slice_reader[1], &slice_reader[2]) { @@ -570,17 +785,14 @@ impl View { let (&ax1, op1) = op_xy.get_index(1).unwrap(); let a = arr_frame.slice(xys).to_owned(); let b = op0.operate(a, ax0 as usize - 3)?; - let c: &Array1 = unsafe { transmute(&b) }; - let d = op1.operate(c.to_owned(), ax1 as usize - 3)?; - let e: &Array0 = unsafe { transmute(&d) }; - e.to_owned().into_dyn() + let c = op1.operate(b.to_owned(), ax1 as usize - 3)?; + c.to_owned().into_dyn() } else if op_xy.contains_key(&Axis::X) || op_xy.contains_key(&Axis::Y) { let xys = slice_info::(&xy)?; let (&ax, op) = op_xy.first().unwrap(); let a = arr_frame.slice(xys).to_owned(); let b = op.operate(a, ax as usize - 3)?; - let c: &Array0 = unsafe { transmute(&b) }; - c.to_owned().into_dyn() + b.to_owned().into_dyn() } else { let xys = slice_info::(&xy)?; arr_frame.slice(xys).to_owned().into_dyn() @@ -592,8 +804,7 @@ impl View { let (&ax, op) = op_xy.first().unwrap(); let a = arr_frame.slice(xys).to_owned(); let b = op.operate(a, ax as usize - 3)?; - let c: &Array1 = unsafe { transmute(&b) }; - c.to_owned().into_dyn() + b.to_owned().into_dyn() } else { let xys = slice_info::(&xy)?; arr_frame.slice(xys).to_owned().into_dyn() @@ -643,15 +854,21 @@ impl View { } } let mut out = Some(array); - let ax_out: HashMap = ax_out + let mut ax_out: HashMap = ax_out .into_iter() .enumerate() .map(|(i, a)| (a, i)) .collect(); for (ax, op) in self.operations.iter().skip(op_xy.len() + op_czt.len()) { - if let Some(&idx) = ax_out.get(ax) { + if let Some(idx) = ax_out.remove(ax) { + for (_, i) in ax_out.iter_mut() { + if *i > idx { + *i -= 1; + } + } let arr = out.take().unwrap(); - let _ = out.insert(unsafe { transmute_copy(&op.operate(arr, idx)?) }); + let a = op.operate(arr, idx)?; + let _ = out.insert(a); } } let mut n = 1; @@ -671,23 +888,59 @@ impl View { Ok(array) } - /// retrieve a single frame at czt, sliced accordingly - pub fn get_frame(&self, c: isize, z: isize, t: isize) -> Result> + /// turn the view into a 1d array + pub fn flatten(&self) -> Result> where T: Number, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, { + Ok(Array1::from_iter(self.as_array()?.iter().cloned())) + } + + /// turn the data into a byte vector + pub fn to_bytes(&self) -> Result> + where + T: Number + ToBytesVec, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, + { + Ok(self + .as_array()? + .iter() + .flat_map(|i| i.to_bytes_vec()) + .collect()) + } + + /// retrieve a single frame at czt, sliced accordingly + pub fn get_frame(&self, c: N, z: N, t: N) -> Result> + where + T: Number, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, + N: Display + ToPrimitive, + { + let c = c + .to_isize() + .ok_or_else(|| anyhow!("cannot convert {} into isize", c))?; + let z = z + .to_isize() + .ok_or_else(|| anyhow!("cannot convert {} into isize", z))?; + let t = t + .to_isize() + .ok_or_else(|| anyhow!("cannot convert {} into isize", t))?; self.slice_cztyx(s![c, z, t, .., ..])?.as_array() } fn get_stat(&self, operation: Operation) -> Result where T: Number + Sum, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, { let arr: ArrayD = self.as_array_dyn()?; Ok(match operation { @@ -715,9 +968,9 @@ impl View { pub fn max(&self) -> Result where T: Number + Sum, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, { self.get_stat(Operation::Max) } @@ -726,9 +979,9 @@ impl View { pub fn min(&self) -> Result where T: Number + Sum, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, { self.get_stat(Operation::Min) } @@ -737,9 +990,9 @@ impl View { pub fn sum(&self) -> Result where T: Number + Sum, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, { self.get_stat(Operation::Sum) } @@ -748,12 +1001,35 @@ impl View { pub fn mean(&self) -> Result where T: Number + Sum, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, { self.get_stat(Operation::Mean) } + + /// gives a helpful summary of the recorded experiment + pub fn summary(&self) -> Result { + let mut s = "".to_string(); + s.push_str(&format!("path/filename: {}\n", self.path.display())); + s.push_str(&format!("series/pos: {}\n", self.series)); + s.push_str(&format!("dtype: {:?}\n", self.pixel_type)); + let axes = self + .axes() + .into_iter() + .map(|ax| format!("{}", ax)) + .join("") + .to_lowercase(); + let shape = self + .shape() + .into_iter() + .map(|s| format!("{}", s)) + .join(" x "); + let space = " ".repeat(6usize.saturating_sub(axes.len())); + s.push_str(&format!("shape ({}):{}{}\n", axes, space, shape)); + s.push_str(&self.get_ome()?.summary()?); + Ok(s) + } } impl Deref for View { @@ -768,9 +1044,9 @@ impl TryFrom> for Array where T: Number, D: Dimension, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, { type Error = Error; @@ -783,9 +1059,9 @@ impl TryFrom<&View> for Array where T: Number, D: Dimension, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, { type Error = Error; @@ -799,18 +1075,37 @@ pub trait Item { fn item(&self) -> Result where T: Number, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax; + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>; +} + +impl View { + pub fn from_path

(path: P, series: usize) -> Result + where + P: AsRef, + { + let mut path = path.as_ref().to_path_buf(); + if path.is_dir() { + for file in path.read_dir()?.flatten() { + let p = file.path(); + if file.path().is_file() && (p.extension() == Some("tif".as_ref())) { + path = p; + break; + } + } + } + Ok(Reader::new(path, series)?.view()) + } } impl Item for View { fn item(&self) -> Result where T: Number, - ArrayD: MinMax, - Array1: MinMax, - Array2: MinMax, + ArrayD: MinMax>, + Array1: MinMax>, + Array2: MinMax>, { Ok(self .as_array()? @@ -819,3 +1114,36 @@ impl Item for View { .clone()) } } + +impl Display for View { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + if let Ok(summary) = self.summary() { + write!(f, "{}", summary) + } else { + write!(f, "{}", self.path.display()) + } + } +} + +/// trait to convert numbers to bytes +pub trait ToBytesVec { + fn to_bytes_vec(&self) -> Vec; +} + +macro_rules! to_bytes_vec_impl { + ($($t:ty $(,)?)*) => { + $( + impl ToBytesVec for $t { + + #[inline] + fn to_bytes_vec(&self) -> Vec { + self.to_ne_bytes().to_vec() + } + } + )* + }; +} + +to_bytes_vec_impl!( + u8, u16, u32, u64, u128, usize, i8, i16, i32, i64, i128, isize, f32, f64 +);