- more arithmetics for Transform

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<D: Dimension> Mul<Transform<D>> for &Transform<D> {
 impl<D: Dimension> Mul<&Transform<D>> for &Transform<D> {
     type Output = Transform<D>;
 
-    #[allow(clippy::suspicious_arithmetic_impl)]
+    #[expect(clippy::suspicious_arithmetic_impl)]
     fn mul(self, rhs: &Transform<D>) -> 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<D: Dimension> Mul<&Transform<D>> for &Transform<D> {
     }
 }
 
+impl<D: Dimension> Add for Transform<D> {
+    type Output = Transform<D>;
+
+    fn add(self, rhs: Self) -> Self::Output {
+        &self + &rhs
+    }
+}
+
+impl<D: Dimension> Add<&Transform<D>> for Transform<D> {
+    type Output = Transform<D>;
+
+    fn add(self, rhs: &Transform<D>) -> Self::Output {
+        &self + rhs
+    }
+}
+
+impl<D: Dimension> Add<Transform<D>> for &Transform<D> {
+    type Output = Transform<D>;
+
+    fn add(self, rhs: Transform<D>) -> Self::Output {
+        self + &rhs
+    }
+}
+
+impl<D: Dimension> Add<&Transform<D>> for &Transform<D> {
+    type Output = Transform<D>;
+
+    fn add(self, rhs: &Transform<D>) -> 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<D: Dimension> Sub for Transform<D> {
+    type Output = Transform<D>;
+
+    fn sub(self, rhs: Self) -> Self::Output {
+        &self - &rhs
+    }
+}
+
+impl<D: Dimension> Sub<&Transform<D>> for Transform<D> {
+    type Output = Transform<D>;
+
+    fn sub(self, rhs: &Transform<D>) -> Self::Output {
+        &self - rhs
+    }
+}
+
+impl<D: Dimension> Sub<Transform<D>> for &Transform<D> {
+    type Output = Transform<D>;
+
+    fn sub(self, rhs: Transform<D>) -> Self::Output {
+        self - &rhs
+    }
+}
+
+impl<D: Dimension> Sub<&Transform<D>> for &Transform<D> {
+    type Output = Transform<D>;
+
+    fn sub(self, rhs: &Transform<D>) -> 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<D: Dimension> Mul<Transform<D>> for f64 {
+    type Output = Transform<D>;
+
+    fn mul(self, rhs: Transform<D>) -> Self::Output {
+        rhs * self
+    }
+}
+
+impl<D: Dimension> Mul<&Transform<D>> for f64 {
+    type Output = Transform<D>;
+
+    fn mul(self, rhs: &Transform<D>) -> Self::Output {
+        rhs * self
+    }
+}
+
+impl<D: Dimension> Mul<f64> for Transform<D> {
+    type Output = Transform<D>;
+
+    fn mul(self, rhs: f64) -> Self::Output {
+        &self * rhs
+    }
+}
+
+impl<D: Dimension> Mul<f64> for &Transform<D> {
+    type Output = Transform<D>;
+
+    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<D: Dimension> Div<f64> for Transform<D> {
+    type Output = Transform<D>;
+
+    fn div(self, rhs: f64) -> Self::Output {
+        &self / rhs
+    }
+}
+
+impl<D: Dimension> Div<f64> for &Transform<D> {
+    type Output = Transform<D>;
+
+    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<D: Dimension> Eq for Transform<D> {}
 
 impl<D: Dimension> Default for Transform<D> {
@@ -481,6 +677,44 @@ impl<D: Dimension> Transform<D> {
         matrix
     }
 
+    pub fn with_matrix(mut self, matrix: ArrayView2<f64>) -> 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<f64>) {
+        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<f64>) -> 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<f64>) {
+        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<Self, Error> {
         let matrix = self.matrix();
@@ -540,6 +774,48 @@ impl<D: Dimension> Transform<D> {
         let bspline = BSpline::<B, _>::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<D2: Dimension>(self) -> Result<Transform<D2>, 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<IxDyn> {
+        Transform {
+            parameters: self.parameters,
+            dparameters: self.dparameters,
+            center: self.center,
+            shape: self.shape,
+            ndim: self.ndim,
+            dimension: Default::default(),
+        }
+    }
 }
 
 impl Transform<Ix2> {