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

- default value argument for interpolation

Browse Source
This commit is contained in:
Wim Pomp committed 2026-10-25 10:17:31 +01:00
1 parent b12982da48
commit 9257dd1fcf
4 files changed
+157 -35
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

src/transform.rs
+31 -28
View File
@@ -1,4 +1,4 @@
use crate::bspline::{BSpline, BSplineTrait};
use crate::bspline::{BSpline, BSplineTrait, DefaultValue};
use crate::error::Error;
use crate::filter::{fft, ifft};
use crate::metric::FixedMu;
@@ -811,7 +811,7 @@ impl<D: Dimension> Transform<D> {
})
}
pub fn interpolate<'a, const B: usize, T, A>(&self, image: A) -> Result<Array<f64, D>, Error>
pub fn interpolate<'a, const B: usize, T, A>(&self, image: A, default_value: &DefaultValue) -> Result<Array<f64, D>, Error>
where
A: AsArray<'a, T, D>,
D: Dimension,
@@ -820,12 +820,13 @@ impl<D: Dimension> Transform<D> {
{
let image = image.into();
let bspline = BSpline::<B, _>::new(&image);
bspline.interpolate(self)
bspline.interpolate(self, default_value)
}
pub fn interpolate_par<'a, const B: usize, T, A>(
&self,
image: A,
default_value: &DefaultValue,
) -> Result<Array<f64, D>, Error>
where
A: AsArray<'a, T, D>,
@@ -835,13 +836,14 @@ impl<D: Dimension> Transform<D> {
{
let image = image.into();
let bspline = BSpline::<B, _>::new(&image);
bspline.interpolate_par(self)
bspline.interpolate_par(self, default_value)
}
pub fn interpolate_order<'a, T, A>(
&self,
image: A,
order: usize,
default_value: &DefaultValue,
) -> Result<Array<f64, D>, Error>
where
A: AsArray<'a, T, D>,
@@ -849,12 +851,12 @@ impl<D: Dimension> Transform<D> {
T: 'a + Clone + AsPrimitive<f64>,
{
match order {
0 => Ok(self.interpolate::<0, _, _>(image)?),
1 => Ok(self.interpolate::<1, _, _>(image)?),
2 => Ok(self.interpolate::<2, _, _>(image)?),
3 => Ok(self.interpolate::<3, _, _>(image)?),
4 => Ok(self.interpolate::<4, _, _>(image)?),
5 => Ok(self.interpolate::<5, _, _>(image)?),
0 => Ok(self.interpolate::<0, _, _>(image, default_value)?),
1 => Ok(self.interpolate::<1, _, _>(image, default_value)?),
2 => Ok(self.interpolate::<2, _, _>(image, default_value)?),
3 => Ok(self.interpolate::<3, _, _>(image, default_value)?),
4 => Ok(self.interpolate::<4, _, _>(image, default_value)?),
5 => Ok(self.interpolate::<5, _, _>(image, default_value)?),
_ => Err(Error::InterpolationNotImplemented(order)),
}
}
@@ -863,6 +865,7 @@ impl<D: Dimension> Transform<D> {
&self,
image: A,
order: usize,
default_value: &DefaultValue,
) -> Result<Array<f64, D>, Error>
where
A: AsArray<'a, T, D>,
@@ -870,12 +873,12 @@ impl<D: Dimension> Transform<D> {
T: 'a + Clone + AsPrimitive<f64>,
{
match order {
0 => Ok(self.interpolate_par::<0, _, _>(image)?),
1 => Ok(self.interpolate_par::<1, _, _>(image)?),
2 => Ok(self.interpolate_par::<2, _, _>(image)?),
3 => Ok(self.interpolate_par::<3, _, _>(image)?),
4 => Ok(self.interpolate_par::<4, _, _>(image)?),
5 => Ok(self.interpolate_par::<5, _, _>(image)?),
0 => Ok(self.interpolate_par::<0, _, _>(image, default_value)?),
1 => Ok(self.interpolate_par::<1, _, _>(image, default_value)?),
2 => Ok(self.interpolate_par::<2, _, _>(image, default_value)?),
3 => Ok(self.interpolate_par::<3, _, _>(image, default_value)?),
4 => Ok(self.interpolate_par::<4, _, _>(image, default_value)?),
5 => Ok(self.interpolate_par::<5, _, _>(image, default_value)?),
_ => Err(Error::InterpolationNotImplemented(order)),
}
}
@@ -1045,8 +1048,8 @@ mod tests {
&[-0.8, 0.156],
);
let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
let image_b = transform.interpolate_par::<3, _, _>(image_a.view())?;
let image_c = transform.interpolate::<3, _, _>(image_a.view())?;
let image_b = transform.interpolate_par::<3, _, _>(image_a.view(), &Default::default())?;
let image_c = transform.interpolate::<3, _, _>(image_a.view(), &Default::default())?;
assert!(
image_b
@@ -1066,7 +1069,7 @@ mod tests {
&[-0.8, 0.156],
);
let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
let image_b = transform.interpolate::<3, _, _>(image_a.view())?;
let image_b = transform.interpolate::<3, _, _>(image_a.view(), &Default::default())?;
assert_ne!(image_b[[4000, 3000]], 0.0);
Ok(())
}
@@ -1080,7 +1083,7 @@ mod tests {
&[-0.8, 0.156],
);
let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
let image_b = transform.interpolate_par::<3, _, _>(image_a.view())?;
let image_b = transform.interpolate_par::<3, _, _>(image_a.view(), &Default::default())?;
assert_ne!(image_b[[4000, 3000]], 0.0);
Ok(())
}
@@ -1099,7 +1102,7 @@ mod tests {
.slice(s![.., 0])
.mapv(|i| i as f64);
let q = vec![0.85, 4.0];
let im_b = Transform::new(q).interpolate::<1, _, _>(&im_a)?;
let im_b = Transform::new(q).interpolate::<1, _, _>(&im_a, &Default::default())?;
// Test with ALL points (like grad_check uses)
let metric_all = MattesMetric::<ndarray::Ix1>::new(
@@ -1183,7 +1186,7 @@ mod tests {
)
.slice(s![.., 0])
.mapv(|i| i as f64);
let im_b = Transform::new(vec![0.85, 4.0]).interpolate::<1, _, _>(&im_a)?;
let im_b = Transform::new(vec![0.85, 4.0]).interpolate::<1, _, _>(&im_a, &Default::default())?;
// Smooth both images with sigma=8 like level 0
let sigma = Sigma::Absolute(vec![8.0]);
@@ -1292,7 +1295,7 @@ mod tests {
];
let q = vec![0.85, 2.0];
let moving = Transform::new(q.clone()).interpolate::<3, _, _>(&fixed)?;
let moving = Transform::new(q.clone()).interpolate::<3, _, _>(&fixed, &Default::default())?;
// The registration finds T such that im_b(T(x)) = im_a(x), which is the inverse of q
let q_inv = Transform::<ndarray::Ix1>::new(q.clone()).inverse()?;
@@ -1348,7 +1351,7 @@ mod tests {
.slice(s![.., 0])
.mapv(|i| i as f64);
let q = vec![0.85, 2.0];
let im_b = Transform::new(q.clone()).interpolate::<1, _, _>(&im_a)?;
let im_b = Transform::new(q.clone()).interpolate::<1, _, _>(&im_a, &Default::default())?;
let f = Sigma::Absolute(vec![1.0]).smooth(im_a.view())?;
let m = Sigma::Absolute(vec![1.0]).smooth(im_b.view())?;
@@ -1468,7 +1471,7 @@ mod tests {
let params = vec![c * sx, -s * sy + shx, s * sx + shy, c * sy, tx, ty];
let transform = Transform::<Ix2>::new_around_center(params.clone(), &center);
let moving = transform.interpolate::<3, _, _>(fixed.view())?;
let moving = transform.interpolate::<3, _, _>(fixed.view(), &Default::default())?;
let q_inv = transform.inverse()?.parameters;
@@ -1510,7 +1513,7 @@ mod tests {
)
.mapv(|i| i as f64);
let rotation = Transform::<Ix2>::from_rotation(f64::PI() / 4.0, &center);
let im_b = rotation.interpolate::<3, _, _>(im_a.view())?;
let im_b = rotation.interpolate::<3, _, _>(im_a.view(), &Default::default())?;
let q_inv = rotation.inverse()?.parameters;
let p = rotation.parameters.clone();
@@ -1592,7 +1595,7 @@ mod tests {
// Full resolution, sigma=4.0 — matching elastix FixedSmoothingImagePyramid (no downsampling)
for angle_deg in [5.0f64, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0] {
let rotation = Transform::<Ix2>::from_rotation(angle_deg.to_radians(), &center);
let im_b = rotation.interpolate::<3, _, _>(im_a.view())?;
let im_b = rotation.interpolate::<3, _, _>(im_a.view(), &Default::default())?;
let q_inv = rotation.inverse()?.parameters;
let p = rotation.parameters.clone();
let f = gaussian_smooth(im_a.view(), &[4.0, 4.0])?;
@@ -1727,7 +1730,7 @@ mod tests {
fn register_real_images2() -> Result<(), Box<dyn std::error::Error>> {
let fixed = read_tiff("test_files/fixed.tif")?;
let e = Transform::<Ix2>::new(vec![0.8, 0.0, 0.0, 1.0, 0.0, 0.0]).inverse()?;
let moving = e.interpolate::<3, _, _>(fixed.view())?;
let moving = e.interpolate::<3, _, _>(fixed.view(), &Default::default())?;
let t = Transform::<Ix2>::register(
fixed.view(),