From 0f315e66cabf996f33a43f9058838cb6c48a2241 Mon Sep 17 00:00:00 2001 From: "w.pomp" Date: Tue, 28 Jul 2026 11:51:04 +0200 Subject: [PATCH] - all tests passing --- src/filter.rs | 8 +--- src/metric.rs | 89 +++++------------------------------------ src/par_indexed_iter.rs | 2 +- src/register.rs | 7 ++-- src/transform.rs | 6 +-- 5 files changed, 20 insertions(+), 92 deletions(-) diff --git a/src/filter.rs b/src/filter.rs index 8f4a5d8..292301c 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -112,7 +112,7 @@ where mod tests { use super::*; use crate::julia_image; - use ndarray::{Ix1, Ix2, array}; + use ndarray::{Ix1, array}; use tiffwrite::IJTiffFile; #[test] @@ -142,12 +142,6 @@ mod tests { Ok(()) } - #[test] - fn gaussian_kernel_test2() -> Result<(), Box> { - let k = gaussian_kernel::(&[60, 80], &[16.0, 16.0])?; - ndarray_npy::write_npy("/home/wim/tmp/kernel.npy", &k)?; - Ok(()) - } #[test] fn fft_test() -> Result<(), Box> { diff --git a/src/metric.rs b/src/metric.rs index 1e98e91..e1d3f5f 100644 --- a/src/metric.rs +++ b/src/metric.rs @@ -167,15 +167,13 @@ impl Sampling { fn index(&self) -> Vec> { match self { Self::Fixed(index) => index.clone(), - Self::Random((n_samples, shape, _cached)) => { - rand::rng() - .random_iter::() - .take(n_samples * shape.len()) - .chunks(shape.len()) - .into_iter() - .map(|c| c.zip_eq(shape.iter()).map(|(i, s)| i * s - 0.5).collect()) - .collect() - } + Self::Random((n_samples, shape, _cached)) => rand::rng() + .random_iter::() + .take(n_samples * shape.len()) + .chunks(shape.len()) + .into_iter() + .map(|c| c.zip_eq(shape.iter()).map(|(i, s)| i * s - 0.5).collect()) + .collect(), } } } @@ -208,7 +206,7 @@ where shape: Vec, center: Vec, fixed: BSpline<0, D>, - moving: BSpline<3, D>, + moving: BSpline<1, D>, fixed_mu: FixedMu, minmax: [f64; 2], sampling: Sampling, @@ -226,7 +224,7 @@ where { pub fn new( fixed: BSpline<0, D>, - moving: BSpline<3, D>, + moving: BSpline<1, D>, sampling: SamplingArg, n_bins: usize, edge: f64, @@ -523,8 +521,6 @@ mod tests { use algos::optimization::{bfgs_minimize, gradient_descent_minimize}; use ndarray::{MeshIndex, array, meshgrid, stack}; use ndarray_npy::NpzWriter; - use num::integer::Roots; - use num::traits::FloatConst; use std::fs::File; #[test] @@ -546,69 +542,6 @@ mod tests { Ok(()) } - #[test] - fn grad() -> Result<(), Box> { - let mut rng = rand::rng(); - let u1 = (&mut rng) - .random_iter::() - .take(1000) - .collect::>(); - let u2 = (&mut rng) - .random_iter::() - .take(1000) - .collect::>(); - let f = u1 - .iter() - .zip_eq(u2.iter()) - .map(|(i, j)| (-2.0 * i.ln()).sqrt() * (2.0 * f64::PI() * j).cos() - 1.0) - .collect::>(); - let m = u1 - .iter() - .zip_eq(u2.iter()) - .map(|(i, j)| (-2.0 * i.ln()).sqrt() * (2.0 * f64::PI() * j).sin() + 1.0) - .collect::>(); - - let s = f.len().sqrt(); - let n_bins = s.max(11); - let e = 2.0 / (n_bins - 3) as f64; - let first_bin = -5.0 - e; - let last_bin = 5.0 + e; - // println!("f: {:?}, m: {:?}", f, m); - // println!("first_bin: {}, last_bin: {}, n_bins: {}, e: {}", first_bin, last_bin, n_bins, e); - let w = vec![1.0; f.len()]; - let dw = vec![vec![0.0]; f.len()]; - let (jpdf, d_jpdf_m) = parzen( - 1.0, - &[1.0], - &f, - &m, - &[], - &w, - &dw, - first_bin, - last_bin, - n_bins, - ); - let pdf_f = jpdf.rows().into_iter().map(|i| i.sum()).collect::>(); - let pdf_m = jpdf - .columns() - .into_iter() - .map(|i| i.sum()) - .collect::>(); - let d_pdf_m = d_jpdf_m - .columns() - .into_iter() - .map(|i| i.sum()) - .collect::>(); - println!( - "pdf_f = {:?}\npdf_m = {:?}\nd_pdf_m = {:?}", - pdf_f.to_vec(), - pdf_m.to_vec(), - d_pdf_m.to_vec() - ); - Ok(()) - } - #[test] fn metric() -> Result<(), Box> { let im_a = array![0.0, 0.0, 1.0, 0.0, 1.0, 2.0, 1.0, 0.0, 1.0, 0.0, 0.0]; @@ -837,7 +770,7 @@ mod tests { let a = array![0.0, 0.0, 1.0, 0.0, 1.0, 2.0, 1.0, 0.0, 1.0, 0.0, 0.0]; let b = array![0.0, 0.0, 1.0, 0.0, 1.0, 2.0, 1.0, 0.0, 1.0, 0.0, 0.0]; let fixed = BSpline::<0, _>::new(a.view()); - let moving = BSpline::<3, _>::new(b.view()); + let moving = BSpline::<1, _>::new(b.view()); let mus = vec![0.0, 0.001]; let mut npz = NpzWriter::new(File::create( @@ -918,7 +851,7 @@ mod tests { let a = array![0.0, 0.0, 1.0, 0.0, 1.0, 2.0, 1.0, 0.0, 1.0, 0.0, 0.0]; let b = array![0.0, 0.0, 1.0, 0.0, 1.0, 2.0, 1.0, 0.0, 1.0, 0.0, 0.0]; let fixed = BSpline::<0, _>::new(a.view()); - let moving = BSpline::<3, _>::new(b.view()); + let moving = BSpline::<1, _>::new(b.view()); let mus = Array1::linspace(-10.0, 10.0, 500); let s = a.len() as f64; diff --git a/src/par_indexed_iter.rs b/src/par_indexed_iter.rs index cfacf28..24a4b24 100644 --- a/src/par_indexed_iter.rs +++ b/src/par_indexed_iter.rs @@ -208,7 +208,7 @@ impl<'a, T> BaseIndexedIterMut<'a, T> { } } -/// An parallel iterator over array indices and mutable values. +/// A parallel iterator over array indices and mutable values. pub struct IndexedIterMut<'a, T>(BaseIndexedIterMut<'a, T>); /// An iterator over array indices and mutable values. diff --git a/src/register.rs b/src/register.rs index 3d7b7d4..7c7e0c1 100644 --- a/src/register.rs +++ b/src/register.rs @@ -64,6 +64,7 @@ impl RegistrationStep { /// **NOT downsampled**. See: `itkMultiResolutionGaussianSmoothingPyramidImageFilter.hxx` /// /// Schedule [8, 4, 2, 1] means σ = [4.0, 2.0, 1.0, 0.5] (spacing=1). + /// MaximumNumberOfIterations: 256 (matching SimpleElastix default affine). pub fn default_steps(ndim: usize, n: usize) -> Vec { let sigma_schedule: Vec = vec![4.0, 2.0, 1.0, 0.5]; let nlevels = sigma_schedule.len(); @@ -81,7 +82,7 @@ impl RegistrationStep { 32, 1e-6, 0.05, - 512, + 256, 1.0, )); } @@ -309,7 +310,7 @@ impl Registration { continue; } let bf = BSpline::<0, _>::new(f.view()); - let bm = BSpline::<3, _>::new(m.view()); + let bm = BSpline::<1, _>::new(m.view()); let metric = MattesMetric::new(bf, bm, samples, n_bins, edge)? .with_fixed_mu(self.fixed_mu.clone()); @@ -409,7 +410,7 @@ impl Registration { } let bf = BSpline::<0, _>::new(f.view()); - let bm = BSpline::<3, _>::new(m.view()); + let bm = BSpline::<1, _>::new(m.view()); let n_samples = match &samples { SamplingArg::Fixed(n) => *n, SamplingArg::Random(n) => *n, diff --git a/src/transform.rs b/src/transform.rs index 6d7f76f..ed0864d 100644 --- a/src/transform.rs +++ b/src/transform.rs @@ -879,7 +879,7 @@ mod tests { let f = gaussian_smooth(im_a.view(), &[*sigma_val; 2])?; let m = gaussian_smooth(im_b.view(), &[*sigma_val; 2])?; let bf = BSpline::<0, _>::new(f.view()); - let bm = BSpline::<3, _>::new(m.view()); + let bm = BSpline::<1, _>::new(m.view()); let metric = MattesMetric::new(bf, bm, SamplingArg::Random(3000), 128, edge)?; let mi_id = metric.evaluate(&identity); @@ -907,7 +907,7 @@ mod tests { let f4 = gaussian_smooth(im_a.view(), &[4.0, 4.0])?; let m4 = gaussian_smooth(im_b.view(), &[4.0, 4.0])?; let bf4 = BSpline::<0, _>::new(f4.view()); - let bm4 = BSpline::<3, _>::new(m4.view()); + let bm4 = BSpline::<1, _>::new(m4.view()); let metric4 = MattesMetric::new(bf4, bm4, SamplingArg::Random(5000), 128, edge)?; let mi_p_coarse = metric4.evaluate(&p); let mi_qinv_coarse = metric4.evaluate(&q_inv); @@ -951,7 +951,7 @@ mod tests { let m = gaussian_smooth(im_b.view(), &[4.0, 4.0])?; let metric = MattesMetric::new( BSpline::<0, _>::new(f.view()), - BSpline::<3, _>::new(m.view()), + BSpline::<1, _>::new(m.view()), SamplingArg::Random(3000), 128, edge,