diff --git a/README.md b/README.md index ff74007..022e7a4 100644 --- a/README.md +++ b/README.md @@ -336,7 +336,7 @@ fn main() -> Result<(), Error> { | Module | Contents | |----------------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------| -| `image_registration::transform` | `Transform`, the affine transform type: construction (`new`, `new_with_center`, `from_translation`, `from_scaling`, `from_rotation`), algebra (`inverse`, `with_translation`, `with_scaling`, `with_rotation`, `with_rotation_around`, `Mul`), coordinate transforms (`transform_point`, `transform_points`, `matrix`, `dmatrix`), resampling (`interpolate`, `interpolate_par`), registration (`register`, `register_debug`, `register_affine`, `register_translation`), serialization (`to_file`, `from_file`), plus the free function `transform_point` and a blas-free `matrix_inverse`. | +| `image_registration::transform` | `Transform`, the affine transform type: construction (`new`, `new_with_center`, `from_translation`, `from_scaling`, `from_rotation`), algebra (`inverse`, `add_translation`, `add_scaling`, `with_rotation`, `with_rotation_around`, `Mul`), coordinate transforms (`transform_point`, `transform_points`, `matrix`, `dmatrix`), resampling (`interpolate`, `interpolate_par`), registration (`register`, `register_debug`, `register_affine`, `register_translation`), serialization (`to_file`, `from_file`), plus the free function `transform_point` and a blas-free `matrix_inverse`. | | `image_registration::register` | `Registration`, `RegistrationStep`, `RegistrationResult`, `Optimizer` (LBFGS/ASGD). | | `image_registration::metric` | `MattesMetric` (Mattes mutual information with analytic gradient), `SamplingArg` (Fixed/FixedAt/Random), `Sigma` (None/Absolute/Relative), `FixedMu`. | | `image_registration::bspline` | `BSpline` B-spline interpolation of arbitrary order, `BSplineTrait`, `BSplineMem`, Parzen kernel helpers. | diff --git a/src/error.rs b/src/error.rs index d1f66f7..6796bc4 100644 --- a/src/error.rs +++ b/src/error.rs @@ -17,4 +17,6 @@ pub enum Error { FixedMuLength(usize, usize), #[error("expected initial guess length of {0}, got {1}")] InitialGuessLength(usize, usize), + #[error("interpolation order {0} is not implemented")] + InterpolationNotImplemented(usize), } diff --git a/src/transform.rs b/src/transform.rs index 4a62396..d377109 100644 --- a/src/transform.rs +++ b/src/transform.rs @@ -8,6 +8,7 @@ use num::cast::AsPrimitive; use serde::{Deserialize, Serialize}; use serde_yaml::{from_reader, to_writer}; use std::fs::File; +use std::hash::{Hash, Hasher}; use std::marker::PhantomData; use std::ops::{Add, Div, Mul, Sub}; use std::path::PathBuf; @@ -110,13 +111,11 @@ struct TransformS { impl From for Transform { fn from(t: TransformS) -> Self { - let ndim = t.shape.len(); Self { parameters: t.parameters, dparameters: t.dparameters, center: t.center, shape: t.shape, - ndim, dimension: PhantomData, } } @@ -136,8 +135,6 @@ pub struct Transform { /// the shape of images for which this transform is meant pub shape: Vec, #[serde(skip)] - ndim: usize, - #[serde(skip)] dimension: PhantomData, } @@ -174,22 +171,21 @@ impl Mul<&Transform> for &Transform { let dm = self.dmatrix().dot(&rhs.matrix()) + self.matrix().dot(&rhs.dmatrix()); Self::Output { parameters: m - .slice(s![..self.ndim, ..self.ndim]) + .slice(s![..self.ndim(), ..self.ndim()]) .flatten() .iter() - .chain(m.slice(s![..self.ndim, self.ndim]).iter()) + .chain(m.slice(s![..self.ndim(), self.ndim()]).iter()) .cloned() .collect(), dparameters: dm - .slice(s![..self.ndim, ..]) + .slice(s![..self.ndim(), ..]) .flatten() .iter() - .chain(dm.slice(s![..self.ndim, ..]).iter()) + .chain(dm.slice(s![..self.ndim(), ..]).iter()) .cloned() .collect(), center: self.center.clone(), shape: self.shape.clone(), - ndim: self.ndim, dimension: self.dimension, } } @@ -227,22 +223,21 @@ impl Add<&Transform> for &Transform { let dm = (self.dmatrix().powi(2) + rhs.matrix().powi(2)).sqrt(); Self::Output { parameters: m - .slice(s![..self.ndim, ..self.ndim]) + .slice(s![..self.ndim(), ..self.ndim()]) .flatten() .iter() - .chain(m.slice(s![..self.ndim, self.ndim]).iter()) + .chain(m.slice(s![..self.ndim(), self.ndim()]).iter()) .cloned() .collect(), dparameters: dm - .slice(s![..self.ndim, ..]) + .slice(s![..self.ndim(), ..]) .flatten() .iter() - .chain(dm.slice(s![..self.ndim, ..]).iter()) + .chain(dm.slice(s![..self.ndim(), ..]).iter()) .cloned() .collect(), center: self.center.clone(), shape: self.shape.clone(), - ndim: self.ndim, dimension: self.dimension, } } @@ -280,22 +275,21 @@ impl Sub<&Transform> for &Transform { let dm = (self.dmatrix().powi(2) + rhs.matrix().powi(2)).sqrt(); Self::Output { parameters: m - .slice(s![..self.ndim, ..self.ndim]) + .slice(s![..self.ndim(), ..self.ndim()]) .flatten() .iter() - .chain(m.slice(s![..self.ndim, self.ndim]).iter()) + .chain(m.slice(s![..self.ndim(), self.ndim()]).iter()) .cloned() .collect(), dparameters: dm - .slice(s![..self.ndim, ..]) + .slice(s![..self.ndim(), ..]) .flatten() .iter() - .chain(dm.slice(s![..self.ndim, ..]).iter()) + .chain(dm.slice(s![..self.ndim(), ..]).iter()) .cloned() .collect(), center: self.center.clone(), shape: self.shape.clone(), - ndim: self.ndim, dimension: self.dimension, } } @@ -333,22 +327,21 @@ impl Mul for &Transform { let dm = self.dmatrix() * rhs; Self::Output { parameters: m - .slice(s![..self.ndim, ..self.ndim]) + .slice(s![..self.ndim(), ..self.ndim()]) .flatten() .iter() - .chain(m.slice(s![..self.ndim, self.ndim]).iter()) + .chain(m.slice(s![..self.ndim(), self.ndim()]).iter()) .cloned() .collect(), dparameters: dm - .slice(s![..self.ndim, ..]) + .slice(s![..self.ndim(), ..]) .flatten() .iter() - .chain(dm.slice(s![..self.ndim, ..]).iter()) + .chain(dm.slice(s![..self.ndim(), ..]).iter()) .cloned() .collect(), center: self.center.clone(), shape: self.shape.clone(), - ndim: self.ndim, dimension: self.dimension, } } @@ -370,22 +363,21 @@ impl Div for &Transform { let dm = self.dmatrix() / rhs; Self::Output { parameters: m - .slice(s![..self.ndim, ..self.ndim]) + .slice(s![..self.ndim(), ..self.ndim()]) .flatten() .iter() - .chain(m.slice(s![..self.ndim, self.ndim]).iter()) + .chain(m.slice(s![..self.ndim(), self.ndim()]).iter()) .cloned() .collect(), dparameters: dm - .slice(s![..self.ndim, ..]) + .slice(s![..self.ndim(), ..]) .flatten() .iter() - .chain(dm.slice(s![..self.ndim, ..]).iter()) + .chain(dm.slice(s![..self.ndim(), ..]).iter()) .cloned() .collect(), center: self.center.clone(), shape: self.shape.clone(), - ndim: self.ndim, dimension: self.dimension, } } @@ -393,6 +385,19 @@ impl Div for &Transform { impl Eq for Transform {} +impl Hash for Transform { + fn hash(&self, state: &mut H) { + for values in [&self.parameters, &self.dparameters, &self.center] { + values.len().hash(state); + for v in values { + let v = if *v == 0.0 { 0.0 } else { *v }; + v.to_bits().hash(state); + } + } + self.shape.hash(state); + } +} + impl Default for Transform { /// the unit transform fn default() -> Self { @@ -407,7 +412,6 @@ impl Default for Transform { dparameters: vec![0f64; ndim * ndim + ndim], center: vec![0f64; ndim], shape: vec![0usize; ndim], - ndim, dimension: PhantomData, } } @@ -429,7 +433,6 @@ impl Transform { dparameters: vec![0.0; ndim * ndim + ndim], center, shape, - ndim, dimension: PhantomData, } } @@ -448,7 +451,6 @@ impl Transform { dparameters: vec![0.0; ndim * ndim + ndim], center, shape, - ndim, dimension: PhantomData, } } @@ -548,7 +550,6 @@ impl Transform { dparameters: vec![0f64; ndim * ndim + ndim], center: vec![0f64; ndim], shape: vec![0usize; ndim], - ndim, dimension: PhantomData, } } @@ -581,19 +582,79 @@ impl Transform { dparameters: vec![0f64; ndim * ndim + ndim], center: vec![0f64; ndim], shape: vec![0usize; ndim], - ndim, dimension: PhantomData, } } + /// replace the transform by a scaling + pub fn as_scaling(&self, scaling: &[f64]) -> Self { + let ndim = if let Some(ndim) = D::NDIM { + ndim + } else { + scaling.len() + }; + debug_assert_eq!(scaling.len(), ndim); + let mut parameters = vec![0.0; ndim * ndim + ndim]; + parameters[ndim * ndim + ndim - 1] = 1.0; + for (p, s) in parameters + .iter_mut() + .step_by(ndim + 1) + .zip_eq(scaling.iter()) + { + *p = *s; + } + + Self { + parameters, + dparameters: vec![0f64; ndim * ndim + ndim], + center: self.center.clone(), + shape: self.shape.clone(), + dimension: PhantomData, + } + } + + /// scale self + pub fn add_scaling(&self, scaling: &[f64]) -> Self { + Self::from_scaling(scaling) * self + } + /// add a translation to self pub fn with_translation(&self, translation: &[f64]) -> Self { Self::from_translation(translation) * self } - /// scale self - pub fn with_scaling(&self, scaling: &[f64]) -> Self { - Self::from_scaling(scaling) * self + /// replace the transform by a translation + pub fn as_translation(&self, translation: &[f64]) -> Self { + let ndim = if let Some(ndim) = D::NDIM { + ndim + } else { + translation.len() + }; + debug_assert_eq!(translation.len(), ndim); + let mut parameters = vec![0.0; ndim * ndim + ndim]; + for i in 0..ndim { + parameters[(ndim + 1) * i] = 1.0; + } + for (p, t) in parameters + .iter_mut() + .skip(ndim * ndim) + .zip_eq(translation.iter()) + { + *p = *t; + } + + Self { + parameters, + dparameters: vec![0f64; ndim * ndim + ndim], + center: self.center.clone(), + shape: self.shape.clone(), + dimension: PhantomData, + } + } + + /// add a translation to self + pub fn add_translation(&self, translation: &[f64]) -> Self { + Self::from_translation(translation) * self } pub fn with_center(mut self, center: Vec) -> Self { @@ -620,8 +681,8 @@ impl Transform { /// true if transform does nothing pub fn is_unity(&self) -> bool { - let n = self.ndim * self.ndim; - let m = self.ndim + 1; + let n = self.ndim() * self.ndim(); + let m = self.ndim() + 1; self.parameters .iter() .take(n) @@ -667,17 +728,17 @@ impl Transform { /// get the matrix defining the transform pub fn matrix(&self) -> Array2 { - let n = self.ndim * self.ndim; - let mut matrix = Array2::eye(self.ndim + 1); + let n = self.ndim() * self.ndim(); + let mut matrix = Array2::eye(self.ndim() + 1); for (m, p) in matrix - .slice_mut(s![..self.ndim, ..self.ndim]) + .slice_mut(s![..self.ndim(), ..self.ndim()]) .iter_mut() .zip_eq(self.parameters.iter().take(n)) { *m = *p; } for (m, p) in matrix - .slice_mut(s![..self.ndim, self.ndim]) + .slice_mut(s![..self.ndim(), self.ndim()]) .iter_mut() .zip_eq(self.parameters.iter().skip(n)) { @@ -688,17 +749,17 @@ impl Transform { /// get the matrix describing the error of the transform pub fn dmatrix(&self) -> Array2 { - let n = self.ndim * self.ndim; - let mut matrix = Array2::zeros([self.ndim + 1, self.ndim + 1]); + let n = self.ndim() * self.ndim(); + let mut matrix = Array2::zeros([self.ndim() + 1, self.ndim() + 1]); for (m, p) in matrix - .slice_mut(s![..self.ndim, ..self.ndim]) + .slice_mut(s![..self.ndim(), ..self.ndim()]) .iter_mut() .zip_eq(self.dparameters.iter().take(n)) { *m = *p; } for (m, p) in matrix - .slice_mut(s![..self.ndim, self.ndim]) + .slice_mut(s![..self.ndim(), self.ndim()]) .iter_mut() .zip_eq(self.dparameters.iter().skip(n)) { @@ -719,39 +780,39 @@ impl Transform { pub fn with_matrix(mut self, matrix: ArrayView2) -> Self { let mut parameters = matrix - .slice(s![..self.ndim, ..self.ndim]) + .slice(s![..self.ndim(), ..self.ndim()]) .flatten() .to_vec(); - parameters.extend(matrix.slice(s![..self.ndim, self.ndim]).to_vec()); + parameters.extend(matrix.slice(s![..self.ndim(), self.ndim()]).to_vec()); self.parameters = parameters; self } pub fn set_matrix(&mut self, matrix: ArrayView2) { let mut parameters = matrix - .slice(s![..self.ndim, ..self.ndim]) + .slice(s![..self.ndim(), ..self.ndim()]) .flatten() .to_vec(); - parameters.extend(matrix.slice(s![..self.ndim, self.ndim]).to_vec()); + parameters.extend(matrix.slice(s![..self.ndim(), self.ndim()]).to_vec()); self.parameters = parameters; } pub fn with_dmatrix(mut self, dmatrix: ArrayView2) -> Self { let mut dparameters = dmatrix - .slice(s![..self.ndim, ..self.ndim]) + .slice(s![..self.ndim(), ..self.ndim()]) .flatten() .to_vec(); - dparameters.extend(dmatrix.slice(s![..self.ndim, self.ndim]).to_vec()); + dparameters.extend(dmatrix.slice(s![..self.ndim(), self.ndim()]).to_vec()); self.dparameters = dparameters; self } pub fn set_dmatrix(&mut self, dmatrix: ArrayView2) { let mut dparameters = dmatrix - .slice(s![..self.ndim, ..self.ndim]) + .slice(s![..self.ndim(), ..self.ndim()]) .flatten() .to_vec(); - dparameters.extend(dmatrix.slice(s![..self.ndim, self.ndim]).to_vec()); + dparameters.extend(dmatrix.slice(s![..self.ndim(), self.ndim()]).to_vec()); self.dparameters = dparameters; } @@ -760,18 +821,17 @@ impl Transform { let matrix = self.matrix(); let inverse = matrix_inverse(matrix); let parameters = inverse - .slice(s![..self.ndim, ..self.ndim]) + .slice(s![..self.ndim(), ..self.ndim()]) .iter() - .chain(inverse.slice(s![..self.ndim, self.ndim])) + .chain(inverse.slice(s![..self.ndim(), self.ndim()])) .cloned() .collect(); Ok(Self { parameters, - dparameters: vec![0f64; self.ndim * self.ndim + self.ndim], + dparameters: vec![0f64; self.ndim() * self.ndim() + self.ndim()], center: self.center.clone(), shape: self.shape.clone(), - ndim: self.ndim, dimension: PhantomData, }) } @@ -815,9 +875,55 @@ impl Transform { bspline.interpolate_par(self) } + pub fn interpolate_order<'a, T, A>( + &self, + image: A, + order: usize, + ) -> Result, Error> + where + A: AsArray<'a, T, D>, + D: Dimension, + T: 'a + Clone + AsPrimitive, + { + match order { + 0 => Ok(self.interpolate::<0, _, _>(image)?), + 1 => Ok(self.interpolate::<1, _, _>(image)?), + 2 => Ok(self.interpolate::<2, _, _>(image)?), + 3 => Ok(self.interpolate::<3, _, _>(image)?), + 4 => Ok(self.interpolate::<4, _, _>(image)?), + 5 => Ok(self.interpolate::<5, _, _>(image)?), + _ => Err(Error::InterpolationNotImplemented(order)), + } + } + + pub fn interpolate_par_order<'a, T, A>( + &self, + image: A, + order: usize, + ) -> Result, Error> + where + A: AsArray<'a, T, D>, + D: Dimension, + T: 'a + Clone + AsPrimitive, + { + match order { + 0 => Ok(self.interpolate_par::<0, _, _>(image)?), + 1 => Ok(self.interpolate_par::<1, _, _>(image)?), + 2 => Ok(self.interpolate_par::<2, _, _>(image)?), + 3 => Ok(self.interpolate_par::<3, _, _>(image)?), + 4 => Ok(self.interpolate_par::<4, _, _>(image)?), + 5 => Ok(self.interpolate_par::<5, _, _>(image)?), + _ => Err(Error::InterpolationNotImplemented(order)), + } + } + /// the number of dimensions in the view pub fn ndim(&self) -> usize { - if let Some(d) = D::NDIM { d } else { self.ndim } + if let Some(d) = D::NDIM { + d + } else { + self.center.len() + } } pub fn into_dimensionality(self) -> Result, Error> { @@ -828,7 +934,6 @@ impl Transform { dparameters: self.dparameters, center: self.center, shape: self.shape, - ndim: self.ndim, dimension: Default::default(), }) } else { @@ -840,7 +945,6 @@ impl Transform { dparameters: self.dparameters, center: self.center, shape: self.shape, - ndim: self.ndim, dimension: Default::default(), }) } @@ -852,7 +956,44 @@ impl Transform { dparameters: self.dparameters, center: self.center, shape: self.shape, - ndim: self.ndim, + dimension: Default::default(), + } + } + + fn parameters_unit(ndim: usize) -> Vec { + let mut p = vec![0.0; ndim * (ndim + 1)]; + for i in 0..ndim { + p[i * (ndim + 1)] = 1.0; + } + p + } + + pub fn add_dimension(&self, axis: usize, size: usize) -> Transform { + let ndim = self.ndim() + 1; + let mut parameters = Self::parameters_unit(ndim); + let mut dparameters = Self::parameters_unit(ndim); + let mut j = 0; + for p in &self.parameters { + if j == axis * ndim { + j += ndim; + } + if j % ndim == axis { + j += 1; + } + parameters[j] = *p; + dparameters[j] = *p; + } + + let mut center = self.center.clone(); + center.insert(axis, ((size as f64) - 1.0) / 2.0); + let mut shape = self.shape.clone(); + shape.insert(axis, size); + + Transform { + parameters, + dparameters, + center, + shape, dimension: Default::default(), } } @@ -873,7 +1014,6 @@ impl Transform { dparameters: vec![0f64; 6], center: center.to_vec(), shape: vec![0usize; 2], - ndim: 2, dimension: PhantomData, } } @@ -893,13 +1033,13 @@ impl Transform { mod tests { use super::matrix_inverse; use crate::julia_image; + use crate::transform::File; use crate::transform::Transform; use itertools::Itertools; use ndarray::{Array2, Ix2, array, s}; use num::traits::FloatConst; - use std::fs::File; use std::path::Path; - use tiff::decoder::{Decoder, DecodingResult}; + use tiff::decoder::{Decoder, DecodingSampleBuffer}; use tiff::tags::Tag; #[test] @@ -1521,23 +1661,42 @@ mod tests { } fn read_tiff>(path: P) -> Result, Box> { - let mut reader = Decoder::new(File::open(path)?)?; + let mut reader = Decoder::open(File::open(path)?)?; reader.seek_to_image(0)?; let bytes = match reader.read_image()? { - DecodingResult::U8(data) => data.into_iter().map(|i| i as f64).collect::>(), - DecodingResult::U16(data) => data.into_iter().map(|i| i as f64).collect::>(), - DecodingResult::U32(data) => data.into_iter().map(|i| i as f64).collect::>(), - DecodingResult::U64(data) => data.into_iter().map(|i| i as f64).collect::>(), - DecodingResult::I8(data) => data.into_iter().map(|i| i as f64).collect::>(), - DecodingResult::I16(data) => data.into_iter().map(|i| i as f64).collect::>(), - DecodingResult::I32(data) => data.into_iter().map(|i| i as f64).collect::>(), - DecodingResult::I64(data) => data.into_iter().map(|i| i as f64).collect::>(), - DecodingResult::F16(data) => data.into_iter().map(f64::from).collect::>(), - DecodingResult::F32(data) => data.into_iter().map(|i| i as f64).collect::>(), - DecodingResult::F64(data) => data, + DecodingSampleBuffer::U8(data) => { + data.into_iter().map(|i| i as f64).collect::>() + } + DecodingSampleBuffer::U16(data) => { + data.into_iter().map(|i| i as f64).collect::>() + } + DecodingSampleBuffer::U32(data) => { + data.into_iter().map(|i| i as f64).collect::>() + } + DecodingSampleBuffer::U64(data) => { + data.into_iter().map(|i| i as f64).collect::>() + } + DecodingSampleBuffer::I8(data) => { + data.into_iter().map(|i| i as f64).collect::>() + } + DecodingSampleBuffer::I16(data) => { + data.into_iter().map(|i| i as f64).collect::>() + } + DecodingSampleBuffer::I32(data) => { + data.into_iter().map(|i| i as f64).collect::>() + } + DecodingSampleBuffer::I64(data) => { + data.into_iter().map(|i| i as f64).collect::>() + } + DecodingSampleBuffer::F16(data) => data.into_iter().map(f64::from).collect::>(), + DecodingSampleBuffer::F32(data) => { + data.into_iter().map(|i| i as f64).collect::>() + } + DecodingSampleBuffer::F64(data) => data, }; - let width = reader.get_tag(Tag::ImageWidth)?.into_u32()? as usize; - let height = reader.get_tag(Tag::ImageLength)?.into_u32()? as usize; + let mut ifd = reader.current_ifd(); + let width = ifd.get_tag(Tag::ImageWidth)?.into_u32()? as usize; + let height = ifd.get_tag(Tag::ImageLength)?.into_u32()? as usize; Ok(Array2::from_shape_vec((width, height), bytes)?) }