diff --git a/Cargo.toml b/Cargo.toml index 05692c0..2433823 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,9 +3,9 @@ resolver = "3" members = ["ndbioimage", "image-tiff", "image-registration"] [workspace.package] -version = "0.2.0" +version = "0.3.0" edition = "2024" -rust-version = "1.94.0" +rust-version = "1.94.1" authors = ["Wim Pomp "] license = "MIT OR Apache-2.0" description = "Read bio image formats using the bio-formats java package." diff --git a/image-registration b/image-registration index 7eb39a2..5c580dc 160000 --- a/image-registration +++ b/image-registration @@ -1 +1 @@ -Subproject commit 7eb39a2fbf1c0ec5aeb118a47a5f641525434ec7 +Subproject commit 5c580dc80f8552b91a5c2dae4fe1f1f8e3142415 diff --git a/image-tiff b/image-tiff index 9b452e9..a6e8ed7 160000 --- a/image-tiff +++ b/image-tiff @@ -1 +1 @@ -Subproject commit 9b452e93b708db6844a4402f4aac76311afbf44b +Subproject commit a6e8ed772de58cedf3ceb9a5720fe52997433f62 diff --git a/ndbioimage/Cargo.toml b/ndbioimage/Cargo.toml index 0483a61..a1635f6 100644 --- a/ndbioimage/Cargo.toml +++ b/ndbioimage/Cargo.toml @@ -37,7 +37,7 @@ ome-metadata = "0.6" ordered-float = { version = "5", optional = true } phf = { version = "0.14", features = ["macros"] } postcard = { version = "1", features = ["use-std"], optional = true } -pyo3 = { version = "0.29", features = ["abi3-py310", "anyhow", "generate-import-lib"], optional = true } +pyo3 = { version = "0.29", features = ["abi3-py310", "anyhow", "generate-import-lib", "multiple-pymethods"], optional = true } pyo3-stub-gen = { version = "0.23", optional = true } rayon = { version = "1", optional = true } regex = "1" diff --git a/ndbioimage/pyproject.toml b/ndbioimage/pyproject.toml index ca28e32..97652a9 100644 --- a/ndbioimage/pyproject.toml +++ b/ndbioimage/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "ndbioimage" -version = "2027.0.5" +version = "2027.0.6" requires-python = ">=3.10" classifiers = [ "License :: OSI Approved :: MIT License", diff --git a/ndbioimage/src/axes.rs b/ndbioimage/src/axes.rs index b28a97a..a374e2c 100644 --- a/ndbioimage/src/axes.rs +++ b/ndbioimage/src/axes.rs @@ -198,9 +198,9 @@ impl<'de> DeserializeAs<'de, SliceInfoElem> for SliceInfoElemDef { #[derive(Clone, Debug)] pub(crate) struct Slice { - start: isize, - end: isize, - step: isize, + pub(crate) start: isize, + pub(crate) end: isize, + pub(crate) step: isize, } impl Slice { diff --git a/ndbioimage/src/cache.rs b/ndbioimage/src/cache.rs index f74f2e0..0ea21b5 100644 --- a/ndbioimage/src/cache.rs +++ b/ndbioimage/src/cache.rs @@ -1,14 +1,22 @@ use crate::axes::{Axis, Operation}; +use crate::readers::Frame; +#[cfg(feature = "transforms")] +use image_registration::transform::Transform; use indexmap::{Equivalent, IndexMap}; +#[cfg(feature = "transforms")] +use ndarray::IxDyn; use ndarray::{ArrayD, SliceInfoElem}; use std::any::Any; +#[cfg(feature = "transforms")] +use std::borrow::Cow; +use std::cmp::PartialEq; use std::fmt::{Debug, Formatter}; use std::hash::{Hash, Hasher}; use std::path::{Path, PathBuf}; use std::sync::{Arc, Mutex, OnceLock}; use std::time::SystemTime; - -use crate::readers::Frame; +#[cfg(feature = "transforms")] +use crate::transforms::Transforms; /// maximum number of frames held in the cache const DEFAULT_FRAME_CACHE_SIZE: usize = 128; @@ -132,7 +140,15 @@ impl ReaderKeyRef<'_> { } } -pub(crate) type FrameKey = (ReaderKey, usize, usize, usize); +#[derive(Debug, PartialEq, Eq, Hash)] +pub(crate) struct FrameKey { + pub(crate) reader_key: ReaderKey, + pub(crate) c: usize, + pub(crate) z: usize, + pub(crate) t: usize, + #[cfg(feature = "transforms")] + pub(crate) transform: Option<(Transform, usize)>, +} /// borrowed view of [`FrameKey`] for cache lookups without allocation. #[derive(Debug, PartialEq, Eq)] @@ -141,6 +157,8 @@ pub(crate) struct FrameKeyRef<'a> { pub(crate) c: usize, pub(crate) z: usize, pub(crate) t: usize, + #[cfg(feature = "transforms")] + pub(crate) transform: Option<(Cow<'a, Transform>, usize)>, } impl Hash for FrameKeyRef<'_> { @@ -149,19 +167,50 @@ impl Hash for FrameKeyRef<'_> { self.c.hash(state); self.z.hash(state); self.t.hash(state); + #[cfg(feature = "transforms")] + self.transform.hash(state); } } impl Equivalent for FrameKeyRef<'_> { fn equivalent(&self, key: &FrameKey) -> bool { - let (rk, c, z, t) = key; - self.reader.equivalent(rk) && self.c == *c && self.z == *z && self.t == *t + let FrameKey { + reader_key, + c, + z, + t, + #[cfg(feature = "transforms")] + transform, + } = key; + #[cfg(feature = "transforms")] + { + self.reader.equivalent(reader_key) + && self.c == *c + && self.z == *z + && self.t == *t + && self.transform.as_ref().map(|(t, i)| (t.as_ref(), i)) + == transform.as_ref().map(|(t, i)| (t, i)) + } + #[cfg(not(feature = "transforms"))] + { + self.reader.equivalent(reader_key) && self.c == *c && self.z == *z && self.t == *t + } } } impl FrameKeyRef<'_> { pub(crate) fn to_owned(&self) -> FrameKey { - (self.reader.to_owned(), self.c, self.z, self.t) + FrameKey { + reader_key: self.reader.to_owned(), + c: self.c, + z: self.z, + t: self.t, + #[cfg(feature = "transforms")] + transform: self + .transform + .as_ref() + .map(|(t, i)| (t.clone().into_owned(), *i)), + } } } @@ -173,6 +222,8 @@ pub(crate) struct ArrayKey { pub(crate) slice: Vec, pub(crate) axes: Vec, pub(crate) operations: Vec<(Axis, Operation)>, + #[cfg(feature = "transforms")] + pub(crate) transforms: Transforms, } /// borrowed view of [`ArrayKey`] for cache lookups without allocation. @@ -183,6 +234,8 @@ pub(crate) struct ArrayKeyRef<'a> { pub(crate) slice: &'a [SliceInfoElem], pub(crate) axes: &'a [Axis], pub(crate) operations: &'a IndexMap, + #[cfg(feature = "transforms")] + pub(crate) transforms: &'a Transforms, } impl Hash for ArrayKeyRef<'_> { @@ -196,21 +249,40 @@ impl Hash for ArrayKeyRef<'_> { ax.hash(state); op.hash(state); } + #[cfg(feature = "transforms")] + self.transforms.hash(state); } } impl Equivalent for ArrayKeyRef<'_> { fn equivalent(&self, key: &ArrayKey) -> bool { - self.reader.equivalent(&key.reader) - && self.dtype == key.dtype - && self.slice == key.slice - && self.axes == key.axes - && self.operations.len() == key.operations.len() - && self + #[cfg(feature = "transforms")] + { + self.reader.equivalent(&key.reader) + && self.dtype == key.dtype + && self.slice == key.slice + && self.axes == key.axes + && self.operations.len() == key.operations.len() + && self .operations .iter() .zip(&key.operations) .all(|((ax, op), (key_ax, key_op))| ax == key_ax && op == key_op) + && self.transforms == key.transforms + } + #[cfg(not(feature = "transforms"))] + { + self.reader.equivalent(&key.reader) + && self.dtype == key.dtype + && self.slice == key.slice + && self.axes == key.axes + && self.operations.len() == key.operations.len() + && self + .operations + .iter() + .zip(&key.operations) + .all(|((ax, op), (key_ax, key_op))| ax == key_ax && op == key_op) + } } } @@ -226,6 +298,8 @@ impl ArrayKeyRef<'_> { .iter() .map(|(ax, op)| (*ax, op.clone())) .collect(), + #[cfg(feature = "transforms")] + transforms: self.transforms.to_owned(), } } } diff --git a/ndbioimage/src/error.rs b/ndbioimage/src/error.rs index ea74b58..5b548c3 100644 --- a/ndbioimage/src/error.rs +++ b/ndbioimage/src/error.rs @@ -179,6 +179,10 @@ pub enum Error { it cannot use. read the file before forking, or start the child with spawn/exec instead." )] ForkedJvm(u32, u32), + #[error("{0}D transformation is not implemented")] + TransformationNotImplemented(usize), + #[error("no channel {0} present")] + NoChannel(usize), } impl Error { diff --git a/ndbioimage/src/lib.rs b/ndbioimage/src/lib.rs index ab9cfcf..b4a182e 100644 --- a/ndbioimage/src/lib.rs +++ b/ndbioimage/src/lib.rs @@ -1,4 +1,4 @@ -#![cfg_attr(docsrs, feature(doc_cfg))] +#![cfg_attr(docsrs, feature(doc_cfg, doc_auto_cfg))] //! The ndbioimage crate exposes (bio) images as a struct that can be sliced like an ndarray Array //! (Rust), but without loading the whole image into memory, reading from the file only when needed. //! Some metadata is read @@ -16,7 +16,7 @@ //! # fn main() -> Result<(), ndbioimage::error::Error> { //! let path = "/path/to/file"; //! let reader = DynReader::new(&path, 0, 0)?; -//! println!("shape: {}", reader.shape()); +//! println!("shape: {}", reader.raw_shape()); //! let frame = reader.get_frame(0, 0, 0)?; //! if let Ok(arr) = >>::try_into(frame) { //! println!("{:?}", arr); @@ -70,7 +70,7 @@ pub mod readers; pub mod tiffwrite; #[cfg(feature = "transforms")] /// image registration and warping images using affine transforms -mod transforms; +pub mod transforms; mod utils; /// main entry point for the application @@ -338,7 +338,7 @@ mod tests { fn read_ser() -> Result<(), Error> { let file = "czi/Experiment-2029.czi"; let reader = open(file)?; - println!("shape: {}", reader.shape()); + println!("shape: {}", reader.raw_shape()); let frame = reader.get_frame(0, 0, 0)?; if let Ok(arr) = >>::try_into(frame) { println!("{:?}", arr); diff --git a/ndbioimage/src/py.rs b/ndbioimage/src/py.rs index 19486b6..481d1e5 100644 --- a/ndbioimage/src/py.rs +++ b/ndbioimage/src/py.rs @@ -67,6 +67,7 @@ fn main() -> PyResult<()> { #[pymodule] #[pyo3(name = "ndbioimage_rs")] mod ndbioimage_rs { + #[cfg(feature = "transforms")] use pyo3::prelude::*; #[pymodule_export] diff --git a/ndbioimage/src/py/imread.rs b/ndbioimage/src/py/imread.rs index 55fb7d0..b25176c 100644 --- a/ndbioimage/src/py/imread.rs +++ b/ndbioimage/src/py/imread.rs @@ -3,8 +3,6 @@ use crate::error::Error; use crate::metadata::Metadata; use crate::movie::MovieOptions; -#[cfg(feature = "transforms")] -use crate::py::transforms::{PyTransform, PyTransforms}; use crate::readers::{DynReader, PixelType, Reader}; use crate::view::{Item, View}; use itertools::Itertools; @@ -27,7 +25,7 @@ pub(crate) use pyo3::{ use pyo3_stub_gen::derive::*; use pyo3_stub_gen::inventory::submit; use serde::{Deserialize, Serialize}; -use std::collections::{HashMap, HashSet}; +use std::collections::HashSet; use std::path::PathBuf; /// class to read image files, while taking good care of important metadata, @@ -61,7 +59,7 @@ use std::path::PathBuf; #[pyclass(subclass, from_py_object, name = "Imread", module = "ndbioimage")] #[derive(Clone, Debug, Serialize, Deserialize)] pub(crate) struct PyView { - view: View, + pub(crate) view: View, dtype: PixelType, #[serde(skip)] ome: Ome, @@ -296,141 +294,6 @@ impl PyView { todo!() } - /// return a new view with transformations applied (channel alignment, drift correction) - /// backward compatibility - #[cfg(feature = "transforms")] - #[pyo3(signature = (channels = true, drift = false, file = None, bead_files = None, main_channel = None, default_transform = None))] - fn with_transform<'py>( - &self, - channels: bool, - drift: bool, - file: Option>, - bead_files: Option>, - main_channel: Option, - default_transform: Option>, - ) -> PyResult { - self.with_transforms( - channels, - drift, - file, - bead_files, - main_channel, - default_transform, - ) - } - - /// return a new view with transformations applied (channel alignment, drift correction) - #[cfg(feature = "transforms")] - #[expect(unused_variables)] - #[pyo3(signature = (channels = true, drift = false, file = None, bead_files = None, main_channel = None, default_transform = None))] - fn with_transforms<'py>( - &self, - channels: bool, - drift: bool, - file: Option>, - bead_files: Option>, - main_channel: Option, - default_transform: Option>, - ) -> PyResult { - todo!() - } - - /// get the transformation - #[cfg(feature = "transforms")] - #[getter] - fn get_transforms(&self) -> PyResult { - Ok(PyTransforms { - inner: self.view.get_transforms().clone(), - }) - } - - /// set the transformation - #[cfg(feature = "transforms")] - #[setter] - fn set_transforms(&mut self, transform: PyTransforms) { - self.view.set_transforms(transform.inner); - } - - #[cfg(feature = "transforms")] - fn load_transforms_from_yaml(&mut self, path: PathBuf) -> PyResult<()> { - self.view.load_transforms_from_yaml(path.as_path())?; - Ok(()) - } - - #[cfg(feature = "transforms")] - fn calculate_channel_transforms_2d( - &self, - py: Python, - main_channel: usize, - ) -> PyResult> { - py.detach(|| { - Ok(self - .view - .calculate_channel_transforms_2d(main_channel)? - .into_iter() - .map(|(n, t)| { - ( - n, - PyTransform { - inner: t.into_dyn(), - }, - ) - }) - .collect()) - }) - } - - #[cfg(feature = "transforms")] - fn calculate_channel_transforms_3d( - &self, - py: Python, - main_channel: usize, - ) -> PyResult> { - py.detach(|| { - Ok(self - .view - .calculate_channel_transforms_3d(main_channel)? - .into_iter() - .map(|(n, t)| { - ( - n, - PyTransform { - inner: t.into_dyn(), - }, - ) - }) - .collect()) - }) - } - - #[cfg(feature = "transforms")] - fn calculate_drift_transforms_2d(&self, py: Python) -> PyResult> { - py.detach(|| { - Ok(self - .view - .calculate_drift_transforms_2d()? - .into_iter() - .map(|t| PyTransform { - inner: t.into_dyn(), - }) - .collect()) - }) - } - - #[cfg(feature = "transforms")] - fn calculate_drift_transforms_3d(&self, py: Python) -> PyResult> { - py.detach(|| { - Ok(self - .view - .calculate_drift_transforms_3d()? - .into_iter() - .map(|t| PyTransform { - inner: t.into_dyn(), - }) - .collect()) - }) - } - #[gen_stub(override_return_type(type_repr="numpy.ndarray | int | float", imports=("numpy")))] fn squeeze<'py>(&self, py: Python<'py>) -> PyResult> { let view = self.view.squeeze()?; @@ -1481,6 +1344,7 @@ impl PyView { ))] axes: Option>>, ) -> PyResult { + // TODO: transpose("xyz"), transpose(*axes) let view = if let Some(axes) = axes { let ax = axes .into_iter() diff --git a/ndbioimage/src/py/transforms.rs b/ndbioimage/src/py/transforms.rs index f488ab3..824b0b0 100644 --- a/ndbioimage/src/py/transforms.rs +++ b/ndbioimage/src/py/transforms.rs @@ -1,4 +1,5 @@ use crate::error::Error; +use crate::py::imread::PyView; use crate::transforms::Transforms; use image_registration::metric::{SamplingArg, Sigma}; use image_registration::register::{Optimizer, RegistrationStep}; @@ -709,6 +710,138 @@ pub(crate) struct PyTransforms { pub(crate) inner: Transforms, } +#[gen_stub_pymethods] +#[pymethods] +impl PyView { + /// return a new view with transformations applied (channel alignment, drift correction) + /// backward compatibility + #[cfg(feature = "transforms")] + #[pyo3(signature = (channels = true, drift = false, file = None, bead_files = None, main_channel = None, default_transform = None))] + fn with_transform<'py>( + &self, + channels: bool, + drift: bool, + file: Option>, + bead_files: Option>, + main_channel: Option, + default_transform: Option>, + ) -> PyResult { + self.with_transforms( + channels, + drift, + file, + bead_files, + main_channel, + default_transform, + ) + } + + /// return a new view with transformations applied (channel alignment, drift correction) + #[cfg(feature = "transforms")] + #[expect(unused_variables)] + #[pyo3(signature = (channels = true, drift = false, file = None, bead_files = None, main_channel = None, default_transform = None))] + fn with_transforms<'py>( + &self, + channels: bool, + drift: bool, + file: Option>, + bead_files: Option>, + main_channel: Option, + default_transform: Option>, + ) -> PyResult { + todo!() + } + + /// get the transformation + #[getter] + fn get_transforms(&self) -> PyResult { + Ok(PyTransforms { + inner: self.view.get_transforms().clone(), + }) + } + + /// set the transformation + #[setter] + fn set_transforms(&mut self, transform: PyTransforms) { + self.view.set_transforms(transform.inner); + } + + fn load_transforms_from_yaml(&mut self, path: PathBuf) -> PyResult<()> { + self.view.load_transforms_from_yaml(path.as_path())?; + Ok(()) + } + + fn calculate_channel_transforms_2d( + &self, + py: Python, + main_channel: usize, + ) -> PyResult> { + py.detach(|| { + Ok(self + .view + .calculate_channel_transforms_2d(main_channel)? + .into_iter() + .map(|(n, t)| { + ( + n, + PyTransform { + inner: t.into_dyn(), + }, + ) + }) + .collect()) + }) + } + + fn calculate_channel_transforms_3d( + &self, + py: Python, + main_channel: usize, + ) -> PyResult> { + py.detach(|| { + Ok(self + .view + .calculate_channel_transforms_3d(main_channel)? + .into_iter() + .map(|(n, t)| { + ( + n, + PyTransform { + inner: t.into_dyn(), + }, + ) + }) + .collect()) + }) + } + + fn calculate_drift_transforms_2d(&self, py: Python) -> PyResult> { + py.detach(|| { + Ok(self + .view + .calculate_drift_transforms_2d()? + .into_iter() + .map(|t| PyTransform { + inner: t.into_dyn(), + }) + .collect()) + }) + } + + fn calculate_drift_transforms_3d(&self, py: Python) -> PyResult> { + py.detach(|| { + Ok(self + .view + .calculate_drift_transforms_3d()? + .into_iter() + .map(|t| PyTransform { + inner: t.into_dyn(), + }) + .collect()) + }) + } +} + #[gen_stub_pymethods] #[pymethods] impl PyTransforms { diff --git a/ndbioimage/src/readers.rs b/ndbioimage/src/readers.rs index 4827f3d..c1e2d19 100644 --- a/ndbioimage/src/readers.rs +++ b/ndbioimage/src/readers.rs @@ -2,7 +2,7 @@ use crate::axes::{Axis, Shape}; use crate::error::Error; use crate::view::View; use ndarray::{Array, Dimension, Ix2, Ix5, s}; -use num::{FromPrimitive, Zero}; +use num::{FromPrimitive, One, Zero}; use ome_metadata::Ome; use regex::Regex; use serde::{Deserialize, Serialize}; @@ -148,6 +148,8 @@ impl PixelType { /// array data with a specific pixel type #[derive(Clone, Debug)] pub enum ArrayT { + /// boolean array + Bool(Array), /// signed 8-bit integer array I8(Array), /// unsigned 8-bit integer array @@ -197,7 +199,7 @@ pub trait Reader: Clone + Sized + Debug + Send + Hash + Into { /// get a sliceable view on the image file fn view(&self) -> View { - let shape = self.shape(); + let shape = self.raw_shape(); let slice = s![0..shape.c, 0..shape.z, 0..shape.t, 0..shape.y, 0..shape.x,]; View::new( self.clone(), @@ -216,7 +218,7 @@ pub trait Reader: Clone + Sized + Debug + Send + Hash + Into { /// the position index fn position(&self) -> usize; /// the shape of the image data - fn shape(&self) -> &Shape; + fn raw_shape(&self) -> &Shape; /// the pixel type of the image data fn pixel_type(&self) -> &PixelType; /// get all available positions for a given series @@ -314,6 +316,7 @@ macro_rules! try_into_array_body { ($self:expr, $map:ident) => {{ let mut err = Ok(()); let arr = match $self { + ArrayT::Bool(v) => v.$map(|x| if x { T::one() } else { T::zero() }), ArrayT::I8(v) => v.$map(|x| { T::from_i8(x).unwrap_or_else(|| { err = Err(Error::Cast(x.to_string(), type_name::().to_string())); @@ -397,7 +400,7 @@ macro_rules! try_into_array_body { impl TryInto> for ArrayT where D: Dimension, - T: FromPrimitive + Zero + 'static, + T: FromPrimitive + Zero + One + 'static, { type Error = Error; @@ -409,7 +412,7 @@ where impl TryInto> for &ArrayT where D: Dimension, - T: Clone + FromPrimitive + Zero + 'static, + T: Clone + FromPrimitive + Zero + One + 'static, { type Error = Error; @@ -577,18 +580,18 @@ impl Reader for DynReader { } } - fn shape(&self) -> &Shape { + fn raw_shape(&self) -> &Shape { match self { #[cfg(feature = "tiff")] - DynReader::Tiff(r) => r.shape(), + DynReader::Tiff(r) => r.raw_shape(), #[cfg(feature = "tiffseq")] - DynReader::TiffSeq(r) => r.shape(), + DynReader::TiffSeq(r) => r.raw_shape(), #[cfg(feature = "czi")] - DynReader::Czi(r) => r.shape(), + DynReader::Czi(r) => r.raw_shape(), #[cfg(feature = "bioformats_rust")] - DynReader::BioFormatsRust(r) => r.shape(), + DynReader::BioFormatsRust(r) => r.raw_shape(), #[cfg(feature = "bioformats_java")] - DynReader::BioFormatsJava(r) => r.shape(), + DynReader::BioFormatsJava(r) => r.raw_shape(), #[expect(unreachable_patterns)] _ => unreachable!(), } diff --git a/ndbioimage/src/readers/bioformats_java.rs b/ndbioimage/src/readers/bioformats_java.rs index b8b47cb..b7fc7a6 100644 --- a/ndbioimage/src/readers/bioformats_java.rs +++ b/ndbioimage/src/readers/bioformats_java.rs @@ -38,7 +38,7 @@ static JVM_PID: AtomicU32 = AtomicU32::new(0); /// else -- verified by `JNI_GetCreatedJavaVMs`, which still reports 1 in the child, so a second /// `JNI_CreateJavaVM` (which is what j4rs' `Jvm::attach_thread` calls) takes the already-exists /// path. An `exec` replaces the address space, and that does reset the count to 0. -fn check_not_forked() -> Result<(), Error> { +pub fn check_not_forked() -> Result<(), Error> { let built_in = JVM_PID.load(Ordering::Relaxed); if built_in != 0 && built_in != std::process::id() { return Err(Error::ForkedJvm(built_in, std::process::id())); @@ -698,7 +698,7 @@ impl Reader for BioFormatsJavaReader { 0 } - fn shape(&self) -> &Shape { + fn raw_shape(&self) -> &Shape { &self.shape } diff --git a/ndbioimage/src/readers/bioformats_rust.rs b/ndbioimage/src/readers/bioformats_rust.rs index 878bee0..a23290d 100644 --- a/ndbioimage/src/readers/bioformats_rust.rs +++ b/ndbioimage/src/readers/bioformats_rust.rs @@ -309,7 +309,7 @@ impl Reader for BioFormatsRustReader { 0 } - fn shape(&self) -> &Shape { + fn raw_shape(&self) -> &Shape { &self.shape } diff --git a/ndbioimage/src/readers/czi.rs b/ndbioimage/src/readers/czi.rs index 13c6467..717e510 100644 --- a/ndbioimage/src/readers/czi.rs +++ b/ndbioimage/src/readers/czi.rs @@ -1149,7 +1149,7 @@ impl Reader for CziReader { self.position } - fn shape(&self) -> &Shape { + fn raw_shape(&self) -> &Shape { &self.shape } diff --git a/ndbioimage/src/readers/tiff.rs b/ndbioimage/src/readers/tiff.rs index 110610c..8cd9797 100644 --- a/ndbioimage/src/readers/tiff.rs +++ b/ndbioimage/src/readers/tiff.rs @@ -459,7 +459,7 @@ impl Reader for TiffReader { fn position(&self) -> usize { self.position } - fn shape(&self) -> &Shape { + fn raw_shape(&self) -> &Shape { &self.shape } fn pixel_type(&self) -> &PixelType { diff --git a/ndbioimage/src/readers/tiffseq.rs b/ndbioimage/src/readers/tiffseq.rs index a392ed3..03faddb 100644 --- a/ndbioimage/src/readers/tiffseq.rs +++ b/ndbioimage/src/readers/tiffseq.rs @@ -572,7 +572,7 @@ impl Reader for TiffSeqReader { self.position } - fn shape(&self) -> &Shape { + fn raw_shape(&self) -> &Shape { &self.shape } diff --git a/ndbioimage/src/transforms.rs b/ndbioimage/src/transforms.rs index cf6e975..97779fb 100644 --- a/ndbioimage/src/transforms.rs +++ b/ndbioimage/src/transforms.rs @@ -1,13 +1,16 @@ use crate::axes::Axis; use crate::error::Error; -use crate::readers::{DynReader, Reader}; +use crate::readers::{DynReader, Frame, PixelType, Reader}; use crate::view::View; pub use image_registration::transform::Transform; -use ndarray::{Array2, Dimension, Ix2, Ix3, IxDyn, s}; +use ndarray::{Array2, Dimension, Ix2, Ix3, IxDyn, array, s}; use serde::{Deserialize, Serialize}; +use std::borrow::Cow; use std::collections::{HashMap, HashSet}; +use std::hash::{Hash, Hasher}; use std::io::Seek; use std::path::Path; +use std::sync::Arc; #[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)] struct TransformOld { @@ -54,6 +57,24 @@ pub struct Transforms { pub drift: Vec>, } +impl Hash for Transforms { + fn hash(&self, state: &mut H) { + let mut keys: Vec<_> = self.channel.keys().collect(); + keys.sort(); + for key in keys { + key.hash(state); + self.channel[key].hash(state); + } + self.drift.hash(state); + } +} + +impl PartialEq for &Transforms { + fn eq(&self, other: &Transforms) -> bool { + **self == *other + } +} + impl Serialize for Transforms { fn serialize(&self, serializer: S) -> Result where @@ -89,7 +110,23 @@ impl Serialize for Transforms { } } +pub enum TransformKind { + None, + Channel, + Drift, + Both, +} + impl Transforms { + pub fn kind(&self) -> TransformKind { + match (self.channel.is_empty(), self.drift.is_empty()) { + (true, true) => TransformKind::None, + (false, true) => TransformKind::Channel, + (true, false) => TransformKind::Drift, + (false, false) => TransformKind::Both, + } + } + pub fn load>(path: P) -> Result { let mut file = std::fs::File::open(&path)?; Ok( @@ -210,16 +247,32 @@ impl Transforms { } Ok(average_transforms) } + + pub fn get_transform( + &self, + channel: &str, + t: usize, + ) -> Result>>, Error> { + match (self.channel.get(channel), self.drift.get(t)) { + (Some(channel), Some(drift)) => match (channel.ndim(), drift.ndim()) { + (2, 2) | (3, 3) => Ok(Some(Cow::Owned(drift * channel))), + (2, 3) => Ok(Some(Cow::Owned(drift.add_dimension(0, 1) * channel))), + (3, 2) => Ok(Some(Cow::Owned(drift * channel.add_dimension(0, 1)))), + (2 | 3, n) | (n, _) => Err(Error::TransformationNotImplemented(n)), + }, + (Some(channel), None) => Ok(Some(Cow::Borrowed(channel))), + (None, Some(drift)) => Ok(Some(Cow::Borrowed(drift))), + (None, None) => Ok(None), + } + } } impl View { - #[cfg(feature = "transforms")] - pub(crate) fn get_transforms(&self) -> &Transforms { + pub fn get_transforms(&self) -> &Transforms { &self.transforms } - #[cfg(feature = "transforms")] - pub(crate) fn set_transforms(&mut self, transforms: Transforms) { + pub fn set_transforms(&mut self, transforms: Transforms) { self.transforms = transforms; } @@ -318,6 +371,181 @@ impl View { pub fn calculate_drift_transforms_3d(&self) -> Result>, Error> { todo!() } + + pub(crate) fn get_transformed_frame( + &self, + c: usize, + z: usize, + t: usize, + transform: &Transform, + order: usize, + ) -> Result, Error> { + match transform.ndim() { + 2 => { + let tr = transform.clone().into_dimensionality::()?; + + macro_rules! interpolate { + ($array:expr ,$var:ident, bool) => { + Ok(Arc::from(Frame::$var( + tr.interpolate_par_order( + &$array.mapv(|i| if i { 1.0 } else { 0.0 }), + order, + )? + .mapv(|i| i > 0.5), + ))) + }; + ($array:expr ,$var:ident, f64) => { + Ok(Arc::from(Frame::$var( + tr.interpolate_par_order($array, order)?, + ))) + }; + ($array:expr ,$var:ident, $t:ty) => { + Ok(Arc::from(Frame::$var( + tr.interpolate_par_order($array, order)?.mapv(|i| i as $t), + ))) + }; + } + + match self.get_cached_frame(c, z, t, None)?.as_ref() { + Frame::Bool(frame) => interpolate!(frame, Bool, bool), + Frame::I8(frame) => interpolate!(frame, I8, i8), + Frame::U8(frame) => interpolate!(frame, U8, u8), + Frame::I16(frame) => interpolate!(frame, I16, i16), + Frame::U16(frame) => interpolate!(frame, U16, u16), + Frame::I32(frame) => interpolate!(frame, I32, i32), + Frame::U32(frame) => interpolate!(frame, U32, u32), + Frame::U64(frame) => interpolate!(frame, U64, u64), + Frame::I64(frame) => interpolate!(frame, I64, i64), + Frame::I128(frame) => interpolate!(frame, I128, i128), + Frame::U128(frame) => interpolate!(frame, U128, u128), + Frame::F32(frame) => interpolate!(frame, F32, f32), + Frame::F64(frame) => interpolate!(frame, F64, f64), + Frame::F128(frame) => interpolate!(frame, F128, f64), + } + } + 3 => { + let tr = transform.clone().into_dimensionality::()?; + let shape = self.shape(); + let zf = z as f64; + let yf = shape.y as f64 - 1.0; + let xf = shape.x as f64 - 1.0; + let corners = array![[zf, 0.0, 0.0], [zf, 0.0, xf], [zf, yf, 0.0], [zf, yf, xf]]; + let extent = tr.transform_points(&corners)?; + let z_extent = extent.slice(s![.., 2]); + let z_min = (z_extent + .iter() + .min_by(|a, b| a.total_cmp(b)) + .unwrap() + .floor() as isize + - 2) + .max(0) as usize; + let z_max = (z_extent + .iter() + .max_by(|a, b| a.total_cmp(b)) + .unwrap() + .ceil() as isize + + 2) + .min(shape.z as isize - 1) + .max(0) as usize; + let frames = (z_min..=z_max) + .map(|zi| self.get_cached_frame(c, zi, t, None)) + .collect::, Error>>()?; + + // TODO; check these translations + let u0 = tr.as_translation(&[(z_min as f64), 0.0, 0.0]); + let u1 = tr.as_translation(&[-(z_min as f64), 0.0, 0.0]); + let tr = u0 * tr * u1; + + match frames.first().map(|f| f.as_ref()) { + Some(first) => { + macro_rules! interpolate { + ($var:ident, bool) => {{ + let views = frames + .iter() + .map(|f| match f.as_ref() { + Frame::$var(a) => Ok(a.view()), + _ => unreachable!(), + }) + .collect::, Error>>()?; + let stack = ndarray::stack(ndarray::Axis(0), &views)?; + let transformed = tr.interpolate_par_order( + &stack.mapv(|i| if i { 1.0 } else { 0.0 }), + order, + )?; + Ok(Arc::from(Frame::$var( + transformed.slice(s![0, .., ..]).mapv(|i| i > 0.5), + ))) + }}; + ($var:ident, f64) => {{ + let views = frames + .iter() + .map(|f| match f.as_ref() { + Frame::$var(a) => Ok(a.view()), + _ => unreachable!(), + }) + .collect::, Error>>()?; + let stack = ndarray::stack(ndarray::Axis(0), &views)?; + let transformed = tr.interpolate_par_order(&stack, order)?; + Ok(Arc::from(Frame::$var( + transformed.slice(s![0, .., ..]).to_owned(), + ))) + }}; + ($var:ident, $t:ty) => {{ + let views = frames + .iter() + .map(|f| match f.as_ref() { + Frame::$var(a) => Ok(a.view()), + _ => unreachable!(), + }) + .collect::, Error>>()?; + let stack = ndarray::stack(ndarray::Axis(0), &views)?; + let transformed = tr.interpolate_par_order(&stack, order)?; + Ok(Arc::from(Frame::$var( + transformed.slice(s![0, .., ..]).mapv(|i| i as $t), + ))) + }}; + } + match first { + Frame::Bool(_) => interpolate!(Bool, bool), + Frame::I8(_) => interpolate!(I8, i8), + Frame::U8(_) => interpolate!(U8, u8), + Frame::I16(_) => interpolate!(I16, i16), + Frame::U16(_) => interpolate!(U16, u16), + Frame::I32(_) => interpolate!(I32, i32), + Frame::U32(_) => interpolate!(U32, u32), + Frame::F32(_) => interpolate!(F32, f32), + Frame::F64(_) => interpolate!(F64, f64), + Frame::I64(_) => interpolate!(I64, i64), + Frame::U64(_) => interpolate!(U64, u64), + Frame::I128(_) => interpolate!(I128, i128), + Frame::U128(_) => interpolate!(U128, u128), + Frame::F128(_) => interpolate!(F128, f64), + } + } + None => { + let shape = self.shape(); + Ok(Arc::new(match self.pixel_type() { + PixelType::Bool => Frame::Bool(Array2::default((shape.y, shape.x))), + PixelType::I8 => Frame::I8(Array2::zeros((shape.y, shape.x))), + PixelType::U8 => Frame::U8(Array2::zeros((shape.y, shape.x))), + PixelType::I16 => Frame::I16(Array2::zeros((shape.y, shape.x))), + PixelType::U16 => Frame::U16(Array2::zeros((shape.y, shape.x))), + PixelType::I32 => Frame::I32(Array2::zeros((shape.y, shape.x))), + PixelType::U32 => Frame::U32(Array2::zeros((shape.y, shape.x))), + PixelType::F32 => Frame::F32(Array2::zeros((shape.y, shape.x))), + PixelType::F64 => Frame::F64(Array2::zeros((shape.y, shape.x))), + PixelType::I64 => Frame::I64(Array2::zeros((shape.y, shape.x))), + PixelType::U64 => Frame::U64(Array2::zeros((shape.y, shape.x))), + PixelType::I128 => Frame::I128(Array2::zeros((shape.y, shape.x))), + PixelType::U128 => Frame::U128(Array2::zeros((shape.y, shape.x))), + PixelType::F128 => Frame::F128(Array2::zeros((shape.y, shape.x))), + })) + } + } + } + d => Err(Error::TransformationNotImplemented(d)), + } + } } #[cfg(test)] diff --git a/ndbioimage/src/view.rs b/ndbioimage/src/view.rs index eb0ebb3..b7ddb74 100644 --- a/ndbioimage/src/view.rs +++ b/ndbioimage/src/view.rs @@ -5,18 +5,23 @@ use crate::metadata::Metadata; use crate::readers::{Dimensions, DynReader, Frame, Reader}; use crate::stats::MinMax; #[cfg(feature = "transforms")] -use crate::transforms::Transforms; +use crate::transforms::{TransformKind, Transforms}; +#[cfg(feature = "transforms")] +use image_registration::transform::Transform; use indexmap::IndexMap; use itertools::{Itertools, iproduct}; use ndarray::{ Array, Array0, Array1, Array2, ArrayD, Dimension, IntoDimension, Ix0, Ix1, Ix2, Ix5, IxDyn, SliceArg, SliceInfoElem, s, }; -use num::traits::ToBytes; -use num::{Bounded, FromPrimitive, ToPrimitive, Zero}; +use num::traits::{AsPrimitive, ToBytes}; +use num::{Bounded, FromPrimitive, One, ToPrimitive, Zero}; use serde::{Deserialize, Serialize}; use serde_with::serde_as; use std::any::type_name; +#[cfg(feature = "transforms")] +use std::borrow::Cow; +use std::cell::OnceCell; use std::collections::{HashMap, HashSet}; use std::fmt::{Debug, Display, Formatter}; use std::hash::{Hash, Hasher}; @@ -57,27 +62,43 @@ pub trait Number: + Sync + AddAssign + Bounded - + Clone + + Copy + Div + FromPrimitive + PartialOrd + Zero + + One + + AsPrimitive { + const ISFLOAT: bool; } -impl Number for T where - T: 'static - + Send - + Sync - + AddAssign - + Bounded - + Clone - + Div - + FromPrimitive - + PartialOrd - + Zero -{ + +macro_rules! impl_int { + ($($t:tt),+ $(,)?) => { + $( + impl Number for $t { + const ISFLOAT: bool = false; + } + )* + }; } +impl_int!( + i8, i16, i32, i64, i128, isize, u8, u16, u32, u64, u128, usize +); + +macro_rules! impl_float { + ($($t:tt),+ $(,)?) => { + $( + impl Number for $t { + const ISFLOAT: bool = true; + } + )* + }; +} + +impl_float!(f32, f64); + /// sliceable view on an image file #[serde_as] #[derive(Clone, Debug, PartialEq, Eq, Serialize, Deserialize)] @@ -89,6 +110,8 @@ pub struct View { /// always has all of cztyx with possibly some new axes added axes: Vec, operations: IndexMap, + #[serde(skip)] + channel_names: OnceCell>, dimensionality: PhantomData, #[cfg(feature = "transforms")] pub(crate) transforms: Transforms, @@ -117,6 +140,7 @@ impl View { slice, axes, operations: IndexMap::new(), + channel_names: OnceCell::new(), dimensionality: PhantomData, #[cfg(feature = "transforms")] transforms: Transforms::default(), @@ -126,7 +150,7 @@ impl View { #[expect(dead_code)] pub(crate) fn new_with_axes(reader: R, axes: Vec) -> Result { let mut slice = Vec::new(); - let shape = reader.shape(); + let shape = reader.raw_shape(); for axis in axes.iter() { match axis { Axis::C => slice.push(SliceInfoElem::Slice { @@ -183,6 +207,7 @@ impl View { axes, operations: IndexMap::new(), dimensionality: PhantomData, + channel_names: OnceCell::new(), #[cfg(feature = "transforms")] transforms: Transforms::default(), }) @@ -198,6 +223,26 @@ impl View { self.reader.series() } + /// the name of a particular channel + pub fn channel_name(&self, channel: usize) -> Result<&str, Error> { + if self.channel_names.get().is_none() { + let names = self.metadata()?.image[0] + .pixels + .channel + .iter() + .enumerate() + .map(|(i, c)| (i, c.name.clone().unwrap_or_else(|| format!("{}", i)))) + .collect::>(); + let _ = self.channel_names.set(names); + } + + self.channel_names + .get() + .and_then(|names| names.get(&channel)) + .map(String::as_str) + .ok_or(Error::NoChannel(channel)) + } + /// set the maximum number of frames the process-wide cache may hold, /// evicting the least recently used frames pub fn with_cache_capacity(self, capacity: usize) -> Self { @@ -234,6 +279,7 @@ impl View { slice: self.slice, axes: self.axes, operations: self.operations, + channel_names: self.channel_names, dimensionality: PhantomData, #[cfg(feature = "transforms")] transforms: Transforms::default(), @@ -249,6 +295,7 @@ impl View { slice: self.slice, axes: self.axes, operations: self.operations, + channel_names: self.channel_names, dimensionality: PhantomData, #[cfg(feature = "transforms")] transforms: self.transforms, @@ -262,6 +309,7 @@ impl View { slice: self.slice, axes: self.axes, operations: self.operations, + channel_names: self.channel_names, dimensionality: PhantomData, #[cfg(feature = "transforms")] transforms: self.transforms, @@ -671,7 +719,7 @@ impl View { 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()) + Ok(*arr.first().unwrap()) } /// collect the view into an ndarray @@ -714,7 +762,7 @@ impl View { IndexMap::new() }; let mut shape_out = Vec::new(); - let mut slice = Vec::new(); + let mut slice_out = Vec::new(); let mut ax_out = Vec::new(); for (s, a) in self.slice.iter().zip(&self.axes) { match s { @@ -722,7 +770,7 @@ impl View { let end = end.expect("slice has no end"); if !op_xy.contains_key(a) && !op_czt.contains_key(a) { shape_out.push(((end - start).max(0) / step) as usize); - slice.push(SliceInfoElem::Slice { + slice_out.push(SliceInfoElem::Slice { start: 0, end: None, step: 1, @@ -733,59 +781,54 @@ impl View { SliceInfoElem::Index(_) => {} SliceInfoElem::NewAxis => { shape_out.push(1); - slice.push(SliceInfoElem::Index(0)); + slice_out.push(SliceInfoElem::Index(0)); ax_out.push(*a); } } } let mut slice_reader = vec![Slice::empty(); 5]; let mut xy_dim = 0usize; - let shape_reader = self.reader.shape(); + let shape_reader = self.reader.raw_shape(); 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 - && !op_xy.contains_key(&axis) - { - xy_dim += 1; - } - slice_reader[axis as usize] = Slice::new( - idx_bnd(*start, shape_reader[axis] as isize)?, - slc_bnd(end.unwrap(), shape_reader[axis] as isize)?, - *step, - ); + match (axis, s) { + (Axis::New, _) => {} + (_, SliceInfoElem::Slice { start, end, step }) => { + if let Axis::X | Axis::Y = axis + && !op_xy.contains_key(&axis) + { + xy_dim += 1; } - SliceInfoElem::Index(j) => { - slice_reader[axis as usize] = Slice::new( - idx_bnd(*j, shape_reader[axis] as isize)?, - slc_bnd(*j + 1, shape_reader[axis] as isize)?, - 1, - ); - } - SliceInfoElem::NewAxis => panic!("axis cannot be a new axis"), - }, + slice_reader[axis as usize] = Slice::new( + idx_bnd(*start, shape_reader[axis] as isize)?, + slc_bnd(end.unwrap(), shape_reader[axis] as isize)?, + *step, + ); + } + (_, SliceInfoElem::Index(j)) => { + slice_reader[axis as usize] = Slice::new( + idx_bnd(*j, shape_reader[axis] as isize)?, + slc_bnd(*j + 1, shape_reader[axis] as isize)?, + 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()), + let mut array = match op_czt.first() { + Some((_, Operation::Max)) => { + ArrayD::::from_elem(shape_out.into_dimension(), T::min_value()) } - } else { - ArrayD::::zeros(shape_out.into_dimension()) + Some((_, Operation::Min)) => { + ArrayD::::from_elem(shape_out.into_dimension(), T::max_value()) + } + _ => ArrayD::::zeros(shape_out.into_dimension()), }; - let shape = self.reader.shape(); + + let shape = self.reader.raw_shape(); let mut axes_out_idx = [None; 5]; for (i, ax) in ax_out.iter().enumerate() { if *ax < Axis::New { @@ -793,22 +836,56 @@ impl View { } } + #[cfg(feature = "transforms")] + let kind = self.transforms.kind(); + #[cfg(feature = "transforms")] + let order = if T::ISFLOAT { 3 } else { 0 }; + #[cfg(feature = "transforms")] + let mut transform_done = HashSet::new(); + 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) + slice_out[i] = SliceInfoElem::Index(c) }; if let Some(i) = axes_out_idx[1] { - slice[i] = SliceInfoElem::Index(z) + slice_out[i] = SliceInfoElem::Index(z) }; if let Some(i) = axes_out_idx[2] { - slice[i] = SliceInfoElem::Index(t) + slice_out[i] = SliceInfoElem::Index(t) }; + #[cfg(feature = "transforms")] + let transform = self + .transforms + .get_transform(self.channel_name(c as usize)?, t as usize)?; + #[cfg(feature = "transforms")] + let transform_now = if let (TransformKind::Channel, Some((&Axis::C, _)), _) + | (TransformKind::Drift, Some((&Axis::T, _)), _) + | (TransformKind::Both, Some((&Axis::C | &Axis::T, _)), _) + | ( + TransformKind::Channel | TransformKind::Drift | TransformKind::Both, + _, + false, + ) = (&kind, op_czt.first(), xy_dim == 2) + { + transform_done.insert((c, z, t)); + transform.map(|transform| (transform, order)) + } else { + None + }; + #[cfg(feature = "transforms")] + let frame = self.get_cached_frame( + (c % shape.c as isize) as usize, + (z % shape.z as isize) as usize, + (t % shape.t as isize) as usize, + transform_now, + )?; + #[cfg(not(feature = "transforms"))] let frame = self.get_cached_frame( (c % shape.c as isize) as usize, (z % shape.z as isize) as usize, (t % shape.t as isize) as usize, )?; - + // 0, 1 or 2d frame depending on x or y operations let arr_frame: Array2 = frame.as_ref().try_into()?; let arr_frame = match xy_dim { 0 => { @@ -855,35 +932,28 @@ impl View { unreachable!("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()); - } + // assign arr_frame to output array, apply first of c, z or t operation + match op_czt.first() { + Some((_, Operation::Max)) => { + array + .slice_mut(slice_out.as_slice()) + .zip_mut_with(&arr_frame, |x, y| { + *x = if *x >= *y { *x } else { *y } + }); } - } else { - array.slice_mut(slice.as_slice()).assign(&arr_frame) + Some((_, Operation::Min)) => { + array + .slice_mut(slice_out.as_slice()) + .zip_mut_with(&arr_frame, |x, y| { + *x = if *x < *y { *x } else { *y } + }); + } + Some((_, Operation::Mean | Operation::Sum)) => { + array + .slice_mut(slice_out.as_slice()) + .zip_mut_with(&arr_frame, |x, y| *x += *y); + } + _ => array.slice_mut(slice_out.as_slice()).assign(&arr_frame), } } let mut out = Some(array); @@ -894,23 +964,61 @@ impl View { .collect(); for (ax, op) in self.operations.iter().skip(op_xy.len() + op_czt.len()) { if let Some(idx) = ax_out.remove(ax) { + #[cfg(feature = "transforms")] + let mut arr = out.take().unwrap(); + #[cfg(not(feature = "transforms"))] + let arr = out.take().unwrap(); + + #[cfg(feature = "transforms")] + if let (Axis::C, TransformKind::Channel | TransformKind::Both) + | (Axis::T, TransformKind::Drift | TransformKind::Both) = (ax, &kind) + { + let mut slice_transform = vec![ + SliceInfoElem::Slice { + start: 0, + end: None, + step: 1 + }; + ax_out.len() + ]; + for (c, z, t) in iproduct!(&slice_reader[0], &slice_reader[1], &slice_reader[2]) + { + if let Some(transform) = self + .transforms + .get_transform(self.channel_name(c as usize)?, t as usize)? + && !transform_done.contains(&(c, z, t)) + { + for (axis, i) in [&Axis::C, &Axis::Z, &Axis::T].iter().zip([c, z, t]) { + if let Some(axis_pos) = ax_out.get(axis) { + slice_transform[*axis_pos] = SliceInfoElem::Index(i); + } + } + let untransformed = arr.slice(slice_transform.as_slice()); + let transformed = transform + .interpolate_par_order(untransformed, order)? + .mapv(|i| T::from_f64(i).unwrap()); + arr.slice_mut(slice_transform.as_slice()) + .assign(&transformed); + } + transform_done.insert((c, z, t)); + } + } + + let a = op.operate(arr, idx)?; + let _ = out.insert(a); + for i in ax_out.values_mut() { if *i > idx { *i -= 1; } } - let arr = out.take().unwrap(); - let a = op.operate(arr, idx)?; - let _ = out.insert(a); } } - let n = if let Some((&ax, op)) = op_czt.first() - && *op == Operation::Mean - { + let n = if let Some((ax, Operation::Mean)) = op_czt.first() { self.axes .iter() .zip(self.slice.iter()) - .find(|(a, _)| **a == ax) + .find(|(a, _)| *a == ax) .and_then(|(_, s)| match s { SliceInfoElem::Slice { start, end, step } => { end.map(|e| (((e - start).max(0) / step) as usize).max(1)) @@ -921,12 +1029,29 @@ impl View { } else { 1 }; + #[cfg(not(feature = "transforms"))] 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()) + out.take().unwrap().mapv(|x| x / m) }; + #[cfg(feature = "transforms")] + let mut 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) + }; + #[cfg(feature = "transforms")] + if transform_done.is_empty() { + let c = slice_reader[0].start; + let t = slice_reader[2].start; + if let Some(transform) = self.transforms.get_transform(self.channel_name(c as usize)?, t as usize)? { + array = transform.interpolate_par_order(&array, order)?.mapv(|i| T::from_f64(i).unwrap()); + } + } + ArrayCache::global().insert(key.to_owned(), array.clone()); Ok(array) } @@ -957,7 +1082,13 @@ impl View { .collect()) } - fn frame_key(&self, c: usize, z: usize, t: usize) -> FrameKeyRef<'_> { + fn frame_key<'a>( + &'a self, + c: usize, + z: usize, + t: usize, + #[cfg(feature = "transforms")] transform: Option<(Cow<'a, Transform>, usize)>, + ) -> FrameKeyRef<'a> { FrameKeyRef { reader: ReaderKeyRef { name: self.reader.reader_name(), @@ -969,6 +1100,8 @@ impl View { c, z, t, + #[cfg(feature = "transforms")] + transform, } } @@ -985,14 +1118,37 @@ impl View { slice: &self.slice, axes: &self.axes, operations: &self.operations, + #[cfg(feature = "transforms")] + transforms: &self.transforms, } } - fn get_cached_frame(&self, c: usize, z: usize, t: usize) -> Result, Error> { - let key = self.frame_key(c, z, t); + pub(crate) fn get_cached_frame( + &self, + c: usize, + z: usize, + t: usize, + #[cfg(feature = "transforms")] transform: Option<(Cow>, usize)>, + ) -> Result, Error> { + let key = self.frame_key( + c, + z, + t, + #[cfg(feature = "transforms")] + transform + .as_ref() + .map(|(transform, order)| (Cow::Borrowed(&**transform), *order)), + ); if let Some(frame) = FrameCache::global().get(&key) { return Ok(frame); } + #[cfg(feature = "transforms")] + let frame = if let Some((transform, order)) = transform.as_ref() { + self.get_transformed_frame(c, z, t, transform, *order)? + } else { + Arc::new(self.reader.get_frame(c, z, t)?) + }; + #[cfg(not(feature = "transforms"))] let frame = Arc::new(self.reader.get_frame(c, z, t)?); FrameCache::global().insert(key.to_owned(), frame.clone()); Ok(frame) @@ -1187,7 +1343,7 @@ impl Item for View { Array1: MinMax>, Array2: MinMax>, { - Ok(self.as_array()?.first().ok_or(Error::EmptyView)?.clone()) + Ok(*self.as_array()?.first().ok_or(Error::EmptyView)?) } } @@ -1229,7 +1385,7 @@ to_bytes_vec_impl!( mod tests { use crate::axes::{Axis, Operation}; use crate::cache::{ - ArrayCache, ArrayKey, ArrayKeyRef, FileStamp, FrameCache, ReaderKey, ReaderKeyRef, + ArrayCache, ArrayKey, ArrayKeyRef, FileStamp, FrameCache, FrameKey, ReaderKey, ReaderKeyRef, }; use crate::error::Error; use crate::readers::{DynReader, Frame, Reader}; @@ -1241,6 +1397,8 @@ mod tests { use postcard::{from_bytes, to_stdvec}; use std::path::{Path, PathBuf}; use std::sync::Arc; + #[cfg(feature = "transforms")] + use crate::transforms::Transforms; fn open(file: &str) -> Result { let path = std::env::current_dir()? @@ -1512,19 +1670,19 @@ mod tests { #[test] fn frame_cache() -> Result<(), Error> { let cache = FrameCache::new(2); - let key = |c: usize, z: usize, t: usize| { - ( - ReaderKey { - name: "test".to_string(), - path: PathBuf::from("test.tif"), - series: 0, - position: 0, - stamp: FileStamp::default(), - }, - c, - z, - t, - ) + let key = |c: usize, z: usize, t: usize| FrameKey { + reader_key: ReaderKey { + name: "test".to_string(), + path: PathBuf::from("test.tif"), + series: 0, + position: 0, + stamp: FileStamp::default(), + }, + c, + z, + t, + #[cfg(feature = "transforms")] + transform: None, }; cache.insert( key(0, 0, 0), @@ -1561,6 +1719,8 @@ mod tests { slice: vec![SliceInfoElem::Index(index)], axes: axes.to_vec(), operations: vec![], + #[cfg(feature = "transforms")] + transforms: Transforms::default(), }; let k0 = key("u16", 0, &[Axis::T]); let k1 = key("u16", 1, &[Axis::T]); @@ -1587,6 +1747,8 @@ mod tests { slice: &[SliceInfoElem::Index(0)], axes: &[Axis::T], operations: &ops, + #[cfg(feature = "transforms")] + transforms: &Transforms::default(), }; assert!(cache.get::(&borrowed).is_none()); assert!(cache.get::(&borrowed).is_some()); @@ -1602,12 +1764,37 @@ mod tests { let (c, z, t) = (shape[0] - 1, shape[1] - 1, shape[2] - 1); let cache = FrameCache::global(); let a: Array5 = view.clone().try_into()?; - assert!(cache.get(&view.frame_key(c, z, t)).is_some()); + assert!( + cache + .get(&view.frame_key( + c, + z, + t, + #[cfg(feature = "transforms")] + None + )) + .is_some() + ); let b: Array5 = view.clone().try_into()?; assert_eq!(a, b); // a fresh reader on the same file shares the same cache key let view2 = open(file)?.view(); - assert_eq!(view.frame_key(c, z, t), view2.frame_key(c, z, t)); + assert_eq!( + view.frame_key( + c, + z, + t, + #[cfg(feature = "transforms")] + None + ), + view2.frame_key( + c, + z, + t, + #[cfg(feature = "transforms")] + None + ) + ); let c2: Array5 = view2.clone().try_into()?; assert_eq!(a, c2); Ok(())