From 22e8ab1dd2e86daf0616fdc1c6637dad252ff42e Mon Sep 17 00:00:00 2001 From: Wim Pomp Date: Sat, 24 Oct 2026 11:58:06 +0200 Subject: [PATCH] - some more macro use - implement reshape for python - make Reader and View Sync - register_translation_fft - calculate_transforms_from and functions it needs --- image-registration | 2 +- .../py/ndbioimage/ndbioimage_rs/__init__.pyi | 7 +- .../ndbioimage/ndbioimage_rs/transforms.pyi | 4 + ndbioimage/src/py/imread.rs | 601 ++++++++---------- ndbioimage/src/py/transforms.rs | 23 + ndbioimage/src/readers.rs | 2 +- ndbioimage/src/transforms.rs | 161 ++++- ndbioimage/src/view.rs | 52 +- 8 files changed, 479 insertions(+), 373 deletions(-) diff --git a/image-registration b/image-registration index aac61b9..b12982d 160000 --- a/image-registration +++ b/image-registration @@ -1 +1 @@ -Subproject commit aac61b981f539bbdce5a005706c6b114da80c5a6 +Subproject commit b12982da48e422087030a90b2850992b2312bcaa diff --git a/ndbioimage/py/ndbioimage/ndbioimage_rs/__init__.pyi b/ndbioimage/py/ndbioimage/ndbioimage_rs/__init__.pyi index 11f57f8..f19db1f 100644 --- a/ndbioimage/py/ndbioimage/ndbioimage_rs/__init__.pyi +++ b/ndbioimage/py/ndbioimage/ndbioimage_rs/__init__.pyi @@ -301,7 +301,12 @@ class Imread: r""" only remains for backwards compatibility """ - def reshape(self, order: builtins.str, copy: builtins.bool) -> typing.Any: + def reshape( + self, + *shape: typing.Any, + order: builtins.str = "C", + copy: typing.Optional[builtins.bool] = None, + ) -> typing.Any: r""" reshape the view with a new axis order """ diff --git a/ndbioimage/py/ndbioimage/ndbioimage_rs/transforms.pyi b/ndbioimage/py/ndbioimage/ndbioimage_rs/transforms.pyi index 9df8502..b80e77f 100644 --- a/ndbioimage/py/ndbioimage/ndbioimage_rs/transforms.pyi +++ b/ndbioimage/py/ndbioimage/ndbioimage_rs/transforms.pyi @@ -193,6 +193,10 @@ class Transform: def register_translation( fixed: numpy.typing.ArrayLike, moving: numpy.typing.ArrayLike ) -> Transform: ... + @staticmethod + def register_translation_fft( + fixed: numpy.typing.ArrayLike, moving: numpy.typing.ArrayLike + ) -> Transform: ... class Transforms: @property diff --git a/ndbioimage/src/py/imread.rs b/ndbioimage/src/py/imread.rs index 3088925..b390817 100644 --- a/ndbioimage/src/py/imread.rs +++ b/ndbioimage/src/py/imread.rs @@ -6,7 +6,7 @@ use crate::movie::MovieOptions; use crate::readers::{DynReader, PixelType, Reader}; use crate::view::{Item, View}; use itertools::Itertools; -use ndarray::{Ix0, Ix1, IxDyn, SliceInfo, SliceInfoElem}; +use ndarray::{ArrayD, Ix0, Ix1, IxDyn, SliceInfo, SliceInfoElem}; use numpy::{ AllowTypeChange, IntoPyArray, PyArray, PyArrayDescr, PyArrayLike0, PyArrayLike1, PyArrayMethods, dtype, @@ -71,73 +71,30 @@ unsafe impl Sync for PyView {} impl PyView { fn item(py: Python, view: View, dtype: PixelType) -> PyResult> { + macro_rules! item { + ($t:tt) => { + view.into_dimensionality::()? + .item::<$t>()? + .into_pyobject(py)? + .into_any() + }; + } + Ok(match dtype { PixelType::Bool => todo!(), - 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(), + PixelType::I8 => item!(i8), + PixelType::U8 => item!(u8), + PixelType::I16 => item!(i16), + PixelType::U16 => item!(u16), + PixelType::I32 => item!(i32), + PixelType::U32 => item!(u32), + PixelType::F32 => item!(f32), + PixelType::F64 => item!(f64), + PixelType::I64 => item!(i64), + PixelType::U64 => item!(u64), + PixelType::I128 => item!(i128), + PixelType::U128 => item!(u128), + PixelType::F128 => item!(f64), }) } } @@ -299,9 +256,45 @@ impl PyView { } /// reshape the view with a new axis order - #[expect(unused_variables)] - fn reshape<'py>(&self, order: &str, copy: bool) -> PyResult> { - todo!() + #[pyo3(signature = (*shape, order="C", copy=None))] + fn reshape<'py>( + &self, + py: Python<'py>, + shape: Vec, + order: &str, + copy: Option, + ) -> PyResult> { + // TODO: check shape argument + if order != "C" { + return Err(PyNotImplementedError::new_err( + "order except C is not implemented", + )); + } + if let Some(copy) = copy + && !copy + { + return Err(PyNotImplementedError::new_err( + "copy = False is not implemented", + )); + } + let n = self.view.shape().iter().product::(); + let m = shape + .iter() + .fold(1, |i, j| if *j < 0 { i } else { i * (*j as usize) }); + let shape = shape + .iter() + .map(|&s| if s < 0 { n / m } else { s as usize }) + .collect::>(); + if shape.iter().zip(self.view.shape()).all(|(&i, j)| i == j) { + self.clone().into_bound_py_any(py) + } else { + let vec = + py.detach(|| Ok::<_, Error>(self.view.as_array::()?.flatten().to_vec()))?; + Ok(ArrayD::from_shape_vec(shape, vec) + .map_err(Error::from)? + .into_pyarray(py) + .into_any()) + } } #[gen_stub(override_return_type(type_repr="Imread | float | int", imports=("numpy", "numpy.typing")))] @@ -351,73 +344,30 @@ impl PyView { self.view.squeeze()? }; if view.ndim() == 0 { + macro_rules! squeeze { + ($t:tt) => { + view.into_dimensionality::()? + .item::<$t>()? + .into_pyobject(py)? + .into_any() + }; + } + Ok(match self.dtype { PixelType::Bool => todo!(), - PixelType::I8 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), - PixelType::U8 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), - PixelType::I16 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), - PixelType::U16 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), - PixelType::I32 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), - PixelType::U32 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), - PixelType::I64 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), - PixelType::U64 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), - PixelType::I128 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), - PixelType::U128 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), - PixelType::F32 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), - PixelType::F64 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), - PixelType::F128 => view - .into_dimensionality::()? - .item::()? - .into_pyobject(py)? - .into_any(), + PixelType::I8 => squeeze!(i8), + PixelType::U8 => squeeze!(u8), + PixelType::I16 => squeeze!(i16), + PixelType::U16 => squeeze!(u16), + PixelType::I32 => squeeze!(i32), + PixelType::U32 => squeeze!(u32), + PixelType::I64 => squeeze!(i64), + PixelType::U64 => squeeze!(u64), + PixelType::I128 => squeeze!(i128), + PixelType::U128 => squeeze!(u128), + PixelType::F32 => squeeze!(f32), + PixelType::F64 => squeeze!(f64), + PixelType::F128 => squeeze!(f64), }) } else { PyView { @@ -588,7 +538,6 @@ impl PyView { } /// convert to a numpy array, optionally with a different dtype - #[expect(unused_variables)] #[pyo3(signature = (dtype = None, copy = None))] #[gen_stub(override_return_type(type_repr="numpy.typing.NDArray", imports=("numpy", "numpy.typing")))] fn __array__<'py>( @@ -602,6 +551,13 @@ impl PyView { dtype: Option>, copy: Option, ) -> PyResult> { + if let Some(copy) = copy + && !copy + { + return Err(PyNotImplementedError::new_err( + "copy = False is not implemented", + )); + } if let Some(dtype) = dtype { self.as_type(py, dtype)?.as_array(py) } else { @@ -1232,114 +1188,83 @@ impl PyView { z: isize, t: isize, ) -> PyResult> { + macro_rules! get_frame { + ($t:tt) => { + self.view + .get_frame::<$t, _>(c, z, t)? + .into_pyarray(py) + .into_any() + }; + } + Ok(match self.dtype { PixelType::Bool => todo!(), - 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(), + PixelType::I8 => get_frame!(i8), + PixelType::U8 => get_frame!(u8), + PixelType::I16 => get_frame!(i16), + PixelType::U16 => get_frame!(u16), + PixelType::I32 => get_frame!(i32), + PixelType::U32 => get_frame!(u32), + PixelType::F32 => get_frame!(f32), + PixelType::F64 => get_frame!(f64), + PixelType::I64 => get_frame!(i64), + PixelType::U64 => get_frame!(u64), + PixelType::I128 => get_frame!(i64), + PixelType::U128 => get_frame!(u64), + PixelType::F128 => get_frame!(f64), }) } /// flatten the view into a 1D numpy array #[gen_stub(override_return_type(type_repr="numpy.typing.NDArray | float", imports=("numpy", "numpy.typing")))] fn flatten<'py>(&self, py: Python<'py>) -> PyResult> { + macro_rules! flatten { + ($t:tt) => { + self.view.flatten::<$t>()?.into_pyarray(py).into_any() + }; + } + Ok(match self.dtype { PixelType::Bool => todo!(), - PixelType::I8 => self.view.flatten::()?.into_pyarray(py).into_any(), - PixelType::U8 => self.view.flatten::()?.into_pyarray(py).into_any(), - PixelType::I16 => self.view.flatten::()?.into_pyarray(py).into_any(), - PixelType::U16 => self.view.flatten::()?.into_pyarray(py).into_any(), - PixelType::I32 => self.view.flatten::()?.into_pyarray(py).into_any(), - PixelType::U32 => self.view.flatten::()?.into_pyarray(py).into_any(), - PixelType::F32 => self.view.flatten::()?.into_pyarray(py).into_any(), - PixelType::F64 => self.view.flatten::()?.into_pyarray(py).into_any(), - PixelType::I64 => self.view.flatten::()?.into_pyarray(py).into_any(), - PixelType::U64 => self.view.flatten::()?.into_pyarray(py).into_any(), - PixelType::I128 => self.view.flatten::()?.into_pyarray(py).into_any(), - PixelType::U128 => self.view.flatten::()?.into_pyarray(py).into_any(), - PixelType::F128 => self.view.flatten::()?.into_pyarray(py).into_any(), + PixelType::I8 => flatten!(i8), + PixelType::U8 => flatten!(u8), + PixelType::I16 => flatten!(i16), + PixelType::U16 => flatten!(u16), + PixelType::I32 => flatten!(i32), + PixelType::U32 => flatten!(u32), + PixelType::F32 => flatten!(f32), + PixelType::F64 => flatten!(f64), + PixelType::I64 => flatten!(i64), + PixelType::U64 => flatten!(u64), + PixelType::I128 => flatten!(i64), + PixelType::U128 => flatten!(u64), + PixelType::F128 => flatten!(f64), }) } /// convert the view to bytes fn to_bytes(&self) -> PyResult> { + macro_rules! to_bytes { + ($t:tt) => { + self.view.to_bytes::<$t>()? + }; + } + Ok(match self.dtype { PixelType::Bool => todo!(), - PixelType::I8 => self.view.to_bytes::()?, - PixelType::U8 => self.view.to_bytes::()?, - PixelType::I16 => self.view.to_bytes::()?, - PixelType::U16 => self.view.to_bytes::()?, - PixelType::I32 => self.view.to_bytes::()?, - PixelType::U32 => self.view.to_bytes::()?, - PixelType::F32 => self.view.to_bytes::()?, - PixelType::F64 => self.view.to_bytes::()?, - PixelType::I64 => self.view.to_bytes::()?, - PixelType::U64 => self.view.to_bytes::()?, - PixelType::I128 => self.view.to_bytes::()?, - PixelType::U128 => self.view.to_bytes::()?, - PixelType::F128 => self.view.to_bytes::()?, + PixelType::I8 => to_bytes!(i8), + PixelType::U8 => to_bytes!(u8), + PixelType::I16 => to_bytes!(i16), + PixelType::U16 => to_bytes!(u16), + PixelType::I32 => to_bytes!(i32), + PixelType::U32 => to_bytes!(u32), + PixelType::F32 => to_bytes!(f32), + PixelType::F64 => to_bytes!(f64), + PixelType::I64 => to_bytes!(i64), + PixelType::U64 => to_bytes!(u64), + PixelType::I128 => to_bytes!(i64), + PixelType::U128 => to_bytes!(u64), + PixelType::F128 => to_bytes!(f64), }) } @@ -1488,60 +1413,29 @@ impl PyView { /// collect data into a numpy array #[gen_stub(override_return_type(type_repr="numpy.typing.NDArray | float", imports=("numpy", "numpy.typing")))] fn as_array<'py>(&self, py: Python<'py>) -> PyResult> { + macro_rules! as_array { + ($t:tt) => { + py.detach(|| self.view.as_array_dyn::<$t>())? + .into_pyarray(py) + .into_any() + }; + } + Ok(match self.dtype { PixelType::Bool => todo!(), - 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(), + PixelType::I8 => as_array!(i8), + PixelType::U8 => as_array!(u8), + PixelType::I16 => as_array!(i16), + PixelType::U16 => as_array!(u16), + PixelType::I32 => as_array!(i32), + PixelType::U32 => as_array!(u32), + PixelType::F32 => as_array!(f32), + PixelType::F64 => as_array!(f64), + PixelType::I64 => as_array!(i64), + PixelType::U64 => as_array!(u64), + PixelType::I128 => as_array!(i64), + PixelType::U128 => as_array!(u64), + PixelType::F128 => as_array!(f64), }) } @@ -1549,25 +1443,33 @@ impl PyView { #[gen_stub(override_return_type(type_repr = "numpy.dtype", imports=("numpy")))] #[getter] fn get_dtype<'py>(&self, py: Python<'py>) -> PyResult> { + macro_rules! get_dtype { + ($t:tt) => { + Ok(dtype::<$t>(py)) + }; + } + match self.dtype { - PixelType::Bool => todo!(), - PixelType::I8 => Ok(dtype::(py)), - PixelType::U8 => Ok(dtype::(py)), - PixelType::I16 => Ok(dtype::(py)), - PixelType::U16 => Ok(dtype::(py)), - PixelType::I32 => Ok(dtype::(py)), - PixelType::U32 => Ok(dtype::(py)), - PixelType::F32 => Ok(dtype::(py)), - PixelType::F64 => Ok(dtype::(py)), - PixelType::I64 => Ok(dtype::(py)), - PixelType::U64 => Ok(dtype::(py)), + PixelType::Bool => get_dtype!(bool), + PixelType::I8 => get_dtype!(i8), + PixelType::U8 => get_dtype!(u8), + PixelType::I16 => get_dtype!(i16), + PixelType::U16 => get_dtype!(u16), + PixelType::I32 => get_dtype!(i32), + PixelType::U32 => get_dtype!(u32), + PixelType::F32 => get_dtype!(f32), + PixelType::F64 => get_dtype!(f64), + PixelType::I64 => get_dtype!(i64), + PixelType::U64 => get_dtype!(u64), PixelType::I128 => Err(PyTypeError::new_err( "type is i128, but this cannot be represented by numpy.dtype", )), PixelType::U128 => Err(PyTypeError::new_err( "type is u128, but this cannot be represented by numpy.dtype", )), - PixelType::F128 => Ok(dtype::(py)), + PixelType::F128 => Err(PyTypeError::new_err( + "type is f128, but this is not supported by rust", + )), } } @@ -1619,21 +1521,27 @@ impl PyView { } .into_bound_py_any(py) } else { + macro_rules! max { + ($t:tt) => { + self.view.max::<$t>()?.into_pyobject(py)?.into_any() + }; + } + Ok(match self.dtype { PixelType::Bool => todo!(), - 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(), + PixelType::I8 => max!(i8), + PixelType::U8 => max!(u8), + PixelType::I16 => max!(i16), + PixelType::U16 => max!(u16), + PixelType::I32 => max!(i32), + PixelType::U32 => max!(u32), + PixelType::F32 => max!(f32), + PixelType::F64 => max!(f64), + PixelType::I64 => max!(i64), + PixelType::U64 => max!(u64), + PixelType::I128 => max!(i64), + PixelType::U128 => max!(u64), + PixelType::F128 => max!(f64), }) } } @@ -1674,21 +1582,27 @@ impl PyView { } .into_bound_py_any(py) } else { + macro_rules! min { + ($t:tt) => { + self.view.min::<$t>()?.into_pyobject(py)?.into_any() + }; + } + Ok(match self.dtype { PixelType::Bool => todo!(), - 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(), + PixelType::I8 => min!(i8), + PixelType::U8 => min!(u8), + PixelType::I16 => min!(i16), + PixelType::U16 => min!(u16), + PixelType::I32 => min!(i32), + PixelType::U32 => min!(u32), + PixelType::F32 => min!(f32), + PixelType::F64 => min!(f64), + PixelType::I64 => min!(i64), + PixelType::U64 => min!(u64), + PixelType::I128 => min!(i64), + PixelType::U128 => min!(u64), + PixelType::F128 => min!(f64), }) } } @@ -1760,7 +1674,7 @@ impl PyView { ))?; } let dtype = match self.dtype { - PixelType::Bool => todo!(), + PixelType::Bool => PixelType::U8, PixelType::I8 => PixelType::I16, PixelType::U8 => PixelType::U16, PixelType::I16 => PixelType::I32, @@ -1784,15 +1698,21 @@ impl PyView { } .into_bound_py_any(py) } else { + macro_rules! sum { + ($t:tt) => { + self.view.sum::<$t>()?.into_pyobject(py)?.into_any() + }; + } + 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(), + PixelType::F32 => sum!(f32), + PixelType::F64 => sum!(f64), + PixelType::I64 => sum!(i64), + PixelType::U64 => sum!(u64), + PixelType::I128 => sum!(i128), + PixelType::U128 => sum!(u128), + PixelType::F128 => sum!(f64), + _ => sum!(i64), }) } } @@ -2021,20 +1941,27 @@ impl PyView { }; let options = crate::tiffwrite::TiffOptions::new(bar, None, colors.unwrap_or_default(), overwrite)?; + + macro_rules! save_as_tiff { + ($t:tt) => { + self.view.save_as_tiff_with_type::<$t, _>(file, &options) + }; + } + py.detach(|| match self.dtype { PixelType::Bool => todo!(), - PixelType::I8 => self.view.save_as_tiff_with_type::(file, &options), - PixelType::U8 => self.view.save_as_tiff_with_type::(file, &options), - PixelType::I16 => self.view.save_as_tiff_with_type::(file, &options), - PixelType::U16 => self.view.save_as_tiff_with_type::(file, &options), - PixelType::I32 => self.view.save_as_tiff_with_type::(file, &options), - PixelType::U32 => self.view.save_as_tiff_with_type::(file, &options), - PixelType::I64 => self.view.save_as_tiff_with_type::(file, &options), - PixelType::U64 => self.view.save_as_tiff_with_type::(file, &options), - PixelType::I128 => self.view.save_as_tiff_with_type::(file, &options), - PixelType::U128 => self.view.save_as_tiff_with_type::(file, &options), - PixelType::F32 => self.view.save_as_tiff_with_type::(file, &options), - PixelType::F64 => self.view.save_as_tiff_with_type::(file, &options), + PixelType::I8 => save_as_tiff!(i8), + PixelType::U8 => save_as_tiff!(u8), + PixelType::I16 => save_as_tiff!(i16), + PixelType::U16 => save_as_tiff!(u16), + PixelType::I32 => save_as_tiff!(i32), + PixelType::U32 => save_as_tiff!(u32), + PixelType::I64 => save_as_tiff!(i64), + PixelType::U64 => save_as_tiff!(u64), + PixelType::I128 => save_as_tiff!(i128), + PixelType::U128 => save_as_tiff!(u128), + PixelType::F32 => save_as_tiff!(f32), + PixelType::F64 => save_as_tiff!(f64), PixelType::F128 => Err(Error::NotImplemented( "saving as f128 is not implemented".to_string(), )), diff --git a/ndbioimage/src/py/transforms.rs b/ndbioimage/src/py/transforms.rs index c446183..203542a 100644 --- a/ndbioimage/src/py/transforms.rs +++ b/ndbioimage/src/py/transforms.rs @@ -664,6 +664,29 @@ impl PyTransform { ))), }) } + + #[staticmethod] + fn register_translation_fft<'py>( + py: Python<'py>, + #[gen_stub(override_type(type_repr="numpy.typing.ArrayLike", imports=("numpy.typing")))] + fixed: PyArrayLikeDyn, + #[gen_stub(override_type(type_repr="numpy.typing.ArrayLike", imports=("numpy.typing")))] + moving: PyArrayLikeDyn, + ) -> PyResult { + let fixed = fixed.as_array(); + let moving = moving.as_array(); + if fixed.shape() != moving.shape() { + return Err(PyErr::from(Error::ShapeMismatch( + fixed.shape().to_vec(), + moving.shape().to_vec(), + ))); + } + py.detach(|| { + Ok(PyTransform { + inner: Transform::register_translation_fft(fixed, moving).map_err(Error::from)?, + }) + }) + } } #[gen_stub_pyclass] diff --git a/ndbioimage/src/readers.rs b/ndbioimage/src/readers.rs index d87eb78..36bd46f 100644 --- a/ndbioimage/src/readers.rs +++ b/ndbioimage/src/readers.rs @@ -182,7 +182,7 @@ pub enum ArrayT { pub type Frame = ArrayT; /// trait for reading image files -pub trait Reader: Clone + Sized + Debug + Send + Hash + Into { +pub trait Reader: Clone + Sized + Debug + Sync + Send + Hash + Into { /// create a new reader for the given path, series, and position fn new

(path: P, series: usize, position: usize) -> Result where diff --git a/ndbioimage/src/transforms.rs b/ndbioimage/src/transforms.rs index 4bcf4c6..95c1788 100644 --- a/ndbioimage/src/transforms.rs +++ b/ndbioimage/src/transforms.rs @@ -1,9 +1,11 @@ +use rayon::iter::ParallelIterator; use crate::axes::Axis; use crate::error::Error; use crate::readers::{DynReader, Frame, PixelType, Reader}; use crate::view::View; +use image_registration::filter; pub use image_registration::transform::Transform; -use ndarray::{Array2, Dimension, Ix2, Ix3, IxDyn, array, s}; +use ndarray::{array, s, Array2, Dimension, Ix2, Ix3, IxDyn, concatenate}; use serde::{Deserialize, Serialize}; use std::borrow::Cow; use std::collections::{HashMap, HashSet}; @@ -11,6 +13,8 @@ use std::hash::{Hash, Hasher}; use std::io::Seek; use std::path::Path; use std::sync::Arc; +use rayon::prelude::IntoParallelIterator; +use tiffwrite::IJTiffFile; #[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)] struct TransformOld { @@ -158,22 +162,27 @@ 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>, + pub fn calculate_transforms_from( + transform_file: Option

, + bead_files: Option<&[Q]>, + drift_file: Option, main_channel: Option, default_transform: Option>, - ) -> Result { - if let Some(transform_file) = transform_file { - Self::load(transform_file) + ) -> Result + where + P: AsRef, + Q: AsRef, + R: AsRef, + { + if let Some(transform_file) = transform_file.as_ref() && let Ok(transforms) = Self::load(transform_file) { + Ok(transforms) } else { let default_transform = default_transform.map(Transform::new); let mut transforms = Self::default(); - if let Some(channels_files) = channels_files { + if let Some(bead_files) = bead_files { let main_channel = main_channel.unwrap_or(0); transforms.channel = Self::calculate_channel_transforms_2d( - channels_files, + bead_files, main_channel, default_transform, )? @@ -187,20 +196,58 @@ impl Transforms { .map(|t| t.into_dyn()) .collect(); }; + if let Some(transform_file) = transform_file { + let transform_file = transform_file.as_ref(); + transforms.save(transform_file)?; + if let Some(bead_files) = bead_files { + transforms.save_channel_transform_tiff(transform_file.with_extension("tif"), bead_files)?; + } + + } Ok(transforms) } } - pub fn calculate_drift_transforms_2d(file: &Path) -> Result>, Error> { + /// save a tiff file with the uncorrected and corrected frames of the bead files side by side + pub fn save_channel_transform_tiff(&self, file: P, bead_files: &[Q]) -> Result<(), Error> + where + P: AsRef, + Q: AsRef, + { + let mut transforms = self.clone(); + transforms.drift = Vec::new(); + let mut tif = IJTiffFile::new(file)?; + for (t, bead_file) in bead_files.iter().enumerate() { + let mut view = View::<_, DynReader>::from_path(bead_file)?.into_dyn(); + let axes = view.axes(); + if axes.contains(&Axis::Z) { + view = view.max_proj(Axis::Z)?; + } + if axes.contains(&Axis::T) { + view = view.max_proj(Axis::T)?; + } + let view_u = view.into_dimensionality::()?; + let view_t = view_u.with_transforms(transforms.clone()); + for c in 0..view_u.shape().c { + let a = view_u.slice_cztyx(s![c, 0, 0, .., ..])?.as_array::()?; + let b = view_t.slice_cztyx(s![c, 0, 0, .., ..])?.as_array::()?; + let frame = concatenate(ndarray::Axis(1), &[a.view(), b.view()])?; + tif.save(frame, c, 0, t)?; + } + } + Ok(()) + } + + pub fn calculate_drift_transforms_2d>(file: P) -> Result>, Error> { View::<_, DynReader>::from_path(file)?.calculate_drift_transforms_2d() } - pub fn calculate_drift_transforms_3d(file: &Path) -> Result>, Error> { + pub fn calculate_drift_transforms_3d>(file: P) -> Result>, Error> { View::<_, DynReader>::from_path(file)?.calculate_drift_transforms_3d() } - pub fn calculate_channel_transforms_2d( - bead_files: &[&Path], + pub fn calculate_channel_transforms_2d>( + bead_files: &[P], main_channel: usize, default_transform: Option>, ) -> Result>, Error> { @@ -249,8 +296,8 @@ impl Transforms { Ok(average_transforms) } - pub fn calculate_channel_transforms_3d( - bead_files: &[&Path], + pub fn calculate_channel_transforms_3d>( + bead_files: &[P], main_channel: usize, default_transform: Option>, ) -> Result>, Error> { @@ -401,7 +448,85 @@ impl View { } pub fn calculate_drift_transforms_2d(&self) -> Result>, Error> { - todo!() + let mut view = self.clone().into_dyn(); + let axes = view.axes(); + if axes.contains(&Axis::Z) { + view = view.max_proj(Axis::Z)?; + } + if axes.contains(&Axis::C) { + view = view.max_proj(Axis::C)?; + } + view = view.permute_axes(&[Axis::T, Axis::Y, Axis::X])?; + let n = view.shape().t; + let m = (n as f64).sqrt().round() as usize; + let t_groups = (0..=m) + .map(|i| (((i * n) as f64) / (m as f64)).round() as usize) + .collect::>(); + let t_keys = t_groups + .array_windows() + .map(|[start, end]| (start + end) / 2) + .collect::>(); + let t_pairs = t_keys + .iter() + .zip(t_groups.array_windows()) + .flat_map(|(t_key, [start, end])| (*start..*end).map(|i| (*t_key, i))) + .chain(t_keys.array_windows().map(|[start, end]| (*start, *end))) + .collect::>(); + let fmaxz_keys = t_keys + .iter() + .map(|t| { + Ok(( + t, + filter::gaussian_smooth( + &view.slice(s![*t, .., ..])?.as_array::()?, + &[5.0; 2], + )?, + )) + }) + .collect::, Error>>()?; + let shifts = t_pairs + .into_par_iter() + .map(|(t_key, t)| { + if t_key == t { + Ok(Transform::::from_translation(&vec![0.0; view.ndim()])) + } else { + if let Some(fmaxz) = fmaxz_keys.get(&t_key) { + Ok(Transform::::register_translation( + fmaxz_keys[&t_key].view(), + fmaxz, + )?) + } else { + let fmaxz = filter::gaussian_smooth( + &view.slice(s![t, .., ..])?.as_array::()?, + &[5.0; 2], + )?; + Ok(Transform::::register_translation( + fmaxz_keys[&t_key].view(), + fmaxz.view(), + )?) + } + } + }) + .collect::, Error>>()?; + let mut t = Transform::::from_translation(&[0.0; 2]); + let mut key_shifts = vec![t.clone()]; + key_shifts.extend(shifts.iter().skip(n).map(|i| { + t = i * &t; + t.clone() + })); + Ok(key_shifts + .iter() + .zip(t_groups.array_windows()) + .flat_map(|(s, [start, end])| { + let mut t = s.clone(); + (*start..*end) + .map(|i| { + t = &shifts[i] * &t; + t.clone() + }) + .collect::>() + }) + .collect()) } pub fn calculate_drift_transforms_3d(&self) -> Result>, Error> { @@ -489,7 +614,7 @@ impl View { .collect::, Error>>()?; // TODO; check these translations - let u0 = tr.as_translation(&[(z_min as f64), 0.0, 0.0]); + 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; diff --git a/ndbioimage/src/view.rs b/ndbioimage/src/view.rs index 46c435b..643ee48 100644 --- a/ndbioimage/src/view.rs +++ b/ndbioimage/src/view.rs @@ -11,8 +11,8 @@ 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, + Array, Array0, Array1, Array2, ArrayD, Dimension, IntoDimension, Ix0, Ix1, + Ix2, Ix5, IxDyn, SliceArg, SliceInfoElem, s, }; use num::traits::{AsPrimitive, ToBytes}; use num::{Bounded, FromPrimitive, One, ToPrimitive, Zero}; @@ -21,7 +21,6 @@ use serde_with::serde_as; 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}; @@ -30,7 +29,7 @@ use std::iter::Sum; use std::marker::PhantomData; use std::ops::{AddAssign, Deref, Div}; use std::path::Path; -use std::sync::Arc; +use std::sync::{Arc, OnceLock}; fn idx_bnd(idx: isize, bnd: isize) -> Result { if idx < -bnd { @@ -112,7 +111,7 @@ pub struct View { axes: Vec, operations: IndexMap, #[serde(skip)] - channel_names: OnceCell>, + channel_names: OnceLock>, #[cfg(feature = "transforms")] pub(crate) transforms: Transforms, #[serde(skip)] @@ -830,8 +829,20 @@ impl View { } else { Default::default() }; + #[cfg(feature = "transforms")] let mut slice_xy = false; let mut xy_dim = 0usize; + + #[cfg(not(feature = "transforms"))] + for &axis in &self.axes { + if let Axis::Y | Axis::X = axis + && !op_xy.contains_key(&axis) + { + xy_dim += 1; + } + } + + #[cfg(feature = "transforms")] 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) { @@ -861,7 +872,7 @@ impl View { #[cfg(feature = "transforms")] let defer_slice_xy = kind != TransformKind::None && xy_dim == 2 && slice_xy; #[cfg(not(feature = "transforms"))] - let defer_slice = false; + let defer_slice_xy = false; let mut shape_out = Vec::new(); let mut slice_out = Vec::new(); @@ -900,7 +911,7 @@ impl View { } } } - let mut slice_reader = vec![Slice::empty(); 5]; + let mut slice_reader = [Slice::empty(); 5]; let shape_reader = self.reader.raw_shape(); for (s, &axis) in self.slice.iter().zip(&self.axes) { match (axis, s) { @@ -1115,14 +1126,6 @@ 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) - }; - #[cfg(feature = "transforms")] let mut array = if n == 1 { out.take().unwrap() } else { @@ -1585,6 +1588,7 @@ mod tests { DynReader::new(&path, 0, 0) } + #[cfg(feature = "tiffseq")] #[test] fn view() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; @@ -1602,6 +1606,7 @@ mod tests { Ok(()) } + #[cfg(feature = "tiffseq")] #[test] fn view_shape() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; @@ -1613,6 +1618,7 @@ mod tests { Ok(()) } + #[cfg(feature = "tiffseq")] #[test] fn view_new_axis() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; @@ -1628,6 +1634,7 @@ mod tests { Ok(()) } + #[cfg(feature = "tiffseq")] #[test] fn view_permute_axes() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; @@ -1655,6 +1662,7 @@ mod tests { macro_rules! test_max { ($($name:ident: $b:expr $(,)?)*) => { $( + #[cfg(feature = "tiffseq")] #[test] fn $name() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; @@ -1683,6 +1691,7 @@ mod tests { macro_rules! test_index { ($($name:ident: $b:expr $(,)?)*) => { $( + #[cfg(feature = "tiffseq")] #[test] fn $name() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; @@ -1729,6 +1738,7 @@ mod tests { index_r: s![0, 0, 0, 0, 0] } + #[cfg(feature = "tiffseq")] #[test] fn dyn_view() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; @@ -1741,6 +1751,7 @@ mod tests { Ok(()) } + #[cfg(feature = "czi")] #[test] fn item() -> Result<(), Error> { let file = "czi/1xp53-01-AP1.czi"; @@ -1753,6 +1764,7 @@ mod tests { Ok(()) } + #[cfg(feature = "czi")] #[test] fn slice_cztyx() -> Result<(), Error> { let file = "czi/1xp53-01-AP1.czi"; @@ -1770,6 +1782,7 @@ mod tests { Ok(()) } + #[cfg(feature = "czi")] #[test] fn reset_axes() -> Result<(), Error> { let file = "czi/1xp53-01-AP1.czi"; @@ -1782,6 +1795,7 @@ mod tests { Ok(()) } + #[cfg(feature = "czi")] #[test] fn reset_axes2() -> Result<(), Error> { let file = "czi/Experiment-2029.czi"; @@ -1792,6 +1806,7 @@ mod tests { Ok(()) } + #[cfg(feature = "czi")] #[test] fn reset_axes3() -> Result<(), Error> { let file = "czi/Experiment-2029.czi"; @@ -1820,6 +1835,7 @@ mod tests { Ok(()) } + #[cfg(feature = "czi")] #[test] fn max() -> Result<(), Error> { let file = "czi/Experiment-2029.czi"; @@ -1932,6 +1948,7 @@ mod tests { Ok(()) } + #[cfg(feature = "tiffseq")] #[test] fn as_array_caches_frames() -> Result<(), Error> { let file = "tiffseq/YTL1841B2-2-1_1hr_DMSO_galinduction_1"; @@ -1977,6 +1994,11 @@ mod tests { Ok(()) } + #[cfg(any( + feature = "tiff", + feature = "bioformats_java", + feature = "bioformats_rust" + ))] #[test] fn serde() -> Result<(), Error> { let file = "tiff/20251014_20-Pos_000_000_max.tif";