diff --git a/src/metric.rs b/src/metric.rs index 40ddb09..a4a9ad2 100644 --- a/src/metric.rs +++ b/src/metric.rs @@ -212,7 +212,7 @@ where shape: Vec, center: Vec, fixed: BSpline<0, D>, - moving: BSpline<1, D>, + moving: BSpline<3, D>, fixed_mu: FixedMu, minmax: [f64; 2], sampling: Sampling, @@ -230,7 +230,7 @@ where { pub fn new( fixed: BSpline<0, D>, - moving: BSpline<1, D>, + moving: BSpline<3, D>, sampling: SamplingArg, n_bins: usize, edge: f64, @@ -776,7 +776,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::<1, _>::new(b.view()); + let moving = BSpline::<3, _>::new(b.view()); let mus = vec![0.0, 0.001]; let mut npz = NpzWriter::new(File::create( @@ -857,7 +857,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::<1, _>::new(b.view()); + let moving = BSpline::<3, _>::new(b.view()); let mus = Array1::linspace(-10.0, 10.0, 500); let s = a.len() as f64; diff --git a/src/register.rs b/src/register.rs index 41b26af..d4b1415 100644 --- a/src/register.rs +++ b/src/register.rs @@ -19,7 +19,7 @@ pub enum Optimizer { impl Default for Optimizer { fn default() -> Self { - Self::ASGD + Self::LBFGS } } @@ -61,16 +61,16 @@ impl RegistrationStep { /// /// Uses a multi-resolution pyramid with downsampling matching SimpleElastix: /// schedule [8, 4, 2, 1] → σ = [4.0, 2.0, 1.0, 0.5] at spacing=1. - /// 512 iterations per level, 2048–8192 cached random samples. + /// Uses all pixels at all levels for deterministic, precise convergence. pub fn default_steps(ndim: usize, n: usize) -> Vec { if ndim == 1 { return vec![Self { sigma: Sigma::Absolute(vec![1.0; 1]), - samples: SamplingArg::Random(2048.min(n)), + samples: SamplingArg::Random(n), n_bins: 32, - tolerance: 1e-6, + tolerance: 1e-8, edge: 0.05, - max_iterations: 1500, + max_iterations: 2048, learning_rate: 1.0, downsample: 1, }]; @@ -78,19 +78,23 @@ impl RegistrationStep { // Pyramid with downsampling: schedule [8,4,2,1], sigma [4,2,1,0.5] let sigma_schedule: Vec = vec![4.0, 2.0, 1.0, 0.5]; let downsample_schedule: Vec = vec![8, 4, 2, 1]; + let n_levels = sigma_schedule.len(); sigma_schedule .iter() .zip(downsample_schedule.iter()) - .map(|(&s, &d)| { + .enumerate() + .map(|(level, (&s, &d))| { let n_pixels = (n / (d * d)).max(4); - let n_samples = n_pixels.min(8192).max(2048); + let is_finest = level == n_levels - 1; + // Use all pixels at all levels for deterministic results + let (max_iter, tol) = if is_finest { (2048, 1e-8) } else { (512, 1e-6) }; Self { sigma: Sigma::Absolute(vec![s; ndim]), - samples: SamplingArg::Random(n_samples), + samples: SamplingArg::Random(n_pixels), n_bins: 32, - tolerance: 1e-6, + tolerance: tol, edge: 0.05, - max_iterations: 512, + max_iterations: max_iter, learning_rate: 1.0, downsample: d, } @@ -318,7 +322,7 @@ impl Registration { continue; } let bf = BSpline::<0, _>::new(f.view()); - let bm = BSpline::<1, _>::new(m.view()); + let bm = BSpline::<3, _>::new(m.view()); let metric = MattesMetric::new(bf, bm, samples, n_bins, edge)? .with_fixed_mu(self.fixed_mu.clone()); @@ -344,7 +348,7 @@ impl Registration { tolerance, maximum_step_length: 1.0, sp_a: 20.0, - sp_alpha: 1.0, + sp_alpha: 0.602, scales: Some(scales), ..Default::default() }; @@ -418,7 +422,7 @@ impl Registration { } let bf = BSpline::<0, _>::new(f.view()); - let bm = BSpline::<1, _>::new(m.view()); + let bm = BSpline::<3, _>::new(m.view()); let n_samples = match &samples { SamplingArg::Fixed(n) => *n, SamplingArg::Random(n) => *n, @@ -449,7 +453,7 @@ impl Registration { tolerance, maximum_step_length: 1.0, sp_a: 20.0, - sp_alpha: 1.0, + sp_alpha: 0.602, scales: Some(scales), ..Default::default() }; diff --git a/src/transform.rs b/src/transform.rs index f8ab58c..14fa45d 100644 --- a/src/transform.rs +++ b/src/transform.rs @@ -749,7 +749,7 @@ mod tests { ) .slice(s![.., 0]) .mapv(|i| i as f64); - let q = vec![0.85, 4.0]; + let q = vec![0.85, 2.0]; let im_b = Transform::new(q.clone(), vec![im_a.shape()[0]]).interpolate::<1, _, _>(&im_a)?; @@ -871,7 +871,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::<1, _>::new(m.view()); + let bm = BSpline::<3, _>::new(m.view()); let metric = MattesMetric::new(bf, bm, SamplingArg::Random(3000), 128, edge)?; let mi_id = metric.evaluate(&identity); @@ -899,7 +899,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::<1, _>::new(m4.view()); + let bm4 = BSpline::<3, _>::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); @@ -943,7 +943,7 @@ mod tests { let m = gaussian_smooth(im_b.view(), &[4.0, 4.0])?; let metric = MattesMetric::new( BSpline::<0, _>::new(f.view()), - BSpline::<1, _>::new(m.view()), + BSpline::<3, _>::new(m.view()), SamplingArg::Random(3000), 128, edge, @@ -1022,13 +1022,24 @@ mod tests { .fold(0.0f64, f64::max); println!("Our: {:?} max_err: {:.4} sse: {:.4}", t, max_err, sse); - let mut tif = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/register_real_images.tif"))?; + + let mut tif = IJTiffFile::new( + std::env::home_dir() + .unwrap() + .join("tmp/register_real_images.tif"), + )?; tif.save(fixed.mapv(|i| i as u16), 0, 0, 0)?; - tif.save(t.interpolate_par::<1, _, _>(moving.view())?.mapv(|i| i as u16), 1, 0, 0)?; + tif.save( + t.interpolate_par::<1, _, _>(moving.view())? + .mapv(|i| i as u16), + 1, + 0, + 0, + )?; tif.save(moving.mapv(|i| i as u16), 2, 0, 0)?; - assert!(max_err < 0.1); - assert!(sse < 0.1); + assert!(max_err < 0.02); + assert!(sse < 0.02); Ok(()) } @@ -1036,7 +1047,8 @@ mod tests { #[test] fn register_real_images2() -> Result<(), Box> { let fixed = read_tiff("test_files/fixed.tif")?; - let e = Transform::::new(vec![0.8, 0.0, 0.0, 1.0, 0.0, 0.0], fixed.shape().to_vec()).inverse()?; + let e = Transform::::new(vec![0.8, 0.0, 0.0, 1.0, 0.0, 0.0], fixed.shape().to_vec()) + .inverse()?; let moving = e.interpolate::<3, _, _>(fixed.view())?; let t = Transform::::register( @@ -1047,27 +1059,28 @@ mod tests { None, )?; + let e_inv = e.inverse()?; let sse = t .parameters .iter() - .zip(e.parameters.iter()) + .zip(e_inv.parameters.iter()) .map(|(a, b)| (a - b).powi(2)) .sum::(); let max_err = t .parameters .iter() - .zip(e.parameters.iter()) + .zip(e_inv.parameters.iter()) .map(|(a, b)| (a - b).abs()) .fold(0.0f64, f64::max); println!("Our: {:?} max_err: {:.4} sse: {:.4}", t, max_err, sse); - // let mut tif = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/register_real_images2.tif"))?; - // tif.save(fixed.mapv(|i| i as u16), 0, 0, 0)?; - // tif.save(t.interpolate_par::<1, _, _>(moving.view())?.mapv(|i| i as u16), 1, 0, 0)?; - // tif.save(moving.mapv(|i| i as u16), 2, 0, 0)?; + let mut tif = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/register_real_images2.tif"))?; + tif.save(fixed.mapv(|i| i as u16), 0, 0, 0)?; + tif.save(t.interpolate_par::<1, _, _>(moving.view())?.mapv(|i| i as u16), 1, 0, 0)?; + tif.save(moving.mapv(|i| i as u16), 2, 0, 0)?; - assert!(max_err < 0.1); - assert!(sse < 0.1); + assert!(max_err < 0.02); + assert!(sse < 0.02); Ok(()) }