diff --git a/README.md b/README.md index 022e7a4..c77b00c 100644 --- a/README.md +++ b/README.md @@ -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` B-spline interpolation of arbitrary order, `BSplineTrait`, `BSplineMem`, Parzen kernel helpers. | diff --git a/src/bspline.rs b/src/bspline.rs index 0e36211..7d333e2 100644 --- a/src/bspline.rs +++ b/src/bspline.rs @@ -77,16 +77,17 @@ where pub fn interpolate(&self, transform: &Transform) -> Result, Error> { let mut mem = self.get_mem(); let shape_f = self.shape_f(); + let center = shape_f.iter().map(|i| i / 2.0).collect::>(); let mut new: Array = Array::::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, ¢er, &transform.parameters) } else { let k = j.to_vec(); - transform_point(&k, &transform.center, &transform.parameters) + transform_point(&k, ¢er, &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) -> Result, Error> { let mem = self.get_mem(); let shape_f = self.shape_f(); + let center = shape_f.iter().map(|i| i / 2.0).collect::>(); let mut new: Array = Array::::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, ¢er, &transform.parameters); if shape_f .iter() .zip_eq(v.iter()) diff --git a/src/lib.rs b/src/lib.rs index 986f0e9..453fe9d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -53,11 +53,7 @@ mod tests { #[test] fn serialization() -> Result<(), Error> { let file = NamedTempFile::new()?; - let t = Transform::::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::::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> { - let transform = Transform::::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::::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::>() .try_into() .unwrap(); - let transform = - Transform::::new(transform_j.to_vec(), shape.to_vec()).with_rotation(1.0); + let transform = Transform::::new(transform_j.to_vec()).with_rotation(1.0, ¢er); let transform_j = transform.parameters.clone().try_into().unwrap(); let c = [-0.8, 0.156]; let j = julia_image(&shape, &transform_j, ¢er, &c).mapv(|x| x as f64); @@ -128,8 +119,7 @@ mod tests { .collect::>() .try_into() .unwrap(); - let transform = - Transform::::new(transform_j.to_vec(), shape.to_vec()).with_rotation(1.0); + let transform = Transform::::new(transform_j.to_vec()).with_rotation(1.0, ¢er); let transform_j = transform.parameters.clone().try_into().unwrap(); let c = [-0.8, 0.156]; let j = julia_image(&shape, &transform_j, ¢er, &c).mapv(|x| x as f64); diff --git a/src/register.rs b/src/register.rs index 2e29242..babd86c 100644 --- a/src/register.rs +++ b/src/register.rs @@ -330,10 +330,7 @@ impl Registration { p_var = result.optimal_point; } } - Ok(Transform::::new( - self.fixed_mu.combine(p_var.as_slice()), - fixed.shape().to_vec(), - )) + Ok(Transform::::new(self.fixed_mu.combine(p_var.as_slice()))) } fn optimize_metric>( @@ -480,9 +477,6 @@ impl Registration { converged: optimization_result.converged, }); } - Ok(( - Transform::::new(p, fixed.shape().to_vec()), - registration_results, - )) + Ok((Transform::::new(p), registration_results)) } } diff --git a/src/transform.rs b/src/transform.rs index cae4a10..d3aa10c 100644 --- a/src/transform.rs +++ b/src/transform.rs @@ -101,39 +101,14 @@ fn matrix_inverse(matrix: Array2) -> Array2 { inv } -#[derive(Deserialize)] -struct TransformS { - parameters: Vec, - dparameters: Vec, - center: Vec, - shape: Vec, -} - -impl From for Transform { - 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 { +pub struct Transform { /// flattened NxN rotation matrix + N translation parameters pub parameters: Vec, /// error / significance on parameters pub dparameters: Vec, - /// the point about which rotations are performed - pub center: Vec, - /// the shape of images for which this transform is meant - pub shape: Vec, #[serde(skip)] dimension: PhantomData, } @@ -184,8 +159,6 @@ impl Mul<&Transform> for &Transform { .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 Add<&Transform> for &Transform { .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 Sub<&Transform> for &Transform { .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 Mul for &Transform { .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 Div for &Transform { .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 Eq for Transform {} impl Hash for Transform { fn hash(&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 Default for Transform { Self { parameters, dparameters: vec![0f64; ndim * ndim + ndim], - center: vec![-0.5; ndim], - shape: vec![0usize; ndim], dimension: PhantomData, } } } impl Transform { - /// parameters: flat NxN part of matrix + translation; center: center of rotation - pub fn new_with_center(parameters: Vec, center: Vec, shape: Vec) -> 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) -> 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, shape: Vec) -> 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, dparameters: Vec) -> Self { Self { parameters, - dparameters: vec![0.0; ndim * ndim + ndim], - center, - shape, + dparameters, dimension: PhantomData, } } + pub fn new_around_center(parameters: Vec, 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 where @@ -548,17 +496,10 @@ impl Transform { 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) -> 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 Transform { Self { parameters, dparameters: vec![0f64; ndim * ndim + ndim], - center: vec![0f64; ndim], - shape: vec![0usize; ndim], dimension: PhantomData, } } @@ -607,8 +546,6 @@ impl Transform { Self { parameters, dparameters: vec![0f64; ndim * ndim + ndim], - center: self.center.clone(), - shape: self.shape.clone(), dimension: PhantomData, } } @@ -646,8 +583,6 @@ impl Transform { Self { parameters, dparameters: vec![0f64; ndim * ndim + ndim], - center: self.center.clone(), - shape: self.shape.clone(), dimension: PhantomData, } } @@ -657,11 +592,6 @@ impl Transform { Self::from_translation(translation) * self } - pub fn with_center(mut self, center: Vec) -> Self { - self.center = center; - self - } - /// read a transform from a file pub fn from_file(path: PathBuf) -> Result { let file = File::open(path)?; @@ -693,16 +623,20 @@ impl Transform { /// 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 + pub fn transform_point<'a, T>(&'a self, point: &'a [T], center: &'a [f64]) -> Vec where T: AsPrimitive, { - 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, Error> + pub fn transform_points<'a, A, T>( + &'a self, + points: A, + center: &'a [f64], + ) -> Result, Error> where T: AsPrimitive, A: AsArray<'a, T, Ix2>, @@ -717,9 +651,9 @@ impl Transform { .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::>(); @@ -830,24 +764,10 @@ impl Transform { 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, Error> where A: AsArray<'a, T, D>, @@ -922,7 +842,7 @@ impl Transform { if let Some(d) = D::NDIM { d } else { - self.center.len() + self.parameters.len().isqrt() } } @@ -932,9 +852,7 @@ impl Transform { 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 Transform { Ok(Transform { parameters: self.parameters, dparameters: self.dparameters, - center: self.center, - shape: self.shape, - dimension: Default::default(), + dimension: PhantomData, }) } } @@ -954,9 +870,7 @@ impl Transform { Transform { parameters: self.parameters, dparameters: self.dparameters, - center: self.center, - shape: self.shape, - dimension: Default::default(), + dimension: PhantomData, } } @@ -968,7 +882,7 @@ impl Transform { p } - pub fn add_dimension(&self, axis: usize, size: usize) -> Transform { + pub fn add_dimension(&self, axis: usize) -> Transform { let ndim = self.ndim() + 1; let mut parameters = Self::parameters_unit(ndim); let mut dparameters = Self::parameters_unit(ndim); @@ -984,20 +898,13 @@ impl Transform { 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(), } } - + pub fn swap_axes(&self, a: usize, b: usize) -> Self { let mut new = self.clone(); let n = self.ndim(); @@ -1013,10 +920,6 @@ impl Transform { } 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 Transform { impl Transform { /// 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 { 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::::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::::new(q.clone(), vec![fixed.shape()[0]]).inverse()?; + let q_inv = Transform::::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::>(); 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], ¢er); + let e_point = q_inv.transform_point(&[point], ¢er); 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::::new(q.clone(), vec![im_a.shape()[0]]) + let q_inv = Transform::::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::::new_with_center(params.clone(), center.to_vec(), shape.to_vec()); + let transform = Transform::::new_around_center(params.clone(), ¢er); let moving = transform.interpolate::<3, _, _>(fixed.view())?; let q_inv = transform.inverse()?.parameters; @@ -1745,17 +1640,18 @@ mod tests { None, )?; - let expected_transform = Transform::::new(expected.to_vec(), fixed.shape().to_vec()); + let expected_transform = Transform::::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::>(); 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, ¢er); + let e_point = expected_transform.transform_point(&point, ¢er); 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> { 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]).inverse()?; let moving = e.interpolate::<3, _, _>(fixed.view())?; let t = Transform::::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::>(); 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, ¢er); + let e_point = e_inv.transform_point(&point, ¢er); 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;