From 1563352d42c2ef5e5720841d9a894d8fa6ee6daa Mon Sep 17 00:00:00 2001 From: "w.pomp" Date: Mon, 3 Aug 2026 14:43:04 +0200 Subject: [PATCH] - cleanup --- Cargo.toml | 5 +- src/bspline.rs | 247 +-------------------------------- src/error.rs | 2 - src/metric.rs | 367 +------------------------------------------------ 4 files changed, 7 insertions(+), 614 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 43fb81d..ac27cf1 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -12,13 +12,12 @@ documentation = "https://docs.rs/image-registration" readme = "README.md" keywords = ["affine", "transformation", "ndarray"] categories = ["multimedia::images", "science"] -exclude = ["/tests"] +exclude = ["/test_files"] [dependencies] algos = "0.6" itertools = "0.15" ndarray = { version = "0.17", features = ["rayon"] } -ndarray-npy = { version = "0.10.0", features = ["npz"] } ndrustfft = "0.6" num = "0.4" rayon = "1" @@ -26,11 +25,11 @@ serde = { version = "1", features = ["derive"] } serde_yaml = "0.9" rand = "0.10" thiserror = "2" -tiffwrite = "2026.6.0" [dev-dependencies] tempfile = "3" tiff = "0.11" +tiffwrite = "2026.6.0" [profile.release] debug = true diff --git a/src/bspline.rs b/src/bspline.rs index 1a38996..ef926f9 100644 --- a/src/bspline.rs +++ b/src/bspline.rs @@ -1,17 +1,14 @@ use crate::Error; -use crate::filter::{fft, ifft}; use crate::par_indexed_iter::ParallelIndexedIterMut; use crate::transform::{Transform, transform_point}; use itertools::Itertools; use ndarray::{ Array, Array1, Array2, ArrayD, ArrayViewMut1, AsArray, Axis, Dimension, IntoDimension, IxDyn, - SliceInfoElem, }; use num::traits::FloatConst; use num::{Complex, cast::AsPrimitive, complex::ComplexFloat}; use rayon::iter::ParallelIterator; use std::ops::{Deref, MulAssign}; -use tiffwrite::IJTiffFile; fn in_shape(shape_f: &[f64], v: &[f64]) -> bool { shape_f @@ -786,38 +783,6 @@ where } } - pub fn subsample2(&self) -> Result { - let shape = self.coefficients.shape(); - let v = Self::reduction_filter2(shape)?; - let c = fft(self.coefficients.view())?; - let mut t = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/fft.tif")).unwrap(); - let v0 = v.mapv(|i| i.re).into_dimensionality()?; - let c0 = c.mapv(|i| i.re).into_dimensionality()?; - let filtered = ifft((v * c).view())?.into_dyn().mapv(|i| i.re); - t.save(v0.view(), 0, 0, 0).unwrap(); - t.save(c0.view(), 1, 0, 0).unwrap(); - t.save(filtered.view().into_dimensionality()?.view(), 2, 0, 0) - .unwrap(); - - let slice = vec![ - SliceInfoElem::Slice { - start: 0, - end: None, - step: 2 - }; - filtered.ndim() - ]; - let coefficients = filtered - .slice(slice.as_slice()) - .to_owned() - .into_dimensionality()?; - Ok(BSpline { - coefficients, - max_number_interpolation_points: self.max_number_interpolation_points, - points_to_index: self.points_to_index.clone(), - }) - } - pub fn evaluate(&self, index: Array) -> Result, Error> where D: Dimension, @@ -851,7 +816,7 @@ where } /// table III The L2 Polynomial Spline Pyramid, Unser et al. 1993 - pub(crate) fn reduction_filter2(shape: &[usize]) -> Result, D>, Error> { + pub(crate) fn reduction_filter(shape: &[usize]) -> Result, D>, Error> { let mut filter: Array, D> = ArrayD::ones(shape).into_dimensionality()?; for (i, &s) in shape.iter().enumerate() { let si = s as isize; @@ -880,226 +845,22 @@ where } #[inline] -pub fn cubic_bspline(x: f64) -> f64 { +pub(crate) fn cubic_bspline(x: f64) -> f64 { (x + 2.0).max(0.0).powi(3) / 6.0 - (x + 1.0).max(0.0).powi(3) / 1.5 + x.max(0.0).powi(3) - (x - 1.0).max(0.0).powi(3) / 1.5 + (x - 2.0).max(0.0).powi(3) / 6.0 } #[inline] -pub fn square_bspline(x: f64) -> f64 { +pub(crate) fn square_bspline(x: f64) -> f64 { (x + 1.5).max(0.0).powi(2) / 2.0 - (x + 0.5).max(0.0).powi(2) * 1.5 + (x - 0.5).max(0.0).powi(2) * 1.5 - (x - 1.5).max(0.0).powi(2) / 2.0 } #[inline] -pub fn square_bspline_integral(x: f64) -> f64 { +pub(crate) fn square_bspline_integral(x: f64) -> f64 { (x + 1.5).max(0.0).powi(3) / 6.0 - (x + 0.5).max(0.0).powi(3) / 2.0 + (x - 0.5).max(0.0).powi(3) / 2.0 - (x - 1.5).max(0.0).powi(3) / 6.0 } - -pub fn cubic_bspline_parzen_1d( - data: &[f64], - first_bin: f64, - last_bin: f64, - n_bins: usize, - e: f64, -) -> (Array1, Array1, Array1) { - let n_bins_i = n_bins as isize; - let mut hist = Array1::zeros(n_bins); - let mut grad = Array1::zeros(n_bins); - let d = ((n_bins - 1) as f64) / (last_bin - first_bin); - let m = (2.0 * d / e) as isize; - let bins = Array1::linspace(first_bin, last_bin, n_bins); - for x in data { - let n = (d * (x - first_bin)).round() as isize; - for i in (n - m)..=(n + m) { - if (0 <= i) && (i < n_bins_i) { - let j = i as usize; - let k = (x - bins[j]) / e; - hist[j] += cubic_bspline(k); - grad[j] += square_bspline(k - 0.5) - square_bspline(k + 0.5); - } - } - } - let alpha = (data.len() as f64) * e * d; - (bins, hist / alpha, grad / alpha) -} - -pub fn cubic_bspline_joint_parzen( - data_a: &[f64], - data_b: &[f64], - first_bin: f64, - last_bin: f64, - n_bins: usize, - e: f64, -) -> (Array1, Array2, Array2, Array2) { - let n_bins_i = n_bins as isize; - let mut hist = Array2::zeros([n_bins, n_bins]); - let mut grad_a = Array2::zeros([n_bins, n_bins]); - let mut grad_b = Array2::zeros([n_bins, n_bins]); - let d = ((n_bins - 1) as f64) / (last_bin - first_bin); - let m = (2.0 * d / e) as isize; - let bins = Array1::linspace(first_bin, last_bin, n_bins); - for xa in data_a { - for xb in data_b { - let na = (d * (xa - first_bin)).round() as isize; - let nb = (d * (xb - first_bin)).round() as isize; - for ia in (na - m)..=(na + m) { - for ib in (nb - m)..=(nb + m) { - if (0 <= ia) && (ia < n_bins_i) && (0 <= ib) && (ib < n_bins_i) { - let ja = ia as usize; - let jb = ib as usize; - let ka = (xa - bins[ja]) / e; - let kb = (xb - bins[jb]) / e; - let csa = cubic_bspline(ka); - let csb = cubic_bspline(kb); - hist[[ja, jb]] += csa * csb; - grad_a[[ja, jb]] += - (square_bspline(ka - 0.5) - square_bspline(ka + 0.5)) * csb; - grad_b[[ja, jb]] += - (square_bspline(kb - 0.5) - square_bspline(kb + 0.5)) * csa; - } - } - } - } - } - let alpha = (data_a.len() as f64) * e * d * (data_b.len() as f64) * e * d; - (bins, hist / alpha, grad_a / alpha, grad_b / alpha) -} - -pub fn cubic_bspline_joint_parzen_grad_b( - data_a: &[f64], - data_b: &[f64], - first_bin: f64, - last_bin: f64, - n_bins: usize, - e: f64, -) -> (Array1, Array2, Array2) { - let n_bins_i = n_bins as isize; - let mut hist = Array2::zeros([n_bins, n_bins]); - let mut grad_b = Array2::zeros([n_bins, n_bins]); - let d = ((n_bins - 1) as f64) / (last_bin - first_bin); - let m = (2.0 * d / e) as isize; - let bins = Array1::linspace(first_bin, last_bin, n_bins); - for (xa, xb) in data_a.iter().zip_eq(data_b.iter()) { - let na = (d * (xa - first_bin)).round() as isize; - let nb = (d * (xb - first_bin)).round() as isize; - for ia in (na - m)..=(na + m) { - if (0 <= ia) && (ia < n_bins_i) { - let ja = ia as usize; - let ka = (xa - bins[ja]) / e; - let csa = cubic_bspline(ka); - for ib in (nb - m)..=(nb + m) { - if (0 <= ib) && (ib < n_bins_i) { - let jb = ib as usize; - let kb = (xb - bins[jb]) / e; - hist[[ja, jb]] += csa * cubic_bspline(kb); - grad_b[[ja, jb]] += - csa * (square_bspline(kb - 0.5) - square_bspline(kb + 0.5)); - } - } - } - } - } - let alpha = (data_a.len() as f64) * e * d; - // println!("alpha: {}", alpha); - (bins, hist / alpha, grad_b / alpha) -} - -#[cfg(test)] -mod test { - use super::{BSpline, BSplineTrait, cubic_bspline_joint_parzen_grad_b}; - use crate::error::Error; - use crate::julia_image; - use ndarray::{Array1, Array3, Ix1, array, s}; - use ndarray_npy::NpzWriter; - use std::fs::File; - - #[test] - fn bspline_grad() -> Result<(), Box> { - let a = vec![0.0, 0.0, 1.0, 0.0, 1.0, 5.0, 1.0, 0.0, 1.0, 0.0, 0.0]; - let first_bin = 0.0; - let last_bin = 5.0; - let n_bins = 3; - let e = (last_bin - first_bin) / (n_bins - 1) as f64; - let (bins, jpdf, _d_jpdf_m) = - cubic_bspline_joint_parzen_grad_b(&a, &a, first_bin, last_bin, n_bins, e); - println!("bins: {:?}, jpdf: {:?}", bins, jpdf); - println!("sum: {}", jpdf.sum()); - Ok(()) - } - - #[test] - fn reduction_filter() -> Result<(), Error> { - let _v = BSpline::<3, Ix1>::reduction_filter2(&[10]); - Ok(()) - } - - #[test] - fn subsample2() -> Result<(), Box> { - let shape = [600, 800]; - let transform = [1.0, 0.0, 0.0, 1.0, 0.0, 0.0]; - let center = shape - .iter() - .map(|&s| (s - 1) as f64 / 2.0) - .collect::>(); - let center = [center[0], center[1]]; - let k = julia_image(&shape, &transform, ¢er, &[-0.4, 0.6]); - let b = BSpline::<3, _>::new(k.view()); - let b2 = b.subsample2()?; - let shape = k.shape().iter().map(|i| i / 2).collect::>(); - let mut xy = Array3::zeros((2, shape[0], shape[1])); - for (i, mut x) in xy - .slice_mut(s![1, .., ..]) - .columns_mut() - .into_iter() - .enumerate() - { - x.fill(i as f64); - } - for (i, mut y) in xy - .slice_mut(s![0, .., ..]) - .rows_mut() - .into_iter() - .enumerate() - { - y.fill(i as f64); - } - let k2 = b2.evaluate(xy)?.mapv(|i| i.clamp(0.0, 255.0) as u8); - let mut t1 = - tiffwrite::IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/subsample1.tif"))?; - // t1.save(b.coefficients.view(), 0, 0, 0)?; - t1.save(k.view(), 0, 0, 0)?; - - let mut t2 = - tiffwrite::IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/subsample2.tif"))?; - // t2.save(b2.coefficients.view(), 0, 0, 0)?; - t2.save(k2.view(), 0, 0, 0)?; - Ok(()) - } - - #[test] - fn derivative() -> Result<(), Box> { - 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 = BSpline::<3, _>::new(a.view()); - let mut mem = b.get_mem(); - let x = Array1::linspace(-0.5, 10.5, 500); - let mut v = Vec::new(); - let mut d = Vec::new(); - for i in &x { - let (vi, di) = b.evaluate_value_and_derivative_at_continuous_index(&[*i], &mut mem)?; - v.push(vi); - d.push(di[0]); - } - let mut npz = NpzWriter::new(File::create( - std::env::home_dir().unwrap().join("tmp/metric.npz"), - )?); - npz.add_array("x", &x)?; - npz.add_array("v", &Array1::from(v))?; - npz.add_array("d", &Array1::from(d))?; - Ok(()) - } -} diff --git a/src/error.rs b/src/error.rs index bac24cd..189974d 100644 --- a/src/error.rs +++ b/src/error.rs @@ -8,8 +8,6 @@ pub enum Error { SerdeYAML(#[from] serde_yaml::Error), #[error(transparent)] ShapeError(#[from] ndarray::ShapeError), - #[error(transparent)] - NpyError(#[from] ndarray_npy::WriteNpzError), #[error("number of dimensions is not defined")] NumberOfDimensionsNotDefined, } diff --git a/src/metric.rs b/src/metric.rs index 45ec974..813b784 100644 --- a/src/metric.rs +++ b/src/metric.rs @@ -499,13 +499,9 @@ impl Sigma { mod tests { use super::*; use crate::filter::gaussian_smooth; - use crate::julia_image; - use crate::transform::Transform; use algos::OptimizationConfig; use algos::optimization::{bfgs_minimize, gradient_descent_minimize}; - use ndarray::{MeshIndex, array, meshgrid, stack}; - use ndarray_npy::NpzWriter; - use std::fs::File; + use ndarray::array; #[test] fn derivative() -> Result<(), Box> { @@ -526,39 +522,6 @@ mod tests { 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]; - let im_b = array![1.0, 0.0, 1.0, 2.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0]; - // let a = gaussian_smooth(im_a.view(), &[4.0])?; - // let b = gaussian_smooth(im_b.view(), &[4.0])?; - let a = im_a; - let b = im_b; - - // let points = vec![vec![3.0], vec![4.0], vec![5.0], vec![6.0], vec![7.0], vec![8.0]]; - // let m = MattesMetric::new_from_arrays(a.view(), b.view(), SamplingArg::FixedAt(points))? - // .with_fixed_mu(vec![Some(1.0), None]); - let m = MattesMetric::new_from_arrays(a.view(), b.view(), SamplingArg::Fixed(20), 5, 0.05)? - .with_fixed_mu(vec![Some(1.0), None]); - // let m = MattesMetric::new_all_from_arrays(im_a.view(), im_b.view())?.with_fixed_mu(vec![Some(1.0), None]); - let mu = Array1::linspace(-4.0, 2.0, 601); - let mut v = Vec::new(); - let mut d = Vec::new(); - for x in &mu { - v.push(m.evaluate(&[*x])); - d.push(m.gradient(&[*x]).unwrap()[0]); - } - let mut npz = NpzWriter::new(File::create( - std::env::home_dir().unwrap().join("tmp/metric.npz"), - )?); - npz.add_array("mu", &mu)?; - npz.add_array("v", &Array1::from_vec(v))?; - npz.add_array("d", &Array1::from_vec(d))?; - npz.add_array("a", &a)?; - npz.add_array("b", &b)?; - Ok(()) - } - #[test] fn metric2() -> 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]; @@ -597,135 +560,6 @@ mod tests { Ok(()) } - #[test] - fn metric1() -> Result<(), Box> { - let im_a = julia_image( - &[60, 80], - &[1.0, 0.0, 0.0, 1.0, 0.0, 0.0], - &[29.5, 39.5], - &[-0.8, 0.156], - ) - .mapv(|i| i as f64); - let im_b = julia_image( - &[60, 80], - &[1.0, 0.0, 0.0, 1.0, 10.0, -20.0], - &[29.5, 39.5], - &[-0.8, 0.156], - ) - .mapv(|i| i as f64); - let a = gaussian_smooth(im_a.view(), &[30.0, 40.0])?; - let b = gaussian_smooth(im_b.view(), &[30.0, 40.0])?; - let m = MattesMetric::new_from_arrays(a.view(), b.view(), SamplingArg::Fixed(250), 3, 0.2)? - .with_fixed_mu(vec![Some(1.0), Some(0.0), Some(0.0), Some(1.0), None, None]); - let mut npz = NpzWriter::new(File::create( - std::env::home_dir().unwrap().join("tmp/metric.npz"), - )?); - - let s = 200; - let mu = Array1::linspace(-50.0, 50.0, s); - let (mux, muy) = meshgrid((&mu, &mu), MeshIndex::XY); - let mu = stack(Axis(0), &[mux, muy])?; - let mut v = Array2::zeros([s, s]); - let mut d = Array3::zeros([2, s, s]); - for (mui, (vi, mut di)) in mu - .lanes(Axis(0)) - .into_iter() - .zip_eq(v.iter_mut().zip_eq(d.lanes_mut(Axis(0)))) - { - *vi = m.evaluate(mui.to_vec().as_slice()); - di.assign(&Array1::from_vec( - m.gradient(mui.to_vec().as_slice()).unwrap(), - )); - } - npz.add_array("a", &a)?; - npz.add_array("b", &b)?; - npz.add_array("mu", &mu)?; - npz.add_array("v", &v)?; - npz.add_array("d", &d)?; - Ok(()) - } - - #[test] - fn metric3() -> Result<(), Box> { - let im_a = julia_image( - &[100, 1], - &[1.0, 0.0, 0.0, 0.01, 0.0, 0.0], - &[99.5, 0.5], - &[-0.8, 0.156], - ) - .slice(s![.., 0]) - .mapv(|i| i as f64); - let im_b = - Transform::new(vec![1.0, 0.0], vec![im_a.shape()[0]]).interpolate::<1, _, _>(&im_a)?; - let a = gaussian_smooth(im_a.view(), &[25.0])?; - let b = gaussian_smooth(im_b.view(), &[25.0])?; - let m = MattesMetric::new_from_arrays(a.view(), b.view(), SamplingArg::Fixed(100), 6, 0.1)? - .with_fixed_mu(vec![None, None]); - let mut npz = NpzWriter::new(File::create( - std::env::home_dir().unwrap().join("tmp/metric.npz"), - )?); - - let s = 200; - let translate = Array1::linspace(-50.0, 50.0, s); - let scale = Array1::logspace(10.0, -1.0, 1.0, s); - let (mux, muy) = meshgrid((&scale, &translate), MeshIndex::XY); - let mu = stack(Axis(0), &[mux, muy])?; - let mut v = Array2::zeros([s, s]); - let mut d = Array3::zeros([2, s, s]); - for (mui, (vi, mut di)) in mu - .lanes(Axis(0)) - .into_iter() - .zip_eq(v.iter_mut().zip_eq(d.lanes_mut(Axis(0)))) - { - *vi = m.evaluate(mui.to_vec().as_slice()); - di.assign(&Array1::from_vec( - m.gradient(mui.to_vec().as_slice()).unwrap(), - )); - } - npz.add_array("a", &a)?; - npz.add_array("b", &b)?; - npz.add_array("mu", &mu)?; - npz.add_array("v", &v)?; - npz.add_array("d", &d)?; - Ok(()) - } - - #[test] - fn metric4() -> Result<(), Box> { - let im_a = julia_image( - &[100, 1], - &[1.0, 0.0, 0.0, 0.01, 0.0, 0.0], - &[99.5, 0.5], - &[-0.8, 0.156], - ) - .slice(s![.., 0]) - .mapv(|i| i as f64); - let im_b = - Transform::new(vec![1.0, 0.0], vec![im_a.shape()[0]]).interpolate::<1, _, _>(&im_a)?; - let a = gaussian_smooth(im_a.view(), &[50.0])?; - let b = gaussian_smooth(im_b.view(), &[50.0])?; - let m = - MattesMetric::new_from_arrays(a.view(), b.view(), SamplingArg::Fixed(1000), 3, 0.05)? - .with_fixed_mu(vec![None, Some(0.0)]); - let mut npz = NpzWriter::new(File::create( - std::env::home_dir().unwrap().join("tmp/metric.npz"), - )?); - - let mu = Array1::logspace(2.0, -1.0, 1.0, 1000); - let mut v = Vec::new(); - let mut d = Vec::new(); - for mui in &mu { - v.push(m.evaluate(&[*mui])); - d.push(m.gradient(&[*mui]).unwrap()[0]); - } - npz.add_array("a", &a)?; - npz.add_array("b", &b)?; - npz.add_array("mu", &mu)?; - npz.add_array("v", &Array1::from(v))?; - npz.add_array("d", &Array1::from(d))?; - Ok(()) - } - #[test] fn optimize() -> 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]; @@ -749,205 +583,6 @@ mod tests { Ok(()) } - #[test] - fn parzen_test() -> Result<(), Box> { - 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 mus = [0.0, 0.001]; - let mut npz = NpzWriter::new(File::create( - std::env::home_dir().unwrap().join("tmp/metric.npz"), - )?); - let s = a.len() as f64; - // let center = (s - 1.0) / 2.0; - let mut rng = rand::rng(); - let idx = (&mut rng) - .random_iter::() - .take(25) - .map(|i| i * s - 0.5) - .collect::>(); - - for (index, mu) in mus.iter().enumerate() { - let mut mem_f = fixed.get_mem(); - let mut mem_m = moving.get_mem(); - let mut w = Vec::new(); - let mut dw = Vec::new(); - let mut a = Vec::new(); - let mut b = Vec::new(); - let mut db = Vec::new(); - - for u in idx.clone() { - let v = u + mu; - if (-0.5 < v) && (v <= s - 0.5) { - let f = fixed.evaluate_at_continuous_index(&[u], &mut mem_f)?; - let (m, d) = moving - .evaluate_value_and_derivative_at_continuous_index(&[v], &mut mem_m)?; - // let j = image_jacobian(&[v], &[center], &d); - a.push(f); - b.push(m); - db.push(d); - // the metric wouldn't be smooth if we'd drop a point abruptly, - // so apply a 0 <= weight <= 1 to points within 5% of the edge - if 20.0 * (v + 0.5) < s { - let k = 60.0 * (v + 0.5) / s - 1.5; - w.push(square_bspline_integral(k)); - dw.push(vec![square_bspline(k) * 60.0 / s]); - } else if 20.0 * (v + 0.5) > 19.0 * s { - let k = 60.0 * (s - v - 0.5) / s - 1.5; - w.push(square_bspline_integral(k)); - dw.push(vec![-square_bspline(k) * 60.0 / s]); - } else { - w.push(1.0); - dw.push(vec![0.0]); - } - // let dw = image_jacobian(&[v], &[center], &dw); - } - } - - let alpha = w.iter().sum::(); - let mut dalpha = [0.0]; - for i in dw.iter() { - for (j, a) in i.iter().zip_eq(dalpha.iter_mut()) { - *a += j; - } - } - - let (j, dj) = parzen(alpha, &dalpha, &a, &b, &db, &w, &dw, 0.0, 2.0, 5); - let db = db.iter().map(|i| i[0]).collect::>(); - let dw = dw.iter().map(|i| i[0]).collect::>(); - - npz.add_array(format!("idx{}", index), &Array1::from_vec(idx.clone()))?; - npz.add_array(format!("a{}", index), &Array1::from_vec(a))?; - npz.add_array(format!("b{}", index), &Array1::from_vec(b))?; - npz.add_array(format!("db{}", index), &Array1::from_vec(db))?; - npz.add_array(format!("w{}", index), &Array1::from_vec(w))?; - npz.add_array(format!("dw{}", index), &Array1::from_vec(dw))?; - npz.add_array(format!("j{}", index), &j)?; - npz.add_array(format!("dj{}", index), &dj)?; - } - Ok(()) - } - - #[test] - fn metric_test() -> Result<(), Box> { - 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 mus = Array1::linspace(-10.0, 10.0, 500); - let s = a.len() as f64; - let mut rng = rand::rng(); - let idx = (&mut rng) - .random_iter::() - .take(25) - .map(|i| i * s - 0.5) - .collect::>(); - - let mut metrics = Vec::new(); - let mut dmetrics = Vec::new(); - - for mu in &mus { - let mut mem_f = fixed.get_mem(); - let mut mem_m = moving.get_mem(); - let mut w = Vec::new(); - let mut dw = Vec::new(); - let mut a = Vec::new(); - let mut b = Vec::new(); - let mut db = Vec::new(); - - for u in idx.clone() { - let v = u + mu; - if (-0.5 < v) && (v <= s - 0.5) { - let f = fixed.evaluate_at_continuous_index(&[u], &mut mem_f)?; - let (m, d) = moving - .evaluate_value_and_derivative_at_continuous_index(&[v], &mut mem_m)?; - // let j = image_jacobian(&[v], &[center], &d); - a.push(f); - b.push(m); - db.push(d); - // the metric wouldn't be smooth if we'd drop a point abruptly, - // so apply a 0 <= weight <= 1 to points within 5% of the edge - if 20.0 * (v + 0.5) < s { - let k = 60.0 * (v + 0.5) / s - 1.5; - w.push(square_bspline_integral(k)); - dw.push(vec![square_bspline(k) * 60.0 / s]); - } else if 20.0 * (v + 0.5) > 19.0 * s { - let k = 60.0 * (s - v - 0.5) / s - 1.5; - w.push(square_bspline_integral(k)); - dw.push(vec![-square_bspline(k) * 60.0 / s]); - } else { - w.push(1.0); - dw.push(vec![0.0]); - } - // let dw = image_jacobian(&[v], &[center], &dw); - } - } - - let alpha = w.iter().sum::(); - let mut dalpha = [0.0]; - for i in dw.iter() { - for (j, a) in i.iter().zip_eq(dalpha.iter_mut()) { - *a += j; - } - } - - let (jpdf, d_jpdf_m) = parzen(alpha, &dalpha, &a, &b, &db, &w, &dw, 0.0, 2.0, 5); - - let pdf_f = jpdf.axis_iter(Axis(0)).map(|i| i.sum()).collect::>(); - let pdf_m = jpdf.axis_iter(Axis(1)).map(|i| i.sum()).collect::>(); - let d_pdf_m = d_jpdf_m - .axis_iter(Axis(1)) - .map(|i| i.sum()) - .collect::>(); - - let mut metric = 0.0; - let mut dmetric = vec![0.0]; - let mut n; - let mut m; - - for ((row, d_row), k) in jpdf - .axis_iter(Axis(0)) - .zip_eq(d_jpdf_m.axis_iter(Axis(0))) - .zip_eq(pdf_f) - { - debug_assert_eq!(row.len(), d_row.shape()[0]); - debug_assert_eq!(row.len(), pdf_m.len()); - - for ((&j, dj), (i, di)) in row - .into_iter() - .zip_eq(d_row.axis_iter(Axis(0))) - .zip_eq(pdf_m.iter().zip(d_pdf_m.iter())) - { - // check for non-zero bin contribution - if (j > 1e-16) && (i * k > 1e-16) { - metric -= j * (j / (i * k)).ln(); - n = (j / (i * k)).ln() + 1.0; - m = di * j / i; - for (d, g) in dmetric.iter_mut().zip_eq(dj) { - *d += m - n * g; // eq 23 of Thevenaz & Unser paper [3] - } - } - } - } - - metrics.push(metric); - dmetrics.push(dmetric); - } - let mut npz = NpzWriter::new(File::create( - std::env::home_dir().unwrap().join("tmp/metric.npz"), - )?); - npz.add_array("mu", &mus)?; - npz.add_array("metric", &Array1::from_vec(metrics))?; - npz.add_array( - "dmetric", - &Array1::from_vec(dmetrics.iter().map(|d| d[0]).collect()), - )?; - Ok(()) - } - #[test] fn test_square_bspline_integral() -> Result<(), Box> { let x = Array1::linspace(-2.0, 2.0, 500);