diff --git a/src/register.rs b/src/register.rs index babd86c..d6f5cb9 100644 --- a/src/register.rs +++ b/src/register.rs @@ -175,7 +175,7 @@ impl Registration { } 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; diff --git a/src/transform.rs b/src/transform.rs index d3aa10c..aa6add8 100644 --- a/src/transform.rs +++ b/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 Mul<&Transform> for &Transform { .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 Transform { } pub fn new_around_center(parameters: Vec, center: &[f64]) -> Self { - (Self::new(vec![1.0, 0.0, 0.0, 1.0, -center[0], -center[1]]) * Self::new(parameters)) - * Self::new(vec![1.0, 0.0, 0.0, 1.0, center[0], center[1]]) + Self::new(parameters).with_translation(¢er.iter().map(|i| -i).collect::>()) + * Self::from_translation(center) } /// find the affine transform which transforms moving into fixed @@ -420,7 +424,46 @@ impl Transform { M: AsArray<'a, T, D>, T: 'a + Clone + AsPrimitive, { - 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 + where + F: AsArray<'a, T, D>, + M: AsArray<'a, T, D>, + T: 'a + Clone + AsPrimitive, + { + 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::>(); + Ok(Self::from_translation(&j)) } /// find the transform which transforms moving into fixed, using fixed_mu to specify which