wim/image-registration
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

- cleanup

Browse Source
This commit is contained in:
w.pomp
2026-08-03 14:43:04 +02:00
parent 687c37360d
commit 1563352d42
4 changed files with 7 additions and 614 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Cargo.toml
+2 -3
View File
@@ -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
src/bspline.rs
+4 -243
View File
@@ -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<Self, Error> {
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<f64, D::Larger>) -> Result<Array<f64, D>, 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<Array<Complex<f64>, D>, Error> {
pub(crate) fn reduction_filter(shape: &[usize]) -> Result<Array<Complex<f64>, D>, Error> {
let mut filter: Array<Complex<f64>, 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<f64>, Array1<f64>, Array1<f64>) {
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<f64>, Array2<f64>, Array2<f64>, Array2<f64>) {
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<f64>, Array2<f64>, Array2<f64>) {
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<dyn std::error::Error>> {
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<dyn std::error::Error>> {
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::<Vec<_>>();
let center = [center[0], center[1]];
let k = julia_image(&shape, &transform, &center, &[-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::<Vec<_>>();
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<dyn std::error::Error>> {
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(())
}
}
src/error.rs
-2
View File
@@ -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,
}
src/metric.rs
+1 -366
View File
@@ -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<dyn std::error::Error>> {
@@ -526,39 +522,6 @@ mod tests {
Ok(())
}
#[test]
fn metric() -> Result<(), Box<dyn std::error::Error>> {
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<dyn std::error::Error>> {
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<dyn std::error::Error>> {
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<dyn std::error::Error>> {
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<dyn std::error::Error>> {
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<dyn std::error::Error>> {
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<dyn std::error::Error>> {
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::<f64>()
.take(25)
.map(|i| i * s - 0.5)
.collect::<Vec<_>>();
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::<f64>();
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::<Vec<_>>();
let dw = dw.iter().map(|i| i[0]).collect::<Vec<_>>();
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<dyn std::error::Error>> {
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::<f64>()
.take(25)
.map(|i| i * s - 0.5)
.collect::<Vec<_>>();
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::<f64>();
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::<Vec<_>>();
let pdf_m = jpdf.axis_iter(Axis(1)).map(|i| i.sum()).collect::<Vec<_>>();
let d_pdf_m = d_jpdf_m
.axis_iter(Axis(1))
.map(|i| i.sum())
.collect::<Vec<_>>();
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<dyn std::error::Error>> {
let x = Array1::linspace(-2.0, 2.0, 500);