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

- remove center and shape from Transform

Browse Source
This commit is contained in:
w.pomp committed 2026-10-22 11:42:42 +02:00
1 parent dc1d5d8364
commit aac61b981f
5 files changed
+67 -185
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

README.md
+1 -1
View File
@@ -336,7 +336,7 @@ fn main() -> Result<(), Error> {
| Module | Contents |
|----------------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| `image_registration::transform` | `Transform`, the affine transform type: construction (`new`, `new_with_center`, `from_translation`, `from_scaling`, `from_rotation`), algebra (`inverse`, `add_translation`, `add_scaling`, `with_rotation`, `with_rotation_around`, `Mul`), coordinate transforms (`transform_point`, `transform_points`, `matrix`, `dmatrix`), resampling (`interpolate`, `interpolate_par`), registration (`register`, `register_debug`, `register_affine`, `register_translation`), serialization (`to_file`, `from_file`), plus the free function `transform_point` and a blas-free `matrix_inverse`. |
| `image_registration::transform` | `Transform`, the affine transform type: construction (`new`, `new_around_center`, `from_translation`, `from_scaling`, `from_rotation`), algebra (`inverse`, `add_translation`, `add_scaling`, `with_rotation`, `with_rotation_around`, `Mul`), coordinate transforms (`transform_point`, `transform_points`, `matrix`, `dmatrix`), resampling (`interpolate`, `interpolate_par`), registration (`register`, `register_debug`, `register_affine`, `register_translation`), serialization (`to_file`, `from_file`), plus the free function `transform_point` and a blas-free `matrix_inverse`. |
| `image_registration::register` | `Registration`, `RegistrationStep`, `RegistrationResult`, `Optimizer` (LBFGS/ASGD). |
| `image_registration::metric` | `MattesMetric` (Mattes mutual information with analytic gradient), `SamplingArg` (Fixed/FixedAt/Random), `Sigma` (None/Absolute/Relative), `FixedMu`. |
| `image_registration::bspline` | `BSpline<N, D>` B-spline interpolation of arbitrary order, `BSplineTrait`, `BSplineMem`, Parzen kernel helpers. |
src/bspline.rs
+5 -3
View File
@@ -77,16 +77,17 @@ where
pub fn interpolate(&self, transform: &Transform<D>) -> Result<Array<f64, D>, Error> {
let mut mem = self.get_mem();
let shape_f = self.shape_f();
let center = shape_f.iter().map(|i| i / 2.0).collect::<Vec<_>>();
let mut new: Array<f64, D> =
Array::<f64, IxDyn>::zeros(self.shape()).into_dimensionality()?;
for (i, x) in new.indexed_iter_mut() {
let d = i.into_dimension();
let j = d.as_array_view();
let v = if let Some(k) = j.as_slice() {
transform_point(k, &transform.center, &transform.parameters)
transform_point(k, &center, &transform.parameters)
} else {
let k = j.to_vec();
transform_point(&k, &transform.center, &transform.parameters)
transform_point(&k, &center, &transform.parameters)
};
if in_shape(&shape_f, &v) {
*x = self.evaluate_at_continuous_index(&v, &mut mem)?;
@@ -99,13 +100,14 @@ where
pub fn interpolate_par(&self, transform: &Transform<D>) -> Result<Array<f64, D>, Error> {
let mem = self.get_mem();
let shape_f = self.shape_f();
let center = shape_f.iter().map(|i| i / 2.0).collect::<Vec<_>>();
let mut new: Array<f64, D> =
Array::<f64, IxDyn>::zeros(self.shape()).into_dimensionality()?;
new.par_indexed_iter_mut()
.try_fold(
|| mem.clone(),
|mut mem, (i, x)| {
let v = transform_point(&i, &transform.center, &transform.parameters);
let v = transform_point(&i, &center, &transform.parameters);
if shape_f
.iter()
.zip_eq(v.iter())
src/lib.rs
+5 -15
View File
@@ -53,11 +53,7 @@ mod tests {
#[test]
fn serialization() -> Result<(), Error> {
let file = NamedTempFile::new()?;
let t = Transform::<Ix2>::new_with_center(
vec![1.2, 0.3, -0.4, 0.9, 10.2, -9.5],
vec![59.5, 49.5],
vec![120, 100],
);
let t = Transform::<Ix2>::new(vec![1.2, 0.3, -0.4, 0.9, 10.2, -9.5]);
t.to_file(file.path().to_path_buf())?;
let s = Transform::from_file(file.path().to_path_buf())?;
assert_eq!(s, t);
@@ -66,12 +62,8 @@ mod tests {
#[test]
fn transform_point() -> Result<(), Box<dyn std::error::Error>> {
let transform = Transform::<Ix2>::new_with_center(
vec![1.0, 0.0, 0.0, 1.0, -120.0, 10.0],
vec![299.5, 399.5],
vec![600, 800],
);
let p = transform.transform_point(&[300, 400]);
let transform = Transform::<Ix2>::new(vec![1.0, 0.0, 0.0, 1.0, -120.0, 10.0]);
let p = transform.transform_point(&[300, 400], &[299.5, 399.5]);
assert_eq!(p, [180.0, 410.0]);
Ok(())
}
@@ -87,8 +79,7 @@ mod tests {
.collect::<Vec<_>>()
.try_into()
.unwrap();
let transform =
Transform::<Ix2>::new(transform_j.to_vec(), shape.to_vec()).with_rotation(1.0);
let transform = Transform::<Ix2>::new(transform_j.to_vec()).with_rotation(1.0, &center);
let transform_j = transform.parameters.clone().try_into().unwrap();
let c = [-0.8, 0.156];
let j = julia_image(&shape, &transform_j, &center, &c).mapv(|x| x as f64);
@@ -128,8 +119,7 @@ mod tests {
.collect::<Vec<_>>()
.try_into()
.unwrap();
let transform =
Transform::<Ix2>::new(transform_j.to_vec(), shape.to_vec()).with_rotation(1.0);
let transform = Transform::<Ix2>::new(transform_j.to_vec()).with_rotation(1.0, &center);
let transform_j = transform.parameters.clone().try_into().unwrap();
let c = [-0.8, 0.156];
let j = julia_image(&shape, &transform_j, &center, &c).mapv(|x| x as f64);
src/register.rs
+2 -8
View File
@@ -330,10 +330,7 @@ impl<D: Dimension> Registration<D> {
p_var = result.optimal_point;
}
}
Ok(Transform::<D>::new(
self.fixed_mu.combine(p_var.as_slice()),
fixed.shape().to_vec(),
))
Ok(Transform::<D>::new(self.fixed_mu.combine(p_var.as_slice())))
}
fn optimize_metric<M: ObjectiveFunction<f64>>(
@@ -480,9 +477,6 @@ impl<D: Dimension> Registration<D> {
converged: optimization_result.converged,
});
}
Ok((
Transform::<D>::new(p, fixed.shape().to_vec()),
registration_results,
))
Ok((Transform::<D>::new(p), registration_results))
}
}
src/transform.rs
+54 -158
View File
@@ -101,39 +101,14 @@ fn matrix_inverse(matrix: Array2<f64>) -> Array2<f64> {
inv
}
#[derive(Deserialize)]
struct TransformS {
parameters: Vec<f64>,
dparameters: Vec<f64>,
center: Vec<f64>,
shape: Vec<usize>,
}
impl<D: Dimension> From<TransformS> for Transform<D> {
fn from(t: TransformS) -> Self {
Self {
parameters: t.parameters,
dparameters: t.dparameters,
center: t.center,
shape: t.shape,
dimension: PhantomData,
}
}
}
/// a struct describing the transform
/// generic parameter N = # image dimensions
#[derive(Clone, Debug, Deserialize, Serialize, PartialEq)]
#[serde(from = "TransformS")]
pub struct Transform<D: Dimension = Ix2> {
pub struct Transform<D: Dimension = IxDyn> {
/// flattened NxN rotation matrix + N translation parameters
pub parameters: Vec<f64>,
/// error / significance on parameters
pub dparameters: Vec<f64>,
/// the point about which rotations are performed
pub center: Vec<f64>,
/// the shape of images for which this transform is meant
pub shape: Vec<usize>,
#[serde(skip)]
dimension: PhantomData<D>,
}
@@ -184,8 +159,6 @@ impl<D: Dimension> Mul<&Transform<D>> for &Transform<D> {
.chain(dm.slice(s![..self.ndim(), ..]).iter())
.cloned()
.collect(),
center: self.center.clone(),
shape: self.shape.clone(),
dimension: self.dimension,
}
}
@@ -236,8 +209,6 @@ impl<D: Dimension> Add<&Transform<D>> for &Transform<D> {
.chain(dm.slice(s![..self.ndim(), ..]).iter())
.cloned()
.collect(),
center: self.center.clone(),
shape: self.shape.clone(),
dimension: self.dimension,
}
}
@@ -288,8 +259,6 @@ impl<D: Dimension> Sub<&Transform<D>> for &Transform<D> {
.chain(dm.slice(s![..self.ndim(), ..]).iter())
.cloned()
.collect(),
center: self.center.clone(),
shape: self.shape.clone(),
dimension: self.dimension,
}
}
@@ -340,8 +309,6 @@ impl<D: Dimension> Mul<f64> for &Transform<D> {
.chain(dm.slice(s![..self.ndim(), ..]).iter())
.cloned()
.collect(),
center: self.center.clone(),
shape: self.shape.clone(),
dimension: self.dimension,
}
}
@@ -376,8 +343,6 @@ impl<D: Dimension> Div<f64> for &Transform<D> {
.chain(dm.slice(s![..self.ndim(), ..]).iter())
.cloned()
.collect(),
center: self.center.clone(),
shape: self.shape.clone(),
dimension: self.dimension,
}
}
@@ -387,14 +352,13 @@ impl<D: Dimension> Eq for Transform<D> {}
impl<D: Dimension> Hash for Transform<D> {
fn hash<H: Hasher>(&self, state: &mut H) {
for values in [&self.parameters, &self.dparameters, &self.center] {
for values in [&self.parameters, &self.dparameters] {
values.len().hash(state);
for v in values {
let v = if *v == 0.0 { 0.0 } else { *v };
v.to_bits().hash(state);
}
}
self.shape.hash(state);
}
}
@@ -410,51 +374,35 @@ impl<D: Dimension> Default for Transform<D> {
Self {
parameters,
dparameters: vec![0f64; ndim * ndim + ndim],
center: vec![-0.5; ndim],
shape: vec![0usize; ndim],
dimension: PhantomData,
}
}
}
impl<D: Dimension> Transform<D> {
/// parameters: flat NxN part of matrix + translation; center: center of rotation
pub fn new_with_center(parameters: Vec<f64>, center: Vec<f64>, shape: Vec<usize>) -> Self {
let ndim = if let Some(ndim) = D::NDIM {
ndim
} else {
center.len()
};
debug_assert_eq!(parameters.len(), ndim * ndim + ndim);
debug_assert_eq!(center.len(), ndim);
debug_assert_eq!(shape.len(), ndim);
/// parameters: flat NxN part of matrix
pub fn new(parameters: Vec<f64>) -> Self {
let dparameters = vec![0.0; parameters.len()];
Self {
parameters,
dparameters: vec![0.0; ndim * ndim + ndim],
center,
shape,
dparameters,
dimension: PhantomData,
}
}
pub fn new(parameters: Vec<f64>, shape: Vec<usize>) -> Self {
let ndim = if let Some(ndim) = D::NDIM {
ndim
} else {
shape.len()
};
debug_assert_eq!(parameters.len(), ndim * ndim + ndim);
debug_assert_eq!(shape.len(), ndim);
let center = shape.iter().map(|s| ((*s as f64) - 1.0) / 2.0).collect();
pub fn new_dparameters(parameters: Vec<f64>, dparameters: Vec<f64>) -> Self {
Self {
parameters,
dparameters: vec![0.0; ndim * ndim + ndim],
center,
shape,
dparameters,
dimension: PhantomData,
}
}
pub fn new_around_center(parameters: Vec<f64>, center: &[f64]) -> Self {
(Self::new(vec![1.0, 0.0, 0.0, 1.0, -center[0], -center[1]]) * Self::new(parameters))
* Self::new(vec![1.0, 0.0, 0.0, 1.0, center[0], center[1]])
}
/// find the affine transform which transforms moving into fixed
pub fn register_affine<'a, F, M, T>(fixed: F, moving: M) -> Result<Self, Error>
where
@@ -548,17 +496,10 @@ impl<D: Dimension> Transform<D> {
Self {
parameters,
dparameters: vec![0f64; ndim * ndim + ndim],
center: vec![0f64; ndim],
shape: vec![0usize; ndim],
dimension: PhantomData,
}
}
pub fn with_shape(mut self, shape: Vec<usize>) -> Self {
self.shape = shape;
self
}
/// create a transform from a scaling
pub fn from_scaling(scaling: &[f64]) -> Self {
let ndim = if let Some(ndim) = D::NDIM {
@@ -580,8 +521,6 @@ impl<D: Dimension> Transform<D> {
Self {
parameters,
dparameters: vec![0f64; ndim * ndim + ndim],
center: vec![0f64; ndim],
shape: vec![0usize; ndim],
dimension: PhantomData,
}
}
@@ -607,8 +546,6 @@ impl<D: Dimension> Transform<D> {
Self {
parameters,
dparameters: vec![0f64; ndim * ndim + ndim],
center: self.center.clone(),
shape: self.shape.clone(),
dimension: PhantomData,
}
}
@@ -646,8 +583,6 @@ impl<D: Dimension> Transform<D> {
Self {
parameters,
dparameters: vec![0f64; ndim * ndim + ndim],
center: self.center.clone(),
shape: self.shape.clone(),
dimension: PhantomData,
}
}
@@ -657,11 +592,6 @@ impl<D: Dimension> Transform<D> {
Self::from_translation(translation) * self
}
pub fn with_center(mut self, center: Vec<f64>) -> Self {
self.center = center;
self
}
/// read a transform from a file
pub fn from_file(path: PathBuf) -> Result<Self, Error> {
let file = File::open(path)?;
@@ -693,16 +623,20 @@ impl<D: Dimension> Transform<D> {
/// get coordinates resulting from transforming input coordinates, coordinate must have N
/// columns: x, y, z, ...
pub fn transform_point<'a, T>(&'a self, point: &'a [T]) -> Vec<f64>
pub fn transform_point<'a, T>(&'a self, point: &'a [T], center: &'a [f64]) -> Vec<f64>
where
T: AsPrimitive<f64>,
{
transform_point(point, &self.center, &self.parameters)
transform_point(point, center, &self.parameters)
}
/// get coordinates resulting from transforming input coordinates, coordinates must have N
/// columns: x, y, z, ...
pub fn transform_points<'a, A, T>(&'a self, points: A) -> Result<Array2<f64>, Error>
pub fn transform_points<'a, A, T>(
&'a self,
points: A,
center: &'a [f64],
) -> Result<Array2<f64>, Error>
where
T: AsPrimitive<f64>,
A: AsArray<'a, T, Ix2>,
@@ -717,9 +651,9 @@ impl<D: Dimension> Transform<D> {
.into_iter()
.flat_map(|p| {
if let Some(q) = p.as_slice() {
transform_point(q, &self.center, &self.parameters)
transform_point(q, center, &self.parameters)
} else {
transform_point(&p.to_vec(), &self.center, &self.parameters)
transform_point(&p.to_vec(), center, &self.parameters)
}
})
.collect::<Vec<_>>();
@@ -830,24 +764,10 @@ impl<D: Dimension> Transform<D> {
Ok(Self {
parameters,
dparameters: vec![0f64; self.ndim() * self.ndim() + self.ndim()],
center: self.center.clone(),
shape: self.shape.clone(),
dimension: PhantomData,
})
}
/// adapt the transform to a new center and shape
pub fn adapt(&mut self, center: &[f64], shape: &[usize]) {
self.center = self
.shape
.iter()
.zip_eq(shape.iter())
.zip_eq(center)
.map(|((a, b), o)| o + (((a - b) as f64) / 2.0))
.collect();
self.shape = shape.to_vec();
}
pub fn interpolate<'a, const B: usize, T, A>(&self, image: A) -> Result<Array<f64, D>, Error>
where
A: AsArray<'a, T, D>,
@@ -922,7 +842,7 @@ impl<D: Dimension> Transform<D> {
if let Some(d) = D::NDIM {
d
} else {
self.center.len()
self.parameters.len().isqrt()
}
}
@@ -932,9 +852,7 @@ impl<D: Dimension> Transform<D> {
Ok(Transform {
parameters: self.parameters,
dparameters: self.dparameters,
center: self.center,
shape: self.shape,
dimension: Default::default(),
dimension: PhantomData,
})
} else {
Err(Error::DimensionalityMismatch(d, self.ndim()))
@@ -943,9 +861,7 @@ impl<D: Dimension> Transform<D> {
Ok(Transform {
parameters: self.parameters,
dparameters: self.dparameters,
center: self.center,
shape: self.shape,
dimension: Default::default(),
dimension: PhantomData,
})
}
}
@@ -954,9 +870,7 @@ impl<D: Dimension> Transform<D> {
Transform {
parameters: self.parameters,
dparameters: self.dparameters,
center: self.center,
shape: self.shape,
dimension: Default::default(),
dimension: PhantomData,
}
}
@@ -968,7 +882,7 @@ impl<D: Dimension> Transform<D> {
p
}
pub fn add_dimension(&self, axis: usize, size: usize) -> Transform<D::Larger> {
pub fn add_dimension(&self, axis: usize) -> Transform<D::Larger> {
let ndim = self.ndim() + 1;
let mut parameters = Self::parameters_unit(ndim);
let mut dparameters = Self::parameters_unit(ndim);
@@ -984,16 +898,9 @@ impl<D: Dimension> Transform<D> {
dparameters[j] = *p;
}
let mut center = self.center.clone();
center.insert(axis, ((size as f64) - 1.0) / 2.0);
let mut shape = self.shape.clone();
shape.insert(axis, size);
Transform {
parameters,
dparameters,
center,
shape,
dimension: Default::default(),
}
}
@@ -1013,10 +920,6 @@ impl<D: Dimension> Transform<D> {
}
new.parameters[n + a] = self.parameters[n + b];
new.dparameters[n + a] = self.dparameters[n + b];
new.center[a] = self.center[b];
new.center[b] = self.center[a];
new.shape[a] = self.shape[b];
new.shape[b] = self.shape[a];
new
}
}
@@ -1024,8 +927,8 @@ impl<D: Dimension> Transform<D> {
impl Transform<Ix2> {
/// create a transform from a rotation in radians
pub fn from_rotation(theta: f64, center: &[f64]) -> Self {
Self {
parameters: vec![
Self::new_around_center(
vec![
theta.cos(),
-theta.sin(),
theta.sin(),
@@ -1033,20 +936,12 @@ impl Transform<Ix2> {
0.0,
0.0,
],
dparameters: vec![0f64; 6],
center: center.to_vec(),
shape: vec![0usize; 2],
dimension: PhantomData,
}
center,
)
}
/// multiply self with a rotation matrix
pub fn with_rotation(&self, theta: f64) -> Self {
Self::from_rotation(theta, &self.center) * self
}
/// multiply self with a rotation matrix
pub fn with_rotation_around(&self, theta: f64, center: &[f64]) -> Self {
pub fn with_rotation(&self, theta: f64, center: &[f64]) -> Self {
Self::from_rotation(theta, center) * self
}
}
@@ -1161,8 +1056,7 @@ mod tests {
.slice(s![.., 0])
.mapv(|i| i as f64);
let q = vec![0.85, 4.0];
let im_b =
Transform::new(q.clone(), vec![im_a.shape()[0]]).interpolate::<1, _, _>(&im_a)?;
let im_b = Transform::new(q).interpolate::<1, _, _>(&im_a)?;
// Test with ALL points (like grad_check uses)
let metric_all = MattesMetric::<ndarray::Ix1>::new(
@@ -1246,8 +1140,7 @@ mod tests {
)
.slice(s![.., 0])
.mapv(|i| i as f64);
let im_b =
Transform::new(vec![0.85, 4.0], vec![im_a.shape()[0]]).interpolate::<1, _, _>(&im_a)?;
let im_b = Transform::new(vec![0.85, 4.0]).interpolate::<1, _, _>(&im_a)?;
// Smooth both images with sigma=8 like level 0
let sigma = Sigma::Absolute(vec![8.0]);
@@ -1356,11 +1249,10 @@ mod tests {
];
let q = vec![0.85, 2.0];
let moving =
Transform::new(q.clone(), vec![fixed.shape()[0]]).interpolate::<3, _, _>(&fixed)?;
let moving = Transform::new(q.clone()).interpolate::<3, _, _>(&fixed)?;
// 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(), vec![fixed.shape()[0]]).inverse()?;
let q_inv = Transform::<ndarray::Ix1>::new(q.clone()).inverse()?;
// Use default steps with ASGD — the multi-resolution pyramid converges
// reliably for this 1D affine problem.
@@ -1369,9 +1261,14 @@ mod tests {
// Compare transformed coordinates of all pixels
let mut sum_diff = 0.0;
let mut count = 0;
let center = fixed
.shape()
.iter()
.map(|&i| (i as f64) / 2.0)
.collect::<Vec<_>>();
for point in 0..fixed.len() {
let t_point = t.transform_point(&[point]);
let e_point = q_inv.transform_point(&[point]);
let t_point = t.transform_point(&[point], &center);
let e_point = q_inv.transform_point(&[point], &center);
let diff_sq = (t_point[0] - e_point[0]).powi(2);
sum_diff += diff_sq.sqrt();
count += 1;
@@ -1408,8 +1305,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(), vec![im_a.shape()[0]]).interpolate::<1, _, _>(&im_a)?;
let im_b = Transform::new(q.clone()).interpolate::<1, _, _>(&im_a)?;
let f = Sigma::Absolute(vec![1.0]).smooth(im_a.view())?;
let m = Sigma::Absolute(vec![1.0]).smooth(im_b.view())?;
@@ -1423,7 +1319,7 @@ mod tests {
)?
.with_fixed_mu(FixedMu::new_none(1));
let q_inv = Transform::<ndarray::Ix1>::new(q.clone(), vec![im_a.shape()[0]])
let q_inv = Transform::<ndarray::Ix1>::new(q.clone())
.inverse()?
.parameters;
@@ -1528,8 +1424,7 @@ mod tests {
// [m00, m01, m10, m11, tx, ty]
let params = vec![c * sx, -s * sy + shx, s * sx + shy, c * sy, tx, ty];
let transform =
Transform::<Ix2>::new_with_center(params.clone(), center.to_vec(), shape.to_vec());
let transform = Transform::<Ix2>::new_around_center(params.clone(), &center);
let moving = transform.interpolate::<3, _, _>(fixed.view())?;
let q_inv = transform.inverse()?.parameters;
@@ -1745,17 +1640,18 @@ mod tests {
None,
)?;
let expected_transform = Transform::<Ix2>::new(expected.to_vec(), fixed.shape().to_vec());
let expected_transform = Transform::<Ix2>::new(expected.to_vec());
// Compare transformed coordinates of all pixels
let shape = fixed.shape();
let mut sum_diff = 0.0;
let mut count = 0;
let center = shape.iter().map(|&i| (i as f64) / 2.0).collect::<Vec<_>>();
for row in 0..shape[0] {
for col in 0..shape[1] {
let point = [row as f64, col as f64];
let t_point = t.transform_point(&point);
let e_point = expected_transform.transform_point(&point);
let t_point = t.transform_point(&point, &center);
let e_point = expected_transform.transform_point(&point, &center);
let diff_sq = (t_point[0] - e_point[0]).powi(2) + (t_point[1] - e_point[1]).powi(2);
sum_diff += diff_sq.sqrt();
count += 1;
@@ -1787,8 +1683,7 @@ mod tests {
#[test]
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], fixed.shape().to_vec())
.inverse()?;
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 t = Transform::<Ix2>::register(
@@ -1803,13 +1698,14 @@ mod tests {
// Compare transformed coordinates of all pixels
let shape = fixed.shape();
let center = shape.iter().map(|&i| (i as f64) / 2.0).collect::<Vec<_>>();
let mut sum_diff = 0.0;
let mut count = 0;
for row in 0..shape[0] {
for col in 0..shape[1] {
let point = [row as f64, col as f64];
let t_point = t.transform_point(&point);
let e_point = e_inv.transform_point(&point);
let t_point = t.transform_point(&point, &center);
let e_point = e_inv.transform_point(&point, &center);
let diff_sq = (t_point[0] - e_point[0]).powi(2) + (t_point[1] - e_point[1]).powi(2);
sum_diff += diff_sq.sqrt();
count += 1;