wim/ndbioimage
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

- fix Imread.transpose

Browse Source 
- implement with_transform
- implement 2d transforms in View.as_array_dyn
This commit is contained in:
w.pomp committed 2026-10-21 20:08:42 +02:00
1 parent 453237cae9
commit e1caa1c0b0
12 files changed
+662 -187
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

image-registration
+1 -1
Submodule image-registration updated: 5c580dc80f...dc1d5d8364.
ndbioimage/py/ndbioimage/ndbioimage_rs/__init__.pyi
+80 -48
View File
@@ -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"""
ndbioimage/py/ndbioimage/ndbioimage_rs/transforms.pyi
+20
View File
@@ -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]]: ...
ndbioimage/src/axes.rs
+2 -2
View File
@@ -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
}
}
ndbioimage/src/cache.rs
+2 -2
View File
@@ -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;
ndbioimage/src/error.rs
+6 -5
View File
@@ -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 {
ndbioimage/src/py/imread.rs
+72 -19
View File
@@ -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<Vec<Bound<PyAny>>>,
axes: Option<Bound<PyAny>>,
) -> PyResult<Self> {
// TODO: transpose("xyz"), transpose(*axes)
let view = if let Some(axes) = axes {
let ax = axes
.into_iter()
.map(|a| self.get_ax(a))
.collect::<Result<Vec<_>, _>>()?;
.try_iter()?
.map(|a| self.get_ax(a?))
.collect::<PyResult<Vec<_>>>()?;
self.view.permute_axes(&ax)?
} else {
self.view.transpose()?
@@ -1373,19 +1382,58 @@ impl PyView {
fn as_array<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyAny>> {
Ok(match self.dtype {
PixelType::Bool => todo!(),
PixelType::I8 => self.view.as_array_dyn::<i8>()?.into_pyarray(py).into_any(),
PixelType::U8 => self.view.as_array_dyn::<u8>()?.into_pyarray(py).into_any(),
PixelType::I16 => self.view.as_array_dyn::<i16>()?.into_pyarray(py).into_any(),
PixelType::U16 => self.view.as_array_dyn::<u16>()?.into_pyarray(py).into_any(),
PixelType::I32 => self.view.as_array_dyn::<i32>()?.into_pyarray(py).into_any(),
PixelType::U32 => self.view.as_array_dyn::<u32>()?.into_pyarray(py).into_any(),
PixelType::F32 => self.view.as_array_dyn::<f32>()?.into_pyarray(py).into_any(),
PixelType::F64 => self.view.as_array_dyn::<f64>()?.into_pyarray(py).into_any(),
PixelType::I64 => self.view.as_array_dyn::<i64>()?.into_pyarray(py).into_any(),
PixelType::U64 => self.view.as_array_dyn::<u64>()?.into_pyarray(py).into_any(),
PixelType::I128 => self.view.as_array_dyn::<i64>()?.into_pyarray(py).into_any(),
PixelType::U128 => self.view.as_array_dyn::<u64>()?.into_pyarray(py).into_any(),
PixelType::F128 => self.view.as_array_dyn::<f64>()?.into_pyarray(py).into_any(),
PixelType::I8 => py
.detach(|| self.view.as_array_dyn::<i8>())?
.into_pyarray(py)
.into_any(),
PixelType::U8 => py
.detach(|| self.view.as_array_dyn::<u8>())?
.into_pyarray(py)
.into_any(),
PixelType::I16 => py
.detach(|| self.view.as_array_dyn::<i16>())?
.into_pyarray(py)
.into_any(),
PixelType::U16 => py
.detach(|| self.view.as_array_dyn::<u16>())?
.into_pyarray(py)
.into_any(),
PixelType::I32 => py
.detach(|| self.view.as_array_dyn::<i32>())?
.into_pyarray(py)
.into_any(),
PixelType::U32 => py
.detach(|| self.view.as_array_dyn::<u32>())?
.into_pyarray(py)
.into_any(),
PixelType::F32 => py
.detach(|| self.view.as_array_dyn::<f32>())?
.into_pyarray(py)
.into_any(),
PixelType::F64 => py
.detach(|| self.view.as_array_dyn::<f64>())?
.into_pyarray(py)
.into_any(),
PixelType::I64 => py
.detach(|| self.view.as_array_dyn::<i64>())?
.into_pyarray(py)
.into_any(),
PixelType::U64 => py
.detach(|| self.view.as_array_dyn::<u64>())?
.into_pyarray(py)
.into_any(),
PixelType::I128 => py
.detach(|| self.view.as_array_dyn::<i64>())?
.into_pyarray(py)
.into_any(),
PixelType::U128 => py
.detach(|| self.view.as_array_dyn::<u64>())?
.into_pyarray(py)
.into_any(),
PixelType::F128 => py
.detach(|| self.view.as_array_dyn::<f64>())?
.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<String> {
Ok(self.view.channel_name(channel).map(|s| s.to_owned())?)
}
/// gives a helpful summary of the recorded experiment
fn summary(&self) -> PyResult<String> {
Ok(self.view.summary()?)
ndbioimage/src/py/transforms.rs
+96 -11
View File
@@ -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<usize> {
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<usize>) {
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<usize>) -> 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<Bound<'py, PyAny>>,
bead_files: Option<Bound<'py, PyAny>>,
file: Option<PathBuf>,
bead_files: Option<Vec<PathBuf>>,
main_channel: Option<usize>,
default_transform: Option<Vec<f64>>,
) -> PyResult<PyView> {
@@ -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<Bound<'py, PyAny>>,
bead_files: Option<Bound<'py, PyAny>>,
file: Option<PathBuf>,
bead_files: Option<Vec<PathBuf>>,
main_channel: Option<usize>,
default_transform: Option<Vec<f64>>,
) -> PyResult<PyView> {
todo!()
let mut new = self.clone();
let bead_files = bead_files
.as_ref()
.map(|i| i.iter().map(|j| j.as_path()).collect::<Vec<_>>());
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<String, PyTransform>,
drift: Vec<PyTransform>,
interpolation_order: Option<usize>,
) -> PyResult<PyTransforms> {
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<usize> {
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<usize>) {
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<usize>) -> Self {
let mut new = self.clone();
new.inner.interpolation_order = order;
new
}
#[getter]
fn get_channel_transforms(&self) -> HashMap<String, PyTransform> {
self.inner
ndbioimage/src/readers.rs
+1 -1
View File
@@ -201,7 +201,7 @@ pub trait Reader: Clone + Sized + Debug + Send + Hash + Into<DynReader> {
fn view(&self) -> View<Ix5, Self> {
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],
ndbioimage/src/readers/czi.rs
+2 -10
View File
@@ -122,14 +122,6 @@ impl From<Option<String>> 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<P>(path: P, series: usize, position: usize) -> Result<Self, Error>
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)
}
ndbioimage/src/transforms.rs
+143
View File
@@ -29,6 +29,7 @@ impl TryFrom<HashMap<String, TransformOld>> for Transforms {
fn try_from(t: HashMap<String, TransformOld>) -> Result<Self, Error> {
Ok(Self {
interpolation_order: None,
channel: t
.into_iter()
.map(|(n, i)| {
@@ -51,6 +52,8 @@ impl TryFrom<HashMap<String, TransformOld>> for Transforms {
#[derive(Clone, Debug, Default, Eq, PartialEq, Deserialize)]
pub struct Transforms {
#[serde(default)]
pub(crate) interpolation_order: Option<usize>,
#[serde(default)]
pub channel: HashMap<String, Transform<IxDyn>>,
#[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<usize>,
default_transform: Option<Vec<f64>>,
) -> Result<Self, Error> {
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<Vec<Transform<Ix2>>, Error> {
View::<_, DynReader>::from_path(file)?.calculate_drift_transforms_2d()
}
pub fn calculate_drift_transforms_3d(file: &Path) -> Result<Vec<Transform<Ix3>>, 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<Transform<Ix2>>,
) -> Result<HashMap<String, Transform<Ix2>>, 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<D: Dimension, R: Reader> View<D, R> {
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, Error> {
self.set_transforms(Transforms::load(path)?);
Ok(self)
@@ -550,7 +603,21 @@ impl<D: Dimension, R: Reader> View<D, R> {
#[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<DynReader, Error> {
let path = std::env::current_dir()?
.join("tests")
.join("files")
.join(file);
DynReader::new(&path, 0, 0)
}
#[test]
fn test_transforms() -> Result<(), Box<dyn std::error::Error>> {
@@ -576,4 +643,80 @@ mod tests {
assert_eq!(t.channel["TV2"].ndim(), 2);
Ok(())
}
#[test]
fn transforms_as_array() -> Result<(), Box<dyn std::error::Error>> {
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::<u16>()?;
let b = s_b.as_array::<u16>()?;
assert_ne!(a, b);
// let a = view.slice_cztyx(s![1, 0, 0, .., ..])?.as_array::<u16>()?;
// let b = view_t.slice_cztyx(s![1, 0, 0, .., ..])?.as_array::<u16>()?;
// assert_ne!(a, b);
//
// let a = view
// .max_proj(Axis::Z)?
// .slice(s![1, 0, .., ..])?
// .as_array::<u16>()?;
// let b = view_t
// .max_proj(Axis::Z)?
// .slice(s![1, 0, .., ..])?
// .as_array::<u16>()?;
// assert_ne!(a, b);
Ok(())
}
}
ndbioimage/src/view.rs
+237 -88
View File
@@ -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<D: Dimension, R: Reader = DynReader> {
operations: IndexMap<Axis, Operation>,
#[serde(skip)]
channel_names: OnceCell<HashMap<usize, String>>,
dimensionality: PhantomData<D>,
#[cfg(feature = "transforms")]
pub(crate) transforms: Transforms,
#[serde(skip)]
dimensionality: PhantomData<D>,
}
impl<D, R> Hash for View<D, R>
@@ -134,16 +136,20 @@ where
}
impl<D: Dimension, R: Reader> View<D, R> {
pub(crate) fn new(reader: R, slice: Vec<SliceInfoElem>, axes: Vec<Axis>) -> Self {
pub(crate) fn new_with_slice_and_axes(
reader: R,
slice: Vec<SliceInfoElem>,
axes: Vec<Axis>,
) -> 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<D: Dimension, R: Reader> View<D, R> {
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<D2: Dimension>(
&self,
slice: Vec<SliceInfoElem>,
axes: Vec<Axis>,
) -> View<D2, R> {
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<D2: Dimension>(
&self,
slice: Vec<SliceInfoElem>,
axes: Vec<Axis>,
operations: IndexMap<Axis, Operation>,
) -> View<D2, R> {
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<D: Dimension, R: Reader> View<D, R> {
ArrayCache::global().capacity()
}
fn with_operations(mut self, operations: IndexMap<Axis, Operation>) -> Self {
self.operations = operations;
self
}
/// change the dimension into a dynamic dimension
pub fn into_dyn(self) -> View<IxDyn, R> {
View {
@@ -282,7 +318,7 @@ impl<D: Dimension, R: Reader> View<D, R> {
channel_names: self.channel_names,
dimensionality: PhantomData,
#[cfg(feature = "transforms")]
transforms: Transforms::default(),
transforms: Default::default(),
}
}
@@ -438,7 +474,7 @@ impl<D: Dimension, R: Reader> View<D, R> {
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<D: Dimension, R: Reader> View<D, R> {
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<D: Dimension, R: Reader> View<D, R> {
}
}
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<Self, Error> {
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<D: Dimension, R: Reader> View<D, R> {
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<D: Dimension, R: Reader> View<D, R> {
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<D: Dimension, R: Reader> View<D, R> {
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,12 +791,56 @@ impl<D: Dimension, R: Reader> View<D, R> {
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) {
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");
@@ -786,18 +862,13 @@ impl<D: Dimension, R: Reader> View<D, R> {
}
}
}
}
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<D: Dimension, R: Reader> View<D, R> {
(_, SliceInfoElem::NewAxis) => panic!("axis cannot be a new axis"),
}
}
let xy = [
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::<T>::from_elem(shape_out.into_dimension(), T::min_value())
@@ -836,13 +917,6 @@ impl<D: Dimension, R: Reader> View<D, R> {
}
}
#[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<D: Dimension, R: Reader> View<D, R> {
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<D: Dimension, R: Reader> View<D, R> {
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<D: Dimension, R: Reader> View<D, R> {
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<T, D2>(
&self,
transform_done: &mut HashSet<(isize, isize, isize)>,
array: &mut Array<T, D2>,
slice: &[Slice],
ax: &HashMap<Axis, usize>,
defer_slice_xy: bool,
order: usize,
) -> Result<(), Error>
where
T: Number,
D2: Dimension,
[SliceInfoElem]: SliceArg<D2>,
{
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::<Vec<_>>();
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<T>(&self) -> Result<Array1<T>, 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<DynReader, Error> {
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::<u16, _>(&k0).is_some());
assert!(cache.get::<u16, _>(&k1).is_none());
assert_eq!(cache.len(), 2);
let ops: IndexMap<Axis, Operation> = IndexMap::new();
let ops: IndexMap<Axis, Operation> = 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::<f64, _>(&borrowed).is_none());
assert!(cache.get::<u16, _>(&borrowed).is_some());