From 5195ccfcb557126362d784367545222c5ca63658 Mon Sep 17 00:00:00 2001 From: Wim Pomp Date: Sun, 27 Apr 2025 20:07:49 +0200 Subject: [PATCH] - implement sliced views, including min, max, sum and mean operations --- Cargo.toml | 21 +- README.md | 32 +- build.rs | 2 +- py/ndbioimage/__init__.py | 1414 +++++++++++++++++++++-------------- py/ndbioimage/transforms.py | 296 +++++--- src/axes.rs | 218 ++++++ src/bioformats.rs | 18 +- src/lib.rs | 488 +++++------- src/py.rs | 605 ++++++++++++++- src/reader.rs | 442 +++++++++++ src/stats.rs | 201 +++++ src/view.rs | 821 ++++++++++++++++++++ tests/test_open.py | 24 +- tests/test_slicing.py | 15 +- tests/test_ufuncs.py | 37 +- 15 files changed, 3566 insertions(+), 1068 deletions(-) create mode 100644 src/axes.rs create mode 100644 src/reader.rs create mode 100644 src/stats.rs create mode 100644 src/view.rs diff --git a/Cargo.toml b/Cargo.toml index 3eb797f..95e9133 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "ndbioimage" -version = "2025.2.3" +version = "2025.4.0" edition = "2021" rust-version = "1.78.0" authors = ["Wim Pomp "] @@ -19,15 +19,20 @@ name = "ndbioimage" crate-type = ["cdylib", "rlib"] [dependencies] -anyhow = "1.0.95" +anyhow = "1.0.98" +itertools = "0.14.0" +indexmap = { version = "2.9.0", features = ["serde"] } j4rs = "0.22.0" -ndarray = "0.16.1" +ndarray = { version = "0.16.1", features = ["serde"] } num = "0.4.3" -numpy = { version = "0.23.0", optional = true } +numpy = { version = "0.24.0", optional = true } +serde = { version = "1.0.219", features = ["rc"] } +serde_json = { version = "1.0.140", optional = true } +serde_with = "3.12.0" thread_local = "1.1.8" [dependencies.pyo3] -version = "0.23.4" +version = "0.24.2" features = ["extension-module", "abi3-py310", "generate-import-lib", "anyhow"] optional = true @@ -35,10 +40,10 @@ optional = true rayon = "1.10.0" [build-dependencies] -anyhow = "1.0.95" +anyhow = "1.0.98" j4rs = "0.22.0" -retry = "2.0.0" +retry = "2.1.0" [features] -python = ["dep:pyo3", "dep:numpy"] +python = ["dep:pyo3", "dep:numpy", "dep:serde_json"] gpl-formats = [] \ No newline at end of file diff --git a/README.md b/README.md index 78d0171..f8080f0 100644 --- a/README.md +++ b/README.md @@ -74,7 +74,7 @@ use ndarray::Array2; use ndbioimage::Reader; let path = "/path/to/file"; -let reader = Reader::new(&path, 0).unwrap(); +let reader = Reader::new(&path, 0)?; println!("size: {}, {}", reader.size_y, reader.size_y); let frame = reader.get_frame(0, 0, 0).unwrap(); if let Ok(arr) = >>::try_into(frame) { @@ -86,30 +86,22 @@ let xml = reader.get_ome_xml().unwrap(); println!("{}", xml); ``` +``` +use ndarray::Array2; +use ndbioimage::Reader; + +let path = "/path/to/file"; +let reader = Reader::new(&path, 0)?; +let view = reader.view(); +let view = view.max_proj(3)?; +let array = view.as_array::()? +``` + ### Command line ```ndbioimage --help```: show help ```ndbioimage image```: show metadata about image ```ndbioimage image -w {name}.tif -r```: copy image into image.tif (replacing {name} with image), while registering channels ```ndbioimage image -w image.mp4 -C cyan lime red``` copy image into image.mp4 (z will be max projected), make channel colors cyan lime and red -## Adding more formats -Readers for image formats subclass AbstractReader. When an image reader is imported, Imread will -automatically recognize it and use it to open the appropriate file format. Image readers -are required to implement the following methods: - -- staticmethod _can_open(path): return True if path can be opened by this reader -- \_\_frame__(self, c, z, t): return the frame at channel=c, z-slice=z, time=t from the file - -Optional methods: -- get_ome: reads metadata from file and adds them to an OME object imported - from the ome-types library -- open(self): maybe open some file handle -- close(self): close any file handles - -Optional fields: -- priority (int): Imread will try readers with a lower number first, default: 99 -- do_not_pickle (strings): any attributes that should not be included when the object is pickled, - for example: any file handles - # TODO - more image formats diff --git a/build.rs b/build.rs index e1495e7..99f9ae9 100644 --- a/build.rs +++ b/build.rs @@ -8,7 +8,7 @@ use retry::{delay, delay::Exponential, retry}; use j4rs::Jvm; fn main() -> anyhow::Result<()> { - println!("cargo:rerun-if-changed=build.rs"); + println!("cargo::rerun-if-changed=build.rs"); #[cfg(not(feature = "python"))] if std::env::var("DOCS_RS").is_err() { diff --git a/py/ndbioimage/__init__.py b/py/ndbioimage/__init__.py index fac8e54..0f20518 100755 --- a/py/ndbioimage/__init__.py +++ b/py/ndbioimage/__init__.py @@ -12,7 +12,7 @@ from importlib.metadata import version from itertools import product from operator import truediv from pathlib import Path -from typing import Any, Callable, Generator, Iterable, Optional, Sequence, TypeVar +from typing import Any, Callable, Optional, Sequence, TypeVar import numpy as np import yaml @@ -22,25 +22,28 @@ from pint import set_application_registry from tiffwrite import FrameInfo, IJTiffParallel from tqdm.auto import tqdm -from .ndbioimage_rs import Reader, download_bioformats # noqa +from . import ndbioimage_rs as rs # noqa from .transforms import Transform, Transforms # noqa: F401 try: __version__ = version(Path(__file__).parent.name) except Exception: # noqa - __version__ = 'unknown' + __version__ = "unknown" try: - with open(Path(__file__).parent.parent / '.git' / 'HEAD') as g: - head = g.read().split(':')[1].strip() - with open(Path(__file__).parent.parent / '.git' / head) as h: - __git_commit_hash__ = h.read().rstrip('\n') + with open(Path(__file__).parent.parent / ".git" / "HEAD") as g: + head = g.read().split(":")[1].strip() + with open(Path(__file__).parent.parent / ".git" / head) as h: + __git_commit_hash__ = h.read().rstrip("\n") except Exception: # noqa - __git_commit_hash__ = 'unknown' + __git_commit_hash__ = "unknown" -ureg.formatter.default_format = '~P' +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') +warnings.filterwarnings("ignore", "Reference to unknown ID") Number = int | float | np.integer | np.floating @@ -48,17 +51,18 @@ for dep_file in (Path(__file__).parent / "deps").glob("*_"): dep_file.rename(str(dep_file)[:-1]) if not list((Path(__file__).parent / "jassets").glob("bioformats*.jar")): - download_bioformats(True) + rs.download_bioformats(True) class TransformTiff(IJTiffParallel): - """ transform frames in a parallel process to speed up saving """ + """transform frames in a parallel process to speed up saving""" + def __init__(self, image: Imread, *args: Any, **kwargs: Any) -> None: self.image = image super().__init__(*args, **kwargs) def parallel(self, frame: tuple[int, int, int]) -> tuple[FrameInfo]: - return (np.asarray(self.image(*frame)), 0, 0, 0), + return ((np.asarray(self.image(*frame)), 0, 0, 0),) class DequeDict(OrderedDict): @@ -87,9 +91,10 @@ def find(obj: Sequence[Any], **kwargs: Any) -> Any: return item except AttributeError: pass + return None -R = TypeVar('R') +R = TypeVar("R") def try_default(fun: Callable[..., R], default: Any, *args: Any, **kwargs: Any) -> R: @@ -100,7 +105,7 @@ def try_default(fun: Callable[..., R], default: Any, *args: Any, **kwargs: Any) class Shape(tuple): - def __new__(cls, shape: Sequence[int] | Shape, axes: str = 'yxczt') -> Shape: + def __new__(cls, shape: Sequence[int] | Shape, axes: str = "yxczt") -> Shape: if isinstance(shape, Shape): axes = shape.axes # type: ignore new = super().__new__(cls, shape) @@ -117,11 +122,11 @@ class Shape(tuple): @cached_property def yxczt(self) -> tuple[int, int, int, int, int]: - return tuple(self[i] for i in 'yxczt') # type: ignore + return tuple(self[i] for i in "yxczt") # type: ignore class OmeCache(DequeDict): - """ prevent (potentially expensive) rereading of ome data by caching """ + """prevent (potentially expensive) rereading of ome data by caching""" instance = None @@ -155,10 +160,19 @@ class OmeCache(DequeDict): return super().__contains__(self.path_and_lstat(path)) @staticmethod - def path_and_lstat(path: str | Path) -> tuple[Path, Optional[os.stat_result], Optional[os.stat_result]]: + def path_and_lstat( + path: str | Path, + ) -> tuple[Path, Optional[os.stat_result], Optional[os.stat_result]]: path = Path(path) - return (path, (path.lstat() if path.exists() else None), - (path.with_suffix('.ome.xml').lstat() if path.with_suffix('.ome.xml').exists() else None)) + return ( + path, + (path.lstat() if path.exists() else None), + ( + path.with_suffix(".ome.xml").lstat() + if path.with_suffix(".ome.xml").exists() + else None + ), + ) def get_positions(path: str | Path) -> Optional[list[int]]: # noqa @@ -167,42 +181,35 @@ def get_positions(path: str | Path) -> Optional[list[int]]: # noqa class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): - """ class to read image files, while taking good care of important metadata, - currently optimized for .czi files, but can open anything that bioformats can handle - path: path to the image file - optional: - axes: order of axes, default: cztyx, but omitting any axes with lenght 1 - dtype: datatype to be used when returning frames + """class to read image files, while taking good care of important metadata, + currently optimized for .czi files, but can open anything that bioformats can handle + path: path to the image file + optional: + axes: order of axes, default: cztyx, but omitting any axes with lenght 1 + dtype: datatype to be used when returning frames - modify images on the fly with a decorator function: - define a function which takes an instance of this object, one image frame, - and the coordinates c, z, t as arguments, and one image frame as return - >> imread.frame_decorator = fun - then use imread as usually + Examples: + >> im = Imread('/path/to/file.image', axes='czt) + >> im + << shows summary + >> im.shape + << (15, 26, 1000, 1000) + >> im.axes + << 'ztyx' + >> plt.imshow(im[1, 0]) + << plots frame at position z=1, t=0 (python type indexing) + >> plt.imshow(im[:, 0].max('z')) + << plots max-z projection at t=0 + >> im.pxsize_um + << 0.09708737864077668 image-plane pixel size in um + >> im.laserwavelengths + << [642, 488] + >> im.laserpowers + << [0.02, 0.0005] in % - Examples: - >> im = Imread('/path/to/file.image', axes='czt) - >> im - << shows summary - >> im.shape - << (15, 26, 1000, 1000) - >> im.axes - << 'ztyx' - >> plt.imshow(im[1, 0]) - << plots frame at position z=1, t=0 (python type indexing) - >> plt.imshow(im[:, 0].max('z')) - << plots max-z projection at t=0 - >> im.pxsize_um - << 0.09708737864077668 image-plane pixel size in um - >> im.laserwavelengths - << [642, 488] - >> im.laserpowers - << [0.02, 0.0005] in % - - See __init__ and other functions for more ideas. + See __init__ and other functions for more ideas. """ - do_not_pickle = 'cache', 'reader' ureg = ureg def close(self) -> None: @@ -213,79 +220,94 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): # TODO return None - def __init__(self, path: Path | str | Imread | Any = None, dtype: DTypeLike = None, axes: str = None) -> None: + def __init__( + self, path: Path | str = None, dtype: DTypeLike = None, axes: str = None + ) -> None: super().__init__() - if isinstance(path, Imread): - self.base = path.base or path - else: - self.base = self - self.isclosed = False - if isinstance(path, str): - path = Path(path) - self.path, self.series = self.split_path_series(path) - self.reader = Reader(str(self.path), int(self.series)) - self.frame_decorator = None - self.transform = Transforms() + 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') + 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.acquisitiondate = "now" self.pcf = None self.powermode = None self.collimator = None self.tirfangle = None - self.cyllens = ['None', 'None'] - self.duolink = 'None' + self.cyllens = ["None", "None"] + self.duolink = "None" self.detector = [0, 1] self.track = [0] - self.cache = DequeDict(16) self.frameoffset = 0, 0 # how far apart the centers of frame and sensor are # extract some metadata from ome instrument = self.ome.instruments[0] if self.ome.instruments else None image = self.ome.images[self.series if len(self.ome.images) > 1 else 0] pixels = image.pixels - self._shape = Shape((pixels.size_y, pixels.size_x, pixels.size_c, pixels.size_z, pixels.size_t), - 'yxczt') - self.base_shape = Shape((pixels.size_y, pixels.size_x, pixels.size_c, pixels.size_z, pixels.size_t), - 'yxczt') - self.dtype = pixels.type.value if dtype is None else dtype + self.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'])) + 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 + 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) + 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 + 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')] + 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.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()] + optovars = [ + objective + for objective in instrument.objectives + if "tubelens" in objective.id.lower() + ] except AttributeError: optovars = [] if len(optovars) == 0: @@ -294,7 +316,10 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): self.tubelens = optovars[0] if self.objective: if self.tubelens: - self.magnification = self.objective.nominal_magnification * self.tubelens.nominal_magnification + self.magnification = ( + self.objective.nominal_magnification + * self.tubelens.nominal_magnification + ) else: self.magnification = self.objective.nominal_magnification self.NA = self.objective.lens_na @@ -302,61 +327,91 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): 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.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] + 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' + self.axes = "".join(i for i in "cztyx" if self.shape[i] > 1) + elif axes.lower() == "full": + self.axes = "cztyx" else: self.axes = axes - self.slice = [np.arange(s, dtype=int) for s in self.shape.yxczt] m = self.extrametadata if m is not None: try: - self.cyllens = m['CylLens'] - self.duolink = m['DLFilterSet'].split(' & ')[m['DLFilterChannel']] - if 'FeedbackChannels' in m: - self.feedback = m['FeedbackChannels'] + self.cyllens = m["CylLens"] + self.duolink = m["DLFilterSet"].split(" & ")[m["DLFilterChannel"]] + if "FeedbackChannels" in m: + self.feedback = m["FeedbackChannels"] else: - self.feedback = m['FeedbackChannel'] + self.feedback = m["FeedbackChannel"] except Exception: # noqa - self.cyllens = ['None', 'None'] - self.duolink = 'None' + 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() + 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 + 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.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'] + self.sigma = [2] * self.shape["c"] if not self.NA: self.immersionN = 1 elif 1.5 < self.NA: @@ -368,98 +423,67 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): else: self.immersionN = 1 - p = re.compile(r'(\d+):(\d+)$') + 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'])]) + 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 made with slice() or np.s_ """ - return self[{k: slice(v) if v is None else v for k, v in dict(c=c, z=z, t=t, x=x, y=y).items()}] + def __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: - def unique_yield(a: Iterable[Any], b: Iterable[Any]) -> Generator[Any, None, None]: - for k in a: - yield k - for k in b: - if k not in a: - yield k - - for idx in unique_yield([key[:3] for key in self.cache.keys()], - product(range(self.shape['c']), range(self.shape['z']), range(self.shape['t']))): - yxczt = (slice(None), slice(None)) + idx - in_idx = tuple(yxczt['yxczt'.find(i)] for i in self.axes) - if item in np.asarray(self[in_idx]): - return True - return False + raise NotImplementedError def __enter__(self) -> Imread: return self def __exit__(self, *args: Any, **kwargs: Any) -> None: - if not self.isclosed: - self.isclosed = True - if hasattr(self, 'close'): - self.close() + 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 | np.ndarray: - """ slice like a numpy array but return an Imread instance """ - if self.isclosed: - raise OSError('file is closed') - if isinstance(n, (slice, Number)): # None = : - n = (n,) - elif isinstance(n, type(Ellipsis)): - n = (None,) * len(self.axes) - elif isinstance(n, dict): # allow im[dict(z=0)] etc. - n = [n.get(i, slice(None)) for i in self.axes] - n = list(n) - - # deal with ... - ell = [i for i, e in enumerate(n) if isinstance(e, type(Ellipsis))] - if len(ell) > 1: - raise IndexError('an index can only have a single ellipsis (...)') - if len(ell): - if len(n) > self.ndim: - n.remove(Ellipsis) - else: - n[ell[0]] = None - while len(n) < self.ndim: - n.insert(ell[0], None) - while len(n) < self.ndim: - n.append(None) - - axes_idx = [self.shape.axes.find(i) for i in 'yxczt'] - n = [n[j] if 0 <= j < len(n) else None for j in axes_idx] # reorder n - - new_slice = [] - for s, e in zip(self.slice, n): - if e is None: - new_slice.append(s) - else: - new_slice.append(s[e]) - - # TODO: check output dimensionality when requested shape in some dimension is 1 - if all([isinstance(s, Number) or a < 0 and s.size == 1 for s, a in zip(new_slice, axes_idx)]): - return self.block(*new_slice).item() # type: ignore - else: - new = View(self) - new.slice = new_slice - new._shape = Shape([1 if isinstance(s, Number) else len(s) for s in new_slice]) - new.axes = ''.join(j for j in self.axes if j in [i for i, s in zip('yxczt', new_slice) - if not isinstance(s, Number)]) + 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 - - def __getstate__(self) -> dict[str: Any]: - return ({key: value for key, value in self.__dict__.items() if key not in self.do_not_pickle} | - {'cache_size': self.cache.maxlen}) + else: + return item def __len__(self) -> int: return self.shape[0] @@ -467,15 +491,18 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): def __repr__(self) -> str: return self.summary - def __setstate__(self, state: dict[str, Any]) -> None: - """ What happens during unpickling """ - self.__dict__.update({key: value for key, value in state.items() if key != 'cache_size'}) - self.reader = Reader(str(self.path), int(self.series)) - self.cache = DequeDict(state.get('cache_size', 16)) - def __str__(self) -> str: return str(self.path) + @staticmethod + def as_axis(axis): + if axis is None: + return None + elif isinstance(axis, int): + return axis + 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()) @@ -483,30 +510,28 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): def tobytes(self) -> bytes: return self.flatten().tobytes() - def __array__(self, dtype: DTypeLike = None, copy: bool = None) -> np.ndarray: - if copy is False: - raise ValueError("`copy=False` isn't supported. A copy is always created.") - block = self.block(*self.slice) - axes_idx = [self.shape.axes.find(i) for i in 'yxczt'] - axes_squeeze = tuple({i for i, j in enumerate(axes_idx) if j == -1}.union( - {i for i, j in enumerate(self.slice) if isinstance(j, Number)})) - block = block.squeeze(axes_squeeze) - if dtype is not None: - block = block.astype(dtype) - if block.ndim == 0: - return block.item() # type: ignore - axes = ''.join(j for i, j in enumerate('yxczt') if i not in axes_squeeze) - return block.transpose([axes.find(i) for i in self.shape.axes if i in axes]) + def __array__(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 """ + 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'])): + 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) + 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] @@ -514,15 +539,21 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): 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) + 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) + arg = np.ravel_multi_index( + [arg[axes.find(i)] for i in self.axes], self.shape + ) if out is None: return arg else: @@ -534,29 +565,37 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): 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}.') + 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'])): + 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) + 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'])): + 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) + 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: + 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: @@ -566,13 +605,20 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): 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') + 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] @@ -581,21 +627,29 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): 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'])): + 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) + 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)] + 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)) + 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: @@ -610,7 +664,7 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): 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}.') + raise IndexError(f"Axis {axis_str} not in {self.axes}.") out_shape = list(self.shape) out_axes = list(self.axes) if not keepdims: @@ -618,63 +672,105 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): 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' + 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'])): + 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)) + 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'])): + 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) + 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],) + 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 + 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]) + 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)) + 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 self.shape.axes + return "".join(self.reader.axes).lower() @axes.setter def axes(self, value: str) -> None: - shape = self.shape[value] - if isinstance(shape, Number): - shape = (shape,) - self._shape = Shape(shape, value) + 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 self._dtype + return np.dtype(self.reader.dtype.lower()) @dtype.setter def dtype(self, value: DTypeLike) -> None: - self._dtype = np.dtype(value) + 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') + 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: @@ -685,65 +781,84 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): @property def ndim(self) -> int: - return len(self.shape) + return self.reader.ndim @property def size(self) -> int: - return np.prod(self.shape) # type: ignore + return self.reader.size @property def shape(self) -> Shape: - return self._shape - - @shape.setter - def shape(self, value: Shape | tuple[int, ...]) -> None: - if isinstance(value, Shape): - self._shape = value - else: - self._shape = Shape([value['yxczt'.find(i.lower())] for i in self.axes], self.axes) + 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)}")) + """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') + 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') + 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') + 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') + s.append(f"time interval: {self.timeinterval:.3f} s") if self.binning: - s.append('binning: {}x{}'.format(*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') + 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]) + ' %') + 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}') + s.append(f"objective: {self.objective.model}") if self.magnification: - s.append(f'magnification: {self.magnification}x') + s.append(f"magnification: {self.magnification}x") if self.tubelens and self.tubelens.model: - s.append(f'tubelens: {self.tubelens.model}') + s.append(f"tubelens: {self.tubelens.model}") if self.filter: - s.append(f'filterset: {self.filter}') + s.append(f"filterset: {self.filter}") if self.powermode: - s.append(f'powermode: {self.powermode}') + s.append(f"powermode: {self.powermode}") if self.collimator: - s.append('collimator: ' + (' {}' * len(self.collimator)).format(*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)) + s.append( + "TIRF angle: " + + (" {:.2f}°" * len(self.tirfangle)).format(*self.tirfangle) + ) if self.gain: - s.append('gain: ' + (' {:.0f}' * len(self.gain)).format(*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) + s.append("pcf: " + (" {:.2f}" * len(self.pcf)).format(*self.pcf)) + return "\n".join(s) @property def T(self) -> Imread: # noqa @@ -751,11 +866,11 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): @cached_property def timeseries(self) -> bool: - return self.shape['t'] > 1 + return self.shape["t"] > 1 @cached_property def zstack(self) -> bool: - return self.shape['z'] > 1 + return self.shape["z"] > 1 @wraps(np.ndarray.argmax) def argmax(self, *args, **kwargs): @@ -766,36 +881,61 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): return self.__array_arg_fun__(np.argmin, *args, **kwargs) @wraps(np.ndarray.max) - def max(self, axis=None, out=None, keepdims=False, initial=None, where=True, **_): - return self.__array_fun__([np.max], axis, None, out, keepdims, [initial], where) + 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, dtype=None, out=None, keepdims=False, *, where=True, **_): - dtype = dtype or float - n = np.prod(self.shape) if axis is None else self.shape[axis] - - def sfun(frame: ArrayLike) -> np.ndarray: - return np.asarray(frame).astype(float) - - def cfun(res: np.ndarray) -> np.ndarray: - return res / n - - return self.__array_fun__([np.sum], axis, dtype, out, keepdims, None, where, [sfun], cfun) + 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, out=None, keepdims=False, initial=None, where=True, **_): - return self.__array_fun__([np.min], axis, None, out, keepdims, [initial], where) + 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') + 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) + 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, **_): + def nanmean( + self, axis=None, dtype=None, out=None, keepdims=False, *, where=True, **_ + ): dtype = dtype or float def sfun(frame): @@ -804,22 +944,59 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): 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) + 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) + 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) + 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): + 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): + def nanvar( + self, + axis=None, + dtype=None, + out=None, + ddof=0, + keepdims=None, + *, + where=True, + std=False, + ): dtype = dtype or float def sfun(frame): @@ -832,20 +1009,32 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): return np.invert(np.isnan(frame)) if std: + def cfun(s, s2, n): - return np.sqrt((s2 - s ** 2 / n) / (n - ddof)) + 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) + 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): + def flatten(self, *args, **kwargs) -> np.ndarray: return np.asarray(self).flatten(*args, **kwargs) @wraps(np.ndarray.reshape) - def reshape(self, *args, **kwargs): + def reshape(self, *args, **kwargs) -> np.ndarray: return np.asarray(self).reshape(*args, **kwargs) # noqa @wraps(np.ndarray.squeeze) @@ -856,43 +1045,49 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): elif isinstance(axes, Number): axes = (axes,) else: - axes = tuple(new.axes.find(ax) if isinstance(ax, str) else ax for ax in axes) + 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) + 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): + def std( + self, axis=None, dtype=None, out=None, ddof=0, keepdims=None, *, where=True + ): return self.var(axis, dtype, out, ddof, keepdims, where=where, std=True) # type: ignore - @wraps(np.ndarray.sum) - def sum(self, axis=None, dtype=None, out=None, keepdims=False, initial=None, where=True, **_): - return self.__array_fun__([np.sum], axis, dtype, out, keepdims, [initial], where) - @wraps(np.ndarray.swapaxes) def swapaxes(self, axis1, axis2): - new = self.copy() - axes = new.axes - if isinstance(axis1, str): - axis1 = axes.find(axis1) - if isinstance(axis2, str): - axis2 = axes.find(axis2) - axes[axis1], axes[axis2] = axes[axis2], axes[axis1] - new.axes = axes + 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): - new = self.copy() - if not axes: - new.axes = new.axes[::-1] - else: - new.axes = ''.join(ax if isinstance(ax, str) else new.axes[ax] for ax in axes) + 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): + 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] @@ -903,12 +1098,25 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): return np.asarray(frame).astype(float) ** 2 if std: + def cfun(s, s2): - return np.sqrt((s2 - s ** 2 / n) / (n - ddof)) + 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) + return (s2 - s**2 / n) / (n - ddof) + + return self.__array_fun__( + [np.sum, np.sum], + axis, + dtype, + out, + keepdims, + None, + where, + (sfun, s2fun), + cfun, + ) def asarray(self) -> np.ndarray: return self.__array__() @@ -919,78 +1127,57 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): new.dtype = dtype return new - def block(self, y: int | Sequence[int] = None, x: int | Sequence[int] = None, - c: int | Sequence[int] = None, z: int | Sequence[int] = None, - t: int | Sequence[int] = None) -> np.ndarray: - """ returns 5D block of frames """ - y, x, c, z, t = (np.arange(self.shape[i]) if e is None else np.array(e, ndmin=1) - for i, e in zip('yxczt', (y, x, c, z, t))) - d = np.empty((len(y), len(x), len(c), len(z), len(t)), self.dtype) - for (ci, cj), (zi, zj), (ti, tj) in product(enumerate(c), enumerate(z), enumerate(t)): - d[:, :, ci, zi, ti] = self.frame(cj, zj, tj)[y][:, x] # type: ignore - return d + def copy(self) -> Imread: + new = self.new() + new.reader = self.reader.copy() + return new - def copy(self) -> View: - return View(self) - - def data(self, c: int | Sequence[int] = 0, z: int | Sequence[int] = 0, t: int | Sequence[int] = 0) -> np.ndarray: - """ returns 3D stack of frames """ - c, z, t = (np.arange(self.shape[i]) if e is None else np.array(e, ndmin=1) for i, e in zip('czt', (c, z, t))) - return np.dstack([self.frame(ci, zi, ti) for ci, zi, ti in product(c, z, t)]) - - def frame(self, c: int = 0, z: int = 0, t: int = 0) -> np.ndarray: - """ returns single 2D frame """ - c = self.get_channel(c) - c %= self.base.shape['c'] - z %= self.base.shape['z'] - t %= self.base.shape['t'] - - # cache last n (default 16) frames in memory for speed (~250x faster) - key = (c, z, t, self.transform, self.frame_decorator) - if self.cache.maxlen and key in self.cache: - self.cache.move_to_end(key) - f = self.cache[key] - else: - f = self.transform[self.channel_names[c], t].frame(self.reader.get_frame(int(c), int(z), int(t))) - if self.frame_decorator is not None: - f = self.frame_decorator(self, f, c, z, t) - if self.cache.maxlen: - self.cache[key] = f - if self.dtype is not None: - return f.copy().astype(self.dtype) - else: - return f.copy() + def 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}' + c = [ + i + for i, c in enumerate(self.channel_names) + if c.lower().startswith(channel_name.lower()) + ] + assert len(c) > 0, f"Channel {c} not found in {self.channel_names}" + assert len(c) < 2, f"Channel {c} not unique in {self.channel_names}" return c[0] @staticmethod def get_config(file: Path | str) -> Any: - """ Open a yml config file """ + """Open a yml config file""" loader = yaml.SafeLoader loader.add_implicit_resolver( - r'tag:yaml.org,2002:float', - re.compile(r'''^(?: + 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.')) + |\.(?: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]]: + 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)): + for i, n in zip("czt", (c, z, t)): if n is None: czt.append(list(range(self.shape[i]))) elif isinstance(n, range): @@ -1012,27 +1199,42 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): # fix ome if necessary for image in ome.images: try: - if image.pixels.physical_size_z is None and len(set([plane.the_z - for plane in image.pixels.planes])) > 1: - z = np.array([(plane.position_z * ureg.Quantity(plane.position_z_unit.value).to(ureg.m).magnitude, - plane.the_z) - for plane in image.pixels.planes if plane.the_c == 0 and plane.the_t == 0]) + if ( + image.pixels.physical_size_z is None + and len(set([plane.the_z for plane in image.pixels.planes])) > 1 + ): + z = np.array( + [ + ( + plane.position_z + * ureg.Quantity(plane.position_z_unit.value) + .to(ureg.m) + .magnitude, + plane.the_z, + ) + for plane in image.pixels.planes + if plane.the_c == 0 and plane.the_t == 0 + ] + ) i = np.argsort(z[:, 1]) - image.pixels.physical_size_z = np.nanmean(np.true_divide(*np.diff(z[i], axis=0).T)) * 1e6 - image.pixels.physical_size_z_unit = 'µm' # type: ignore - except Exception: # noqa + image.pixels.physical_size_z = ( + np.nanmean(np.true_divide(*np.diff(z[i], axis=0).T)) * 1e6 + ) + image.pixels.physical_size_z_unit = "µm" # type: ignore + except Exception: # noqa pass return ome @staticmethod - def read_ome(path: [str, Path]) -> Optional[OME]: - path = Path(path) - if path.with_suffix('.ome.xml').exists(): - return OME.from_xml(path.with_suffix('.ome.xml')) + def 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 None def get_ome(self) -> OME: - """ OME metadata structure """ - return from_xml(self.reader.get_ome_xml(), parser='lxml') + """OME metadata structure""" + return from_xml(self.reader.get_ome_xml(), validate=False) @cached_property def ome(self) -> OME: @@ -1045,38 +1247,42 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): return cache[self.path] def is_noise(self, volume: ArrayLike = None) -> bool: - """ True if volume only has noise """ + """True if volume only has noise""" if volume is None: volume = self fft = np.fft.fftn(volume) - corr = np.fft.fftshift(np.fft.ifftn(fft * fft.conj()).real / np.sum(volume ** 2)) + corr = np.fft.fftshift(np.fft.ifftn(fft * fft.conj()).real / np.sum(volume**2)) return 1 - corr[tuple([0] * corr.ndim)] < 0.0067 @staticmethod def kill_vm() -> None: pass - def new(self, *args: Any, **kwargs: Any) -> View: - warnings.warn('Imread.new has been deprecated, use Imread.view instead.', DeprecationWarning, 2) - return self.view(*args, **kwargs) - - def save_as_movie(self, fname: Path | str = None, - c: int | Sequence[int] = None, z: int | Sequence[int] = None, # noqa - t: str | int | Sequence[int] = None, # noqa - colors: tuple[str] = None, brightnesses: tuple[float] = None, - scale: int = None, bar: bool = True) -> None: - """ saves the image as a mp4 or mkv file """ + def save_as_movie( + self, + fname: Path | str = None, + c: int | Sequence[int] = None, # noqa + z: int | Sequence[int] = None, # noqa + t: str | int | Sequence[int] = None, # noqa + colors: tuple[str] = None, + brightnesses: tuple[float] = None, + scale: int = None, + bar: bool = True, + ) -> None: + """saves the image as a mp4 or mkv file""" from matplotlib.colors import to_rgb from skvideo.io import FFmpegWriter if t is None: - t = np.arange(self.shape['t']) + t = np.arange(self.shape["t"]) elif isinstance(t, str): t = eval(f"np.arange(self.shape['t'])[{t}]") elif np.isscalar(t): t = (t,) - def get_ab(tyx: Imread, p: tuple[float, float] = (1, 99)) -> tuple[float, float]: + def get_ab( + tyx: Imread, p: tuple[float, float] = (1, 99) + ) -> tuple[float, float]: s = tyx.flatten() s = s[s > 0] a, b = np.percentile(s, p) @@ -1086,173 +1292,253 @@ class Imread(np.lib.mixins.NDArrayOperatorsMixin, ABC): a, b = 0, 1 return a, b - def cframe(frame: ArrayLike, color: str, a: float, b: float, scale: float = 1) -> np.ndarray: # noqa + def cframe( + frame: ArrayLike, color: str, a: float, b: float, scale: float = 1 # noqa + ) -> np.ndarray: color = to_rgb(color) frame = (frame - a) / (b - a) frame = np.dstack([255 * frame * i for i in color]) - return np.clip(np.round(frame), 0, 255).astype('uint8') + return np.clip(np.round(frame), 0, 255).astype("uint8") - ab = list(zip(*[get_ab(i) for i in self.transpose('cztyx')])) # type: ignore - colors = colors or ('r', 'g', 'b')[:self.shape['c']] + max(0, self.shape['c'] - 3) * ('w',) - brightnesses = brightnesses or (1,) * self.shape['c'] + ab = list(zip(*[get_ab(i) for i in self.transpose("cztyx")])) # type: ignore + colors = colors or ("r", "g", "b")[: self.shape["c"]] + max( + 0, self.shape["c"] - 3 + ) * ("w",) + brightnesses = brightnesses or (1,) * self.shape["c"] scale = scale or 1 - shape_x = 2 * ((self.shape['x'] * scale + 1) // 2) - shape_y = 2 * ((self.shape['y'] * scale + 1) // 2) + shape_x = 2 * ((self.shape["x"] * scale + 1) // 2) + shape_y = 2 * ((self.shape["y"] * scale + 1) // 2) with FFmpegWriter( - str(fname).format(name=self.path.stem, path=str(self.path.parent)), - outputdict={'-vcodec': 'libx264', '-preset': 'veryslow', '-pix_fmt': 'yuv420p', '-r': '7', - '-vf': f'setpts={25 / 7}*PTS,scale={shape_x}:{shape_y}:flags=neighbor'} + str(fname).format(name=self.path.stem, path=str(self.path.parent)), + outputdict={ + "-vcodec": "libx264", + "-preset": "veryslow", + "-pix_fmt": "yuv420p", + "-r": "7", + "-vf": f"setpts={25 / 7}*PTS,scale={shape_x}:{shape_y}:flags=neighbor", + }, ) as movie: - im = self.transpose('tzcyx') # type: ignore - for ti in tqdm(t, desc='Saving movie', disable=not bar): - movie.writeFrame(np.max([cframe(yx, c, a, b / s, scale) - for yx, a, b, c, s in zip(im[ti].max('z'), *ab, colors, brightnesses)], 0)) + im = self.transpose("tzcyx") # type: ignore + for ti in tqdm(t, desc="Saving movie", disable=not bar): + movie.writeFrame( + np.max( + [ + cframe(yx, c, a, b / s, scale) + for yx, a, b, c, s in zip( + im[ti].max("z"), *ab, colors, brightnesses + ) + ], + 0, + ) + ) - def save_as_tiff(self, fname: Path | str = None, c: int | Sequence[int] = None, z: int | Sequence[int] = None, - t: int | Sequence[int] = None, split: bool = False, bar: bool = True, pixel_type: str = 'uint16', - **kwargs: Any) -> None: - """ saves the image as a tif file - split: split channels into different files """ + def save_as_tiff( + self, + fname: Path | str = None, + c: int | Sequence[int] = None, + z: int | Sequence[int] = None, + t: int | Sequence[int] = None, + split: bool = False, + bar: bool = True, + pixel_type: str = "uint16", + **kwargs: Any, + ) -> None: + """saves the image as a tif file + split: split channels into different files""" fname = Path(str(fname).format(name=self.path.stem, path=str(self.path.parent))) if fname is None: - fname = self.path.with_suffix('.tif') + fname = self.path.with_suffix(".tif") if fname == self.path: - raise FileExistsError(f'File {fname} exists already.') + raise FileExistsError(f"File {fname} exists already.") if not isinstance(fname, Path): fname = Path(fname) if split: - for i in range(self.shape['c']): + for i in range(self.shape["c"]): if self.timeseries: - self.save_as_tiff(fname.with_name(f'{fname.stem}_C{i:01d}').with_suffix('.tif'), i, 0, None, False, - bar, pixel_type) + self.save_as_tiff( + fname.with_name(f"{fname.stem}_C{i:01d}").with_suffix(".tif"), + i, + 0, + None, + False, + bar, + pixel_type, + ) else: - self.save_as_tiff(fname.with_name(f'{fname.stem}_C{i:01d}').with_suffix('.tif'), i, None, 0, False, - bar, pixel_type) + self.save_as_tiff( + fname.with_name(f"{fname.stem}_C{i:01d}").with_suffix(".tif"), + i, + None, + 0, + False, + bar, + pixel_type, + ) else: n = [c, z, t] - for i, ax in enumerate('czt'): + for i, ax in enumerate("czt"): if n[i] is None: n[i] = range(self.shape[ax]) elif not isinstance(n[i], (tuple, list)): n[i] = (n[i],) shape = [len(i) for i in n] - with TransformTiff(self, fname.with_suffix('.tif'), dtype=pixel_type, - pxsize=self.pxsize_um, deltaz=self.deltaz_um, **kwargs) as tif: - for i, m in tqdm(zip(product(*[range(s) for s in shape]), product(*n)), # noqa - total=np.prod(shape), desc='Saving tiff', disable=not bar): - tif.save(m, *i) + with TransformTiff( + self, + fname.with_suffix(".tif"), + dtype=pixel_type, + pxsize=self.pxsize_um, + deltaz=self.deltaz_um, + **kwargs, + ) as tif: + for i, m in tqdm( # noqa + zip(product(*[range(s) for s in shape]), product(*n)), # noqa + total=np.prod(shape), + desc="Saving tiff", + disable=not bar, + ): + tif.save(m, *i) # type: ignore - def with_transform(self, channels: bool = True, drift: bool = False, file: Path | str = None, - bead_files: Sequence[Path | str] = ()) -> View: - """ returns a view where channels and/or frames are registered with an affine transformation - channels: True/False register channels using bead_files - drift: True/False register frames to correct drift - file: load registration from file with name file, default: transform.yml in self.path.parent - bead_files: files used to register channels, default: files in self.path.parent, - with names starting with 'beads' - """ - view = self.view() - if file is None: - file = Path(view.path.parent) / 'transform.yml' - else: - file = Path(file) - if not bead_files: - try: - bead_files = Transforms.get_bead_files(view.path.parent) - except Exception: # noqa - if not file.exists(): - raise Exception('No transform file and no bead file found.') - bead_files = () - - if channels: - try: - view.transform = Transforms.from_file(file, T=drift) - except Exception: # noqa - view.transform = Transforms().with_beads(view.cyllens, bead_files) - if drift: - view.transform = view.transform.with_drift(view) - view.transform.save(file.with_suffix('.yml')) - view.transform.save_channel_transform_tiff(bead_files, file.with_suffix('.tif')) - elif drift: - try: - view.transform = Transforms.from_file(file, C=False) - except Exception: # noqa - view.transform = Transforms().with_drift(self) - view.transform.adapt(view.frameoffset, view.shape.yxczt, view.channel_names) - return view - - def set_cache_size(self, cache_size: int) -> None: - assert isinstance(cache_size, int) and cache_size >= 0 - self.cache.maxlen = cache_size - self.cache.truncate() + def with_transform( + self, + channels: bool = True, + drift: bool = False, + file: Path | str = None, + bead_files: Sequence[Path | str] = (), + ) -> Imread: + """returns a view where channels and/or frames are registered with an affine transformation + channels: True/False register channels using bead_files + drift: True/False register frames to correct drift + file: load registration from file with name file, default: transform.yml in self.path.parent + bead_files: files used to register channels, default: files in self.path.parent, + with names starting with 'beads' + """ + raise NotImplementedError("transforms are not yet implemented") + # view = self.copy() + # if file is None: + # file = Path(view.path.parent) / 'transform.yml' + # else: + # file = Path(file) + # if not bead_files: + # try: + # bead_files = Transforms.get_bead_files(view.path.parent) + # except Exception: # noqa + # if not file.exists(): + # raise Exception('No transform file and no bead file found.') + # bead_files = () + # + # if channels: + # try: + # view.transform = Transforms.from_file(file, T=drift) + # except Exception: # noqa + # view.transform = Transforms().with_beads(view.cyllens, bead_files) + # if drift: + # view.transform = view.transform.with_drift(view) + # view.transform.save(file.with_suffix('.yml')) + # view.transform.save_channel_transform_tiff(bead_files, file.with_suffix('.tif')) + # elif drift: + # try: + # view.transform = Transforms.from_file(file, C=False) + # except Exception: # noqa + # view.transform = Transforms().with_drift(self) + # view.transform.adapt(view.frameoffset, view.shape.yxczt, view.channel_names) + # return view @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')) + if ( + isinstance(path, Path) + and path.name.startswith("Pos") + and path.name.lstrip("Pos").isdigit() + ): + return path.parent, int(path.name.lstrip("Pos")) return path, 0 - def view(self, *args: Any, **kwargs: Any) -> View: - return View(self, *args, **kwargs) - - -class View(Imread, ABC): - def __init__(self, base: Imread, dtype: DTypeLike = None) -> None: - super().__init__(base) - self.dtype = dtype or base.dtype - self.slice = base.slice - self._shape = Shape(base.shape) - self.frame_decorator = base.frame_decorator - self.transform = base.transform - self.flags = dict(C_CONTIGUOUS=False, F_CONTIGUOUS=False, OWNDATA=False, WRITEABLE=False, - ALIGNED=False, WRITEBACKIFCOPY=False, UPDATEIFCOPY=False) - - def __getattr__(self, item: str) -> Any: - if not hasattr(self.base, item): - raise AttributeError(f'{self.__class__} object has no attribute {item}') - return self.base.__getattribute__(item) - def main() -> None: - parser = ArgumentParser(description='Display info and save as tif') - parser.add_argument('-v', '--version', action='version', version=__version__) - parser.add_argument('file', help='image_file', type=str, nargs='*') - parser.add_argument('-w', '--write', help='path to tif/movie out, {folder}, {name} and {ext} take this from file in', - type=str, default=None) - parser.add_argument('-o', '--extract_ome', help='extract ome to xml file', action='store_true') - parser.add_argument('-r', '--register', help='register channels', action='store_true') - parser.add_argument('-c', '--channel', help='channel', type=int, default=None) - parser.add_argument('-z', '--zslice', help='z-slice', type=int, default=None) - parser.add_argument('-t', '--time', help='time (frames) in python slicing notation', type=str, default=None) - parser.add_argument('-s', '--split', help='split channels', action='store_true') - parser.add_argument('-f', '--force', help='force overwrite', action='store_true') - parser.add_argument('-C', '--movie-colors', help='colors for channels in movie', type=str, nargs='*') - parser.add_argument('-B', '--movie-brightnesses', help='scale brightness of each channel', - type=float, nargs='*') - parser.add_argument('-S', '--movie-scale', help='upscale movie xy size, int', type=float) + parser = ArgumentParser(description="Display info and save as tif") + parser.add_argument("-v", "--version", action="version", version=__version__) + parser.add_argument("file", help="image_file", type=str, nargs="*") + parser.add_argument( + "-w", + "--write", + help="path to tif/movie out, {folder}, {name} and {ext} take this from file in", + type=str, + default=None, + ) + parser.add_argument( + "-o", "--extract_ome", help="extract ome to xml file", action="store_true" + ) + parser.add_argument( + "-r", "--register", help="register channels", action="store_true" + ) + parser.add_argument("-c", "--channel", help="channel", type=int, default=None) + parser.add_argument("-z", "--zslice", help="z-slice", type=int, default=None) + parser.add_argument( + "-t", + "--time", + help="time (frames) in python slicing notation", + type=str, + default=None, + ) + parser.add_argument("-s", "--split", help="split channels", action="store_true") + parser.add_argument("-f", "--force", help="force overwrite", action="store_true") + parser.add_argument( + "-C", "--movie-colors", help="colors for channels in movie", type=str, nargs="*" + ) + parser.add_argument( + "-B", + "--movie-brightnesses", + help="scale brightness of each channel", + type=float, + nargs="*", + ) + parser.add_argument( + "-S", "--movie-scale", help="upscale movie xy size, int", type=float + ) args = parser.parse_args() - for file in tqdm(args.file, desc='operating on files', disable=len(args.file) == 1): + for file in tqdm(args.file, desc="operating on files", disable=len(args.file) == 1): file = Path(file) with Imread(file) as im: # noqa if args.register: im = im.with_transform() # noqa if args.write: - write = Path(args.write.format(folder=str(file.parent), name=file.stem, ext=file.suffix)).absolute() # noqa + write = Path( + args.write.format( + folder=str(file.parent), name=file.stem, ext=file.suffix + ) + ).absolute() # noqa write.parent.mkdir(parents=True, exist_ok=True) if write.exists() and not args.force: - print(f'File {args.write} exists already, add the -f flag if you want to overwrite it.') - elif write.suffix in ('.mkv', '.mp4'): - im.save_as_movie(write, args.channel, args.zslice, args.time, args.movie_colors, - args.movie_brightnesses, args.movie_scale, bar=len(args.file) == 1) + print( + f"File {args.write} exists already, add the -f flag if you want to overwrite it." + ) + elif write.suffix in (".mkv", ".mp4"): + im.save_as_movie( + write, + args.channel, + args.zslice, + args.time, + args.movie_colors, + args.movie_brightnesses, + args.movie_scale, + bar=len(args.file) == 1, + ) else: - im.save_as_tiff(write, args.channel, args.zslice, args.time, args.split, bar=len(args.file) == 1) + im.save_as_tiff( + write, + args.channel, + args.zslice, + args.time, + args.split, + bar=len(args.file) == 1, + ) if args.extract_ome: - with open(im.path.with_suffix('.ome.xml'), 'w') as f: + with open(im.path.with_suffix(".ome.xml"), "w") as f: f.write(im.ome.to_xml()) if len(args.file) == 1: print(im.summary) diff --git a/py/ndbioimage/transforms.py b/py/ndbioimage/transforms.py index 4d73a9a..38a0cff 100644 --- a/py/ndbioimage/transforms.py +++ b/py/ndbioimage/transforms.py @@ -21,7 +21,7 @@ except ImportError: DataFrame, Series, concat = None, None, None -if hasattr(yaml, 'full_load'): +if hasattr(yaml, "full_load"): yamlload = yaml.full_load else: yamlload = yaml.load @@ -34,7 +34,7 @@ class Transforms(dict): @classmethod def from_file(cls, file, C=True, T=True): - with open(Path(file).with_suffix('.yml')) as f: + with open(Path(file).with_suffix(".yml")) as f: return cls.from_dict(yamlload(f), C, T) @classmethod @@ -42,7 +42,9 @@ class Transforms(dict): new = cls() for key, value in d.items(): if isinstance(key, str) and C: - new[key.replace(r'\:', ':').replace('\\\\', '\\')] = Transform.from_dict(value) + new[key.replace(r"\:", ":").replace("\\\\", "\\")] = ( + Transform.from_dict(value) + ) elif T: new[key] = Transform.from_dict(value) return new @@ -69,11 +71,19 @@ class Transforms(dict): return new def asdict(self): - return {key.replace('\\', '\\\\').replace(':', r'\:') if isinstance(key, str) else key: value.asdict() - for key, value in self.items()} + return { + key.replace("\\", "\\\\").replace(":", r"\:") + if isinstance(key, str) + else key: value.asdict() + for key, value in self.items() + } def __getitem__(self, item): - return np.prod([self[i] for i in item[::-1]]) if isinstance(item, tuple) else super().__getitem__(item) + return ( + np.prod([self[i] for i in item[::-1]]) + if isinstance(item, tuple) + else super().__getitem__(item) + ) def __missing__(self, key): return self.default @@ -88,7 +98,7 @@ class Transforms(dict): return hash(frozenset((*self.__dict__.items(), *self.items()))) def save(self, file): - with open(Path(file).with_suffix('.yml'), 'w') as f: + with open(Path(file).with_suffix(".yml"), "w") as f: yaml.safe_dump(self.asdict(), f, default_flow_style=None) def copy(self): @@ -109,8 +119,10 @@ class Transforms(dict): transform_channels = {key for key in self.keys() if isinstance(key, str)} if set(channel_names) - transform_channels: mapping = key_map(channel_names, transform_channels) - warnings.warn(f'The image file and the transform do not have the same channels,' - f' creating a mapping: {mapping}') + warnings.warn( + f"The image file and the transform do not have the same channels," + f" creating a mapping: {mapping}" + ) for key_im, key_t in mapping.items(): self[key_im] = self[key_t] @@ -124,37 +136,54 @@ class Transforms(dict): def coords_pandas(self, array, channel_names, columns=None): if isinstance(array, DataFrame): - return concat([self.coords_pandas(row, channel_names, columns) for _, row in array.iterrows()], axis=1).T + return concat( + [ + self.coords_pandas(row, channel_names, columns) + for _, row in array.iterrows() + ], + axis=1, + ).T elif isinstance(array, Series): key = [] - if 'C' in array: - key.append(channel_names[int(array['C'])]) - if 'T' in array: - key.append(int(array['T'])) + if "C" in array: + key.append(channel_names[int(array["C"])]) + if "T" in array: + key.append(int(array["T"])) return self[tuple(key)].coords(array, columns) else: - raise TypeError('Not a pandas DataFrame or Series.') + raise TypeError("Not a pandas DataFrame or Series.") def with_beads(self, cyllens, bead_files): - assert len(bead_files) > 0, 'At least one file is needed to calculate the registration.' - transforms = [self.calculate_channel_transforms(file, cyllens) for file in bead_files] + assert len(bead_files) > 0, ( + "At least one file is needed to calculate the registration." + ) + transforms = [ + self.calculate_channel_transforms(file, cyllens) for file in bead_files + ] for key in {key for transform in transforms for key in transform.keys()}: - new_transforms = [transform[key] for transform in transforms if key in transform] + new_transforms = [ + transform[key] for transform in transforms if key in transform + ] if len(new_transforms) == 1: self[key] = new_transforms[0] else: self[key] = Transform() - self[key].parameters = np.mean([t.parameters for t in new_transforms], 0) - self[key].dparameters = (np.std([t.parameters for t in new_transforms], 0) / - np.sqrt(len(new_transforms))).tolist() + self[key].parameters = np.mean( + [t.parameters for t in new_transforms], 0 + ) + self[key].dparameters = ( + np.std([t.parameters for t in new_transforms], 0) + / np.sqrt(len(new_transforms)) + ).tolist() return self @staticmethod def get_bead_files(path): from . import Imread + files = [] for file in path.iterdir(): - if file.name.lower().startswith('beads'): + if file.name.lower().startswith("beads"): try: with Imread(file): files.append(file) @@ -162,32 +191,36 @@ class Transforms(dict): pass files = sorted(files) if not files: - raise Exception('No bead file found!') + raise Exception("No bead file found!") checked_files = [] for file in files: try: if file.is_dir(): - file /= 'Pos0' + file /= "Pos0" with Imread(file): # check for errors opening the file checked_files.append(file) except (Exception,): continue if not checked_files: - raise Exception('No bead file found!') + raise Exception("No bead file found!") return checked_files @staticmethod def calculate_channel_transforms(bead_file, cyllens): - """ When no channel is not transformed by a cylindrical lens, assume that the image is scaled by a factor 1.162 - in the horizontal direction """ + """When no channel is not transformed by a cylindrical lens, assume that the image is scaled by a factor 1.162 + in the horizontal direction""" from . import Imread - with Imread(bead_file, axes='zcyx') as im: # noqa - max_ims = im.max('z') + with Imread(bead_file, axes="zcyx") as im: # noqa + max_ims = im.max("z") goodch = [c for c, max_im in enumerate(max_ims) if not im.is_noise(max_im)] if not goodch: goodch = list(range(len(max_ims))) - untransformed = [c for c in range(im.shape['c']) if cyllens[im.detector[c]].lower() == 'none'] + untransformed = [ + c + for c in range(im.shape["c"]) + if cyllens[im.detector[c]].lower() == "none" + ] good_and_untrans = sorted(set(goodch) & set(untransformed)) if good_and_untrans: @@ -200,54 +233,81 @@ class Transforms(dict): matrix[0, 0] = 0.86 transform.matrix = matrix transforms = Transforms() - for c in tqdm(goodch, desc='Calculating channel transforms'): # noqa + for c in tqdm(goodch, desc="Calculating channel transforms"): # noqa if c == masterch: transforms[im.channel_names[c]] = transform else: - transforms[im.channel_names[c]] = Transform.register(max_ims[masterch], max_ims[c]) * transform + transforms[im.channel_names[c]] = ( + Transform.register(max_ims[masterch], max_ims[c]) * transform + ) return transforms @staticmethod def save_channel_transform_tiff(bead_files, tiffile): from . import Imread + n_channels = 0 for file in bead_files: with Imread(file) as im: - n_channels = max(n_channels, im.shape['c']) + n_channels = max(n_channels, im.shape["c"]) with IJTiffFile(tiffile) as tif: for t, file in enumerate(bead_files): with Imread(file) as im: with Imread(file).with_transform() as jm: - for c in range(im.shape['c']): - tif.save(np.hstack((im(c=c, t=0).max('z'), jm(c=c, t=0).max('z'))), c, 0, t) + for c in range(im.shape["c"]): + tif.save( + np.hstack( + (im(c=c, t=0).max("z"), jm(c=c, t=0).max("z")) + ), + c, + 0, + t, + ) def with_drift(self, im): - """ Calculate shifts relative to the first frame - divide the sequence into groups, - compare each frame to the frame in the middle of the group and compare these middle frames to each other + """Calculate shifts relative to the first frame + divide the sequence into groups, + compare each frame to the frame in the middle of the group and compare these middle frames to each other """ - im = im.transpose('tzycx') - t_groups = [list(chunk) for chunk in Chunks(range(im.shape['t']), size=round(np.sqrt(im.shape['t'])))] + im = im.transpose("tzycx") + t_groups = [ + list(chunk) + for chunk in Chunks( + range(im.shape["t"]), size=round(np.sqrt(im.shape["t"])) + ) + ] t_keys = [int(np.round(np.mean(t_group))) for t_group in t_groups] - t_pairs = [(int(np.round(np.mean(t_group))), frame) for t_group in t_groups for frame in t_group] + t_pairs = [ + (int(np.round(np.mean(t_group))), frame) + for t_group in t_groups + for frame in t_group + ] t_pairs.extend(zip(t_keys, t_keys[1:])) - fmaxz_keys = {t_key: filters.gaussian(im[t_key].max('z'), 5) for t_key in t_keys} + fmaxz_keys = { + t_key: filters.gaussian(im[t_key].max("z"), 5) for t_key in t_keys + } def fun(t_key_t, im, fmaxz_keys): t_key, t = t_key_t if t_key == t: return 0, 0 else: - fmaxz = filters.gaussian(im[t].max('z'), 5) - return Transform.register(fmaxz_keys[t_key], fmaxz, 'translation').parameters[4:] + fmaxz = filters.gaussian(im[t].max("z"), 5) + return Transform.register( + fmaxz_keys[t_key], fmaxz, "translation" + ).parameters[4:] - shifts = np.array(pmap(fun, t_pairs, (im, fmaxz_keys), desc='Calculating image shifts.')) + shifts = np.array( + pmap(fun, t_pairs, (im, fmaxz_keys), desc="Calculating image shifts.") + ) shift_keys_cum = np.zeros(2) - for shift_keys, t_group in zip(np.vstack((-shifts[0], shifts[im.shape['t']:])), t_groups): + for shift_keys, t_group in zip( + np.vstack((-shifts[0], shifts[im.shape["t"] :])), t_groups + ): shift_keys_cum += shift_keys shifts[t_group] += shift_keys_cum - for i, shift in enumerate(shifts[:im.shape['t']]): + for i, shift in enumerate(shifts[: im.shape["t"]]): self[i] = Transform.from_shift(shift) return self @@ -257,9 +317,11 @@ class Transform: if sitk is None: self.transform = None else: - self.transform = sitk.ReadTransform(str(Path(__file__).parent / 'transform.txt')) - self.dparameters = [0., 0., 0., 0., 0., 0.] - self.shape = [512., 512.] + self.transform = sitk.ReadTransform( + str(Path(__file__).parent / "transform.txt") + ) + self.dparameters = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + self.shape = [512.0, 512.0] self.origin = [255.5, 255.5] self._last, self._inverse = None, None @@ -274,12 +336,14 @@ class Transform: @classmethod def register(cls, fix, mov, kind=None): - """ kind: 'affine', 'translation', 'rigid' """ + """kind: 'affine', 'translation', 'rigid'""" if sitk is None: - raise ImportError('SimpleElastix is not installed: ' - 'https://simpleelastix.readthedocs.io/GettingStarted.html') + raise ImportError( + "SimpleElastix is not installed: " + "https://simpleelastix.readthedocs.io/GettingStarted.html" + ) new = cls() - kind = kind or 'affine' + kind = kind or "affine" new.shape = fix.shape fix, mov = new.cast_image(fix), new.cast_image(mov) # TODO: implement RigidTransform @@ -290,16 +354,18 @@ class Transform: tfilter.SetParameterMap(sitk.GetDefaultParameterMap(kind)) tfilter.Execute() transform = tfilter.GetTransformParameterMap()[0] - if kind == 'affine': - new.parameters = [float(t) for t in transform['TransformParameters']] - new.shape = [float(t) for t in transform['Size']] - new.origin = [float(t) for t in transform['CenterOfRotationPoint']] - elif kind == 'translation': - new.parameters = [1.0, 0.0, 0.0, 1.0] + [float(t) for t in transform['TransformParameters']] - new.shape = [float(t) for t in transform['Size']] + if kind == "affine": + new.parameters = [float(t) for t in transform["TransformParameters"]] + new.shape = [float(t) for t in transform["Size"]] + new.origin = [float(t) for t in transform["CenterOfRotationPoint"]] + elif kind == "translation": + new.parameters = [1.0, 0.0, 0.0, 1.0] + [ + float(t) for t in transform["TransformParameters"] + ] + new.shape = [float(t) for t in transform["Size"]] new.origin = [(t - 1) / 2 for t in new.shape] else: - raise NotImplementedError(f'{kind} tranforms not implemented (yet)') + raise NotImplementedError(f"{kind} tranforms not implemented (yet)") new.dparameters = 6 * [np.nan] return new @@ -315,18 +381,35 @@ class Transform: @classmethod def from_file(cls, file): - with open(Path(file).with_suffix('.yml')) as f: + with open(Path(file).with_suffix(".yml")) as f: return cls.from_dict(yamlload(f)) @classmethod def from_dict(cls, d): new = cls() - new.origin = None if d['CenterOfRotationPoint'] is None else [float(i) for i in d['CenterOfRotationPoint']] - new.parameters = ((1., 0., 0., 1., 0., 0.) if d['TransformParameters'] is None else - [float(i) for i in d['TransformParameters']]) - new.dparameters = ([(0., 0., 0., 0., 0., 0.) if i is None else float(i) for i in d['dTransformParameters']] - if 'dTransformParameters' in d else 6 * [np.nan] and d['dTransformParameters'] is not None) - new.shape = None if d['Size'] is None else [None if i is None else float(i) for i in d['Size']] + new.origin = ( + None + if d["CenterOfRotationPoint"] is None + else [float(i) for i in d["CenterOfRotationPoint"]] + ) + new.parameters = ( + (1.0, 0.0, 0.0, 1.0, 0.0, 0.0) + if d["TransformParameters"] is None + else [float(i) for i in d["TransformParameters"]] + ) + new.dparameters = ( + [ + (0.0, 0.0, 0.0, 0.0, 0.0, 0.0) if i is None else float(i) + for i in d["dTransformParameters"] + ] + if "dTransformParameters" in d + else 6 * [np.nan] and d["dTransformParameters"] is not None + ) + new.shape = ( + None + if d["Size"] is None + else [None if i is None else float(i) for i in d["Size"]] + ) return new def __mul__(self, other): # TODO: take care of dmatrix @@ -359,9 +442,13 @@ class Transform: @property def matrix(self): - return np.array(((*self.parameters[:2], self.parameters[4]), - (*self.parameters[2:4], self.parameters[5]), - (0, 0, 1))) + return np.array( + ( + (*self.parameters[:2], self.parameters[4]), + (*self.parameters[2:4], self.parameters[5]), + (0, 0, 1), + ) + ) @matrix.setter def matrix(self, value): @@ -370,9 +457,13 @@ class Transform: @property def dmatrix(self): - return np.array(((*self.dparameters[:2], self.dparameters[4]), - (*self.dparameters[2:4], self.dparameters[5]), - (0, 0, 0))) + return np.array( + ( + (*self.dparameters[:2], self.dparameters[4]), + (*self.dparameters[2:4], self.dparameters[5]), + (0, 0, 0), + ) + ) @dmatrix.setter def dmatrix(self, value): @@ -384,7 +475,7 @@ class Transform: if self.transform is not None: return list(self.transform.GetParameters()) else: - return [1., 0., 0., 1., 0., 0.] + return [1.0, 0.0, 0.0, 1.0, 0.0, 0.0] @parameters.setter def parameters(self, value): @@ -420,43 +511,62 @@ class Transform: self.shape = shape[:2] def asdict(self): - return {'CenterOfRotationPoint': self.origin, 'Size': self.shape, 'TransformParameters': self.parameters, - 'dTransformParameters': np.nan_to_num(self.dparameters, nan=1e99).tolist()} + return { + "CenterOfRotationPoint": self.origin, + "Size": self.shape, + "TransformParameters": self.parameters, + "dTransformParameters": np.nan_to_num(self.dparameters, nan=1e99).tolist(), + } def frame(self, im, default=0): if self.is_unity(): return im else: if sitk is None: - raise ImportError('SimpleElastix is not installed: ' - 'https://simpleelastix.readthedocs.io/GettingStarted.html') + raise ImportError( + "SimpleElastix is not installed: " + "https://simpleelastix.readthedocs.io/GettingStarted.html" + ) dtype = im.dtype - im = im.astype('float') - intp = sitk.sitkBSpline if np.issubdtype(dtype, np.floating) else sitk.sitkNearestNeighbor - return self.cast_array(sitk.Resample(self.cast_image(im), self.transform, intp, default)).astype(dtype) + im = im.astype("float") + intp = ( + sitk.sitkBSpline + if np.issubdtype(dtype, np.floating) + else sitk.sitkNearestNeighbor + ) + return self.cast_array( + sitk.Resample(self.cast_image(im), self.transform, intp, default) + ).astype(dtype) def coords(self, array, columns=None): - """ Transform coordinates in 2 column numpy array, - or in pandas DataFrame or Series objects in columns ['x', 'y'] + """Transform coordinates in 2 column numpy array, + or in pandas DataFrame or Series objects in columns ['x', 'y'] """ if self.is_unity(): return array.copy() elif DataFrame is not None and isinstance(array, (DataFrame, Series)): - columns = columns or ['x', 'y'] + columns = columns or ["x", "y"] array = array.copy() if isinstance(array, DataFrame): array[columns] = self.coords(np.atleast_2d(array[columns].to_numpy())) elif isinstance(array, Series): - array[columns] = self.coords(np.atleast_2d(array[columns].to_numpy()))[0] + array[columns] = self.coords(np.atleast_2d(array[columns].to_numpy()))[ + 0 + ] return array else: # somehow we need to use the inverse here to get the same effect as when using self.frame - return np.array([self.inverse.transform.TransformPoint(i.tolist()) for i in np.asarray(array)]) + return np.array( + [ + self.inverse.transform.TransformPoint(i.tolist()) + for i in np.asarray(array) + ] + ) def save(self, file): - """ save the parameters of the transform calculated - with affine_registration to a yaml file + """save the parameters of the transform calculated + with affine_registration to a yaml file """ - if not file[-3:] == 'yml': - file += '.yml' - with open(file, 'w') as f: + if not file[-3:] == "yml": + file += ".yml" + with open(file, "w") as f: yaml.safe_dump(self.asdict(), f, default_flow_style=None) diff --git a/src/axes.rs b/src/axes.rs new file mode 100644 index 0000000..03e2575 --- /dev/null +++ b/src/axes.rs @@ -0,0 +1,218 @@ +use crate::stats::MinMax; +use anyhow::{anyhow, Error, Result}; +use ndarray::{Array, Dimension, Ix2, SliceInfo, SliceInfoElem}; +use serde::{Deserialize, Deserializer, Serialize, Serializer}; +use serde_with::{DeserializeAs, SerializeAs}; +use std::hash::{Hash, Hasher}; +use std::str::FromStr; + +/// a trait to find axis indices from any object +pub trait Ax { + /// C: 0, Z: 1, T: 2, Y: 3, X: 4 + fn n(&self) -> usize; + + /// the indices of axes in self.axes, which always has all of CZTYX + fn pos(&self, axes: &[Axis], slice: &[SliceInfoElem]) -> Result; + + /// the indices of axes in self.axes, which always has all of CZTYX, but skip axes with an operation + fn pos_op(&self, axes: &[Axis], slice: &[SliceInfoElem], op_axes: &[Axis]) -> Result; +} + +/// Enum for CZTYX axes or a new axis +#[derive(Clone, Copy, Debug, Eq, Ord, PartialOrd, Serialize, Deserialize)] +pub enum Axis { + C, + Z, + T, + Y, + X, + New, +} + +impl Hash for Axis { + fn hash(&self, state: &mut H) { + (*self as usize).hash(state); + } +} + +impl FromStr for Axis { + type Err = Error; + + fn from_str(s: &str) -> std::result::Result { + match s.to_uppercase().as_str() { + "C" => Ok(Axis::C), + "Z" => Ok(Axis::Z), + "T" => Ok(Axis::T), + "Y" => Ok(Axis::Y), + "X" => Ok(Axis::X), + "NEW" => Ok(Axis::New), + _ => Err(anyhow!("invalid axis: {}", s)), + } + } +} + +impl Ax for Axis { + fn n(&self) -> usize { + *self as usize + } + + fn pos(&self, axes: &[Axis], _slice: &[SliceInfoElem]) -> Result { + if let Some(pos) = axes.iter().position(|a| a == self) { + Ok(pos) + } else { + Err(Error::msg("Axis not found in axes")) + } + } + + fn pos_op(&self, axes: &[Axis], _slice: &[SliceInfoElem], _op_axes: &[Axis]) -> Result { + self.pos(axes, _slice) + } +} + +impl Ax for usize { + fn n(&self) -> usize { + *self + } + + fn pos(&self, _axes: &[Axis], slice: &[SliceInfoElem]) -> Result { + let idx: Vec<_> = slice + .iter() + .enumerate() + .filter_map(|(i, s)| if s.is_index() { None } else { Some(i) }) + .collect(); + Ok(idx[*self]) + } + + fn pos_op(&self, axes: &[Axis], slice: &[SliceInfoElem], op_axes: &[Axis]) -> Result { + let idx: Vec<_> = axes + .iter() + .zip(slice.iter()) + .enumerate() + .filter_map(|(i, (ax, s))| { + if s.is_index() | op_axes.contains(ax) { + None + } else { + Some(i) + } + }) + .collect(); + assert!(*self < idx.len(), "self: {}, idx: {:?}", self, idx); + Ok(idx[*self]) + } +} + +#[derive(Clone, Debug, Serialize, Deserialize)] +pub(crate) enum Operation { + Max, + Min, + Sum, + Mean, +} + +impl Operation { + pub(crate) fn operate( + &self, + array: Array, + axis: usize, + ) -> Result< as MinMax>::Output> + where + D: Dimension, + Array: MinMax, + { + match self { + Operation::Max => array.max(axis), + Operation::Min => array.min(axis), + Operation::Sum => array.sum(axis), + Operation::Mean => array.mean(axis), + } + } +} + +impl PartialEq for Axis { + fn eq(&self, other: &Self) -> bool { + (*self as u8) == (*other as u8) + } +} + +pub(crate) fn slice_info( + info: &[SliceInfoElem], +) -> Result> { + Ok(info.try_into()?) +} + +#[derive(Serialize, Deserialize)] +#[serde(remote = "SliceInfoElem")] +pub(crate) enum SliceInfoElemDef { + Slice { + start: isize, + end: Option, + step: isize, + }, + Index(isize), + NewAxis, +} + +impl SerializeAs for SliceInfoElemDef { + fn serialize_as( + source: &SliceInfoElem, + serializer: S, + ) -> std::result::Result + where + S: Serializer, + { + SliceInfoElemDef::serialize(source, serializer) + } +} + +impl<'de> DeserializeAs<'de, SliceInfoElem> for SliceInfoElemDef { + fn deserialize_as(deserializer: D) -> std::result::Result + where + D: Deserializer<'de>, + { + SliceInfoElemDef::deserialize(deserializer) + } +} + +#[derive(Clone, Debug)] +pub(crate) struct Slice { + start: isize, + end: isize, + step: isize, +} + +impl Slice { + pub(crate) fn new(start: isize, end: isize, step: isize) -> Self { + Self { start, end, step } + } + + pub(crate) fn empty() -> Self { + Self { + start: 0, + end: 0, + step: 1, + } + } +} + +impl Iterator for Slice { + type Item = isize; + + fn next(&mut self) -> Option { + if self.end - self.start >= self.step { + let r = self.start; + self.start += self.step; + Some(r) + } else { + None + } + } +} + +impl IntoIterator for &Slice { + type Item = isize; + type IntoIter = Slice; + + fn into_iter(self) -> Self::IntoIter { + self.clone() + } +} diff --git a/src/bioformats.rs b/src/bioformats.rs index 8ccb391..b83e096 100644 --- a/src/bioformats.rs +++ b/src/bioformats.rs @@ -96,6 +96,18 @@ macro_rules! method { }; } +fn transmute_vec(vec: Vec) -> Vec { + unsafe { + // Ensure the original vector is not dropped. + let mut v_clone = std::mem::ManuallyDrop::new(vec); + Vec::from_raw_parts( + v_clone.as_mut_ptr() as *mut U, + v_clone.len(), + v_clone.capacity(), + ) + } +} + /// Wrapper around bioformats java class loci.common.DebugTools pub struct DebugTools; @@ -125,8 +137,7 @@ impl ChannelSeparator { } pub(crate) fn open_bytes(&self, index: i32) -> Result> { - let bi8 = self.open_bi8(index)?; - Ok(unsafe { std::mem::transmute::, Vec>(bi8) }) + Ok(transmute_vec(self.open_bi8(index)?)) } method!(open_bi8, "openBytes", [index: i32|p] => Vec|c); @@ -149,8 +160,7 @@ impl ImageReader { } pub(crate) fn open_bytes(&self, index: i32) -> Result> { - let bi8 = self.open_bi8(index)?; - Ok(unsafe { std::mem::transmute::, Vec>(bi8) }) + Ok(transmute_vec(self.open_bi8(index)?)) } pub(crate) fn ome_xml(&self) -> Result { diff --git a/src/lib.rs b/src/lib.rs index 2863de5..ba289c0 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,333 +1,23 @@ mod bioformats; +mod axes; #[cfg(feature = "python")] mod py; - -use anyhow::{anyhow, Result}; -use bioformats::{DebugTools, ImageReader, MetadataTools}; -use ndarray::Array2; -use num::{FromPrimitive, Zero}; -use std::any::type_name; -use std::fmt::Debug; -use std::ops::Deref; -use std::path::{Path, PathBuf}; -use thread_local::ThreadLocal; - -/// Pixel types (u)int(8/16/32) or float(32/64) -#[derive(Clone, Debug)] -pub enum PixelType { - INT8 = 0, - UINT8 = 1, - INT16 = 2, - UINT16 = 3, - INT32 = 4, - UINT32 = 5, - FLOAT = 6, - DOUBLE = 7, -} - -impl TryFrom for PixelType { - type Error = anyhow::Error; - - fn try_from(value: i32) -> Result { - match value { - 0 => Ok(PixelType::INT8), - 1 => Ok(PixelType::UINT8), - 2 => Ok(PixelType::INT16), - 3 => Ok(PixelType::UINT16), - 4 => Ok(PixelType::INT32), - 5 => Ok(PixelType::UINT32), - 6 => Ok(PixelType::FLOAT), - 7 => Ok(PixelType::DOUBLE), - _ => Err(anyhow::anyhow!("Unknown pixel type {}", value)), - } - } -} - -/// Struct containing frame data in one of eight pixel types. Cast to `Array2` using try_into. -#[derive(Clone, Debug)] -pub enum Frame { - INT8(Array2), - UINT8(Array2), - INT16(Array2), - UINT16(Array2), - INT32(Array2), - UINT32(Array2), - FLOAT(Array2), - DOUBLE(Array2), -} - -macro_rules! impl_frame_cast { - ($t:tt, $s:ident) => { - impl From> for Frame { - fn from(value: Array2<$t>) -> Self { - Frame::$s(value) - } - } - }; -} - -impl_frame_cast!(i8, INT8); -impl_frame_cast!(u8, UINT8); -impl_frame_cast!(i16, INT16); -impl_frame_cast!(u16, UINT16); -impl_frame_cast!(i32, INT32); -impl_frame_cast!(u32, UINT32); -impl_frame_cast!(f32, FLOAT); -impl_frame_cast!(f64, DOUBLE); - -impl TryInto> for Frame -where - T: FromPrimitive + Zero + 'static, -{ - type Error = anyhow::Error; - - fn try_into(self) -> std::result::Result, Self::Error> { - let mut err = Ok(()); - let arr = match self { - Frame::INT8(v) => v.mapv_into_any(|x| { - T::from_i8(x).unwrap_or_else(|| { - err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); - T::zero() - }) - }), - Frame::UINT8(v) => v.mapv_into_any(|x| { - T::from_u8(x).unwrap_or_else(|| { - err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); - T::zero() - }) - }), - Frame::INT16(v) => v.mapv_into_any(|x| { - T::from_i16(x).unwrap_or_else(|| { - err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); - T::zero() - }) - }), - Frame::UINT16(v) => v.mapv_into_any(|x| { - T::from_u16(x).unwrap_or_else(|| { - err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); - T::zero() - }) - }), - Frame::INT32(v) => v.mapv_into_any(|x| { - T::from_i32(x).unwrap_or_else(|| { - err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); - T::zero() - }) - }), - Frame::UINT32(v) => v.mapv_into_any(|x| { - T::from_u32(x).unwrap_or_else(|| { - err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); - T::zero() - }) - }), - Frame::FLOAT(v) => v.mapv_into_any(|x| { - T::from_f32(x).unwrap_or_else(|| { - err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); - T::zero() - }) - }), - Frame::DOUBLE(v) => v.mapv_into_any(|x| { - T::from_f64(x).unwrap_or_else(|| { - err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); - T::zero() - }) - }), - }; - match err { - Err(err) => Err(err), - Ok(()) => Ok(arr), - } - } -} - -/// Reader interface to file. Use get_frame to get data. -pub struct Reader { - image_reader: ThreadLocal, - /// path to file - pub path: PathBuf, - /// which (if more than 1) of the series in the file to open - pub series: i32, - /// size x (horizontal) - pub size_x: usize, - /// size y (vertical) - pub size_y: usize, - /// size c (# channels) - pub size_c: usize, - /// size z (# slices) - pub size_z: usize, - /// size t (# time/frames) - pub size_t: usize, - /// pixel type ((u)int(8/16/32) or float(32/64)) - pub pixel_type: PixelType, - little_endian: bool, -} - -impl Deref for Reader { - type Target = ImageReader; - - fn deref(&self) -> &Self::Target { - self.image_reader.get_or(|| { - let reader = ImageReader::new().unwrap(); - let meta_data_tools = MetadataTools::new().unwrap(); - 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 - }) - } -} - -impl Clone for Reader { - fn clone(&self) -> Self { - Reader::new(&self.path, self.series).unwrap() - } -} - -impl Debug for Reader { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - f.debug_struct("Reader") - .field("path", &self.path) - .field("series", &self.series) - .field("size_x", &self.size_x) - .field("size_y", &self.size_y) - .field("size_c", &self.size_c) - .field("size_z", &self.size_z) - .field("size_t", &self.size_t) - .field("pixel_type", &self.pixel_type) - .field("little_endian", &self.little_endian) - .finish() - } -} - -impl Reader { - /// Create new reader for image file at path. - pub fn new(path: &Path, series: i32) -> Result { - DebugTools::set_root_level("ERROR")?; - let mut reader = Reader { - image_reader: ThreadLocal::new(), - path: PathBuf::from(path), - series, - size_x: 0, - size_y: 0, - size_c: 0, - size_z: 0, - size_t: 0, - pixel_type: PixelType::INT8, - little_endian: false, - }; - reader.size_x = reader.get_size_x()? as usize; - reader.size_y = reader.get_size_y()? as usize; - reader.size_c = reader.get_size_c()? as usize; - reader.size_z = reader.get_size_z()? as usize; - reader.size_t = reader.get_size_t()? as usize; - reader.pixel_type = PixelType::try_from(reader.get_pixel_type()?)?; - reader.little_endian = reader.is_little_endian()?; - Ok(reader) - } - - /// Get ome metadata as xml string - pub fn get_ome_xml(&self) -> Result { - self.ome_xml() - } - - fn deinterleave(&self, bytes: Vec, channel: usize) -> Result> { - let chunk_size = match self.pixel_type { - PixelType::INT8 => 1, - PixelType::UINT8 => 1, - PixelType::INT16 => 2, - PixelType::UINT16 => 2, - PixelType::INT32 => 4, - PixelType::UINT32 => 4, - PixelType::FLOAT => 4, - PixelType::DOUBLE => 8, - }; - Ok(bytes - .chunks(chunk_size) - .skip(channel) - .step_by(self.size_c) - .flat_map(|a| a.to_vec()) - .collect()) - } - - /// 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 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 { - let channel_separator = bioformats::ChannelSeparator::new(self)?; - let index = channel_separator.get_index(z as i32, c as i32, t as i32)?; - channel_separator.open_bytes(index)? - } else if self.is_indexed()? { - let index = self.get_index(z as i32, 0, t as i32)?; - self.open_bytes(index)? - // TODO: apply LUT - // let _bytes_lut = match self.pixel_type { - // PixelType::INT8 | PixelType::UINT8 => { - // let _lut = self.image_reader.get_8bit_lookup_table()?; - // } - // PixelType::INT16 | PixelType::UINT16 => { - // let _lut = self.image_reader.get_16bit_lookup_table()?; - // } - // _ => {} - // }; - } else { - let index = self.get_index(z as i32, c as i32, t as i32)?; - self.open_bytes(index)? - }; - self.bytes_to_frame(bytes) - } - - fn bytes_to_frame(&self, bytes: Vec) -> Result { - macro_rules! get_frame { - ($t:tt, <$n:expr) => { - Ok(Frame::from(Array2::from_shape_vec( - (self.size_y, self.size_x), - bytes - .chunks($n) - .map(|x| $t::from_le_bytes(x.try_into().unwrap())) - .collect(), - )?)) - }; - ($t:tt, >$n:expr) => { - Ok(Frame::from(Array2::from_shape_vec( - (self.size_y, self.size_x), - bytes - .chunks($n) - .map(|x| $t::from_be_bytes(x.try_into().unwrap())) - .collect(), - )?)) - }; - } - - match (&self.pixel_type, self.little_endian) { - (PixelType::INT8, true) => get_frame!(i8, <1), - (PixelType::UINT8, true) => get_frame!(u8, <1), - (PixelType::INT16, true) => get_frame!(i16, <2), - (PixelType::UINT16, true) => get_frame!(u16, <2), - (PixelType::INT32, true) => get_frame!(i32, <4), - (PixelType::UINT32, true) => get_frame!(u32, <4), - (PixelType::FLOAT, true) => get_frame!(f32, <4), - (PixelType::DOUBLE, true) => get_frame!(f64, <8), - (PixelType::INT8, false) => get_frame!(i8, >1), - (PixelType::UINT8, false) => get_frame!(u8, >1), - (PixelType::INT16, false) => get_frame!(i16, >2), - (PixelType::UINT16, false) => get_frame!(u16, >2), - (PixelType::INT32, false) => get_frame!(i32, >4), - (PixelType::UINT32, false) => get_frame!(u32, >4), - (PixelType::FLOAT, false) => get_frame!(f32, >4), - (PixelType::DOUBLE, false) => get_frame!(f64, >8), - } - } -} +mod reader; +mod stats; +mod view; #[cfg(test)] mod tests { - use super::*; + use crate::stats::MinMax; + use ndarray::{Array, Array4, Array5, NewAxis}; use rayon::prelude::*; + use crate::axes::Axis; + use crate::reader::{Frame, Reader}; + use anyhow::Result; + use ndarray::{s, Array2}; + fn open(file: &str) -> Result { let path = std::env::current_dir()? .join("tests") @@ -413,4 +103,160 @@ mod tests { println!("{}", xml); Ok(()) } + + #[test] + fn view() -> Result<()> { + let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif"; + let reader = open(file)?; + let view = reader.view(); + let a = view.slice(s![0, 5, 0, .., ..])?; + let b = reader.get_frame(0, 5, 0)?; + let c: Array2 = a.try_into()?; + let d: Array2 = b.try_into()?; + assert_eq!(c, d); + Ok(()) + } + + #[test] + fn view_shape() -> Result<()> { + let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif"; + let reader = open(file)?; + let view = reader.view(); + let a = view.slice(s![0, ..5, 0, .., 100..200])?; + let shape = a.shape(); + assert_eq!(shape, vec![5, 1024, 100]); + Ok(()) + } + + #[test] + fn view_new_axis() -> Result<()> { + let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif"; + let reader = open(file)?; + let view = reader.view(); + let a = Array5::::zeros((1, 9, 1, 1024, 1024)); + 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()); + Ok(()) + } + + #[test] + fn view_permute_axes() -> Result<()> { + let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif"; + let reader = open(file)?; + let view = reader.view(); + let s = view.shape(); + let mut a = Array5::::zeros((s[0], s[1], s[2], s[3], s[4])); + assert_eq!(view.shape(), a.shape()); + let b: Array5 = view.clone().try_into()?; + assert_eq!(b.shape(), a.shape()); + + let view = view.swap_axes(Axis::C, Axis::Z)?; + a.swap_axes(0, 1); + 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()?; + assert_eq!(b.shape(), a.shape()); + Ok(()) + } + + macro_rules! test_max { + ($($name:ident: $b:expr $(,)?)*) => { + $( + #[test] + fn $name() -> Result<()> { + let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif"; + let reader = open(file)?; + let view = reader.view(); + let array: Array5 = view.clone().try_into()?; + let view = view.max_proj($b)?; + let a: Array4 = view.clone().try_into()?; + let b = array.max($b)?; + assert_eq!(a.shape(), b.shape()); + assert_eq!(a, b); + Ok(()) + } + )* + }; + } + + test_max! { + max_c: 0 + max_z: 1 + max_t: 2 + max_y: 3 + max_x: 4 + } + + macro_rules! test_index { + ($($name:ident: $b:expr $(,)?)*) => { + $( + #[test] + fn $name() -> Result<()> { + let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif"; + let reader = open(file)?; + let view = reader.view(); + let v4: Array = view.slice($b)?.try_into()?; + let a5: Array5 = reader.view().try_into()?; + let a4 = a5.slice($b).to_owned(); + assert_eq!(a4, v4); + Ok(()) + } + )* + }; + } + + test_index! { + index_0: s![.., .., .., .., ..] + index_1: s![0, .., .., .., ..] + index_2: s![.., 0, .., .., ..] + index_3: s![.., .., 0, .., ..] + index_4: s![.., .., .., 0, ..] + index_5: s![.., .., .., .., 0] + index_6: s![0, 0, .., .., ..] + index_7: s![0, .., 0, .., ..] + index_8: s![0, .., .., 0, ..] + index_9: s![0, .., .., .., 0] + index_a: s![.., 0, 0, .., ..] + index_b: s![.., 0, .., 0, ..] + index_c: s![.., 0, .., .., 0] + index_d: s![.., .., 0, 0, ..] + index_e: s![.., .., 0, .., 0] + index_f: s![.., .., .., 0, 0] + index_g: s![0, 0, 0, .., ..] + index_h: s![0, 0, .., 0, ..] + index_i: s![0, 0, .., .., 0] + index_j: s![0, .., 0, 0, ..] + index_k: s![0, .., 0, .., 0] + index_l: s![0, .., .., 0, 0] + index_m: s![0, 0, 0, 0, ..] + index_n: s![0, 0, 0, .., 0] + index_o: s![0, 0, .., 0, 0] + index_p: s![0, .., 0, 0, 0] + index_q: s![.., 0, 0, 0, 0] + index_r: s![0, 0, 0, 0, 0] + } + + #[test] + fn dyn_view() -> Result<()> { + let file = "YTL1841B2-2-1_1hr_DMSO_galinduction_1/Pos0/img_000000000_mScarlet_GFP-mSc-filter_004.tif"; + let reader = open(file)?; + let a = reader.view().into_dyn(); + let b = a.max_proj(1)?; + let c = b.slice(s![0, 0, .., ..])?; + let d = c.as_array::()?; + assert_eq!(d.shape(), [1024, 1024]); + Ok(()) + } } diff --git a/src/py.rs b/src/py.rs index 16f64e2..cb64a6b 100644 --- a/src/py.rs +++ b/src/py.rs @@ -1,64 +1,600 @@ +use crate::axes::Axis; use crate::bioformats::download_bioformats; -use crate::{Frame, Reader}; -use numpy::ToPyArray; +use crate::reader::{PixelType, Reader}; +use crate::view::{Item, View}; +use anyhow::{anyhow, Result}; +use ndarray::{Ix0, IxDyn, SliceInfoElem}; +use numpy::IntoPyArray; 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; -#[pyclass(subclass)] -#[pyo3(name = "Reader")] -#[derive(Debug)] -struct PyReader { - reader: Reader, +#[pyclass(module = "ndbioimage.ndbioimage_rs")] +struct ViewConstructor; + +#[pymethods] +impl ViewConstructor { + #[new] + fn new() -> Self { + Self + } + + fn __getstate__(&self) -> (u8,) { + (0,) + } + + fn __setstate__(&self, _state: (u8,)) {} + + #[staticmethod] + fn __call__(state: String) -> PyResult { + if let Ok(new) = from_str(&state) { + Ok(new) + } else { + Err(anyhow!("cannot parse state").into()) + } + } +} + +#[pyclass(subclass, module = "ndbioimage.ndbioimage_rs")] +#[pyo3(name = "View")] +#[derive(Debug, Serialize, Deserialize)] +struct PyView { + view: View, + dtype: PixelType, } #[pymethods] -impl PyReader { +impl PyView { #[new] - fn new(path: &str, series: usize) -> PyResult { + #[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()? { - if let Ok(f) = file { - let p = f.path(); - if f.path().is_file() & (p.extension() == Some("tif".as_ref())) { - path = p; - break; - } + 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(PyReader { - reader: Reader::new(&path, series as i32)?, + 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 + fn close(&self) -> PyResult<()> { + Ok(()) + } + + fn copy(&self) -> PyView { + PyView { + view: self.view.clone(), + dtype: self.dtype.clone(), + } + } + + /// slice the view and return a new view or a single number + fn __getitem__<'py>( + &self, + py: Python<'py>, + n: Bound<'py, PyAny>, + ) -> PyResult> { + let slice: Vec<_> = if n.is_instance_of::() { + n.downcast_into::()?.into_iter().collect() + } else if n.is_instance_of::() { + n.downcast_into::()?.into_iter().collect() + } else { + vec![n] + }; + let mut new_slice = Vec::new(); + let mut ellipsis = None; + let shape = self.view.shape(); + for (i, (s, t)) in slice.iter().zip(shape.iter()).enumerate() { + if s.is_instance_of::() { + new_slice.push(SliceInfoElem::Index( + s.downcast::()?.extract::()?, + )); + } else if s.is_instance_of::() { + let u = s.downcast::()?.indices(*t as isize)?; + new_slice.push(SliceInfoElem::Slice { + start: u.start, + end: Some(u.stop), + step: u.step, + }); + } else if s.is_instance_of::() { + if ellipsis.is_some() { + return Err(anyhow!("cannot have more than one ellipsis").into()); + } + let _ = ellipsis.insert(i); + } else { + return Err(anyhow!("cannot convert {:?} to slice", s).into()); + } + } + if new_slice.len() > shape.len() { + return Err(anyhow!( + "got more indices ({}) than dimensions ({})", + new_slice.len(), + shape.len() + ) + .into()); + } + while new_slice.len() < shape.len() { + if let Some(i) = ellipsis { + new_slice.insert( + i, + SliceInfoElem::Slice { + start: 0, + end: None, + step: 1, + }, + ) + } else { + new_slice.push(SliceInfoElem::Slice { + start: 0, + end: None, + step: 1, + }) + } + } + let view = self.view.slice(new_slice.as_slice())?; + 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::F32 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + PixelType::F64 => 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::F128 => view + .into_dimensionality::()? + .item::()? + .into_pyobject(py)? + .into_any(), + }) + } else { + PyView { + view, + dtype: self.dtype.clone(), + } + .into_bound_py_any(py) + } + } + + fn __reduce__(&self) -> PyResult<(ViewConstructor, (String,))> { + if let Ok(s) = to_string(self) { + Ok((ViewConstructor, (s,))) + } else { + Err(anyhow!("cannot get state").into()) + } + } + + /// retrieve a single frame at czt, sliced accordingly fn get_frame<'py>( &self, py: Python<'py>, - c: usize, - z: usize, - t: usize, + c: isize, + z: isize, + t: isize, ) -> PyResult> { - Ok(match self.reader.get_frame(c, z, t)? { - Frame::INT8(arr) => arr.to_pyarray(py).into_any(), - Frame::UINT8(arr) => arr.to_pyarray(py).into_any(), - Frame::INT16(arr) => arr.to_pyarray(py).into_any(), - Frame::UINT16(arr) => arr.to_pyarray(py).into_any(), - Frame::INT32(arr) => arr.to_pyarray(py).into_any(), - Frame::UINT32(arr) => arr.to_pyarray(py).into_any(), - Frame::FLOAT(arr) => arr.to_pyarray(py).into_any(), - Frame::DOUBLE(arr) => arr.to_pyarray(py).into_any(), + Ok(match self.dtype { + PixelType::I8 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), + PixelType::U8 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), + PixelType::I16 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), + PixelType::U16 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), + PixelType::I32 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), + PixelType::U32 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), + PixelType::F32 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), + PixelType::F64 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), + PixelType::I64 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), + PixelType::U64 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), + PixelType::I128 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), + PixelType::U128 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), + PixelType::F128 => self + .view + .get_frame::(c, z, t)? + .into_pyarray(py) + .into_any(), }) } + /// retrieve the ome metadata as an xml string fn get_ome_xml(&self) -> PyResult { - Ok(self.reader.get_ome_xml()?) + Ok(self.view.get_ome_xml()?) } - fn close(&mut self) -> PyResult<()> { - self.reader.close()?; + /// the file path + #[getter] + fn path(&self) -> PyResult { + Ok(self.view.path.display().to_string()) + } + + /// the series in the file + #[getter] + 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() + } + + /// the shape of the view + #[getter] + fn shape(&self) -> Vec { + self.view.shape() + } + + #[getter] + fn slice(&self) -> PyResult> { + Ok(self + .view + .get_slice() + .iter() + .map(|s| format!("{:#?}", s)) + .collect()) + } + + /// the number of pixels in the view + #[getter] + fn size(&self) -> usize { + self.view.size() + } + + /// the number of dimensions in the view + #[getter] + fn ndim(&self) -> usize { + self.view.ndim() + } + + /// find the position of an axis + #[pyo3(text_signature = "axis: str | int")] + fn get_ax(&self, axis: Bound<'_, PyAny>) -> PyResult { + if axis.is_instance_of::() { + let axis = axis + .downcast_into::()? + .extract::()? + .parse::()?; + Ok(self + .view + .axes() + .iter() + .position(|a| *a == axis) + .ok_or_else(|| anyhow!("cannot find axis {:?}", axis))?) + } else if axis.is_instance_of::() { + Ok(axis.downcast_into::()?.extract::()?) + } else { + Err(anyhow!("cannot convert to axis").into()) + } + } + + /// swap two axes + #[pyo3(text_signature = "ax0: str | int, ax1: str | int")] + fn swap_axes(&self, ax0: Bound<'_, PyAny>, ax1: Bound<'_, PyAny>) -> PyResult { + let ax0 = self.get_ax(ax0)?; + let ax1 = self.get_ax(ax1)?; + let view = self.view.swap_axes(ax0, ax1)?; + Ok(PyView { + view, + dtype: self.dtype.clone(), + }) + } + + /// permute the order of the axes + #[pyo3(signature = (axes = None), text_signature = "axes: list[str | int] = None")] + fn transpose(&self, axes: Option>>) -> PyResult { + let view = if let Some(axes) = axes { + let ax = axes + .into_iter() + .map(|a| self.get_ax(a)) + .collect::, _>>()?; + self.view.permute_axes(&ax)? + } else { + self.view.transpose()? + }; + Ok(PyView { + view, + dtype: self.dtype.clone(), + }) + } + + /// collect data into a numpy array + fn as_array<'py>(&self, py: Python<'py>) -> PyResult> { + Ok(match self.dtype { + PixelType::I8 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + PixelType::U8 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + PixelType::I16 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + PixelType::U16 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + PixelType::I32 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + PixelType::U32 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + PixelType::F32 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + PixelType::F64 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + PixelType::I64 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + PixelType::U64 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + PixelType::I128 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + PixelType::U128 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + PixelType::F128 => self.view.as_array_dyn::()?.into_pyarray(py).into_any(), + }) + } + + /// 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 { + PixelType::I8 => "int8", + PixelType::U8 => "uint8", + PixelType::I16 => "int16", + PixelType::U16 => "uint16", + PixelType::I32 => "int32", + PixelType::U32 => "uint32", + PixelType::F32 => "float32", + PixelType::F64 => "float64", + PixelType::I64 => "int64", + PixelType::U64 => "uint64", + PixelType::I128 => "int128", + PixelType::U128 => "uint128", + PixelType::F128 => "float128", + }) + } + + #[setter] + fn set_dtype(&mut self, dtype: String) -> PyResult<()> { + self.dtype = dtype.parse()?; Ok(()) } + + /// get the maximum overall or along a given axis + #[pyo3(signature = (axis = None), text_signature = "axis: str | int")] + fn max<'py>( + &self, + py: Python<'py>, + axis: Option>, + ) -> PyResult> { + if let Some(axis) = axis { + PyView { + dtype: self.dtype.clone(), + view: self.view.max_proj(self.get_ax(axis)?)?, + } + .into_bound_py_any(py) + } else { + Ok(match self.dtype { + PixelType::I8 => self.view.max::()?.into_pyobject(py)?.into_any(), + PixelType::U8 => self.view.max::()?.into_pyobject(py)?.into_any(), + PixelType::I16 => self.view.max::()?.into_pyobject(py)?.into_any(), + PixelType::U16 => self.view.max::()?.into_pyobject(py)?.into_any(), + PixelType::I32 => self.view.max::()?.into_pyobject(py)?.into_any(), + PixelType::U32 => self.view.max::()?.into_pyobject(py)?.into_any(), + PixelType::F32 => self.view.max::()?.into_pyobject(py)?.into_any(), + PixelType::F64 => self.view.max::()?.into_pyobject(py)?.into_any(), + PixelType::I64 => self.view.max::()?.into_pyobject(py)?.into_any(), + PixelType::U64 => self.view.max::()?.into_pyobject(py)?.into_any(), + PixelType::I128 => self.view.max::()?.into_pyobject(py)?.into_any(), + PixelType::U128 => self.view.max::()?.into_pyobject(py)?.into_any(), + PixelType::F128 => self.view.max::()?.into_pyobject(py)?.into_any(), + }) + } + } + + /// get the minimum overall or along a given axis + #[pyo3(signature = (axis = None), text_signature = "axis: str | int")] + fn min<'py>( + &self, + py: Python<'py>, + axis: Option>, + ) -> PyResult> { + if let Some(axis) = axis { + PyView { + dtype: self.dtype.clone(), + view: self.view.min_proj(self.get_ax(axis)?)?, + } + .into_bound_py_any(py) + } else { + Ok(match self.dtype { + PixelType::I8 => self.view.min::()?.into_pyobject(py)?.into_any(), + PixelType::U8 => self.view.min::()?.into_pyobject(py)?.into_any(), + PixelType::I16 => self.view.min::()?.into_pyobject(py)?.into_any(), + PixelType::U16 => self.view.min::()?.into_pyobject(py)?.into_any(), + PixelType::I32 => self.view.min::()?.into_pyobject(py)?.into_any(), + PixelType::U32 => self.view.min::()?.into_pyobject(py)?.into_any(), + PixelType::F32 => self.view.min::()?.into_pyobject(py)?.into_any(), + PixelType::F64 => self.view.min::()?.into_pyobject(py)?.into_any(), + PixelType::I64 => self.view.min::()?.into_pyobject(py)?.into_any(), + PixelType::U64 => self.view.min::()?.into_pyobject(py)?.into_any(), + PixelType::I128 => self.view.min::()?.into_pyobject(py)?.into_any(), + PixelType::U128 => self.view.min::()?.into_pyobject(py)?.into_any(), + PixelType::F128 => self.view.min::()?.into_pyobject(py)?.into_any(), + }) + } + } + + /// get the mean overall or along a given axis + #[pyo3(signature = (axis = None), text_signature = "axis: str | int")] + fn mean<'py>( + &self, + py: Python<'py>, + axis: Option>, + ) -> PyResult> { + if let Some(axis) = axis { + let dtype = if let PixelType::F32 = self.dtype { + PixelType::F32 + } else { + PixelType::F64 + }; + PyView { + dtype, + view: self.view.mean_proj(self.get_ax(axis)?)?, + } + .into_bound_py_any(py) + } else { + Ok(match self.dtype { + PixelType::F32 => self.view.mean::()?.into_pyobject(py)?.into_any(), + _ => self.view.mean::()?.into_pyobject(py)?.into_any(), + }) + } + } + + /// get the sum overall or along a given axis + #[pyo3(signature = (axis = None), text_signature = "axis: str | int")] + fn sum<'py>( + &self, + py: Python<'py>, + axis: Option>, + ) -> PyResult> { + let dtype = match self.dtype { + PixelType::I8 => PixelType::I16, + PixelType::U8 => PixelType::U16, + PixelType::I16 => PixelType::I32, + PixelType::U16 => PixelType::U32, + PixelType::I32 => PixelType::I64, + PixelType::U32 => PixelType::U64, + PixelType::F32 => PixelType::F32, + PixelType::F64 => PixelType::F64, + PixelType::I64 => PixelType::I128, + PixelType::U64 => PixelType::U128, + PixelType::I128 => PixelType::I128, + PixelType::U128 => PixelType::U128, + PixelType::F128 => PixelType::F128, + }; + if let Some(axis) = axis { + PyView { + dtype, + view: self.view.sum_proj(self.get_ax(axis)?)?, + } + .into_bound_py_any(py) + } else { + Ok(match self.dtype { + PixelType::F32 => self.view.sum::()?.into_pyobject(py)?.into_any(), + PixelType::F64 => self.view.sum::()?.into_pyobject(py)?.into_any(), + PixelType::I64 => self.view.sum::()?.into_pyobject(py)?.into_any(), + PixelType::U64 => self.view.sum::()?.into_pyobject(py)?.into_any(), + PixelType::I128 => self.view.sum::()?.into_pyobject(py)?.into_any(), + PixelType::U128 => self.view.sum::()?.into_pyobject(py)?.into_any(), + PixelType::F128 => self.view.sum::()?.into_pyobject(py)?.into_any(), + _ => self.view.sum::()?.into_pyobject(py)?.into_any(), + }) + } + } } pub(crate) fn ndbioimage_file() -> anyhow::Result { @@ -75,7 +611,8 @@ pub(crate) fn ndbioimage_file() -> anyhow::Result { #[pymodule] #[pyo3(name = "ndbioimage_rs")] fn ndbioimage_rs(m: &Bound<'_, PyModule>) -> PyResult<()> { - m.add_class::()?; + m.add_class::()?; + m.add_class::()?; #[pyfn(m)] #[pyo3(name = "download_bioformats")] diff --git a/src/reader.rs b/src/reader.rs new file mode 100644 index 0000000..c7a6e2b --- /dev/null +++ b/src/reader.rs @@ -0,0 +1,442 @@ +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 num::{FromPrimitive, Zero}; +use serde::{Deserialize, Serialize}; +use std::any::type_name; +use std::fmt::Debug; +use std::ops::Deref; +use std::path::{Path, PathBuf}; +use std::str::FromStr; +use std::sync::Arc; +use thread_local::ThreadLocal; + +/// 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)] +pub enum PixelType { + I8, + U8, + I16, + U16, + I32, + U32, + F32, + F64, + I64, + U64, + I128, + U128, + F128, +} + +impl TryFrom for PixelType { + type Error = Error; + + fn try_from(value: i32) -> Result { + match value { + 0 => Ok(PixelType::I8), + 1 => Ok(PixelType::U8), + 2 => Ok(PixelType::I16), + 3 => Ok(PixelType::U16), + 4 => Ok(PixelType::I32), + 5 => Ok(PixelType::U32), + 6 => Ok(PixelType::F32), + 7 => Ok(PixelType::F64), + 8 => Ok(PixelType::I64), + 9 => Ok(PixelType::U64), + 10 => Ok(PixelType::I128), + 11 => Ok(PixelType::U128), + 12 => Ok(PixelType::F128), + _ => Err(anyhow::anyhow!("Unknown pixel type {}", value)), + } + } +} + +impl FromStr for PixelType { + type Err = Error; + + fn from_str(s: &str) -> Result { + match s.to_lowercase().as_str() { + "int8" | "i8" => Ok(PixelType::I8), + "uint8" | "u8" => Ok(PixelType::U8), + "int16" | "i16" => Ok(PixelType::I16), + "uint16" | "u16" => Ok(PixelType::U16), + "int32" | "i32" => Ok(PixelType::I32), + "uint32" | "u32" => Ok(PixelType::U32), + "float" | "f32" | "float32" => Ok(PixelType::F32), + "double" | "f64" | "float64" => Ok(PixelType::F64), + "int64" | "i64" => Ok(PixelType::I64), + "uint64" | "u64" => Ok(PixelType::U64), + "int128" | "i128" => Ok(PixelType::I128), + "uint128" | "u128" => Ok(PixelType::U128), + "extended" | "f128" => Ok(PixelType::F128), + _ => Err(anyhow::anyhow!("Unknown pixel type {}", s)), + } + } +} + +/// Struct containing frame data in one of eight pixel types. Cast to `Array2` using try_into. +#[allow(clippy::upper_case_acronyms)] +#[derive(Clone, Debug)] +pub enum Frame { + I8(Array2), + U8(Array2), + I16(Array2), + U16(Array2), + I32(Array2), + U32(Array2), + F32(Array2), + F64(Array2), + I64(Array2), + U64(Array2), + I128(Array2), + U128(Array2), + F128(Array2), // f128 is nightly +} + +macro_rules! impl_frame_cast { + ($($t:tt: $s:ident $(,)?)*) => { + $( + impl From> for Frame { + fn from(value: Array2<$t>) -> Self { + Frame::$s(value) + } + } + )* + }; +} + +impl_frame_cast! { + u8: U8 + i8: I8 + i16: I16 + u16: U16 + i32: I32 + u32: U32 + f32: F32 + f64: F64 + i64: I64 + u64: U64 + i128: I128 + u128: U128 +} + +#[cfg(target_pointer_width = "32")] +impl_frame_cast! { + usize: UINT32 + isize: INT32 +} + +impl TryInto> for Frame +where + T: FromPrimitive + Zero + 'static, +{ + type Error = Error; + + fn try_into(self) -> Result, Self::Error> { + let mut err = Ok(()); + let arr = match self { + Frame::I8(v) => v.mapv_into_any(|x| { + T::from_i8(x).unwrap_or_else(|| { + err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); + T::zero() + }) + }), + Frame::U8(v) => v.mapv_into_any(|x| { + T::from_u8(x).unwrap_or_else(|| { + err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); + T::zero() + }) + }), + Frame::I16(v) => v.mapv_into_any(|x| { + T::from_i16(x).unwrap_or_else(|| { + err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); + T::zero() + }) + }), + Frame::U16(v) => v.mapv_into_any(|x| { + T::from_u16(x).unwrap_or_else(|| { + err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); + T::zero() + }) + }), + Frame::I32(v) => v.mapv_into_any(|x| { + T::from_i32(x).unwrap_or_else(|| { + err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); + T::zero() + }) + }), + Frame::U32(v) => v.mapv_into_any(|x| { + T::from_u32(x).unwrap_or_else(|| { + err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); + T::zero() + }) + }), + Frame::F32(v) => v.mapv_into_any(|x| { + T::from_f32(x).unwrap_or_else(|| { + err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); + T::zero() + }) + }), + Frame::F64(v) | Frame::F128(v) => v.mapv_into_any(|x| { + T::from_f64(x).unwrap_or_else(|| { + err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); + T::zero() + }) + }), + Frame::I64(v) => v.mapv_into_any(|x| { + T::from_i64(x).unwrap_or_else(|| { + err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); + T::zero() + }) + }), + Frame::U64(v) => v.mapv_into_any(|x| { + T::from_u64(x).unwrap_or_else(|| { + err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); + T::zero() + }) + }), + Frame::I128(v) => v.mapv_into_any(|x| { + T::from_i128(x).unwrap_or_else(|| { + err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); + T::zero() + }) + }), + Frame::U128(v) => v.mapv_into_any(|x| { + T::from_u128(x).unwrap_or_else(|| { + err = Err(anyhow!("cannot convert {} into {}", x, type_name::())); + T::zero() + }) + }), + }; + match err { + Err(err) => Err(err), + Ok(()) => Ok(arr), + } + } +} + +/// Reader interface to file. Use get_frame to get data. +#[derive(Serialize, Deserialize)] +pub struct Reader { + #[serde(skip)] + image_reader: ThreadLocal, + /// path to file + pub path: PathBuf, + /// which (if more than 1) of the series in the file to open + pub series: i32, + /// size x (horizontal) + pub size_x: usize, + /// size y (vertical) + pub size_y: usize, + /// size c (# channels) + pub size_c: usize, + /// size z (# slices) + pub size_z: usize, + /// size t (# time/frames) + pub size_t: usize, + /// pixel type ((u)int(8/16/32) or float(32/64)) + pub pixel_type: PixelType, + little_endian: bool, +} + +impl Deref for Reader { + type Target = ImageReader; + + fn deref(&self) -> &Self::Target { + self.image_reader.get_or(|| { + let reader = ImageReader::new().unwrap(); + let meta_data_tools = MetadataTools::new().unwrap(); + 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 + }) + } +} + +impl Clone for Reader { + fn clone(&self) -> Self { + Reader::new(&self.path, self.series).unwrap() + } +} + +impl Debug for Reader { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("Reader") + .field("path", &self.path) + .field("series", &self.series) + .field("size_x", &self.size_x) + .field("size_y", &self.size_y) + .field("size_c", &self.size_c) + .field("size_z", &self.size_z) + .field("size_t", &self.size_t) + .field("pixel_type", &self.pixel_type) + .field("little_endian", &self.little_endian) + .finish() + } +} + +impl Reader { + /// Create a new reader for the image file at a path, and open series #. + pub fn new(path: &Path, series: i32) -> Result { + DebugTools::set_root_level("ERROR")?; + let mut reader = Reader { + image_reader: ThreadLocal::default(), + path: PathBuf::from(path), + series, + size_x: 0, + size_y: 0, + size_c: 0, + size_z: 0, + size_t: 0, + pixel_type: PixelType::I8, + little_endian: false, + }; + reader.size_x = reader.get_size_x()? as usize; + reader.size_y = reader.get_size_y()? as usize; + reader.size_c = reader.get_size_c()? as usize; + reader.size_z = reader.get_size_z()? as usize; + reader.size_t = reader.get_size_t()? as usize; + reader.pixel_type = PixelType::try_from(reader.get_pixel_type()?)?; + reader.little_endian = reader.is_little_endian()?; + Ok(reader) + } + + /// Get ome metadata as xml string + pub fn get_ome_xml(&self) -> Result { + self.ome_xml() + } + + fn deinterleave(&self, bytes: Vec, channel: usize) -> Result> { + let chunk_size = match self.pixel_type { + PixelType::I8 => 1, + PixelType::U8 => 1, + PixelType::I16 => 2, + PixelType::U16 => 2, + PixelType::I32 => 4, + PixelType::U32 => 4, + PixelType::F32 => 4, + PixelType::F64 => 8, + PixelType::I64 => 8, + PixelType::U64 => 8, + PixelType::I128 => 16, + PixelType::U128 => 16, + PixelType::F128 => 8, + }; + Ok(bytes + .chunks(chunk_size) + .skip(channel) + .step_by(self.size_c) + .flat_map(|a| a.to_vec()) + .collect()) + } + + /// 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 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 { + let channel_separator = bioformats::ChannelSeparator::new(self)?; + let index = channel_separator.get_index(z as i32, c as i32, t as i32)?; + channel_separator.open_bytes(index)? + } else if self.is_indexed()? { + let index = self.get_index(z as i32, 0, t as i32)?; + self.open_bytes(index)? + // TODO: apply LUT + // let _bytes_lut = match self.pixel_type { + // PixelType::INT8 | PixelType::UINT8 => { + // let _lut = self.image_reader.get_8bit_lookup_table()?; + // } + // PixelType::INT16 | PixelType::UINT16 => { + // let _lut = self.image_reader.get_16bit_lookup_table()?; + // } + // _ => {} + // }; + } else { + let index = self.get_index(z as i32, c as i32, t as i32)?; + self.open_bytes(index)? + }; + self.bytes_to_frame(bytes) + } + + fn bytes_to_frame(&self, bytes: Vec) -> Result { + macro_rules! get_frame { + ($t:tt, <$n:expr) => { + Ok(Frame::from(Array2::from_shape_vec( + (self.size_y, self.size_x), + bytes + .chunks($n) + .map(|x| $t::from_le_bytes(x.try_into().unwrap())) + .collect(), + )?)) + }; + ($t:tt, >$n:expr) => { + Ok(Frame::from(Array2::from_shape_vec( + (self.size_y, self.size_x), + bytes + .chunks($n) + .map(|x| $t::from_be_bytes(x.try_into().unwrap())) + .collect(), + )?)) + }; + } + + match (&self.pixel_type, self.little_endian) { + (PixelType::I8, true) => get_frame!(i8, <1), + (PixelType::U8, true) => get_frame!(u8, <1), + (PixelType::I16, true) => get_frame!(i16, <2), + (PixelType::U16, true) => get_frame!(u16, <2), + (PixelType::I32, true) => get_frame!(i32, <4), + (PixelType::U32, true) => get_frame!(u32, <4), + (PixelType::F32, true) => get_frame!(f32, <4), + (PixelType::F64, true) => get_frame!(f64, <8), + (PixelType::I64, true) => get_frame!(i64, <8), + (PixelType::U64, true) => get_frame!(u64, <8), + (PixelType::I128, true) => get_frame!(i128, <16), + (PixelType::U128, true) => get_frame!(u128, <16), + (PixelType::F128, true) => get_frame!(f64, <8), + (PixelType::I8, false) => get_frame!(i8, >1), + (PixelType::U8, false) => get_frame!(u8, >1), + (PixelType::I16, false) => get_frame!(i16, >2), + (PixelType::U16, false) => get_frame!(u16, >2), + (PixelType::I32, false) => get_frame!(i32, >4), + (PixelType::U32, false) => get_frame!(u32, >4), + (PixelType::F32, false) => get_frame!(f32, >4), + (PixelType::F64, false) => get_frame!(f64, >8), + (PixelType::I64, false) => get_frame!(i64, >8), + (PixelType::U64, false) => get_frame!(u64, >8), + (PixelType::I128, false) => get_frame!(i128, >16), + (PixelType::U128, false) => get_frame!(u128, >16), + (PixelType::F128, false) => get_frame!(f64, >8), + } + } + + /// get a slicable view on the image file + pub fn view(&self) -> View { + let slice = s![ + 0..self.size_c, + 0..self.size_z, + 0..self.size_t, + 0..self.size_y, + 0..self.size_x + ]; + View::new( + Arc::new(self.clone()), + slice.as_ref().to_vec(), + vec![Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X], + ) + } +} + +impl Drop for Reader { + fn drop(&mut self) { + let _ = self.close(); + } +} diff --git a/src/stats.rs b/src/stats.rs new file mode 100644 index 0000000..24d128f --- /dev/null +++ b/src/stats.rs @@ -0,0 +1,201 @@ +use anyhow::{anyhow, Result}; +use ndarray::{Array, ArrayD, ArrayView, Axis, Dimension, RemoveAxis}; + +/// a trait to define the min, max, sum and mean operations along an axis +pub trait MinMax { + type Output; + + fn max(self, axis: usize) -> Result; + fn min(self, axis: usize) -> Result; + fn sum(self, axis: usize) -> Result; + fn mean(self, axis: usize) -> Result; +} + +macro_rules! impl_frame_stats_float_view { + ($($t:tt),+ $(,)?) => { + $( + impl MinMax for ArrayView<'_, $t, D> + where + D: Dimension + RemoveAxis, + { + type Output = Array<$t, D::Smaller>; + + fn max(self, axis: usize) -> Result { + let a: Vec<_> = self + .lanes(Axis(axis)) + .into_iter() + .map(|x| { + x.iter() + .fold($t::NEG_INFINITY, |prev, curr| prev.max(*curr)) + }) + .collect(); + let mut shape = self.shape().to_vec(); + shape.remove(axis); + Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?) + } + + fn min(self, axis: usize) -> Result { + let a: Vec<_> = self + .lanes(Axis(axis)) + .into_iter() + .map(|x| { + x.iter() + .fold($t::INFINITY, |prev, curr| prev.min(*curr)) + }) + .collect(); + let mut shape = self.shape().to_vec(); + shape.remove(axis); + Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?) + } + + fn sum(self, axis: usize) -> Result { + Ok(self.sum_axis(Axis(axis))) + } + + fn mean(self, axis: usize) -> Result { + self.mean_axis(Axis(axis)).ok_or_else(|| anyhow!("no mean")) + } + } + )* + }; +} + +macro_rules! impl_frame_stats_int_view { + ($($t:tt),+ $(,)?) => { + $( + impl MinMax for ArrayView<'_, $t, D> + where + D: Dimension + RemoveAxis, + { + type Output = Array<$t, D::Smaller>; + + fn max(self, axis: usize) -> Result { + let a: Vec<_> = self + .lanes(Axis(axis)) + .into_iter() + .map(|x| *x.iter().max().unwrap()) + .collect(); + let mut shape = self.shape().to_vec(); + shape.remove(axis); + Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?) + } + + fn min(self, axis: usize) -> Result { + let a: Vec<_> = self + .lanes(Axis(axis)) + .into_iter() + .map(|x| *x.iter().min().unwrap()) + .collect(); + let mut shape = self.shape().to_vec(); + shape.remove(axis); + Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?) + } + + fn sum(self, axis: usize) -> Result { + Ok(self.sum_axis(Axis(axis))) + } + + fn mean(self, axis: usize) -> Result { + self.mean_axis(Axis(axis)).ok_or_else(|| anyhow!("no mean")) + } + } + )* + }; +} + +macro_rules! impl_frame_stats_float { + ($($t:tt),+ $(,)?) => { + $( + impl MinMax for Array<$t, D> + where + D: Dimension + RemoveAxis, + { + type Output = Array<$t, D::Smaller>; + + fn max(self, axis: usize) -> Result { + let a: Vec<_> = self + .lanes(Axis(axis)) + .into_iter() + .map(|x| { + x.iter() + .fold($t::NEG_INFINITY, |prev, curr| prev.max(*curr)) + }) + .collect(); + let mut shape = self.shape().to_vec(); + shape.remove(axis); + Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?) + } + + fn min(self, axis: usize) -> Result { + let a: Vec<_> = self + .lanes(Axis(axis)) + .into_iter() + .map(|x| { + x.iter() + .fold($t::INFINITY, |prev, curr| prev.min(*curr)) + }) + .collect(); + let mut shape = self.shape().to_vec(); + shape.remove(axis); + Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?) + } + + fn sum(self, axis: usize) -> Result { + Ok(self.sum_axis(Axis(axis))) + } + + fn mean(self, axis: usize) -> Result { + self.mean_axis(Axis(axis)).ok_or_else(|| anyhow!("no mean")) + } + } + )* + }; +} + +macro_rules! impl_frame_stats_int { + ($($t:tt),+ $(,)?) => { + $( + impl MinMax for Array<$t, D> + where + D: Dimension + RemoveAxis, + { + type Output = Array<$t, D::Smaller>; + + fn max(self, axis: usize) -> Result { + let a: Vec<_> = self + .lanes(Axis(axis)) + .into_iter() + .map(|x| *x.iter().max().unwrap()) + .collect(); + let mut shape = self.shape().to_vec(); + shape.remove(axis); + Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?) + } + + fn min(self, axis: usize) -> Result { + let a: Vec<_> = self + .lanes(Axis(axis)) + .into_iter() + .map(|x| *x.iter().min().unwrap()) + .collect(); + let mut shape = self.shape().to_vec(); + shape.remove(axis); + Ok(ArrayD::from_shape_vec(shape, a)?.into_dimensionality()?) + } + + fn sum(self, axis: usize) -> Result { + Ok(self.sum_axis(Axis(axis))) + } + + fn mean(self, axis: usize) -> Result { + self.mean_axis(Axis(axis)).ok_or_else(|| anyhow!("no mean")) + } + } + )* + }; +} + +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_float!(f32, f64); +impl_frame_stats_int!(u8, i8, u16, i16, u32, i32, u64, i64, u128, i128, usize, isize); diff --git a/src/view.rs b/src/view.rs new file mode 100644 index 0000000..87a77e4 --- /dev/null +++ b/src/view.rs @@ -0,0 +1,821 @@ +use crate::axes::{slice_info, Ax, Axis, Operation, Slice, SliceInfoElemDef}; +use crate::reader::Reader; +use crate::stats::MinMax; +use anyhow::{anyhow, Error, Result}; +use indexmap::IndexMap; +use itertools::iproduct; +use ndarray::{ + s, Array, Array0, Array1, Array2, ArrayD, Dimension, IntoDimension, Ix0, Ix1, Ix2, Ix5, IxDyn, + SliceArg, SliceInfoElem, +}; +use num::{Bounded, FromPrimitive, Zero}; +use serde::{Deserialize, Serialize}; +use serde_with::serde_as; +use std::any::type_name; +use std::collections::HashMap; +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::sync::Arc; + +fn idx_bnd(idx: isize, bnd: isize) -> Result { + if idx < -bnd { + Err(anyhow!("Index {} out of bounds {}", idx, bnd)) + } else if idx < 0 { + Ok(bnd - idx) + } else if idx < bnd { + Ok(idx) + } else { + Err(anyhow!("Index {} out of bounds {}", idx, bnd)) + } +} + +fn slc_bnd(idx: isize, bnd: isize) -> Result { + if idx < -bnd { + Err(anyhow!("Index {} out of bounds {}", idx, bnd)) + } else if idx < 0 { + Ok(bnd - idx) + } else if idx <= bnd { + Ok(idx) + } else { + Err(anyhow!("Index {} out of bounds {}", idx, bnd)) + } +} + +pub trait Number: + 'static + AddAssign + Bounded + Clone + Div + FromPrimitive + PartialOrd + Zero +{ +} +impl Number for T where + T: 'static + + AddAssign + + Bounded + + Clone + + Div + + FromPrimitive + + PartialOrd + + Zero +{ +} + +/// sliceable view on an image file +#[serde_as] +#[derive(Clone, Debug, Serialize, Deserialize)] +pub struct View { + reader: Arc, + #[serde_as(as = "Vec")] + slice: Vec, + axes: Vec, + operations: IndexMap, + dimensionality: PhantomData, +} + +impl View { + pub(crate) fn new(reader: Arc, slice: Vec, axes: Vec) -> Self { + 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 { + self.reader.series + } + + fn with_operations(mut self, operations: IndexMap) -> Self { + self.operations = operations; + self + } + + /// change the dimension into a dynamic dimension + pub fn into_dyn(self) -> View { + View { + reader: self.reader, + slice: self.slice, + axes: self.axes, + operations: self.operations, + dimensionality: PhantomData, + } + } + + /// change the dimension into a concrete dimension + pub fn into_dimensionality(self) -> Result> { + if let Some(d) = D2::NDIM { + if d == self.ndim() { + Ok(View { + reader: self.reader, + slice: self.slice, + axes: self.axes, + operations: self.operations, + dimensionality: PhantomData, + }) + } else { + Err(anyhow!("Dimensionality mismatch: {} != {}", d, self.ndim())) + } + } else { + Ok(View { + reader: self.reader, + slice: self.slice, + axes: self.axes, + operations: self.operations, + dimensionality: PhantomData, + }) + } + } + + /// the order of the axes, including axes sliced out + pub fn get_axes(&self) -> &[Axis] { + &self.axes + } + + /// the slice defining the view + pub fn get_slice(&self) -> &[SliceInfoElem] { + &self.slice + } + + /// the axes in the view + pub fn axes(&self) -> Vec { + self.axes + .iter() + .zip(self.slice.iter()) + .filter_map(|(ax, s)| { + if s.is_index() || self.operations.contains_key(ax) { + None + } else { + Some(*ax) + } + }) + .collect() + } + + pub(crate) fn op_axes(&self) -> Vec { + self.operations.keys().cloned().collect() + } + + /// the number of dimensions in the view + pub fn ndim(&self) -> usize { + if let Some(d) = D::NDIM { + d + } else { + self.shape().len() + } + } + + /// the number of pixels in the view + pub fn size(&self) -> usize { + self.shape().into_iter().product() + } + + /// the shape of the view + pub fn shape(&self) -> Vec { + let mut shape = Vec::::new(); + for (ax, s) in self.axes.iter().zip(self.slice.iter()) { + match s { + SliceInfoElem::Slice { start, end, step } => { + if !self.operations.contains_key(ax) { + if let Some(e) = end { + shape.push(((e - start).max(0) / step) as usize); + } else { + panic!("slice has no end") + } + } + } + SliceInfoElem::Index(_) => {} + SliceInfoElem::NewAxis => { + if !self.operations.contains_key(ax) { + shape.push(1); + } + } + } + } + shape + } + + fn shape_all(&self) -> Vec { + let mut shape = Vec::::new(); + for s in self.slice.iter() { + match s { + SliceInfoElem::Slice { start, end, step } => { + if let Some(e) = end { + shape.push(((e - start).max(0) / step) as usize); + } else { + panic!("slice has no end") + } + } + _ => shape.push(1), + } + } + shape + } + + /// swap two axes + pub fn swap_axes(&self, axis0: A, axis1: A) -> Result { + let idx0 = axis0.pos_op(&self.axes, &self.slice, &self.op_axes())?; + let idx1 = axis1.pos_op(&self.axes, &self.slice, &self.op_axes())?; + let mut slice = self.slice.to_vec(); + slice.swap(idx0, idx1); + let mut axes = self.axes.clone(); + axes.swap(idx0, idx1); + Ok(View::new(self.reader.clone(), slice, axes)) + } + + /// subset of gives axes will be reordered in given order + pub fn permute_axes(&self, axes: &[A]) -> Result { + let idx: Vec = axes + .iter() + .map(|a| a.pos_op(&self.axes, &self.slice, &self.op_axes()).unwrap()) + .collect(); + let mut jdx = idx.clone(); + jdx.sort(); + let mut slice = self.slice.to_vec(); + let mut axes = self.axes.clone(); + for (&i, j) in idx.iter().zip(jdx) { + slice[j] = self.slice[i]; + axes[j] = self.axes[i]; + } + Ok(View::new(self.reader.clone(), slice, axes)) + } + + /// reverse the order of the axes + pub fn transpose(&self) -> Result { + Ok(View::new( + self.reader.clone(), + self.slice.iter().rev().cloned().collect(), + self.axes.iter().rev().cloned().collect(), + )) + } + + 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), + ) + } + + /// maximum along axis + pub fn max_proj(&self, axis: A) -> Result> { + self.operate(axis, Operation::Max) + } + + /// minimum along axis + pub fn min_proj(&self, axis: A) -> Result> { + self.operate(axis, Operation::Min) + } + + /// sum along axis + pub fn sum_proj(&self, axis: A) -> Result> { + self.operate(axis, Operation::Sum) + } + + /// mean along axis + pub fn mean_proj(&self, axis: A) -> Result> { + self.operate(axis, Operation::Mean) + } + + /// created a new sliced view + pub fn slice(&self, info: I) -> Result> + where + I: SliceArg, + { + if self.slice.out_ndim() < info.in_ndim() { + return Err(Error::msg("not enough free dimensions")); + } + let info = info.as_ref(); + let mut n_idx = 0; + let mut r_idx = 0; + 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) { + ( + SliceInfoElem::Slice { + start: info_start, + end: info_end, + step: info_step, + }, + 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 new_step = (step * info_step).abs(); + new_slice.push(SliceInfoElem::Slice { + start: new_start, + end: Some(new_end), + step: new_step, + }); + new_axes.push(*a); + 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)) + } else { + new_slice.push(SliceInfoElem::Index(start + step.abs() * k)); + } + new_axes.push(*a); + n_idx += 1; + r_idx += 1; + } + (SliceInfoElem::NewAxis, SliceInfoElem::NewAxis) => { + new_slice.push(SliceInfoElem::NewAxis); + new_slice.push(SliceInfoElem::NewAxis); + new_axes.push(Axis::New); + new_axes.push(Axis::New); + 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, _) => { + new_slice.push(SliceInfoElem::NewAxis); + new_axes.push(Axis::New); + n_idx += 1; + } + (_, SliceInfoElem::Index(k)) => { + new_slice.push(SliceInfoElem::Index(*k)); + new_axes.push(*a); + r_idx += 1; + } + } + } + } + debug_assert_eq!(r_idx, reader_slice.len()); + while n_idx < info.len() { + debug_assert!(info[n_idx].is_new_axis()); + new_slice.push(SliceInfoElem::NewAxis); + new_axes.push(Axis::New); + n_idx += 1; + } + Ok(View::new(self.reader.clone(), new_slice, new_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() + } + + /// the pixel intensity at a given index + pub fn item_at(&self, index: &[isize]) -> Result + where + T: Number, + 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())?; + let arr = view.as_array()?; + Ok(arr.first().unwrap().clone()) + } + + /// collect the view into an ndarray + pub fn as_array(&self) -> Result> + where + T: Number, + ArrayD: MinMax, + Array1: MinMax, + Array2: MinMax, + { + Ok(self.as_array_dyn()?.into_dimensionality()?) + } + + /// collect the view into a dynamic-dimension ndarray + pub fn as_array_dyn(&self) -> Result> + where + T: Number, + ArrayD: MinMax, + Array1: MinMax, + Array2: MinMax, + { + let mut op_xy = IndexMap::new(); + if let Some((&ax, op)) = self.operations.first() { + if (ax == Axis::X) || (ax == Axis::Y) { + op_xy.insert(ax, op.clone()); + if let Some((&ax2, op2)) = self.operations.get_index(1) { + if (ax2 == Axis::X) || (ax2 == Axis::Y) { + op_xy.insert(ax2, op2.clone()); + } + } + } + } + let op_czt = if let Some((&ax, op)) = self.operations.get_index(op_xy.len()) { + IndexMap::from([(ax, op.clone())]) + } else { + IndexMap::new() + }; + let mut shape_out = Vec::new(); + let mut slice = Vec::new(); + let mut ax_out = Vec::new(); + 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") + } + } + SliceInfoElem::Index(_) => {} + SliceInfoElem::NewAxis => { + shape_out.push(1); + slice.push(SliceInfoElem::Index(0)); + ax_out.push(*a); + } + } + } + let mut slice_reader = vec![Slice::empty(); 5]; + let mut xy_dim = 0usize; + let shape = [ + self.size_c as isize, + self.size_z as isize, + self.size_t as isize, + self.size_y as isize, + self.size_x as isize, + ]; + for (s, &axis) in self.slice.iter().zip(&self.axes) { + match axis { + Axis::New => {} + _ => match s { + SliceInfoElem::Slice { start, end, step } => { + if let Axis::X | Axis::Y = axis { + if !op_xy.contains_key(&axis) { + xy_dim += 1; + } + } + slice_reader[axis as usize] = Slice::new( + idx_bnd(*start, shape[axis as usize])?, + slc_bnd(end.unwrap(), shape[axis as usize])?, + *step, + ); + } + SliceInfoElem::Index(j) => { + slice_reader[axis as usize] = Slice::new( + idx_bnd(*j, shape[axis as usize])?, + slc_bnd(*j + 1, shape[axis as usize])?, + 1, + ); + } + SliceInfoElem::NewAxis => panic!("axis cannot be a new axis"), + }, + } + } + let xy = [ + self.slice[Axis::Y.pos(&self.axes, &self.slice)?], + self.slice[Axis::X.pos(&self.axes, &self.slice)?], + ]; + let mut array = if let Some((_, op)) = op_czt.first() { + match op { + Operation::Max => { + ArrayD::::from_elem(shape_out.into_dimension(), T::min_value()) + } + Operation::Min => { + ArrayD::::from_elem(shape_out.into_dimension(), T::max_value()) + } + _ => ArrayD::::zeros(shape_out.into_dimension()), + } + } else { + ArrayD::::zeros(shape_out.into_dimension()) + }; + let size_c = self.reader.size_c as isize; + let size_z = self.reader.size_z as isize; + let size_t = self.reader.size_t as isize; + + let mut axes_out_idx = [None; 5]; + for (i, ax) in ax_out.iter().enumerate() { + axes_out_idx[*ax as usize] = Some(i); + } + + for (c, z, t) in iproduct!(&slice_reader[0], &slice_reader[1], &slice_reader[2]) { + if let Some(i) = axes_out_idx[0] { + slice[i] = SliceInfoElem::Index(c) + }; + if let Some(i) = axes_out_idx[1] { + slice[i] = SliceInfoElem::Index(z) + }; + if let Some(i) = axes_out_idx[2] { + slice[i] = SliceInfoElem::Index(t) + }; + let frame = self.reader.get_frame( + (c % size_c) as usize, + (z % size_z) as usize, + (t % size_t) as usize, + )?; + + let arr_frame: Array2 = frame.try_into()?; + let arr_frame = match xy_dim { + 0 => { + if op_xy.contains_key(&Axis::X) && op_xy.contains_key(&Axis::Y) { + let xys = slice_info::(&xy)?; + let (&ax0, op0) = op_xy.first().unwrap(); + 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() + } 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() + } else { + let xys = slice_info::(&xy)?; + arr_frame.slice(xys).to_owned().into_dyn() + } + } + 1 => { + 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: &Array1 = unsafe { transmute(&b) }; + c.to_owned().into_dyn() + } else { + let xys = slice_info::(&xy)?; + arr_frame.slice(xys).to_owned().into_dyn() + } + } + 2 => { + let xys = slice_info::(&xy)?; + if axes_out_idx[4] < axes_out_idx[3] { + arr_frame.t().slice(xys).to_owned().into_dyn() + } else { + arr_frame.slice(xys).to_owned().into_dyn() + } + } + _ => { + panic!("xy cannot be 3d or more"); + } + }; + if let Some((_, op)) = op_czt.first() { + match op { + Operation::Max => { + array + .slice_mut(slice.as_slice()) + .zip_mut_with(&arr_frame, |x, y| { + *x = if *x >= *y { x.clone() } else { y.clone() } + }); + } + Operation::Min => { + array + .slice_mut(slice.as_slice()) + .zip_mut_with(&arr_frame, |x, y| { + *x = if *x < *y { x.clone() } else { y.clone() } + }); + } + Operation::Sum => { + array + .slice_mut(slice.as_slice()) + .zip_mut_with(&arr_frame, |x, y| *x += y.clone()); + } + Operation::Mean => { + array + .slice_mut(slice.as_slice()) + .zip_mut_with(&arr_frame, |x, y| *x += y.clone()); + } + } + } else { + array.slice_mut(slice.as_slice()).assign(&arr_frame) + } + } + let mut out = Some(array); + let 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) { + let arr = out.take().unwrap(); + let _ = out.insert(unsafe { transmute_copy(&op.operate(arr, idx)?) }); + } + } + let mut n = 1; + for (&ax, size) in self.axes.iter().zip(self.shape_all().iter()) { + if let Some(Operation::Mean) = self.operations.get(&ax) { + if (ax == Axis::C) || (ax == Axis::Z) || (ax == Axis::T) { + n *= size; + } + } + } + let array = if n == 1 { + out.take().unwrap() + } else { + let m = T::from_usize(n).unwrap_or_else(|| T::zero()); + out.take().unwrap().mapv(|x| x / m.clone()) + }; + Ok(array) + } + + /// retrieve a single frame at czt, sliced accordingly + pub fn get_frame(&self, c: isize, z: isize, t: isize) -> Result> + where + T: Number, + ArrayD: MinMax, + Array1: MinMax, + Array2: MinMax, + { + 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, + { + let arr: ArrayD = self.as_array_dyn()?; + Ok(match operation { + Operation::Max => arr + .flatten() + .into_iter() + .reduce(|a, b| if a > b { a } else { b }) + .unwrap_or_else(|| T::min_value()), + Operation::Min => arr + .flatten() + .into_iter() + .reduce(|a, b| if a < b { a } else { b }) + .unwrap_or_else(|| T::max_value()), + Operation::Sum => arr.flatten().into_iter().sum(), + Operation::Mean => { + arr.flatten().into_iter().sum::() + / T::from_usize(arr.len()).ok_or_else(|| { + anyhow!("cannot convert {} into {}", arr.len(), type_name::()) + })? + } + }) + } + + /// maximum intensity + pub fn max(&self) -> Result + where + T: Number + Sum, + ArrayD: MinMax, + Array1: MinMax, + Array2: MinMax, + { + self.get_stat(Operation::Max) + } + + /// minimum intensity + pub fn min(&self) -> Result + where + T: Number + Sum, + ArrayD: MinMax, + Array1: MinMax, + Array2: MinMax, + { + self.get_stat(Operation::Min) + } + + /// sum intensity + pub fn sum(&self) -> Result + where + T: Number + Sum, + ArrayD: MinMax, + Array1: MinMax, + Array2: MinMax, + { + self.get_stat(Operation::Sum) + } + + /// mean intensity + pub fn mean(&self) -> Result + where + T: Number + Sum, + ArrayD: MinMax, + Array1: MinMax, + Array2: MinMax, + { + self.get_stat(Operation::Mean) + } +} + +impl Deref for View { + type Target = Reader; + + fn deref(&self) -> &Self::Target { + self.reader.as_ref() + } +} + +impl TryFrom> for Array +where + T: Number, + D: Dimension, + ArrayD: MinMax, + Array1: MinMax, + Array2: MinMax, +{ + type Error = Error; + + fn try_from(view: View) -> Result { + view.as_array() + } +} + +impl TryFrom<&View> for Array +where + T: Number, + D: Dimension, + ArrayD: MinMax, + Array1: MinMax, + Array2: MinMax, +{ + type Error = Error; + + fn try_from(view: &View) -> Result { + view.as_array() + } +} + +/// trait to define a function to retrieve the only item in a 0d array +pub trait Item { + fn item(&self) -> Result + where + T: Number, + ArrayD: MinMax, + Array1: MinMax, + Array2: MinMax; +} + +impl Item for View { + fn item(&self) -> Result + where + T: Number, + ArrayD: MinMax, + Array1: MinMax, + Array2: MinMax, + { + Ok(self + .as_array()? + .first() + .ok_or_else(|| anyhow!("Empty view"))? + .clone()) + } +} diff --git a/tests/test_open.py b/tests/test_open.py index 4a161d8..3e01d8e 100644 --- a/tests/test_open.py +++ b/tests/test_open.py @@ -6,16 +6,20 @@ import pytest from ndbioimage import Imread -@pytest.mark.parametrize('file', - [file for file in (Path(__file__).parent / 'files').iterdir() if not file.suffix == '.pzl']) +@pytest.mark.parametrize( + "file", + [ + file + for file in (Path(__file__).parent / "files").iterdir() + if not file.suffix == ".pzl" + ], +) def test_open(file): - with Imread(file) as im: - mean = im[dict(c=0, z=0, t=0)].mean() + with Imread(file, axes="cztyx") as im: + mean = im[0, 0, 0].mean() b = pickle.dumps(im) jm = pickle.loads(b) - assert jm[dict(c=0, z=0, t=0)].mean() == mean - v = im.view() - assert v[dict(c=0, z=0, t=0)].mean() == mean - b = pickle.dumps(v) - w = pickle.loads(b) - assert w[dict(c=0, z=0, t=0)].mean() == mean + assert jm.get_frame(0, 0, 0).mean() == mean + b = pickle.dumps(im) + jm = pickle.loads(b) + assert jm[0, 0, 0].mean() == mean diff --git a/tests/test_slicing.py b/tests/test_slicing.py index fb2845c..6aa04b9 100644 --- a/tests/test_slicing.py +++ b/tests/test_slicing.py @@ -12,19 +12,24 @@ from ndbioimage import Imread @pytest.fixture def array(): - return np.random.randint(0, 255, (64, 64, 2, 3, 4), 'uint8') + return np.random.randint(0, 255, (64, 64, 2, 3, 4), "uint8") + @pytest.fixture() def image(array): with tempfile.TemporaryDirectory() as folder: file = Path(folder) / "test.tif" - tiffwrite(file, array, 'yxczt') - with Imread(file, axes='yxczt') as im: + tiffwrite(file, array, "yxczt") + with Imread(file, axes="yxczt") as im: yield im -@pytest.mark.parametrize('s', combinations_with_replacement( - (0, -1, 1, slice(None), slice(0, 1), slice(-1, 0), slice(1, 1)), 5)) +@pytest.mark.parametrize( + "s", + combinations_with_replacement( + (0, -1, 1, slice(None), slice(0, 1), slice(-1, 0), slice(1, 1)), 5 + ), +) def test_slicing(s, image, array): s_im, s_a = image[s], array[s] if isinstance(s_a, Number): diff --git a/tests/test_ufuncs.py b/tests/test_ufuncs.py index 2fa52b8..a10e2ae 100644 --- a/tests/test_ufuncs.py +++ b/tests/test_ufuncs.py @@ -11,21 +11,42 @@ from ndbioimage import Imread @pytest.fixture def array(): - return np.random.randint(0, 255, (64, 64, 2, 3, 4), 'uint16') + return np.random.randint(0, 255, (64, 64, 2, 3, 4), "uint16") + @pytest.fixture() def image(array): with tempfile.TemporaryDirectory() as folder: file = Path(folder) / "test.tif" - tiffwrite(file, array, 'yxczt') - with Imread(file, axes='yxczt') as im: + tiffwrite(file, array, "yxczt") + with Imread(file, axes="yxczt") as im: yield im -@pytest.mark.parametrize('fun_and_axis', product( - (np.sum, np.nansum, np.min, np.nanmin, np.max, np.nanmax, np.argmin, np.argmax, - np.mean, np.nanmean, np.var, np.nanvar, np.std, np.nanstd), (None, 0, 1, 2, 3, 4))) +@pytest.mark.parametrize( + "fun_and_axis", + product( + ( + np.sum, + np.nansum, + np.min, + np.nanmin, + np.max, + np.nanmax, + np.argmin, + np.argmax, + np.mean, + np.nanmean, + np.var, + np.nanvar, + np.std, + np.nanstd, + ), + (None, 0, 1, 2, 3, 4), + ), +) def test_ufuncs(fun_and_axis, image, array): fun, axis = fun_and_axis - assert np.all(np.isclose(fun(image, axis), fun(array, axis))), \ - f'function {fun.__name__} over axis {axis} does not give the correct result' + assert np.all(np.isclose(np.asarray(fun(image, axis)), fun(array, axis))), ( + f"function {fun.__name__} over axis {axis} does not give the correct result" + )