diff --git a/src/bspline.rs b/src/bspline.rs index ef926f9..0e36211 100644 --- a/src/bspline.rs +++ b/src/bspline.rs @@ -798,6 +798,8 @@ where } /// eq 2.8 The L2 Polynomial Spline Pyramid, Unser et al. 1993 + /// may be used for future improvements + #[expect(dead_code)] fn u<'a, A, B>(m: usize, z: A) -> Array, B> where A: AsArray<'a, Complex, B>, @@ -816,6 +818,8 @@ where } /// table III The L2 Polynomial Spline Pyramid, Unser et al. 1993 + /// may be used for future improvements + #[expect(dead_code)] pub(crate) fn reduction_filter(shape: &[usize]) -> Result, D>, Error> { let mut filter: Array, D> = ArrayD::ones(shape).into_dimensionality()?; for (i, &s) in shape.iter().enumerate() { diff --git a/src/error.rs b/src/error.rs index 189974d..33f198c 100644 --- a/src/error.rs +++ b/src/error.rs @@ -10,4 +10,7 @@ pub enum Error { ShapeError(#[from] ndarray::ShapeError), #[error("number of dimensions is not defined")] NumberOfDimensionsNotDefined, + /// the dimensionality of the data does not match + #[error("dimensionality mismatch: {0} != {0}")] + DimensionalityMismatch(usize, usize), } diff --git a/src/metric.rs b/src/metric.rs index 813b784..6a86825 100644 --- a/src/metric.rs +++ b/src/metric.rs @@ -75,7 +75,7 @@ impl FixedMu { /// first_bin, last_bin: center of bin /// n_bins >= 2 -#[allow(clippy::too_many_arguments)] +#[expect(clippy::too_many_arguments)] fn parzen( alpha: f64, dalpha: &[f64], diff --git a/src/transform.rs b/src/transform.rs index fde35d8..ebb9e74 100644 --- a/src/transform.rs +++ b/src/transform.rs @@ -3,13 +3,13 @@ use crate::error::Error; use crate::metric::FixedMu; use crate::register::{Registration, RegistrationResult, RegistrationStep}; use itertools::Itertools; -use ndarray::{Array, Array2, ArrayD, AsArray, Dimension, Ix2, s}; +use ndarray::{Array, Array2, ArrayD, ArrayView2, AsArray, Dimension, Ix2, IxDyn, s}; use num::cast::AsPrimitive; use serde::{Deserialize, Serialize}; use serde_yaml::{from_reader, to_writer}; use std::fs::File; use std::marker::PhantomData; -use std::ops::Mul; +use std::ops::{Add, Div, Mul, Sub}; use std::path::PathBuf; /// get coordinates resulting from transforming input coordinates, coordinate must have N @@ -143,7 +143,7 @@ impl Mul> for &Transform { impl Mul<&Transform> for &Transform { type Output = Transform; - #[allow(clippy::suspicious_arithmetic_impl)] + #[expect(clippy::suspicious_arithmetic_impl)] fn mul(self, rhs: &Transform) -> Self::Output { let m = self.matrix().dot(&rhs.matrix()); let dm = self.dmatrix().dot(&rhs.matrix()) + self.matrix().dot(&rhs.dmatrix()); @@ -170,6 +170,202 @@ impl Mul<&Transform> for &Transform { } } +impl Add for Transform { + type Output = Transform; + + fn add(self, rhs: Self) -> Self::Output { + &self + &rhs + } +} + +impl Add<&Transform> for Transform { + type Output = Transform; + + fn add(self, rhs: &Transform) -> Self::Output { + &self + rhs + } +} + +impl Add> for &Transform { + type Output = Transform; + + fn add(self, rhs: Transform) -> Self::Output { + self + &rhs + } +} + +impl Add<&Transform> for &Transform { + type Output = Transform; + + fn add(self, rhs: &Transform) -> Self::Output { + let m = self.matrix() + &rhs.matrix(); + let dm = (self.dmatrix().powi(2) + rhs.matrix().powi(2)).sqrt(); + Self::Output { + parameters: m + .slice(s![..self.ndim, ..self.ndim]) + .flatten() + .iter() + .chain(m.slice(s![..self.ndim, self.ndim]).iter()) + .cloned() + .collect(), + dparameters: dm + .slice(s![..self.ndim, ..]) + .flatten() + .iter() + .chain(dm.slice(s![..self.ndim, ..]).iter()) + .cloned() + .collect(), + center: self.center.clone(), + shape: self.shape.clone(), + ndim: self.ndim, + dimension: self.dimension, + } + } +} + +impl Sub for Transform { + type Output = Transform; + + fn sub(self, rhs: Self) -> Self::Output { + &self - &rhs + } +} + +impl Sub<&Transform> for Transform { + type Output = Transform; + + fn sub(self, rhs: &Transform) -> Self::Output { + &self - rhs + } +} + +impl Sub> for &Transform { + type Output = Transform; + + fn sub(self, rhs: Transform) -> Self::Output { + self - &rhs + } +} + +impl Sub<&Transform> for &Transform { + type Output = Transform; + + fn sub(self, rhs: &Transform) -> Self::Output { + let m = self.matrix() - &rhs.matrix(); + let dm = (self.dmatrix().powi(2) + rhs.matrix().powi(2)).sqrt(); + Self::Output { + parameters: m + .slice(s![..self.ndim, ..self.ndim]) + .flatten() + .iter() + .chain(m.slice(s![..self.ndim, self.ndim]).iter()) + .cloned() + .collect(), + dparameters: dm + .slice(s![..self.ndim, ..]) + .flatten() + .iter() + .chain(dm.slice(s![..self.ndim, ..]).iter()) + .cloned() + .collect(), + center: self.center.clone(), + shape: self.shape.clone(), + ndim: self.ndim, + dimension: self.dimension, + } + } +} + +impl Mul> for f64 { + type Output = Transform; + + fn mul(self, rhs: Transform) -> Self::Output { + rhs * self + } +} + +impl Mul<&Transform> for f64 { + type Output = Transform; + + fn mul(self, rhs: &Transform) -> Self::Output { + rhs * self + } +} + +impl Mul for Transform { + type Output = Transform; + + fn mul(self, rhs: f64) -> Self::Output { + &self * rhs + } +} + +impl Mul for &Transform { + type Output = Transform; + + fn mul(self, rhs: f64) -> Self::Output { + let m = self.matrix() * rhs; + let dm = self.dmatrix() * rhs; + Self::Output { + parameters: m + .slice(s![..self.ndim, ..self.ndim]) + .flatten() + .iter() + .chain(m.slice(s![..self.ndim, self.ndim]).iter()) + .cloned() + .collect(), + dparameters: dm + .slice(s![..self.ndim, ..]) + .flatten() + .iter() + .chain(dm.slice(s![..self.ndim, ..]).iter()) + .cloned() + .collect(), + center: self.center.clone(), + shape: self.shape.clone(), + ndim: self.ndim, + dimension: self.dimension, + } + } +} + +impl Div for Transform { + type Output = Transform; + + fn div(self, rhs: f64) -> Self::Output { + &self / rhs + } +} + +impl Div for &Transform { + type Output = Transform; + + fn div(self, rhs: f64) -> Self::Output { + let m = self.matrix() / rhs; + let dm = self.dmatrix() / rhs; + Self::Output { + parameters: m + .slice(s![..self.ndim, ..self.ndim]) + .flatten() + .iter() + .chain(m.slice(s![..self.ndim, self.ndim]).iter()) + .cloned() + .collect(), + dparameters: dm + .slice(s![..self.ndim, ..]) + .flatten() + .iter() + .chain(dm.slice(s![..self.ndim, ..]).iter()) + .cloned() + .collect(), + center: self.center.clone(), + shape: self.shape.clone(), + ndim: self.ndim, + dimension: self.dimension, + } + } +} + impl Eq for Transform {} impl Default for Transform { @@ -481,6 +677,44 @@ impl Transform { matrix } + pub fn with_matrix(mut self, matrix: ArrayView2) -> Self { + let mut parameters = matrix + .slice(s![..self.ndim, ..self.ndim]) + .flatten() + .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]) + .flatten() + .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]) + .flatten() + .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]) + .flatten() + .to_vec(); + dparameters.extend(dmatrix.slice(s![..self.ndim, self.ndim]).to_vec()); + self.dparameters = dparameters; + } + /// get the inverse transform pub fn inverse(&self) -> Result { let matrix = self.matrix(); @@ -540,6 +774,48 @@ impl Transform { let bspline = BSpline::::new(&image); bspline.interpolate_par(self) } + + /// the number of dimensions in the view + pub fn ndim(&self) -> usize { + if let Some(d) = D::NDIM { d } else { self.ndim } + } + + pub fn into_dimensionality(self) -> Result, Error> { + if let Some(d) = D2::NDIM { + if d == self.ndim() { + Ok(Transform { + parameters: self.parameters, + dparameters: self.dparameters, + center: self.center, + shape: self.shape, + ndim: self.ndim, + dimension: Default::default(), + }) + } else { + Err(Error::DimensionalityMismatch(d, self.ndim())) + } + } else { + Ok(Transform { + parameters: self.parameters, + dparameters: self.dparameters, + center: self.center, + shape: self.shape, + ndim: self.ndim, + dimension: Default::default(), + }) + } + } + + pub fn into_dyn(self) -> Transform { + Transform { + parameters: self.parameters, + dparameters: self.dparameters, + center: self.center, + shape: self.shape, + ndim: self.ndim, + dimension: Default::default(), + } + } } impl Transform {