- bugfix in new_translation

- new register_translation_fft function

src/register.rs

@@ -175,7 +175,7 @@ impl<D: Dimension> Registration<D> {
         } else {
             return Err(Error::NumberOfDimensionsNotDefined);
         };
-        let mut fixed_mu = vec![Some(0.0); ndim * ndim];
+        let mut fixed_mu = vec![Some(0.0); ndim * ndim + ndim];
         for i in 0..ndim {
             fixed_mu[i * (ndim + 1)] = Some(1.0);
             fixed_mu[ndim * ndim + i] = None;

src/transform.rs

@@ -1,10 +1,14 @@
 use crate::bspline::{BSpline, BSplineTrait};
 use crate::error::Error;
+use crate::filter::{fft, ifft};
 use crate::metric::FixedMu;
 use crate::register::{Registration, RegistrationResult, RegistrationStep};
 use itertools::Itertools;
-use ndarray::{Array, Array2, ArrayD, ArrayView2, AsArray, Dimension, Ix2, IxDyn, s};
+use ndarray::{
+    Array, Array2, ArrayD, ArrayView2, AsArray, Dimension, IntoDimension, Ix2, IxDyn, s,
+};
 use num::cast::AsPrimitive;
+use num::complex::ComplexFloat;
 use serde::{Deserialize, Serialize};
 use serde_yaml::{from_reader, to_writer};
 use std::fs::File;
@@ -153,10 +157,10 @@ impl<D: Dimension> Mul<&Transform<D>> for &Transform<D> {
                 .cloned()
                 .collect(),
             dparameters: dm
-                .slice(s![..self.ndim(), ..])
+                .slice(s![..self.ndim(), ..self.ndim()])
                 .flatten()
                 .iter()
-                .chain(dm.slice(s![..self.ndim(), ..]).iter())
+                .chain(dm.slice(s![..self.ndim(), self.ndim()]).iter())
                 .cloned()
                 .collect(),
             dimension: self.dimension,
@@ -399,8 +403,8 @@ impl<D: Dimension> Transform<D> {
     }
 
     pub fn new_around_center(parameters: Vec<f64>, center: &[f64]) -> Self {
-        (Self::new(vec![1.0, 0.0, 0.0, 1.0, -center[0], -center[1]]) * Self::new(parameters))
+        Self::new(parameters).with_translation(&center.iter().map(|i| -i).collect::<Vec<_>>())
-            * Self::new(vec![1.0, 0.0, 0.0, 1.0, center[0], center[1]])
+            * Self::from_translation(center)
     }
 
     /// find the affine transform which transforms moving into fixed
@@ -420,7 +424,46 @@ impl<D: Dimension> Transform<D> {
         M: AsArray<'a, T, D>,
         T: 'a + Clone + AsPrimitive<f64>,
     {
-        Registration::new_translation()?.register(fixed, moving)
+        let fixed = fixed.into();
+        let moving = moving.into();
+        let initial_guess = Self::register_affine(&fixed, &moving)?;
+        let initial_guess =
+            initial_guess.parameters[initial_guess.ndim() * initial_guess.ndim()..].to_vec();
+        Registration::new_translation()?
+            .with_initial_guess(initial_guess)?
+            .register(fixed, moving)
+    }
+
+    /// find the translation which transforms moving into fixed
+    /// using the cross-correlation theorem and fft
+    pub fn register_translation_fft<'a, F, M, T>(fixed: F, moving: M) -> Result<Self, Error>
+    where
+        F: AsArray<'a, T, D>,
+        M: AsArray<'a, T, D>,
+        T: 'a + Clone + AsPrimitive<f64>,
+    {
+        let f = fft(fixed.into().mapv(|i| i.as_()).view())?.mapv(|i| i.conj());
+        let m = fft(moving.into().mapv(|i| i.as_()).view())?;
+        let c = ifft(&(f * m))?.mapv(|i| i.re);
+        let (i, _) = c
+            .indexed_iter()
+            .max_by(|(_, a), (_, b)| a.total_cmp(b))
+            .expect("no max because array is empty");
+        let d = i.into_dimension();
+        let j = d
+            .as_array_view()
+            .iter()
+            .zip(c.shape())
+            .map(|(&i, &s)| {
+                let mut fi = i as f64;
+                let fs = s as f64;
+                if 2.0 * fi > fs {
+                    fi -= fs;
+                }
+                fi
+            })
+            .collect::<Vec<_>>();
+        Ok(Self::from_translation(&j))
     }
 
     /// find the transform which transforms moving into fixed, using fixed_mu to specify which