From e1caa1c0b0c459c9ae8c84b53278e8f7449e3f60 Mon Sep 17 00:00:00 2001 From: "w.pomp" Date: Wed, 21 Oct 2026 20:08:42 +0200 Subject: [PATCH] - fix Imread.transpose - implement with_transform - implement 2d transforms in View.as_array_dyn --- image-registration | 2 +- .../py/ndbioimage/ndbioimage_rs/__init__.pyi | 128 +++--- .../ndbioimage/ndbioimage_rs/transforms.pyi | 20 + ndbioimage/src/axes.rs | 4 +- ndbioimage/src/cache.rs | 20 +- ndbioimage/src/error.rs | 11 +- ndbioimage/src/py/imread.rs | 91 ++++- ndbioimage/src/py/transforms.rs | 107 +++++- ndbioimage/src/readers.rs | 2 +- ndbioimage/src/readers/czi.rs | 12 +- ndbioimage/src/transforms.rs | 143 +++++++ ndbioimage/src/view.rs | 363 ++++++++++++------ 12 files changed, 689 insertions(+), 214 deletions(-) diff --git a/image-registration b/image-registration index 5c580dc..dc1d5d8 160000 --- a/image-registration +++ b/image-registration @@ -1 +1 @@ -Subproject commit 5c580dc80f8552b91a5c2dae4fe1f1f8e3142415 +Subproject commit dc1d5d8364f51f7323ae3aaac4490b9eece6305d diff --git a/ndbioimage/py/ndbioimage/ndbioimage_rs/__init__.pyi b/ndbioimage/py/ndbioimage/ndbioimage_rs/__init__.pyi index 0f3ba00..f8b4cfe 100644 --- a/ndbioimage/py/ndbioimage/ndbioimage_rs/__init__.pyi +++ b/ndbioimage/py/ndbioimage/ndbioimage_rs/__init__.pyi @@ -96,16 +96,6 @@ class Imread: the name of the reader used to open the file """ @property - def transforms(self) -> transforms.Transforms: - r""" - get the transformation - """ - @transforms.setter - def transforms(self, value: transforms.Transforms) -> None: - r""" - set the transformation - """ - @property def path(self) -> pathlib.Path: r""" the file path @@ -220,6 +210,28 @@ class Imread: r""" the name of the tube lens """ + @property + def interpolation_order(self) -> typing.Optional[builtins.int]: + r""" + the bspline order used for interpolation when using transforms + None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes + """ + @interpolation_order.setter + def interpolation_order(self, value: typing.Optional[builtins.int]) -> None: + r""" + the bspline order used for interpolation when using transforms + None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes + """ + @property + def transforms(self) -> transforms.Transforms: + r""" + get the transformation + """ + @transforms.setter + def transforms(self, value: transforms.Transforms) -> None: + r""" + set the transformation + """ def __new__( cls, path: str | pathlib.Path | Imread | bytes, @@ -236,6 +248,10 @@ class Imread: reader: use a specific reader to open the file, view available choices with Imread.available_readers """ + def debug(self) -> None: + r""" + some debug info + """ @staticmethod def warm_jvm() -> None: r""" @@ -247,7 +263,10 @@ class Imread: @staticmethod def jvm_is_running() -> builtins.bool: r""" - returns whether the jvm was already created on this thread + returns whether the jvm was already created in this process + + this is false in a process that inherited a jvm across a fork, since such a jvm cannot be + used; reading a bio-formats file there raises instead """ def set_flags( self, @@ -274,42 +293,6 @@ class Imread: r""" reshape the view with a new axis order """ - def with_transform( - self, - channels: builtins.bool = True, - drift: builtins.bool = False, - file: typing.Optional[typing.Any] = None, - bead_files: typing.Optional[typing.Any] = None, - main_channel: typing.Optional[builtins.int] = None, - default_transform: typing.Optional[typing.Sequence[builtins.float]] = None, - ) -> Imread: - r""" - return a new view with transformations applied (channel alignment, drift correction) - backward compatibility - """ - def with_transforms( - self, - channels: builtins.bool = True, - drift: builtins.bool = False, - file: typing.Optional[typing.Any] = None, - bead_files: typing.Optional[typing.Any] = None, - main_channel: typing.Optional[builtins.int] = None, - default_transform: typing.Optional[typing.Sequence[builtins.float]] = None, - ) -> Imread: - r""" - return a new view with transformations applied (channel alignment, drift correction) - """ - def load_transforms_from_yaml( - self, path: builtins.str | os.PathLike | pathlib.Path - ) -> None: ... - def calculate_channel_transforms_2d( - self, main_channel: builtins.int - ) -> builtins.dict[builtins.str, transforms.Transform]: ... - def calculate_channel_transforms_3d( - self, main_channel: builtins.int - ) -> builtins.dict[builtins.str, transforms.Transform]: ... - def calculate_drift_transforms_2d(self) -> builtins.list[transforms.Transform]: ... - def calculate_drift_transforms_3d(self) -> builtins.list[transforms.Transform]: ... def squeeze(self) -> numpy.ndarray | int | float: ... def close(self) -> None: r""" @@ -551,7 +534,7 @@ class Imread: swap two axes """ def transpose( - self, axes: typing.Optional[typing.Sequence[int | str]] = None + self, axes: typing.Optional[typing.Sequence[int | str] | str] = None ) -> Imread: r""" permute the order of the axes @@ -586,6 +569,10 @@ class Imread: r""" the detector gain for a given channel """ + def channel_name(self, channel: builtins.int) -> builtins.str: + r""" + the name of a particular channel + """ def summary(self) -> builtins.str: r""" gives a helpful summary of the recorded experiment @@ -692,6 +679,51 @@ class Imread: r""" Return the sum along a given axis. Arguments beyond axis are not implemented """ + def with_interpolation_order(self, order: typing.Optional[builtins.int]) -> Imread: + r""" + the bspline order used for interpolation when using transforms + None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes + """ + def with_transform( + self, + channels: builtins.bool = True, + drift: builtins.bool = False, + file: typing.Optional[builtins.str | os.PathLike | pathlib.Path] = None, + bead_files: typing.Optional[ + typing.Sequence[builtins.str | os.PathLike | pathlib.Path] + ] = None, + main_channel: typing.Optional[builtins.int] = None, + default_transform: typing.Optional[typing.Sequence[builtins.float]] = None, + ) -> Imread: + r""" + return a new view with transformations applied (channel alignment, drift correction) + backward compatibility + """ + def with_transforms( + self, + channels: builtins.bool = True, + drift: builtins.bool = False, + file: typing.Optional[builtins.str | os.PathLike | pathlib.Path] = None, + bead_files: typing.Optional[ + typing.Sequence[builtins.str | os.PathLike | pathlib.Path] + ] = None, + main_channel: typing.Optional[builtins.int] = None, + default_transform: typing.Optional[typing.Sequence[builtins.float]] = None, + ) -> Imread: + r""" + returns a view where channels and/or frames are registered with an affine transformation + """ + def load_transforms_from_yaml( + self, path: builtins.str | os.PathLike | pathlib.Path + ) -> None: ... + def calculate_channel_transforms_2d( + self, main_channel: builtins.int + ) -> builtins.dict[builtins.str, transforms.Transform]: ... + def calculate_channel_transforms_3d( + self, main_channel: builtins.int + ) -> builtins.dict[builtins.str, transforms.Transform]: ... + def calculate_drift_transforms_2d(self) -> builtins.list[transforms.Transform]: ... + def calculate_drift_transforms_3d(self) -> builtins.list[transforms.Transform]: ... class Shape: r""" diff --git a/ndbioimage/py/ndbioimage/ndbioimage_rs/transforms.pyi b/ndbioimage/py/ndbioimage/ndbioimage_rs/transforms.pyi index 284c30a..706f971 100644 --- a/ndbioimage/py/ndbioimage/ndbioimage_rs/transforms.pyi +++ b/ndbioimage/py/ndbioimage/ndbioimage_rs/transforms.pyi @@ -225,6 +225,18 @@ class Transform: class Transforms: @property + def interpolation_order(self) -> typing.Optional[builtins.int]: + r""" + the bspline order used for interpolation when using transforms, + None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes + """ + @interpolation_order.setter + def interpolation_order(self, value: typing.Optional[builtins.int]) -> None: + r""" + the bspline order used for interpolation when using transforms + None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes + """ + @property def channel_transforms(self) -> builtins.dict[builtins.str, Transform]: ... @channel_transforms.setter def channel_transforms( @@ -239,7 +251,15 @@ class Transforms: cls, channel: typing.Mapping[builtins.str, Transform] = {}, drift: typing.Sequence[Transform] = [], + interpolation_order: typing.Optional[builtins.int] = None, ) -> Transforms: ... + def with_interpolation_order( + self, order: typing.Optional[builtins.int] + ) -> Transforms: + r""" + the bspline order used for interpolation when using transforms + None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes + """ def __getnewargs__( self, ) -> tuple[builtins.dict[builtins.str, Transform], builtins.list[Transform]]: ... diff --git a/ndbioimage/src/axes.rs b/ndbioimage/src/axes.rs index a374e2c..6c9c1be 100644 --- a/ndbioimage/src/axes.rs +++ b/ndbioimage/src/axes.rs @@ -196,7 +196,7 @@ impl<'de> DeserializeAs<'de, SliceInfoElem> for SliceInfoElemDef { } } -#[derive(Clone, Debug)] +#[derive(Clone, Debug, Copy)] pub(crate) struct Slice { pub(crate) start: isize, pub(crate) end: isize, @@ -236,7 +236,7 @@ impl IntoIterator for &Slice { type IntoIter = Slice; fn into_iter(self) -> Self::IntoIter { - self.clone() + *self } } diff --git a/ndbioimage/src/cache.rs b/ndbioimage/src/cache.rs index 0ea21b5..49620aa 100644 --- a/ndbioimage/src/cache.rs +++ b/ndbioimage/src/cache.rs @@ -1,6 +1,8 @@ use crate::axes::{Axis, Operation}; use crate::readers::Frame; #[cfg(feature = "transforms")] +use crate::transforms::Transforms; +#[cfg(feature = "transforms")] use image_registration::transform::Transform; use indexmap::{Equivalent, IndexMap}; #[cfg(feature = "transforms")] @@ -15,8 +17,6 @@ use std::hash::{Hash, Hasher}; use std::path::{Path, PathBuf}; use std::sync::{Arc, Mutex, OnceLock}; use std::time::SystemTime; -#[cfg(feature = "transforms")] -use crate::transforms::Transforms; /// maximum number of frames held in the cache const DEFAULT_FRAME_CACHE_SIZE: usize = 128; @@ -264,10 +264,10 @@ impl Equivalent for ArrayKeyRef<'_> { && 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) + .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"))] @@ -278,10 +278,10 @@ impl Equivalent for ArrayKeyRef<'_> { && 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) + .operations + .iter() + .zip(&key.operations) + .all(|((ax, op), (key_ax, key_op))| ax == key_ax && op == key_op) } } } diff --git a/ndbioimage/src/error.rs b/ndbioimage/src/error.rs index 5b548c3..a40a44b 100644 --- a/ndbioimage/src/error.rs +++ b/ndbioimage/src/error.rs @@ -65,10 +65,6 @@ pub enum Error { #[cfg(feature = "czi")] #[error(transparent)] XmlTreeParse(#[from] xmltree::ParseError), - /// a czi-specific error - #[cfg(feature = "czi")] - #[error(transparent)] - Czi(#[from] crate::readers::czi::CziError), /// an error joining a tokio task #[cfg(feature = "movie")] #[error(transparent)] @@ -173,16 +169,21 @@ pub enum Error { #[error("{0} do not match in {1} and {2}")] FileMismatch(String, String, String), /// a jvm created in another process was inherited across a fork - #[cfg(feature = "bioformats_java")] #[error( "the jvm was created in process {0} but this is process {1}: a forked child inherits a jvm \ it cannot use. read the file before forking, or start the child with spawn/exec instead." )] + /// jvm was forked, so it cannot run in this process ForkedJvm(u32, u32), + /// only 2d and 3d transformations are implemented #[error("{0}D transformation is not implemented")] TransformationNotImplemented(usize), + /// cannot find the name for some channel because it does not exist #[error("no channel {0} present")] NoChannel(usize), + /// no valid blocks in czi file + #[error("czi file has no valid blocks")] + CziNoValidBlocks, } impl Error { diff --git a/ndbioimage/src/py/imread.rs b/ndbioimage/src/py/imread.rs index b25176c..5bae44e 100644 --- a/ndbioimage/src/py/imread.rs +++ b/ndbioimage/src/py/imread.rs @@ -210,6 +210,16 @@ impl PyView { } } + /// some debug info + fn debug(&self) { + println!( + "slice: {:?}\naxes: {:?}\noperations: {:?}", + self.view.get_slice(), + self.view.get_axes(), + self.view.get_operations() + ) + } + /// creates the jvm (if it does not exist yet) so that hotspot installs its signal handlers before /// anything else does: see `crate::readers::bioformats_java::warm_jvm` /// @@ -1340,16 +1350,15 @@ impl PyView { #[pyo3(signature = (axes = None))] fn transpose( &self, - #[gen_stub(override_type(type_repr="typing.Optional[typing.Sequence[int | str]]", imports=("typing") + #[gen_stub(override_type(type_repr="typing.Optional[typing.Sequence[int | str] | str]", imports=("typing") ))] - axes: Option>>, + axes: Option>, ) -> PyResult { - // TODO: transpose("xyz"), transpose(*axes) let view = if let Some(axes) = axes { let ax = axes - .into_iter() - .map(|a| self.get_ax(a)) - .collect::, _>>()?; + .try_iter()? + .map(|a| self.get_ax(a?)) + .collect::>>()?; self.view.permute_axes(&ax)? } else { self.view.transpose()? @@ -1373,19 +1382,58 @@ impl PyView { fn as_array<'py>(&self, py: Python<'py>) -> PyResult> { Ok(match self.dtype { PixelType::Bool => todo!(), - 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(), + PixelType::I8 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), + PixelType::U8 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), + PixelType::I16 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), + PixelType::U16 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), + PixelType::I32 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), + PixelType::U32 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), + PixelType::F32 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), + PixelType::F64 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), + PixelType::I64 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), + PixelType::U64 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), + PixelType::I128 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), + PixelType::U128 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), + PixelType::F128 => py + .detach(|| self.view.as_array_dyn::())? + .into_pyarray(py) + .into_any(), }) } @@ -1773,6 +1821,11 @@ impl PyView { self.ome.gain(channel) } + /// the name of a particular channel + fn channel_name(&self, channel: usize) -> PyResult { + Ok(self.view.channel_name(channel).map(|s| s.to_owned())?) + } + /// gives a helpful summary of the recorded experiment fn summary(&self) -> PyResult { Ok(self.view.summary()?) diff --git a/ndbioimage/src/py/transforms.rs b/ndbioimage/src/py/transforms.rs index 824b0b0..58c9e04 100644 --- a/ndbioimage/src/py/transforms.rs +++ b/ndbioimage/src/py/transforms.rs @@ -1,5 +1,6 @@ use crate::error::Error; use crate::py::imread::PyView; +use crate::readers::Reader; use crate::transforms::Transforms; use image_registration::metric::{SamplingArg, Sigma}; use image_registration::register::{Optimizer, RegistrationStep}; @@ -713,16 +714,38 @@ pub(crate) struct PyTransforms { #[gen_stub_pymethods] #[pymethods] impl PyView { + /// the bspline order used for interpolation when using transforms + /// None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes + #[getter] + fn get_interpolation_order(&self) -> Option { + self.view.transforms.interpolation_order + } + + /// the bspline order used for interpolation when using transforms + /// None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes + #[setter] + fn set_interpolation_order(&mut self, order: Option) { + self.view.transforms.interpolation_order = order; + } + + /// the bspline order used for interpolation when using transforms + /// None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes + fn with_interpolation_order(&self, order: Option) -> Self { + let mut new = self.clone(); + new.view.transforms.interpolation_order = order; + new + } + /// 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>( + fn with_transform( &self, channels: bool, drift: bool, - file: Option>, - bead_files: Option>, + file: Option, + bead_files: Option>, main_channel: Option, default_transform: Option>, ) -> PyResult { @@ -736,20 +759,55 @@ impl PyView { ) } - /// return a new view with transformations applied (channel alignment, drift correction) + /// 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' + // main channel: which channel to transform with default_transform + // default_transform: parameters of the default transform: [a, b, c, d, e, f] where the affine transform matrix is + // | a b e | + // | c d f | + // | 0 0 1 | #[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>( + fn with_transforms( &self, channels: bool, drift: bool, - file: Option>, - bead_files: Option>, + file: Option, + bead_files: Option>, main_channel: Option, default_transform: Option>, ) -> PyResult { - todo!() + let mut new = self.clone(); + let bead_files = bead_files + .as_ref() + .map(|i| i.iter().map(|j| j.as_path()).collect::>()); + new.view.transforms = Transforms::calculate_transforms_from( + file.as_deref(), + if channels { + bead_files.as_deref() + } else { + None + }, + if drift { Some(self.view.path()) } else { None }, + main_channel, + default_transform, + )?; + let new_shape = vec![self.view.raw_shape().y, self.view.raw_shape().x]; + for transform in new.view.transforms.channel.values_mut() { + for (c, (old_s, new_s)) in transform + .center + .iter_mut() + .zip(transform.shape.iter().zip(new_shape.iter())) + { + *c -= (*old_s as f64) / 2.0 - (*new_s as f64) / 2.0; + } + transform.shape = new_shape.clone(); + } + Ok(new) } /// get the transformation @@ -846,18 +904,45 @@ impl PyView { #[pymethods] impl PyTransforms { #[new] - #[pyo3(signature = (channel = HashMap::new(), drift = Vec::new()))] + #[pyo3(signature = (channel = HashMap::new(), drift = Vec::new(), interpolation_order = None))] fn new( channel: HashMap, drift: Vec, + interpolation_order: Option, ) -> PyResult { let channel = channel.into_iter().map(|(n, t)| (n, t.inner)).collect(); let drift = drift.into_iter().map(|t| t.inner).collect(); Ok(PyTransforms { - inner: Transforms { channel, drift }, + inner: Transforms { + interpolation_order, + channel, + drift, + }, }) } + /// the bspline order used for interpolation when using transforms, + /// None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes + #[getter] + fn get_interpolation_order(&self) -> Option { + self.inner.interpolation_order + } + + /// the bspline order used for interpolation when using transforms + /// None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes + #[setter] + fn set_interpolation_order(&mut self, order: Option) { + self.inner.interpolation_order = order; + } + + /// the bspline order used for interpolation when using transforms + /// None -> 0 (nearest neighbor) for int dtypes, 3 (cubic) for float dtypes + fn with_interpolation_order(&self, order: Option) -> Self { + let mut new = self.clone(); + new.inner.interpolation_order = order; + new + } + #[getter] fn get_channel_transforms(&self) -> HashMap { self.inner diff --git a/ndbioimage/src/readers.rs b/ndbioimage/src/readers.rs index c1e2d19..2ae3a65 100644 --- a/ndbioimage/src/readers.rs +++ b/ndbioimage/src/readers.rs @@ -201,7 +201,7 @@ pub trait Reader: Clone + Sized + Debug + Send + Hash + Into { fn view(&self) -> View { 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( + View::new_with_slice_and_axes( self.clone(), slice.as_ref().to_vec(), vec![Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X], diff --git a/ndbioimage/src/readers/czi.rs b/ndbioimage/src/readers/czi.rs index 717e510..de9a0e0 100644 --- a/ndbioimage/src/readers/czi.rs +++ b/ndbioimage/src/readers/czi.rs @@ -122,14 +122,6 @@ impl From> for Version { } } -/// errors specific to czi file reading -#[derive(Debug, thiserror::Error)] -pub enum CziError { - /// czi file has no valid blocks - #[error("czi file has no valid blocks")] - NoValidBlocks, -} - impl Reader for CziReader { fn new

(path: P, series: usize, position: usize) -> Result where @@ -195,7 +187,7 @@ impl Reader for CziReader { let mut max_y = i32::MIN; for &i in Ok(block_map.get(&(0, 0, 0))) .transpose() - .unwrap_or_else(|| block_map.values().next().ok_or(CziError::NoValidBlocks))? + .unwrap_or_else(|| block_map.values().next().ok_or(Error::CziNoValidBlocks))? { if let Ok(info) = reader.try_get_sub_block_info_for_index(i) { let rect = info.get_logical_rect(); @@ -219,7 +211,7 @@ impl Reader for CziReader { new.shape.x = (max_y - min_y) as usize; new.shape.y = (max_x - min_x) as usize; new.block_map = block_map; - new.pixel_type = pixel_type.ok_or(CziError::NoValidBlocks)?; + new.pixel_type = pixel_type.ok_or(Error::CziNoValidBlocks)?; Ok(new) } diff --git a/ndbioimage/src/transforms.rs b/ndbioimage/src/transforms.rs index 97779fb..8c8ae3c 100644 --- a/ndbioimage/src/transforms.rs +++ b/ndbioimage/src/transforms.rs @@ -29,6 +29,7 @@ impl TryFrom> for Transforms { fn try_from(t: HashMap) -> Result { Ok(Self { + interpolation_order: None, channel: t .into_iter() .map(|(n, i)| { @@ -51,6 +52,8 @@ impl TryFrom> for Transforms { #[derive(Clone, Debug, Default, Eq, PartialEq, Deserialize)] pub struct Transforms { + #[serde(default)] + pub(crate) interpolation_order: Option, #[serde(default)] pub channel: HashMap>, #[serde(default)] @@ -110,6 +113,7 @@ impl Serialize for Transforms { } } +#[derive(Clone, Debug, PartialEq, Copy)] pub enum TransformKind { None, Channel, @@ -150,11 +154,54 @@ impl Transforms { Ok(yaml_serde::to_writer(file, self)?) } + pub fn calculate_transforms_from( + transform_file: Option<&Path>, + channels_files: Option<&[&Path]>, + drift_file: Option<&Path>, + main_channel: Option, + default_transform: Option>, + ) -> Result { + if let Some(transform_file) = transform_file { + Self::load(transform_file) + } else { + let default_transform = default_transform + .map(|default_transform| Transform::new(default_transform, vec![512, 0])); + let mut transforms = Self::default(); + if let Some(channels_files) = channels_files { + let main_channel = main_channel.unwrap_or(0); + transforms.channel = Self::calculate_channel_transforms_2d( + channels_files, + main_channel, + default_transform, + )? + .into_iter() + .map(|(k, t)| (k, t.into_dyn())) + .collect(); + }; + if let Some(drift_file) = drift_file { + transforms.drift = Self::calculate_drift_transforms_2d(drift_file)? + .into_iter() + .map(|t| t.into_dyn()) + .collect(); + }; + Ok(transforms) + } + } + + pub fn calculate_drift_transforms_2d(file: &Path) -> Result>, Error> { + View::<_, DynReader>::from_path(file)?.calculate_drift_transforms_2d() + } + + pub fn calculate_drift_transforms_3d(file: &Path) -> Result>, Error> { + View::<_, DynReader>::from_path(file)?.calculate_drift_transforms_3d() + } + pub fn calculate_channel_transforms_2d( bead_files: &[&Path], main_channel: usize, default_transform: Option>, ) -> Result>, Error> { + // TODO: autodetect main channel let mut transforms = Vec::new(); let default_transform = default_transform.unwrap_or_default(); for file in bead_files { @@ -276,6 +323,12 @@ impl View { self.transforms = transforms; } + pub fn with_transforms(&self, transforms: Transforms) -> Self { + let mut new = self.clone(); + new.transforms = transforms; + new + } + pub fn with_transforms_from_yaml(mut self, path: &Path) -> Result { self.set_transforms(Transforms::load(path)?); Ok(self) @@ -550,7 +603,21 @@ impl View { #[cfg(test)] mod tests { + use crate::axes::Axis; + use crate::error::Error; + use crate::readers::{DynReader, Reader}; use crate::transforms::Transforms; + use image_registration::transform::Transform; + use ndarray::s; + use std::collections::HashMap; + + fn open(file: &str) -> Result { + let path = std::env::current_dir()? + .join("tests") + .join("files") + .join(file); + DynReader::new(&path, 0, 0) + } #[test] fn test_transforms() -> Result<(), Box> { @@ -576,4 +643,80 @@ mod tests { assert_eq!(t.channel["TV2"].ndim(), 2); Ok(()) } + + #[test] + fn transforms_as_array() -> Result<(), Box> { + let file = "czi/beads_2020_08_05__11_43_07.czi"; + let reader = open(file)?; + let view = reader.view(); + let mut channel = HashMap::new(); + channel.insert( + "TV1".to_string(), + Transform::new_with_center( + vec![1.0, 0.0, 0.0, 1.0, 0.0, 0.0], + vec![255.5, 255.5], + vec![512, 512], + ), + ); + channel.insert( + "TV2".to_string(), + Transform::new_with_center( + vec![ + 0.9995868322142533, + 0.012140955373934502, + 0.012592501668218145, + 0.8571690967587099, + -5.270859902118885, + -5.81118547158953, + ], + vec![255.5, 255.5], + vec![512, 512], + ), + ); + let view_t = view.with_transforms(Transforms { + interpolation_order: None, + channel, + drift: Vec::new(), + }); + + let s_a = view.max_proj(Axis::Z)?.slice_cztyx(s![1, 0, .., .., ..])?; + + println!( + "s_a: {:?}, {:?}, {:?}", + s_a.get_slice(), + s_a.get_axes(), + s_a.get_operations() + ); + + let s_b = view_t + .max_proj(Axis::Z)? + .slice_cztyx(s![1, 0, .., .., ..])?; + + println!( + "s_b: {:?}, {:?}, {:?}", + s_b.get_slice(), + s_b.get_axes(), + s_b.get_operations() + ); + + let a = s_a.as_array::()?; + let b = s_b.as_array::()?; + assert_ne!(a, b); + + // let a = view.slice_cztyx(s![1, 0, 0, .., ..])?.as_array::()?; + // let b = view_t.slice_cztyx(s![1, 0, 0, .., ..])?.as_array::()?; + // assert_ne!(a, b); + // + // let a = view + // .max_proj(Axis::Z)? + // .slice(s![1, 0, .., ..])? + // .as_array::()?; + // let b = view_t + // .max_proj(Axis::Z)? + // .slice(s![1, 0, .., ..])? + // .as_array::()?; + // assert_ne!(a, b); + + Ok(()) + } } diff --git a/ndbioimage/src/view.rs b/ndbioimage/src/view.rs index b7ddb74..7b543cc 100644 --- a/ndbioimage/src/view.rs +++ b/ndbioimage/src/view.rs @@ -22,6 +22,7 @@ use std::any::type_name; #[cfg(feature = "transforms")] use std::borrow::Cow; use std::cell::OnceCell; +use std::cmp::PartialEq; use std::collections::{HashMap, HashSet}; use std::fmt::{Debug, Display, Formatter}; use std::hash::{Hash, Hasher}; @@ -112,9 +113,10 @@ pub struct View { operations: IndexMap, #[serde(skip)] channel_names: OnceCell>, - dimensionality: PhantomData, #[cfg(feature = "transforms")] pub(crate) transforms: Transforms, + #[serde(skip)] + dimensionality: PhantomData, } impl Hash for View @@ -134,16 +136,20 @@ where } impl View { - pub(crate) fn new(reader: R, slice: Vec, axes: Vec) -> Self { + pub(crate) fn new_with_slice_and_axes( + reader: R, + slice: Vec, + axes: Vec, + ) -> Self { Self { reader, slice, axes, - operations: IndexMap::new(), - channel_names: OnceCell::new(), - dimensionality: PhantomData, + operations: Default::default(), + channel_names: Default::default(), #[cfg(feature = "transforms")] - transforms: Transforms::default(), + transforms: Default::default(), + dimensionality: PhantomData, } } @@ -205,14 +211,49 @@ impl View { reader, slice, axes, - operations: IndexMap::new(), + operations: Default::default(), dimensionality: PhantomData, - channel_names: OnceCell::new(), + channel_names: Default::default(), #[cfg(feature = "transforms")] - transforms: Transforms::default(), + transforms: Default::default(), }) } + pub fn clone_with_slice_axes( + &self, + slice: Vec, + axes: Vec, + ) -> View { + View { + reader: self.reader.clone(), + slice, + axes, + operations: self.operations.clone(), + channel_names: self.channel_names.clone(), + #[cfg(feature = "transforms")] + transforms: self.transforms.clone(), + dimensionality: PhantomData, + } + } + + pub fn clone_with_slice_axes_operations( + &self, + slice: Vec, + axes: Vec, + operations: IndexMap, + ) -> View { + View { + reader: self.reader.clone(), + slice, + axes, + operations, + channel_names: self.channel_names.clone(), + #[cfg(feature = "transforms")] + transforms: self.transforms.clone(), + dimensionality: PhantomData, + } + } + /// the file path pub fn path(&self) -> &Path { self.reader.path() @@ -267,11 +308,6 @@ impl View { ArrayCache::global().capacity() } - 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 { @@ -282,7 +318,7 @@ impl View { channel_names: self.channel_names, dimensionality: PhantomData, #[cfg(feature = "transforms")] - transforms: Transforms::default(), + transforms: Default::default(), } } @@ -438,7 +474,7 @@ impl View { slice.swap(idx0, idx1); let mut axes = self.axes.clone(); axes.swap(idx0, idx1); - Ok(View::new(self.reader.clone(), slice, axes).with_operations(self.operations.clone())) + Ok(self.clone_with_slice_axes(slice, axes)) } /// subset of gives axes will be reordered in given order @@ -455,7 +491,7 @@ impl View { slice[j] = self.slice[i]; axes[j] = self.axes[i]; } - Ok(View::new(self.reader.clone(), slice, axes).with_operations(self.operations.clone())) + Ok(self.clone_with_slice_axes(slice, axes)) } /// subset of gives axes will be reordered in given order, only axes in axes will be returned @@ -485,18 +521,15 @@ impl View { } } - Ok(View::new(self.reader.clone(), new_slice, new_axes) - .with_operations(self.operations.clone())) + Ok(self.clone_with_slice_axes(new_slice, new_axes)) } /// reverse the order of the axes pub fn transpose(&self) -> Result { - Ok(View::new( - self.reader.clone(), + Ok(self.clone_with_slice_axes( self.slice.iter().rev().cloned().collect(), self.axes.iter().rev().cloned().collect(), - ) - .with_operations(self.operations.clone())) + )) } /// reduce an axis with an operation (max, min, sum or mean), returning a view of lower dimensionality @@ -528,7 +561,7 @@ impl View { operations.insert(ax, operation); (self.axes.clone(), self.slice.clone(), operations) }; - Ok(View::new(self.reader.clone(), slice, axes).with_operations(operations)) + Ok(self.clone_with_slice_axes_operations(slice, axes, operations)) } /// maximum along axis @@ -667,8 +700,7 @@ impl View { new_axes.push(Axis::New); n_idx += 1; } - Ok(View::new(self.reader.clone(), new_slice, new_axes) - .with_operations(self.operations.clone())) + Ok(self.clone_with_slice_axes(new_slice, new_axes)) } /// resets axes to cztyx order, with all 5 axes present, @@ -696,7 +728,7 @@ impl View { slice.push(SliceInfoElem::NewAxis) } } - Ok(View::new(self.reader.clone(), slice, axes).with_operations(self.operations.clone())) + Ok(self.clone_with_slice_axes(slice, axes)) } /// slice, but slice elements are in cztyx order, all cztyx must be given, @@ -759,45 +791,84 @@ impl View { let op_czt = if let Some((&ax, op)) = self.operations.get_index(op_xy.len()) { IndexMap::from([(ax, op.clone())]) } else { - IndexMap::new() + Default::default() }; + let mut slice_xy = false; + let mut xy_dim = 0usize; + for (s, &axis) in self.slice.iter().zip(&self.axes) { + if let (Axis::Y | Axis::X, SliceInfoElem::Slice { start, end, step }) = (axis, s) { + if !op_xy.contains_key(&axis) { + xy_dim += 1; + } + if *start != 0 { + slice_xy = true; + } else if let Some(end) = end + && *end != self.raw_shape()[axis] as isize + { + slice_xy = true; + } else if *step != 1 { + slice_xy = true; + } + } + } + + #[cfg(feature = "transforms")] + let kind = self.transforms.kind(); + #[cfg(feature = "transforms")] + let order = self + .transforms + .interpolation_order + .unwrap_or(if T::ISFLOAT { 3 } else { 0 }); + #[cfg(feature = "transforms")] + let mut transform_done = HashSet::new(); + #[cfg(feature = "transforms")] + let defer_slice_xy = kind != TransformKind::None && xy_dim == 2 && slice_xy; + #[cfg(not(feature = "transforms"))] + let defer_slice = false; + let mut shape_out = 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 { - SliceInfoElem::Slice { start, end, step } => { - 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_out.push(SliceInfoElem::Slice { - start: 0, - end: None, - step: 1, - }); + if let Axis::X | Axis::Y = a + && defer_slice_xy + { + shape_out.push(self.raw_shape()[a]); + slice_out.push(SliceInfoElem::Slice { + start: 0, + end: None, + step: 1, + }); + ax_out.push(*a) + } else { + match s { + SliceInfoElem::Slice { start, end, step } => { + 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_out.push(SliceInfoElem::Slice { + start: 0, + end: None, + step: 1, + }); + ax_out.push(*a); + } + } + SliceInfoElem::Index(_) => {} + SliceInfoElem::NewAxis => { + shape_out.push(1); + slice_out.push(SliceInfoElem::Index(0)); ax_out.push(*a); } } - SliceInfoElem::Index(_) => {} - SliceInfoElem::NewAxis => { - shape_out.push(1); - 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.raw_shape(); for (s, &axis) in self.slice.iter().zip(&self.axes) { 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; - } slice_reader[axis as usize] = Slice::new( idx_bnd(*start, shape_reader[axis] as isize)?, slc_bnd(end.unwrap(), shape_reader[axis] as isize)?, @@ -814,10 +885,20 @@ impl View { (_, 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 xy = if defer_slice_xy { + [SliceInfoElem::Slice { + start: 0, + end: None, + step: 1, + }; 2] + } else { + [ + self.slice[Axis::Y.pos(&self.axes, &self.slice)?], + self.slice[Axis::X.pos(&self.axes, &self.slice)?], + ] + }; + let mut array = match op_czt.first() { Some((_, Operation::Max)) => { ArrayD::::from_elem(shape_out.into_dimension(), T::min_value()) @@ -836,13 +917,6 @@ 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_out[i] = SliceInfoElem::Index(c) @@ -937,16 +1011,12 @@ impl View { Some((_, Operation::Max)) => { array .slice_mut(slice_out.as_slice()) - .zip_mut_with(&arr_frame, |x, y| { - *x = if *x >= *y { *x } else { *y } - }); + .zip_mut_with(&arr_frame, |x, y| *x = if *x >= *y { *x } else { *y }); } Some((_, Operation::Min)) => { array .slice_mut(slice_out.as_slice()) - .zip_mut_with(&arr_frame, |x, y| { - *x = if *x < *y { *x } else { *y } - }); + .zip_mut_with(&arr_frame, |x, y| *x = if *x < *y { *x } else { *y }); } Some((_, Operation::Mean | Operation::Sum)) => { array @@ -973,35 +1043,14 @@ impl View { 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)); - } + self.transform_block( + &mut transform_done, + &mut arr, + &slice_reader, + &ax_out, + defer_slice_xy, + order, + )?; } let a = op.operate(arr, idx)?; @@ -1044,18 +1093,120 @@ impl View { 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()); - } + if (kind != TransformKind::None) + && ax_out.contains_key(&Axis::X) + && ax_out.contains_key(&Axis::Y) + { + self.transform_block( + &mut transform_done, + &mut array, + &slice_reader, + &ax_out, + defer_slice_xy, + order, + )?; + } + + if defer_slice_xy { + let mut deferred_slice = vec![ + SliceInfoElem::Slice { + start: 0, + end: None, + step: 1 + }; + ax_out.len() + ]; + deferred_slice[ax_out[&Axis::Y]] = self.slice[Axis::Y.pos(&self.axes, &self.slice)?]; + deferred_slice[ax_out[&Axis::X]] = self.slice[Axis::X.pos(&self.axes, &self.slice)?]; + array = array.slice(deferred_slice.as_slice()).to_owned(); } ArrayCache::global().insert(key.to_owned(), array.clone()); Ok(array) } + #[cfg(feature = "transforms")] + fn transform_block( + &self, + transform_done: &mut HashSet<(isize, isize, isize)>, + array: &mut Array, + slice: &[Slice], + ax: &HashMap, + defer_slice_xy: bool, + order: usize, + ) -> Result<(), Error> + where + T: Number, + D2: Dimension, + [SliceInfoElem]: SliceArg, + { + let mut slice_transform = vec![ + SliceInfoElem::Slice { + start: 0, + end: None, + step: 1 + }; + ax.len() + ]; + + let shape = array.shape(); + let slice_out = [Axis::C, Axis::Z, Axis::T, Axis::Y, Axis::X] + .into_iter() + .map(|a| { + Slice::new( + 0, + ax.get(&a).and_then(|&i| shape.get(i)).copied().unwrap_or(1) as isize, + 1, + ) + }) + .collect::>(); + + let transpose = ax[&Axis::X] < ax[&Axis::Y]; + for (ci, c) in slice_out[0].into_iter().zip(&slice[0]) { + let channel_name = self.channel_name(c as usize)?; + for (ti, t) in slice_out[2].into_iter().zip(&slice[2]) { + if let Some(transform) = self.transforms.get_transform(channel_name, t as usize)? { + let mut transform = transform.into_owned(); + if !defer_slice_xy { + if transform.ndim() == 2 { + transform.center[0] -= slice[3].start as f64; + transform.center[1] -= slice[4].start as f64; + if transpose { + transform = transform.swap_axes(0, 1) + } + } else { + transform.center[0] -= slice[1].start as f64; + transform.center[1] -= slice[3].start as f64; + transform.center[2] -= slice[4].start as f64; + if transpose { + transform = transform.swap_axes(1, 2); + } + } + } + + for (zi, z) in slice_out[1].into_iter().zip(&slice[1]) { + if !transform_done.contains(&(c, z, t)) { + for (axis, i) in [Axis::C, Axis::Z, Axis::T].iter().zip([ci, zi, ti]) { + if let Some(axis_pos) = ax.get(axis) { + slice_transform[*axis_pos] = SliceInfoElem::Index(i); + } + } + let untransformed = array.slice(slice_transform.as_slice()).into_dyn(); + let transformed = transform + .interpolate_par_order(untransformed, order)? + .mapv(|i| T::from_f64(i).unwrap()); + array + .slice_mut(slice_transform.as_slice()) + .assign(&transformed) + } + transform_done.insert((c, z, t)); + } + } + } + } + Ok(()) + } + /// turn the view into a 1d array pub fn flatten(&self) -> Result, Error> where @@ -1397,8 +1548,6 @@ 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()? @@ -1720,7 +1869,7 @@ mod tests { axes: axes.to_vec(), operations: vec![], #[cfg(feature = "transforms")] - transforms: Transforms::default(), + transforms: Default::default(), }; let k0 = key("u16", 0, &[Axis::T]); let k1 = key("u16", 1, &[Axis::T]); @@ -1734,7 +1883,7 @@ mod tests { assert!(cache.get::(&k0).is_some()); assert!(cache.get::(&k1).is_none()); assert_eq!(cache.len(), 2); - let ops: IndexMap = IndexMap::new(); + let ops: IndexMap = Default::default(); let borrowed = ArrayKeyRef { reader: ReaderKeyRef { name: "test", @@ -1748,7 +1897,7 @@ mod tests { axes: &[Axis::T], operations: &ops, #[cfg(feature = "transforms")] - transforms: &Transforms::default(), + transforms: &Default::default(), }; assert!(cache.get::(&borrowed).is_none()); assert!(cache.get::(&borrowed).is_some());