- a comment and a TODO

- implement dparameter estimation

Cargo.toml

@@ -1,13 +1,12 @@
 [package]
 name = "image-registration"
-version = "0.1.0"
+version = "0.1.1"
 edition = "2024"
 rust-version = "1.88.0"
 authors = ["Wim Pomp <w.pomp@nki.nl>"]
 license = "MIT OR Apache-2.0"
 description = "Find and apply affine transformations."
-homepage = "https://git.pomppervova.nl/wim/image-registration/src/branch/rs"
+repository = "https://git.pomppervova.nl/wim/image-registration"
-repository = "https://git.pomppervova.nl/wim/image-registration/src/branch/rs"
 documentation = "https://docs.rs/image-registration"
 readme = "README.md"
 keywords = ["affine", "transformation", "ndarray"]

README.md

@@ -4,7 +4,7 @@ Find an ND affine transformation that registers an ND test image to a fixed imag
 other images.
 
 `image-registration` is a pure-Rust, dependency-light implementation of intensity-based image registration in the spirit
-of [SimpleElastix]([https://elastix.lumc.nl](https://simpleelastix.github.io/)). It finds an affine (rigid + scale + shear + translation) transform
+of [SimpleElastix](https://simpleelastix.github.io/). It finds an affine (rigid + scale + shear + translation) transform
 that aligns a *moving* image to a *fixed* image by maximizing
 the [Mattes mutual-information](https://doi.org/10.1109/TMI.2003.809072) metric with a B-spline Parzen
 window, using a multi-resolution Gaussian image pyramid and a hybrid ASGD + L-BFGS optimization strategy.
@@ -12,6 +12,10 @@ window, using a multi-resolution Gaussian image pyramid and a hybrid ASGD + L-BF
 Everything operates on [`ndarray`](https://docs.rs/ndarray) arrays, works in arbitrary dimensions (1D, 2D, 3D, …), and
 is parallelized with [rayon](https://docs.rs/rayon).
 
+So far, 2D registration seems to work, but 3D and higher dimensions often do not. Also, this code is still a lot slower
+than (Simple)Elastix. The biggest benefit of this crate is that it's Rust and does not require a lengthy compilation of
+C++ code.
+
 ## Features
 
 - **Affine registration** in any number of dimensions — rotation, scale, shear, and translation parameters
@@ -36,13 +40,13 @@ Register a moving image to a fixed image and obtain the affine transform:
 ```rust
 use image_registration::error::Error;
 use image_registration::transform::Transform;
-use ndarray::Array2;
+use ndarray::{Array2, Ix2};
 
-fn register(fixed: &Array2<f64>, moving: &Array2<f64>) -> Result<Transform, Error> {
+fn register(fixed: &Array2<f64>, moving: &Array2<f64>) -> Result<Transform<Ix2>, Error> {
     // All 6 parameters of the 2D affine transform are free (None = optimize).
     // steps=None uses the default multi-resolution pyramid; initial_guess=None
     // initializes by aligning the geometric centers of the two images.
-    Transform::<ndarray::Ix2>::register(
+    Transform::<Ix2>::register(
         fixed.view(),
         moving.view(),
         vec![None; 6],
@@ -55,9 +59,10 @@ fn register(fixed: &Array2<f64>, moving: &Array2<f64>) -> Result<Transform, Erro
 The returned `Transform` maps coordinates from the *fixed* image to the *moving* image, i.e.
 `moving(transform_point(x)) ≈ fixed(x)`. To warp the moving image into the fixed image's frame:
 
-```rust
+```ignore
 let transform = register(&fixed, &moving)?;
-let registered = transform.interpolate_par::<3, _, _>(moving.view())?; // cubic B-spline, parallel
+// out-of-range pixels are filled with the default value (0.0 by default)
+let registered = transform.interpolate_par::<3, _, _>(moving.view(), &Default::default())?; // cubic B-spline, parallel
 ```
 
 ## Concepts
@@ -71,15 +76,17 @@ in total:
 - 2D: `[m00, m01, m10, m11, tx, ty]` (row-major linear part, then translation)
 - 3D: `[m00, m01, m02, m10, m11, m12, m20, m21, m22, tx, ty, tz]`
 
-The transformation is applied about a center saved in ```Transform```:
+The transformation is applied about a center `c` that is supplied where the transform is used:
 
 ```text
 v = A·(p − c) + t + c
 ```
 
 where `A` is the linear part, `t` the translation, `c` the center of the transformation, and `p` the input coordinates.
-`Transform::new(parameters, shape)` sets the center to `(shape − 1) / 2` automatically; `Transform::new_with_center`
+`transform_point`/`transform_points` take the center as an argument, and resampling (`interpolate`, `interpolate_par`)
-lets you override it.
+uses the geometric center of the image, `(shape − 1) / 2`. `Transform::new_around_center(parameters, center)`
+re-expresses the parameters so that the transformation acts about a chosen center (this is how
+`from_rotation`/`with_rotation` rotate about a given point).
 
 ### Registration convention
 
@@ -101,19 +108,24 @@ fn main() -> Result<(), Error> {
     let rot = Transform::<Ix2>::from_rotation(std::f64::consts::FRAC_PI_4, &[99.5, 99.5]);
 
     // Compose: translate, then scale, then rotate
-    let t = rot.with_scaling(&[0.9, 1.1]).with_translation(&[10.0, -5.0]);
+    let t = rot.add_scaling(&[0.9, 1.1]).with_translation(&[10.0, -5.0]);
 
     // Transform a single point (the center is subtracted and re-added)
-    let p = t.transform_point(&[200, 300]);
+    let p = t.transform_point(&[200.0, 300.0], &[99.5, 99.5]);
 
     // Transform many points at once (rows of an N×2 array)
     let points = ndarray::array![[0.0, 0.0], [199.0, 199.0]];
-    let pts = t.transform_points(points.view())?;
+    let pts = t.transform_points(points.view(), &[99.5, 99.5])?;
 
     // Matrix form (3×3 homogeneous) and exact inverse
     let m = t.matrix();
     let inv = t.inverse()?;
-    assert!(m.dot(&inv.matrix()).iter().all(|x| (x - 1.0).abs() < 1e-9));
+    let identity = ndarray::Array2::<f64>::eye(3);
+    assert!(m
+        .dot(&inv.matrix())
+        .iter()
+        .zip(identity.iter())
+        .all(|(x, y)| (x - y).abs() < 1e-9));
     Ok(())
 }
 ```
@@ -132,11 +144,11 @@ fn main() -> Result<(), Error> {
     let image = ndarray::Array2::<f64>::zeros((200, 150));
     let transform = Transform::<ndarray::Ix2>::from_rotation(0.3, &[99.5, 74.5]);
 
-    // Order-1 (linear) interpolation, single-threaded
+    // Order-1 (linear) interpolation, single-threaded (out-of-range pixels → 0.0)
-    let warped_linear = transform.interpolate::<1, _, _>(image.view())?;
+    let warped_linear = transform.interpolate::<1, _, _>(image.view(), &Default::default())?;
 
     // Order-3 (cubic) B-spline interpolation, parallel
-    let warped_cubic = transform.interpolate_par::<3, _, _>(image.view())?;
+    let warped_cubic = transform.interpolate_par::<3, _, _>(image.view(), &Default::default())?;
     Ok(())
 }
 ```
@@ -146,7 +158,7 @@ fn main() -> Result<(), Error> {
 Create a moving image by applying a known transform to a fixed image, register, and check that the recovered transform
 matches the inverse of the ground-truth transform (adapted from the crate's `register2_random_affine` test):
 
-```rust
+```no_run
 use image_registration::error::Error;
 use image_registration::julia_image;
 use image_registration::transform::Transform;
@@ -163,11 +175,10 @@ fn main() -> Result<(), Error> {
     let theta = 10.0f64.to_radians();
     let (s, c) = theta.sin_cos();
     let params = vec![c * 1.05, -s, s, c * 0.95, 12.0, -8.0];
-    let ground_truth =
+    let ground_truth = Transform::<Ix2>::new_around_center(params.clone(), &center);
-        Transform::<Ix2>::new_with_center(params.clone(), center.to_vec(), shape.to_vec());
 
     // Moving image = ground truth applied to the fixed image
-    let moving = ground_truth.interpolate::<3, _, _>(fixed.view())?;
+    let moving = ground_truth.interpolate::<3, _, _>(fixed.view(), &Default::default())?;
 
     // Register (all 6 parameters free, default steps)
     let t = Transform::<Ix2>::register(fixed.view(), moving.view(), vec![None; 6], None, None)?;
@@ -185,9 +196,13 @@ fn main() -> Result<(), Error> {
 }
 ```
 
+> **Note:** registration uses random point sampling on its coarser levels, so whether `sse < 1.0` holds for a given
+> image pair and transform is stochastic; this example may need a few tries (or a simpler transform, e.g. pure
+> translation) to converge. The crate's `register2_random_affine` unit test is a verified case for this path.
+
 ### 4. Registration with inspection of every optimization step
 
-```rust
+```ignore
 use image_registration::error::Error;
 use image_registration::transform::Transform;
 
@@ -212,13 +227,13 @@ fn main() -> Result<(), Error> {
 following uses two levels: a coarse level with strong Gaussian smoothing and the robust ASGD optimizer, then a fine
 level with L-BFGS for accuracy.
 
-```rust
+```ignore
 use image_registration::error::Error;
 use image_registration::metric::{SamplingArg, Sigma};
 use image_registration::register::{Optimizer, Registration, RegistrationStep};
 
 fn main() -> Result<(), Error> {
-    let mut coarse = RegistrationStep::new(
+    let coarse = RegistrationStep::new(
         Sigma::Absolute(vec![4.0, 4.0]), // Gaussian sigma per dimension
         SamplingArg::Fixed(3000),        // 3000 sampled points (cached)
         32,                              // histogram bins
@@ -226,10 +241,10 @@ fn main() -> Result<(), Error> {
         0.05,                            // edge fraction
         512,                             // max iterations
         1.0,                             // learning rate
+        Optimizer::ASGD,                 // robust to noisy gradients
     );
-    coarse.optimizer = Optimizer::ASGD; // robust to noisy gradients
 
-    let mut fine = RegistrationStep::new(
+    let fine = RegistrationStep::new(
         Sigma::Absolute(vec![0.5, 0.5]),
         SamplingArg::Fixed(3000),
         32,
@@ -237,9 +252,10 @@ fn main() -> Result<(), Error> {
         0.05,
         2048,
         1.0,
-    ); // optimizer stays LBFGS (the default)
+        Optimizer::LBFGS, // the default optimizer
+    );
 
-    let mut reg = Registration::new(vec![None; 6]);
+    let mut reg = Registration::new(vec![None; 6])?;
     reg.set_steps(vec![coarse, fine]);
     let t = reg.register(fixed.view(), moving.view())?;
     Ok(())
@@ -254,7 +270,7 @@ fn main() -> Result<(), Error> {
 Pass a `Vec<Option<f64>>` as `fixed_mu`: `None` optimizes a parameter, `Some(v)` keeps it fixed at `v`. For example,
 register for translation only:
 
-```rust
+```ignore
 use image_registration::error::Error;
 use image_registration::transform::Transform;
 
@@ -284,18 +300,24 @@ Build a `MattesMetric` and evaluate the (negative) mutual information and its gr
 - `FixedAt(points)` — an explicit list of sample positions (fully deterministic; recommended for reproducible results).
 
 ```rust
-use algos::ObjectiveFunction;
+use algos::ObjectiveFunction; // direct dependency of image-registration; add it to your crate to name the trait
-use image_registration::bspline::BSpline;
+use image_registration::bspline::{BSpline, BSplineTrait};
 use image_registration::error::Error;
+use image_registration::julia_image;
 use image_registration::metric::{MattesMetric, SamplingArg, Sigma};
+use image_registration::transform::Transform;
 
 fn main() -> Result<(), Error> {
-    // `fixed` and `moving` are Array2<f64>
+    // `fixed` and `moving` are Array2<f64>; moving is a rotated + shifted copy of fixed
-    let fixed = ndarray::Array2::<f64>::zeros((64, 64));
+    let center = [31.5f64, 31.5];
-    let moving = ndarray::Array2::<f64>::zeros((64, 64));
+    let fixed = julia_image(&[64, 64], &[1.0, 0.0, 0.0, 1.0, 0.0, 0.0], &center, &[-0.8, 0.156])
+        .mapv(|i| i as f64);
+    let moving = Transform::<ndarray::Ix2>::from_rotation(10.0f64.to_radians(), &center)
+        .interpolate_par::<3, _, _>(fixed.view(), &Default::default())?;
 
-    let f = Sigma::Absolute(vec![2.0, 2.0]).smooth(fixed.view())?;
+    let sigma = Sigma::Absolute(vec![2.0, 2.0]);
-    let m = Sigma::Absolute(vec![2.0, 2.0]).smooth(moving.view())?;
+    let f = sigma.smooth(fixed.view())?;
+    let m = sigma.smooth(moving.view())?;
 
     let metric = MattesMetric::new(
         BSpline::<0, _>::new(f.view()), // fixed image B-spline (nearest-neighbor)
@@ -336,11 +358,11 @@ 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`, `Mul`), coordinate transforms (`transform_point`, `transform_points`, `matrix`, `dmatrix`), resampling (`interpolate`, `interpolate_par`, `interpolate_order`, `interpolate_par_order`), registration (`register`, `register_debug`, `register_affine`, `register_translation`, `register_translation_fft`), serialization (`to_file`, `from_file`, and serde), plus the free function `transform_point`. |
-| `image_registration::register`         | `Registration`, `RegistrationStep`, `RegistrationResult`, `Optimizer` (LBFGS/ASGD).                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
+| `image_registration::register`         | `Registration` (`new_affine`, `new_translation`, `with_steps`, `with_optimizer`, …), `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::metric`           | `MattesMetric` (Mattes mutual information with analytic gradient; `evaluate`/`gradient` via the `algos::ObjectiveFunction` trait), `new_from_arrays`, `residual_fit`, `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.                                                                                                                                                                                                                                                                                                                                                                                                                                                                              |
+| `image_registration::bspline`          | `BSpline<N, D>` B-spline interpolation of arbitrary order, `BSplineTrait`, `BSplineMem`, `DefaultValue` (fill value for out-of-range pixels), Parzen kernel helpers.                                                                                                                                                                                                                                                                                                                                                                                                                                                                              |
-| `image_registration::filter`           | Gaussian smoothing and FFT helpers: `gaussian_smooth`, `gaussian_kernel`, `fft`, `ifft`, `fft_freq`, `fft_shift`.                                                                                                                                                                                                                                                                                                                                                                                                                                                                            |
+| `image_registration::filter`           | Gaussian smoothing and FFT helpers: `gaussian_smooth`, `gaussian_kernel`, `gaussian`, `fft`, `ifft`, `fft_freq`, `fft_shift`.                                                                                                                                                                                                                                                                                                                                                                                                                                                                            |
 | `image_registration::optimize`         | Optimizers: `lbfgs_minimize`, `asgd_minimize`, `AsgdConfig`.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 |
 | `image_registration::par_indexed_iter` | Parallel iterators over array indices.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       |
 | `image_registration::error`            | `Error` type.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                |

src/bspline.rs

@@ -7,8 +7,13 @@ use ndarray::{
 };
 use num::traits::FloatConst;
 use num::{Complex, cast::AsPrimitive, complex::ComplexFloat};
+use rand::RngExt;
+use rand::distr::StandardUniform;
 use rayon::iter::ParallelIterator;
+use serde::{Deserialize, Serialize};
+use std::hash::Hash;
 use std::ops::{Deref, MulAssign};
+use std::sync::{Arc, Mutex};
 
 fn in_shape(shape_f: &[f64], v: &[f64]) -> bool {
     shape_f
@@ -17,6 +22,123 @@ fn in_shape(shape_f: &[f64], v: &[f64]) -> bool {
         .all(|(s, x)| (-0.5 <= *x) && (x <= s))
 }
 
+// TODO: sample nearby points
+#[derive(Debug, Serialize, Deserialize)]
+pub struct DefaultValue {
+    #[serde(skip, default = "DefaultValue::make_rng")]
+    rng: Arc<Mutex<rand::prelude::StdRng>>,
+    values: Vec<f64>,
+    weights: Vec<f64>,
+}
+
+impl Clone for DefaultValue {
+    fn clone(&self) -> Self {
+        Self {
+            rng: Self::make_rng(),
+            values: self.values.clone(),
+            weights: self.weights.clone(),
+        }
+    }
+}
+
+impl Default for DefaultValue {
+    fn default() -> Self {
+        Self {
+            rng: Self::make_rng(),
+            values: vec![0.0],
+            weights: vec![1.0],
+        }
+    }
+}
+
+impl PartialEq for DefaultValue {
+    fn eq(&self, other: &Self) -> bool {
+        (self.values == other.values) && (self.weights == other.weights)
+    }
+}
+
+impl DefaultValue {
+    pub fn new(mut values: Vec<f64>, mut weights: Vec<f64>) -> Self {
+        if values.is_empty() {
+            values.push(f64::NAN);
+        }
+        if weights.len() != values.len() {
+            weights.resize(values.len(), 1.0);
+        }
+        weights.iter_mut().for_each(|w| {
+            if !w.is_finite() {
+                *w = 1.0;
+            } else if *w < 0.0 {
+                *w = 0.0
+            }
+        });
+        let n = weights.iter().sum::<f64>();
+        weights.iter_mut().fold(0.0, |mut acc, w| {
+            acc += *w / n;
+            *w = acc;
+            acc
+        });
+        Self {
+            rng: Self::make_rng(),
+            values,
+            weights,
+        }
+    }
+
+    pub fn values(&self) -> &[f64] {
+        &self.values
+    }
+
+    pub fn weights(&self) -> &[f64] {
+        &self.weights
+    }
+
+    pub fn sample(&self) -> f64 {
+        if self.values.len() == 1 {
+            self.values[0]
+        } else {
+            let i = self.rng.lock().unwrap().sample::<f64, _>(StandardUniform);
+            let n = self
+                .weights
+                .iter()
+                .position(|&w| i <= w)
+                .unwrap_or(self.weights.len() - 1);
+            self.values[n]
+        }
+    }
+
+    fn make_rng() -> Arc<Mutex<rand::prelude::StdRng>> {
+        Arc::new(Mutex::new(rand::make_rng()))
+    }
+}
+
+impl Eq for DefaultValue {}
+
+impl Hash for DefaultValue {
+    fn hash<H: std::hash::Hasher>(&self, state: &mut H) {
+        self.values.iter().for_each(|v| v.to_bits().hash(state));
+        self.weights.iter().for_each(|v| v.to_bits().hash(state));
+    }
+}
+
+impl From<f64> for DefaultValue {
+    fn from(value: f64) -> Self {
+        Self::new(vec![value], vec![1.0])
+    }
+}
+
+impl From<Vec<f64>> for DefaultValue {
+    fn from(value: Vec<f64>) -> Self {
+        Self::new(value, Vec::new())
+    }
+}
+
+impl From<&[f64]> for DefaultValue {
+    fn from(value: &[f64]) -> Self {
+        Self::new(value.to_vec(), Vec::new())
+    }
+}
+
 #[derive(Clone, Debug)]
 pub struct BSplineMem {
     ndim: usize,
@@ -74,44 +196,58 @@ where
     D: Dimension,
     BSpline<N, D>: BSplineTrait<D>,
 {
-    pub fn interpolate(&self, transform: &Transform<D>) -> Result<Array<f64, D>, Error> {
+    pub fn interpolate(
+        &self,
+        transform: &Transform<D>,
+        default_value: &DefaultValue,
+    ) -> 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)?;
+            } else {
+                *x = default_value.sample();
             }
         }
 
         Ok(new)
     }
 
-    pub fn interpolate_par(&self, transform: &Transform<D>) -> Result<Array<f64, D>, Error> {
+    pub fn interpolate_par(
+        &self,
+        transform: &Transform<D>,
+        default_value: &DefaultValue,
+    ) -> 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())
                         .all(|(s, x)| (-0.5 <= *x) && (x <= s))
                     {
                         *x = self.evaluate_at_continuous_index(&v, &mut mem)?;
+                    } else {
+                        *x = default_value.sample();
                     }
                     Ok(mem)
                 },

src/error.rs

@@ -19,4 +19,6 @@ pub enum Error {
     InitialGuessLength(usize, usize),
     #[error("interpolation order {0} is not implemented")]
     InterpolationNotImplemented(usize),
+    #[error("matrix is singular")]
+    SingularMatrix,
 }

src/filter.rs

@@ -1,3 +1,6 @@
+//! Module containing gaussian smooth functions as well as helper functions to calculate ND fast
+//! fourier transforms
+
 use crate::error::Error;
 use itertools::Itertools;
 use ndarray::{Array, Array1, ArrayD, AsArray, Axis, Dimension, RemoveAxis, concatenate, s};
@@ -113,7 +116,6 @@ mod tests {
     use super::*;
     use crate::julia_image;
     use ndarray::{Ix1, array};
-    use tiffwrite::IJTiffFile;
 
     #[test]
     fn smooth() -> Result<(), Box<dyn std::error::Error>> {
@@ -135,10 +137,10 @@ mod tests {
             &[29.5, 39.5],
             &[-0.8, 0.156],
         );
-        let im_b = gaussian_smooth(im_a.mapv(|i| i as f64).view(), &[4.0, 4.0])?;
+        let _im_b = gaussian_smooth(im_a.mapv(|i| i as f64).view(), &[4.0, 4.0])?;
-        let mut t = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/julia.tif"))?;
+        // let mut t = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/julia.tif"))?;
-        t.save(im_a.mapv(|i| i as u32).view(), 0, 0, 0)?;
+        // t.save(im_a.mapv(|i| i as u32).view(), 0, 0, 0)?;
-        t.save(im_b.mapv(|i| i as u32).view(), 1, 0, 0)?;
+        // t.save(im_b.mapv(|i| i as u32).view(), 1, 0, 0)?;
         Ok(())
     }
 

src/lib.rs

@@ -1,3 +1,5 @@
+#![cfg_attr(docsrs, feature(doc_cfg))]
+#![doc = include_str!("../README.md")]
 use rayon::iter::ParallelIterator;
 pub mod bspline;
 pub mod error;
@@ -48,16 +50,11 @@ mod tests {
     use super::*;
     use crate::transform::Transform;
     use tempfile::NamedTempFile;
-    use tiffwrite::IJTiffFile;
 
     #[test]
     fn serialization() -> Result<(), Error> {
         let file = NamedTempFile::new()?;
-        let t = Transform::<Ix2>::new_with_center(
+        let t = Transform::<Ix2>::new(vec![1.2, 0.3, -0.4, 0.9, 10.2, -9.5]);
-            vec![1.2, 0.3, -0.4, 0.9, 10.2, -9.5],
-            vec![59.5, 49.5],
-            vec![120, 100],
-        );
         t.to_file(file.path().to_path_buf())?;
         let s = Transform::from_file(file.path().to_path_buf())?;
         assert_eq!(s, t);
@@ -66,12 +63,8 @@ mod tests {
 
     #[test]
     fn transform_point() -> Result<(), Box<dyn std::error::Error>> {
-        let transform = Transform::<Ix2>::new_with_center(
+        let transform = Transform::<Ix2>::new(vec![1.0, 0.0, 0.0, 1.0, -120.0, 10.0]);
-            vec![1.0, 0.0, 0.0, 1.0, -120.0, 10.0],
+        let p = transform.transform_point(&[300, 400], &[299.5, 399.5]);
-            vec![299.5, 399.5],
-            vec![600, 800],
-        );
-        let p = transform.transform_point(&[300, 400]);
         assert_eq!(p, [180.0, 410.0]);
         Ok(())
     }
@@ -87,13 +80,12 @@ mod tests {
             .collect::<Vec<_>>()
             .try_into()
             .unwrap();
-        let transform =
+        let transform = Transform::<Ix2>::new(transform_j.to_vec()).with_rotation(1.0, &[0.0, 0.0]);
-            Transform::<Ix2>::new(transform_j.to_vec(), shape.to_vec()).with_rotation(1.0);
         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);
         let k = julia_image(&shape, &transform_k, &center, &c).mapv(|x| x as f64);
-        let n = transform.interpolate::<3, _, _>(k.view())?;
+        let n = transform.interpolate::<3, _, _>(k.view(), &Default::default())?;
         let sj = j.iter().sum::<f64>();
         // let sk = k.iter().sum::<f64>();
         let sn = n.iter().sum::<f64>();
@@ -128,13 +120,12 @@ mod tests {
             .collect::<Vec<_>>()
             .try_into()
             .unwrap();
-        let transform =
+        let transform = Transform::<Ix2>::new(transform_j.to_vec()).with_rotation(1.0, &[0.0, 0.0]);
-            Transform::<Ix2>::new(transform_j.to_vec(), shape.to_vec()).with_rotation(1.0);
         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);
         let k = julia_image(&shape, &transform_k, &center, &c).mapv(|x| x as f64);
-        let n = transform.interpolate_par::<3, _, _>(k.view())?;
+        let n = transform.interpolate_par::<3, _, _>(k.view(), &Default::default())?;
         let sj = j.iter().sum::<f64>();
         let sn = n.iter().sum::<f64>();
 
@@ -158,22 +149,22 @@ mod tests {
             &[3999.5, 2999.5],
             &[-0.8, 0.156],
         );
-        let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
+        let transform = Transform::from_rotation(1.0, &[0.0, 0.0]);
-        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(())
     }
 
     #[test]
     fn julia() -> Result<(), Box<dyn std::error::Error>> {
-        let image_a = julia_image(
+        let _image_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],
         );
-        let mut t = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/julia.tif"))?;
+        // let mut t = IJTiffFile::new(std::env::home_dir().unwrap().join("tmp/julia.tif"))?;
-        t.save(image_a.view(), 0, 0, 0)?;
+        // t.save(image_a.view(), 0, 0, 0)?;
         Ok(())
     }
 }

src/metric.rs

@@ -9,8 +9,7 @@ use itertools::Itertools;
 use ndarray::{Array, Array1, Array2, Array3, AsArray, Axis, Dimension, s};
 use num::cast::AsPrimitive;
 use rand::RngExt;
-use rayon::iter::IntoParallelIterator;
+use rayon::iter::{IntoParallelIterator, ParallelIterator};
-use rayon::iter::ParallelIterator;
 use serde::{Deserialize, Serialize};
 use std::cell::RefCell;
 use std::ops::Deref;
@@ -65,6 +64,22 @@ impl FixedMu {
             .collect()
     }
 
+    /// combine the dparameters of the variable parameters into a complete
+    /// parameter set, with no uncertainty (0) on the fixed parameters
+    pub fn combine_dparameters(&self, dparameters: &[f64]) -> Vec<f64> {
+        debug_assert_eq!(dparameters.len(), self.len_variable);
+        let mut dparameters_iter = dparameters.iter();
+        self.iter()
+            .map(|i| {
+                if i.is_some() {
+                    0.0
+                } else {
+                    *dparameters_iter.next().unwrap()
+                }
+            })
+            .collect()
+    }
+
     /// extract mu from complete parameter set
     pub fn extract_variable(&self, mu: &[f64]) -> Vec<f64> {
         self.iter()
@@ -204,6 +219,24 @@ fn image_jacobian(point: &[f64], center: &[f64], image_derivative: &[f64]) -> Ve
         .collect()
 }
 
+/// Local Gauss-Newton statistics of the intensity fit
+/// `fixed(x) = moving(T(x)) + noise` at a parameter point.
+///
+/// `hessian` is the weighted normal matrix `JᵀWJ` and `variance` is the
+/// weighted residual variance `Σ w r² / (Σ w − q − 1)`, with `q` the number of
+/// variable parameters. The extra degree of freedom accounts for the weighted
+/// mean intensity offset that is removed from the residuals. The parameter
+/// covariance of the fit is `variance * hessian⁻¹`.
+#[derive(Clone, Debug)]
+pub struct ResidualFit {
+    /// weighted normal matrix `JᵀWJ`
+    pub hessian: Vec<Vec<f64>>,
+    /// residual variance (mean squared residual, corrected for the dof)
+    pub variance: f64,
+    /// number of sample points that contributed
+    pub n_samples: usize,
+}
+
 pub struct MattesMetric<D>
 where
     D: Dimension,
@@ -440,6 +473,93 @@ where
         }
         Ok(())
     }
+
+    /// Compute the Gauss-Newton normal matrix and the residual variance of the
+    /// intensity fit at `point`, using the same sample points and edge
+    /// weighting as [`Self::eval`].
+    ///
+    /// Residuals are `fixed(x) - moving(T(x))`, so a well-registered image pair
+    /// gives a small variance. Returns `None` when interpolation fails or when
+    /// too few samples contribute to estimate the residual variance.
+    pub fn residual_fit(&self, point: &[f64]) -> Option<ResidualFit> {
+        let complete_mu = self.fixed_mu.combine(point);
+        let idx = self.sampling.index();
+        let shape = self.shape.as_slice();
+        let shape_f = self.shape.iter().map(|s| s - 0.5).collect::<Vec<_>>();
+        let center = self.center.as_slice();
+        let e = 1.0 / self.edge;
+        let n_params = point.len();
+        let mut hessian = vec![vec![0.0; n_params]; n_params];
+        let mut sum_w = 0.0;
+        let mut sum_wr = 0.0;
+        let mut sum_wr2 = 0.0;
+        let mut n_samples = 0;
+        let mut mem_f = self.fixed.get_mem();
+        let mut mem_m = self.moving.get_mem();
+        for u in idx {
+            let v = transform_point(&u, center, &complete_mu);
+            if !shape_f
+                .iter()
+                .zip_eq(v.iter())
+                .all(|(s, x)| (-0.5 <= *x) && (x <= s))
+            {
+                continue;
+            }
+            let f = self
+                .fixed
+                .evaluate_at_continuous_index(&u, &mut mem_f)
+                .ok()?;
+            let (m, d) = self
+                .moving
+                .evaluate_value_and_derivative_at_continuous_index(&v, &mut mem_m)
+                .ok()?;
+            // same smooth edge weighting as in `eval`
+            let mut w = 1.0;
+            for (s, x) in shape.iter().zip_eq(&v) {
+                if e * (x + 0.5) < *s {
+                    let k = 3.0 * e * (x + 0.5) / s - 1.5;
+                    w *= square_bspline_integral(k);
+                } else if e * (x + 0.5) > (e - 1.0) * s {
+                    let k = 3.0 * e * (s - x - 0.5) / s - 1.5;
+                    w *= square_bspline_integral(k);
+                }
+            }
+            if w <= 0.0 || !w.is_finite() {
+                continue;
+            }
+            let r = f - m;
+            sum_w += w;
+            sum_wr += w * r;
+            sum_wr2 += w * r * r;
+            n_samples += 1;
+            let j = self
+                .fixed_mu
+                .extract_variable(&image_jacobian(&u, center, &d));
+            for (a, ja) in j.iter().enumerate() {
+                for (b, jb) in j.iter().enumerate() {
+                    hessian[a][b] += w * ja * jb;
+                }
+            }
+        }
+        // the residual variance needs at least one degree of freedom
+        if n_samples == 0 || sum_w <= n_params as f64 + 1.0 {
+            return None;
+        }
+        // subtract the weighted mean intensity mismatch (e.g. a constant
+        // offset between the images) before estimating the variance
+        let mean_r = sum_wr / sum_w;
+        let sse = (sum_wr2 - sum_w * mean_r * mean_r).max(0.0);
+        // q parameters plus the mean intensity offset that was removed
+        let variance = sse / (sum_w - n_params as f64 - 1.0);
+        if !variance.is_finite() {
+            return None;
+        }
+        Some(ResidualFit {
+            hessian,
+            variance,
+            n_samples,
+        })
+    }
 }
 
 impl<D> ObjectiveFunction<f64> for MattesMetric<D>
@@ -455,6 +575,39 @@ where
         self.eval(point).unwrap();
         Some(self.metric.borrow().derivative.clone())
     }
+
+    /// The hessian is computed numerically: central finite differences of
+    /// `gradient` with respect to each parameter (2 * n full evaluations).
+    /// `MattesMetric` only accumulates first-order terms, so an analytic hessian
+    /// is not available. The result is symmetrized, which also damps the small
+    /// asymmetry introduced by finite differencing.
+    fn hessian(&self, point: &[f64]) -> Option<Vec<Vec<f64>>> {
+        let n = point.len();
+        let mut hessian = vec![vec![0.0; n]; n];
+        let mut p = point.to_vec();
+        for j in 0..n {
+            let h = 1e-5 * f64::max(1.0, point[j].abs());
+            p[j] += h;
+            self.eval(&p).unwrap();
+            let g_plus = self.metric.borrow().derivative.clone();
+            p[j] -= 2.0 * h;
+            self.eval(&p).unwrap();
+            let g_minus = self.metric.borrow().derivative.clone();
+            p[j] += h;
+            for i in 0..n {
+                hessian[i][j] = (g_plus[i] - g_minus[i]) / (2.0 * h);
+            }
+        }
+        #[expect(clippy::needless_range_loop)]
+        for i in 0..n {
+            for j in (i + 1)..n {
+                let w = (hessian[i][j] + hessian[j][i]) / 2.0;
+                hessian[i][j] = w;
+                hessian[j][i] = w;
+            }
+        }
+        Some(hessian)
+    }
 }
 
 #[derive(Clone, Debug, Serialize, Deserialize, PartialEq)]
@@ -588,6 +741,122 @@ mod tests {
         Ok(())
     }
 
+    #[test]
+    fn residual_fit_1d() -> Result<(), Box<dyn std::error::Error>> {
+        let a = array![1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 4.0, 2.0, 0.0, 1.0];
+        // b is a shifted right by 3 (padded on the left with the first value)
+        let b = array![1.0, 1.0, 1.0, 1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 4.0];
+        let points = (0..11).map(|i| vec![i as f64]).collect::<Vec<_>>();
+        let m = MattesMetric::new_from_arrays(
+            a.view(),
+            b.view(),
+            SamplingArg::FixedAt(points),
+            16,
+            0.05,
+        )?
+        .with_fixed_mu(vec![Some(1.0), None]);
+
+        // fixed(x) = moving(x + 3): no residual
+        let fit = m.residual_fit(&[3.0]).expect("residual fit");
+        assert!(fit.n_samples > 0);
+        assert!(fit.hessian[0][0] > 0.0);
+        assert!(fit.variance < 1e-12);
+
+        // away from the optimum the residual variance must grow
+        let worse = m.residual_fit(&[0.0]).expect("residual fit");
+        assert!(worse.variance > fit.variance);
+        Ok(())
+    }
+
+    #[test]
+    fn hessian_1d() -> 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(), &[2.0])?;
+        let b = gaussian_smooth(im_b.view(), &[2.0])?;
+
+        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),
+            5,
+            0.05,
+        )?
+        .with_fixed_mu(vec![Some(1.0), None]);
+
+        let p = -3.49;
+        let h = m.hessian(&[p]).unwrap();
+        assert_eq!(h.len(), 1);
+        assert_eq!(h[0].len(), 1);
+
+        // cross-check against the second central difference of `evaluate`
+        let step = 1e-4;
+        let fd = (m.evaluate(&[p + step]) - 2.0 * m.evaluate(&[p]) + m.evaluate(&[p - step]))
+            / (step * step);
+        println!("hessian = {:?}\nfinite diff = {:?}", h[0][0], fd);
+        assert!((h[0][0] - fd).abs() < 1e-3 + 1e-2 * fd.abs());
+        Ok(())
+    }
+
+    #[test]
+    fn hessian_2d() -> Result<(), Box<dyn std::error::Error>> {
+        let a = Array2::from_shape_fn((12, 12), |(i, j)| i as f64 + 0.7 * j as f64);
+        let b = Array2::from_shape_fn((12, 12), |(i, j)| 0.5 * i as f64 + 24.0 - 0.7 * j as f64);
+        let a = gaussian_smooth(a.view(), &[2.0, 2.0])?;
+        let b = gaussian_smooth(b.view(), &[2.0, 2.0])?;
+
+        let points = (0..12)
+            .flat_map(|i| (0..12).map(move |j| vec![i as f64, j as f64]))
+            .step_by(2)
+            .collect::<Vec<_>>();
+
+        let m = MattesMetric::new_from_arrays(
+            a.view(),
+            b.view(),
+            SamplingArg::FixedAt(points),
+            7,
+            0.05,
+        )?;
+        let n = m.fixed_mu().len_variable();
+        assert_eq!(n, 2 * 2 + 2);
+
+        let p = vec![0.5, -0.3, 0.2, 0.7, -0.4, 0.3];
+        let h = m.hessian(&p).unwrap();
+        assert_eq!(h.len(), n);
+        assert!(h.iter().all(|row| row.len() == n));
+        // symmetrized by construction
+        for (i, row) in h.iter().enumerate() {
+            for (j, v) in row.iter().enumerate() {
+                assert_eq!(*v, h[j][i]);
+            }
+        }
+
+        // cross-check each diagonal entry against the 2nd central difference of `evaluate`
+        let step = 1e-4;
+        for j in 0..n {
+            let mut pp = p.clone();
+            let mut pm = p.clone();
+            pp[j] += step;
+            pm[j] -= step;
+            let fd = (m.evaluate(&pp) - 2.0 * m.evaluate(&p) + m.evaluate(&pm)) / (step * step);
+            println!("diagonal[{j}] = {:?}\nfinite diff   = {:?}", h[j][j], fd);
+            assert!(
+                (h[j][j] - fd).abs() < 1e-3 + 1e-2 * fd.abs(),
+                "diagonal entry {j}"
+            );
+        }
+        Ok(())
+    }
+
     #[test]
     fn test_square_bspline_integral() -> Result<(), Box<dyn std::error::Error>> {
         let x = Array1::linspace(-2.0, 2.0, 500);

src/register.rs

@@ -1,10 +1,10 @@
 use crate::bspline::{BSpline, BSplineTrait};
 use crate::error::Error;
-use crate::metric::{FixedMu, MattesMetric, SamplingArg, Sigma};
+use crate::metric::{FixedMu, MattesMetric, ResidualFit, SamplingArg, Sigma};
 use crate::optimize::{AsgdConfig, asgd_minimize, lbfgs_minimize};
-use crate::transform::Transform;
+use crate::transform::{Transform, try_matrix_inverse};
 use algos::{ObjectiveFunction, OptimizationConfig};
-use ndarray::{AsArray, Dimension};
+use ndarray::{Array2, AsArray, Dimension};
 use num::cast::AsPrimitive;
 use serde::{Deserialize, Serialize};
 use std::marker::PhantomData;
@@ -175,7 +175,7 @@ impl<D: Dimension> Registration<D> {
         } else {
             return Err(Error::NumberOfDimensionsNotDefined);
         };
-        let mut fixed_mu = vec![Some(0.0); ndim * ndim];
+        let mut fixed_mu = vec![Some(0.0); ndim * ndim + ndim];
         for i in 0..ndim {
             fixed_mu[i * (ndim + 1)] = Some(1.0);
             fixed_mu[ndim * ndim + i] = None;
@@ -286,8 +286,11 @@ impl<D: Dimension> Registration<D> {
             .unwrap_or_else(|| RegistrationStep::default_steps(ndim, n));
 
         let mut p_var = self.initial_guess.as_ref().cloned().unwrap_or_else(|| {
-            Self::compute_initial_transform(fixed.shape(), moving.shape(), ndim)
+            let full = Self::compute_initial_transform(fixed.shape(), moving.shape(), ndim);
+            self.fixed_mu.extract_variable(&full)
         });
+        let mut last_metric: Option<MattesMetric<D>> = None;
+        let mut last_result: Option<algos::OptimizationResult<f64>> = None;
 
         for RegistrationStep {
             sigma,
@@ -324,16 +327,23 @@ impl<D: Dimension> Registration<D> {
                 &scales,
             );
 
-            if let Some(result) = optimization_result
+            last_metric = Some(metric);
+            last_result = optimization_result;
+            if let Some(result) = last_result.as_ref()
                 && result.optimal_point.iter().all(|i| i.is_finite())
             {
-                p_var = result.optimal_point;
+                p_var = result.optimal_point.clone();
             }
         }
-        Ok(Transform::<D>::new(
+        let mut transform = Transform::<D>::new(self.fixed_mu.combine(p_var.as_slice()));
-            self.fixed_mu.combine(p_var.as_slice()),
+        if let (Some(metric), Some(result)) = (&last_metric, last_result.as_ref())
-            fixed.shape().to_vec(),
+            && result.converged
-        ))
+            && let Some(fit) = metric.residual_fit(&p_var)
+        {
+            let dparameters = dparameters_from_residual_fit(&fit);
+            transform = transform.with_dparameters(self.fixed_mu.combine_dparameters(&dparameters));
+        }
+        Ok(transform)
     }
 
     fn optimize_metric<M: ObjectiveFunction<f64>>(
@@ -392,9 +402,12 @@ impl<D: Dimension> Registration<D> {
             .unwrap_or_else(|| RegistrationStep::default_steps(ndim, n));
 
         let mut p_var = self.initial_guess.as_ref().cloned().unwrap_or_else(|| {
-            Self::compute_initial_transform(fixed.shape(), moving.shape(), ndim)
+            let full = Self::compute_initial_transform(fixed.shape(), moving.shape(), ndim);
+            self.fixed_mu.extract_variable(&full)
         });
-        let mut p = vec![0.0; n * n + n];
+        let mut p = vec![0.0; ndim * ndim + ndim];
+        let mut last_metric: Option<MattesMetric<D>> = None;
+        let mut last_result: Option<algos::OptimizationResult<f64>> = None;
 
         let mut registration_results = Vec::new();
         for RegistrationStep {
@@ -464,6 +477,8 @@ impl<D: Dimension> Registration<D> {
                     .fixed_mu()
                     .combine(&optimization_result.optimal_point);
             }
+            last_metric = Some(metric);
+            last_result = Some(optimization_result.clone());
 
             registration_results.push(RegistrationResult {
                 sigma_fixed: sigma.sigma(f_orig_shape.as_slice()),
@@ -480,9 +495,265 @@ impl<D: Dimension> Registration<D> {
                 converged: optimization_result.converged,
             });
         }
-        Ok((
+        let mut transform = Transform::<D>::new(p);
-            Transform::<D>::new(p, fixed.shape().to_vec()),
+        if let (Some(metric), Some(result)) = (&last_metric, last_result.as_ref())
-            registration_results,
+            && result.converged
-        ))
+            && let Some(fit) = metric.residual_fit(&p_var)
+        {
+            let dparameters = dparameters_from_residual_fit(&fit);
+            transform = transform.with_dparameters(self.fixed_mu.combine_dparameters(&dparameters));
+        }
+        Ok((transform, registration_results))
+    }
+}
+
+/// Diagonal of the inverse of a symmetric positive-definite matrix.
+///
+/// Non-finite entries are treated as zero. Returns `None` if the matrix is
+/// singular, because a singular normal matrix means the corresponding
+/// directions are not determined by the data. Diagonal entries of the inverse
+/// that are not finite and positive are reported as 0.
+fn variance_from_hessian(hessian: &[Vec<f64>]) -> Option<Vec<f64>> {
+    let n = hessian.len();
+    if hessian.iter().any(|row| row.len() != n) {
+        return None;
+    }
+    let mut data: Vec<f64> = Vec::with_capacity(n * n);
+    data.extend(
+        hessian
+            .iter()
+            .flatten()
+            .map(|v| if v.is_finite() { *v } else { 0.0 }),
+    );
+    let matrix = Array2::from_shape_vec((n, n), data).ok()?;
+    let inv = try_matrix_inverse(matrix)?;
+    Some(
+        (0..n)
+            .map(|i| {
+                let v = inv[[i, i]];
+                if v.is_finite() && v > 0.0 { v } else { 0.0 }
+            })
+            .collect(),
+    )
+}
+
+/// Compute parameter uncertainty (standard error, the `dparameters` field) from
+/// the intensity fit at the optimum.
+///
+/// The moving image is modelled as `moving(T(x)) = fixed(x) + noise` with
+/// independent residuals. Around the optimum this is a linear least-squares
+/// problem with Jacobian `J`, so the parameter covariance is
+///
+/// ```text
+/// Cov = σ² (JᵀWJ)⁻¹
+/// ```
+///
+/// with `σ²` the residual variance and `JᵀWJ` the weighted Gauss-Newton normal
+/// matrix, both computed by [`MattesMetric::residual_fit`]. This is the usual
+/// standard error of a nonlinear least-squares fit: it is calibrated in
+/// parameter units and grows when the images fit poorly or when the image
+/// gradients are small.
+///
+/// This deliberately does *not* use the inverse Hessian of `-MI`: mutual
+/// information is not a likelihood, so its curvature has an arbitrary scale
+/// (it changes with the number of histogram bins and the Parzen width) and
+/// `H⁻¹` is not a standard error.
+///
+/// Each uncertainty is the square root of the corresponding diagonal. If the
+/// normal matrix is singular or the residuals cannot be estimated, the result
+/// is 0 for the affected parameters, meaning "not estimated".
+pub(crate) fn dparameters_from_residual_fit(fit: &ResidualFit) -> Vec<f64> {
+    let n = fit.hessian.len();
+    if !fit.variance.is_finite() || fit.variance <= 0.0 {
+        return vec![0.0; n];
+    }
+    variance_from_hessian(&fit.hessian)
+        .map(|variance| {
+            variance
+                .into_iter()
+                .map(|v| {
+                    let v = fit.variance * v;
+                    if v.is_finite() && v > 0.0 {
+                        v.sqrt()
+                    } else {
+                        0.0
+                    }
+                })
+                .collect()
+        })
+        .unwrap_or_else(|| vec![0.0; n])
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use ndarray::{Array1, Ix1};
+
+    #[test]
+    fn dparameters_from_residual_fit_1d() {
+        // Cov = variance / (JᵀWJ), so with JᵀWJ = 4 and variance = 1/4 the
+        // standard error is sqrt(1/16) = 1/4.
+        let fit = ResidualFit {
+            hessian: vec![vec![4.0]],
+            variance: 0.25,
+            n_samples: 10,
+        };
+        let d = dparameters_from_residual_fit(&fit);
+        assert!((d[0] - 0.25).abs() < 1e-12);
+
+        // the covariance scales linearly with the residual variance
+        let fit = ResidualFit {
+            variance: 1.0,
+            ..fit
+        };
+        let d = dparameters_from_residual_fit(&fit);
+        assert!((d[0] - 0.5).abs() < 1e-12);
+    }
+
+    #[test]
+    fn dparameters_from_residual_fit_2d() {
+        // inverse of [[2, 1], [1, 2]] has diagonal 2/3
+        let fit = ResidualFit {
+            hessian: vec![vec![2.0, 1.0], vec![1.0, 2.0]],
+            variance: 0.5,
+            n_samples: 20,
+        };
+        let d = dparameters_from_residual_fit(&fit);
+        let expected: f64 = (1.0f64 / 3.0).sqrt();
+        assert!((d[0] - expected).abs() < 1e-9);
+        assert!((d[1] - expected).abs() < 1e-9);
+    }
+
+    #[test]
+    fn dparameters_from_residual_fit_not_estimated() {
+        // a singular normal matrix leaves the parameters undetermined; there
+        // is no diagonal fallback
+        let singular = ResidualFit {
+            hessian: vec![vec![2.0, -2.0], vec![-2.0, 2.0]],
+            variance: 1.0,
+            n_samples: 20,
+        };
+        assert_eq!(dparameters_from_residual_fit(&singular), vec![0.0, 0.0]);
+
+        // a non-positive-definite or non-finite normal matrix is not
+        // informative either
+        for hessian in [vec![vec![0.0]], vec![vec![-4.0]], vec![vec![f64::NAN]]] {
+            let fit = ResidualFit {
+                hessian,
+                variance: 1.0,
+                n_samples: 10,
+            };
+            assert_eq!(dparameters_from_residual_fit(&fit), vec![0.0]);
+        }
+
+        // no residual variance means no standard error is available
+        for variance in [0.0, -1.0, f64::NAN, f64::INFINITY] {
+            let fit = ResidualFit {
+                hessian: vec![vec![4.0]],
+                variance,
+                n_samples: 10,
+            };
+            assert_eq!(dparameters_from_residual_fit(&fit), vec![0.0]);
+        }
+    }
+
+    #[test]
+    fn combine_dparameters() {
+        let fixed_mu: FixedMu = vec![Some(1.0), None, None].into();
+        assert_eq!(
+            fixed_mu.combine_dparameters(&[0.1, 0.2]),
+            vec![0.0, 0.1, 0.2]
+        );
+    }
+
+    #[test]
+    fn register_deterministic_dparameters() -> Result<(), Box<dyn std::error::Error>> {
+        // dense, non-periodic 1D signal; b is a, shifted right by 2: the
+        // optimal translation is +2
+        // incommensurate frequencies keep the signal effectively aperiodic,
+        // so the mutual information peaks at the true shift
+        let n = 256;
+        let a: Array1<f64> = (0..n)
+            .map(|i| {
+                let x = i as f64;
+                1.5 + (0.131 * x).sin()
+                    + 0.7 * (0.397 * x + 1.3).sin()
+                    + 0.3 * (0.883 * x + 0.7).sin()
+            })
+            .collect();
+        let b: Array1<f64> = (0..n)
+            .map(|i| if i < 2 { 0.0 } else { a[[i - 2]] })
+            .collect();
+        let points: Vec<Vec<f64>> = (0..n).map(|i| vec![i as f64]).collect();
+        let steps = vec![
+            RegistrationStep::new(
+                Sigma::Absolute(vec![2.0]),
+                SamplingArg::FixedAt(points.clone()),
+                24,
+                1e-6,
+                0.05,
+                1024,
+                1.0,
+                Optimizer::LBFGS,
+            ),
+            RegistrationStep::new(
+                Sigma::None,
+                SamplingArg::FixedAt(points),
+                48,
+                1e-6,
+                0.05,
+                2048,
+                1.0,
+                Optimizer::LBFGS,
+            ),
+        ];
+        let registration = Registration::<Ix1>::new_translation()?.with_steps(steps);
+        let transform = registration.register(a.view(), b.view())?;
+        assert_eq!(transform.parameters.len(), 2);
+        assert_eq!(transform.dparameters.len(), 2);
+        // the fixed scale parameter has no uncertainty
+        assert_eq!(transform.dparameters[0], 0.0);
+        // the data are a perfect shifted copy, so the residual variance is
+        // (numerically) zero and so is the translation uncertainty
+        assert!(transform.dparameters[1].is_finite());
+        assert!(transform.dparameters[1] < 0.05);
+        // the found translation should be close to the applied one
+        assert!((transform.parameters[1] - 2.0).abs() < 1.0);
+        Ok(())
+    }
+
+    #[test]
+    fn noisy_fit_has_positive_dparameters() -> Result<(), Box<dyn std::error::Error>> {
+        let n = 128;
+        let a: Array1<f64> = (0..n)
+            .map(|i| {
+                let x = i as f64;
+                1.5 + (0.131 * x).sin()
+                    + 0.7 * (0.397 * x + 1.3).sin()
+                    + 0.3 * (0.883 * x + 0.7).sin()
+            })
+            .collect();
+        let b: Array1<f64> = (0..n)
+            .map(|i| {
+                let base = if i < 2 { 0.0 } else { a[[i - 2]] };
+                // deterministic broadband perturbation
+                base + 0.05 * ((i * 37 % 19) as f64 / 18.0 - 0.5)
+            })
+            .collect();
+        let points: Vec<Vec<f64>> = (0..n).map(|i| vec![i as f64]).collect();
+        let metric = MattesMetric::new_from_arrays(
+            a.view(),
+            b.view(),
+            SamplingArg::FixedAt(points),
+            16,
+            0.05,
+        )?
+        .with_fixed_mu(vec![Some(1.0), None]);
+
+        let fit = metric.residual_fit(&[2.0]).expect("residual fit");
+        assert!(fit.variance > 0.0);
+        let d = dparameters_from_residual_fit(&fit);
+        assert!(d[0] > 0.0);
+        Ok(())
     }
 }

src/transform.rs

@@ -1,10 +1,14 @@
-use crate::bspline::{BSpline, BSplineTrait};
+use crate::bspline::{BSpline, BSplineTrait, DefaultValue};
 use crate::error::Error;
+use crate::filter::{fft, ifft};
 use crate::metric::FixedMu;
 use crate::register::{Registration, RegistrationResult, RegistrationStep};
 use itertools::Itertools;
-use ndarray::{Array, Array2, ArrayD, ArrayView2, AsArray, Dimension, Ix2, IxDyn, s};
+use ndarray::{
+    Array, Array2, ArrayD, ArrayView2, AsArray, Dimension, IntoDimension, Ix2, IxDyn, s,
+};
 use num::cast::AsPrimitive;
+use num::complex::ComplexFloat;
 use serde::{Deserialize, Serialize};
 use serde_yaml::{from_reader, to_writer};
 use std::fs::File;
@@ -44,11 +48,18 @@ where
         .collect()
 }
 
-/// Calculate the inverse of a square matrix using LU factorization with partial pivoting.
+/// Calculate the inverse of a square matrix using LU factorization with partial
-fn matrix_inverse(matrix: Array2<f64>) -> Array2<f64> {
+/// pivoting.
+///
+/// Returns `None` if the matrix is singular, or numerically rank-deficient.
+pub(crate) fn try_matrix_inverse(matrix: Array2<f64>) -> Option<Array2<f64>> {
     let n = matrix.nrows();
     debug_assert_eq!(n, matrix.ncols(), "matrix must be square");
 
+    // pivot magnitudes must exceed this relative threshold to count as significant
+    let max_abs = matrix.iter().fold(0.0f64, |m, &x| m.max(x.abs()));
+    let threshold = max_abs * (n as f64) * f64::EPSILON;
+
     // PA = LU factorization (Doolittle with partial pivoting)
     let mut lu = matrix;
     let mut perm: Vec<usize> = (0..n).collect();
@@ -60,7 +71,9 @@ fn matrix_inverse(matrix: Array2<f64>) -> Array2<f64> {
                 pivot_row = i;
             }
         }
-        assert_ne!(lu[[pivot_row, k]], 0.0, "matrix is singular");
+        if lu[[pivot_row, k]].abs() <= threshold || !lu[[pivot_row, k]].is_finite() {
+            return None;
+        }
         if pivot_row != k {
             for j in 0..n {
                 lu.swap([k, j], [pivot_row, j]);
@@ -98,42 +111,17 @@ fn matrix_inverse(matrix: Array2<f64>) -> Array2<f64> {
         }
     }
 
-    inv
+    Some(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 = IxDyn> {
-pub struct Transform<D: Dimension = Ix2> {
     /// 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>,
 }
@@ -178,14 +166,12 @@ impl<D: Dimension> Mul<&Transform<D>> for &Transform<D> {
                 .cloned()
                 .collect(),
             dparameters: dm
-                .slice(s![..self.ndim(), ..])
+                .slice(s![..self.ndim(), ..self.ndim()])
                 .flatten()
                 .iter()
-                .chain(dm.slice(s![..self.ndim(), ..]).iter())
+                .chain(dm.slice(s![..self.ndim(), self.ndim()]).iter())
                 .cloned()
                 .collect(),
-            center: self.center.clone(),
-            shape: self.shape.clone(),
             dimension: self.dimension,
         }
     }
@@ -236,8 +222,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 +272,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 +322,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 +356,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 +365,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 +387,35 @@ impl<D: Dimension> Default for Transform<D> {
         Self {
             parameters,
             dparameters: vec![0f64; ndim * ndim + ndim],
-            center: vec![0f64; ndim],
-            shape: vec![0usize; ndim],
             dimension: PhantomData,
         }
     }
 }
 
 impl<D: Dimension> Transform<D> {
-    /// parameters: flat NxN part of matrix + translation; center: center of rotation
+    /// parameters: flat NxN part of matrix
-    pub fn new_with_center(parameters: Vec<f64>, center: Vec<f64>, shape: Vec<usize>) -> Self {
+    pub fn new(parameters: Vec<f64>) -> Self {
-        let ndim = if let Some(ndim) = D::NDIM {
+        let dparameters = vec![0.0; parameters.len()];
-            ndim
-        } else {
-            center.len()
-        };
-        debug_assert_eq!(parameters.len(), ndim * ndim + ndim);
-        debug_assert_eq!(center.len(), ndim);
-        debug_assert_eq!(shape.len(), ndim);
         Self {
             parameters,
-            dparameters: vec![0.0; ndim * ndim + ndim],
+            dparameters,
-            center,
-            shape,
             dimension: PhantomData,
         }
     }
 
-    pub fn new(parameters: Vec<f64>, shape: Vec<usize>) -> Self {
+    pub fn new_dparameters(parameters: Vec<f64>, dparameters: Vec<f64>) -> 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();
         Self {
             parameters,
-            dparameters: vec![0.0; ndim * ndim + ndim],
+            dparameters,
-            center,
-            shape,
             dimension: PhantomData,
         }
     }
 
+    pub fn new_around_center(parameters: Vec<f64>, center: &[f64]) -> Self {
+        Self::new(parameters).with_translation(center)
+            * Self::from_translation(&center.iter().map(|i| -i).collect::<Vec<_>>())
+    }
+
     /// 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
@@ -472,7 +433,46 @@ impl<D: Dimension> Transform<D> {
         M: AsArray<'a, T, D>,
         T: 'a + Clone + AsPrimitive<f64>,
     {
-        Registration::new_translation()?.register(fixed, moving)
+        let fixed = fixed.into();
+        let moving = moving.into();
+        let initial_guess = Self::register_affine(&fixed, &moving)?;
+        let initial_guess =
+            initial_guess.parameters[initial_guess.ndim() * initial_guess.ndim()..].to_vec();
+        Registration::new_translation()?
+            .with_initial_guess(initial_guess)?
+            .register(fixed, moving)
+    }
+
+    /// find the translation which transforms moving into fixed
+    /// using the cross-correlation theorem and fft
+    pub fn register_translation_fft<'a, F, M, T>(fixed: F, moving: M) -> Result<Self, Error>
+    where
+        F: AsArray<'a, T, D>,
+        M: AsArray<'a, T, D>,
+        T: 'a + Clone + AsPrimitive<f64>,
+    {
+        let f = fft(fixed.into().mapv(|i| i.as_()).view())?.mapv(|i| i.conj());
+        let m = fft(moving.into().mapv(|i| i.as_()).view())?;
+        let c = ifft(&(f * m))?.mapv(|i| i.re);
+        let (i, _) = c
+            .indexed_iter()
+            .max_by(|(_, a), (_, b)| a.total_cmp(b))
+            .expect("no max because array is empty");
+        let d = i.into_dimension();
+        let j = d
+            .as_array_view()
+            .iter()
+            .zip(c.shape())
+            .map(|(&i, &s)| {
+                let mut fi = i as f64;
+                let fs = s as f64;
+                if 2.0 * fi > fs {
+                    fi -= fs;
+                }
+                fi
+            })
+            .collect::<Vec<_>>();
+        Ok(Self::from_translation(&j))
     }
 
     /// find the transform which transforms moving into fixed, using fixed_mu to specify which
@@ -548,17 +548,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 +573,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 +598,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 +635,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 +644,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 +675,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 +703,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<_>>();
@@ -819,7 +805,7 @@ impl<D: Dimension> Transform<D> {
     /// get the inverse transform
     pub fn inverse(&self) -> Result<Self, Error> {
         let matrix = self.matrix();
-        let inverse = matrix_inverse(matrix);
+        let inverse = try_matrix_inverse(matrix).ok_or(Error::SingularMatrix)?;
         let parameters = inverse
             .slice(s![..self.ndim(), ..self.ndim()])
             .iter()
@@ -830,39 +816,14 @@ 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 interpolate<'a, const B: usize, T, A>(
-    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>,
-        D: Dimension,
-        T: 'a + Clone + AsPrimitive<f64>,
-        BSpline<B, D>: BSplineTrait<D>,
-    {
-        let image = image.into();
-        let bspline = BSpline::<B, _>::new(&image);
-        bspline.interpolate(self)
-    }
-
-    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>,
@@ -872,13 +833,30 @@ impl<D: Dimension> Transform<D> {
     {
         let image = image.into();
         let bspline = BSpline::<B, _>::new(&image);
-        bspline.interpolate_par(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>,
+        D: Dimension,
+        T: 'a + Clone + AsPrimitive<f64>,
+        BSpline<B, D>: BSplineTrait<D>,
+    {
+        let image = image.into();
+        let bspline = BSpline::<B, _>::new(&image);
+        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>,
@@ -886,12 +864,12 @@ impl<D: Dimension> Transform<D> {
         T: 'a + Clone + AsPrimitive<f64>,
     {
         match order {
-            0 => Ok(self.interpolate::<0, _, _>(image)?),
+            0 => Ok(self.interpolate::<0, _, _>(image, default_value)?),
-            1 => Ok(self.interpolate::<1, _, _>(image)?),
+            1 => Ok(self.interpolate::<1, _, _>(image, default_value)?),
-            2 => Ok(self.interpolate::<2, _, _>(image)?),
+            2 => Ok(self.interpolate::<2, _, _>(image, default_value)?),
-            3 => Ok(self.interpolate::<3, _, _>(image)?),
+            3 => Ok(self.interpolate::<3, _, _>(image, default_value)?),
-            4 => Ok(self.interpolate::<4, _, _>(image)?),
+            4 => Ok(self.interpolate::<4, _, _>(image, default_value)?),
-            5 => Ok(self.interpolate::<5, _, _>(image)?),
+            5 => Ok(self.interpolate::<5, _, _>(image, default_value)?),
             _ => Err(Error::InterpolationNotImplemented(order)),
         }
     }
@@ -900,6 +878,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>,
@@ -907,12 +886,12 @@ impl<D: Dimension> Transform<D> {
         T: 'a + Clone + AsPrimitive<f64>,
     {
         match order {
-            0 => Ok(self.interpolate_par::<0, _, _>(image)?),
+            0 => Ok(self.interpolate_par::<0, _, _>(image, default_value)?),
-            1 => Ok(self.interpolate_par::<1, _, _>(image)?),
+            1 => Ok(self.interpolate_par::<1, _, _>(image, default_value)?),
-            2 => Ok(self.interpolate_par::<2, _, _>(image)?),
+            2 => Ok(self.interpolate_par::<2, _, _>(image, default_value)?),
-            3 => Ok(self.interpolate_par::<3, _, _>(image)?),
+            3 => Ok(self.interpolate_par::<3, _, _>(image, default_value)?),
-            4 => Ok(self.interpolate_par::<4, _, _>(image)?),
+            4 => Ok(self.interpolate_par::<4, _, _>(image, default_value)?),
-            5 => Ok(self.interpolate_par::<5, _, _>(image)?),
+            5 => Ok(self.interpolate_par::<5, _, _>(image, default_value)?),
             _ => Err(Error::InterpolationNotImplemented(order)),
         }
     }
@@ -922,7 +901,7 @@ impl<D: Dimension> Transform<D> {
         if let Some(d) = D::NDIM {
             d
         } else {
-            self.center.len()
+            self.parameters.len().isqrt()
         }
     }
 
@@ -932,9 +911,7 @@ impl<D: Dimension> Transform<D> {
                 Ok(Transform {
                     parameters: self.parameters,
                     dparameters: self.dparameters,
-                    center: self.center,
+                    dimension: PhantomData,
-                    shape: self.shape,
-                    dimension: Default::default(),
                 })
             } else {
                 Err(Error::DimensionalityMismatch(d, self.ndim()))
@@ -943,9 +920,7 @@ impl<D: Dimension> Transform<D> {
             Ok(Transform {
                 parameters: self.parameters,
                 dparameters: self.dparameters,
-                center: self.center,
+                dimension: PhantomData,
-                shape: self.shape,
-                dimension: Default::default(),
             })
         }
     }
@@ -954,9 +929,7 @@ impl<D: Dimension> Transform<D> {
         Transform {
             parameters: self.parameters,
             dparameters: self.dparameters,
-            center: self.center,
+            dimension: PhantomData,
-            shape: self.shape,
-            dimension: Default::default(),
         }
     }
 
@@ -968,7 +941,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,26 +957,37 @@ 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(),
         }
     }
+
+    pub fn swap_axes(&self, a: usize, b: usize) -> Self {
+        let mut new = self.clone();
+        let n = self.ndim();
+        let mut m = (0..n).collect::<Vec<_>>();
+        (m[a], m[b]) = (m[b], m[a]);
+        for (i, j) in (0..n).zip(m) {
+            if i != j {
+                new.parameters[i + n * j] = self.parameters[j + n * i];
+                new.dparameters[i + n * j] = self.dparameters[j + n * i];
+                new.parameters[j + n * i] = self.parameters[i + n * j];
+                new.dparameters[j + n * i] = self.dparameters[i + n * j];
+            }
+        }
+        new.parameters[n + a] = self.parameters[n + b];
+        new.dparameters[n + a] = self.dparameters[n + b];
+        new
+    }
 }
 
 impl Transform<Ix2> {
     /// create a transform from a rotation in radians
     pub fn from_rotation(theta: f64, center: &[f64]) -> Self {
-        Self {
+        Self::new_around_center(
-            parameters: vec![
+            vec![
                 theta.cos(),
                 -theta.sin(),
                 theta.sin(),
@@ -1011,27 +995,19 @@ impl Transform<Ix2> {
                 0.0,
                 0.0,
             ],
-            dparameters: vec![0f64; 6],
+            center,
-            center: center.to_vec(),
+        )
-            shape: vec![0usize; 2],
-            dimension: PhantomData,
-        }
     }
 
     /// multiply self with a rotation matrix
-    pub fn with_rotation(&self, theta: f64) -> Self {
+    pub fn with_rotation(&self, theta: f64, center: &[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 {
         Self::from_rotation(theta, center) * self
     }
 }
 
 #[cfg(test)]
 mod tests {
-    use super::matrix_inverse;
+    use super::try_matrix_inverse;
     use crate::julia_image;
     use crate::transform::File;
     use crate::transform::Transform;
@@ -1048,7 +1024,7 @@ mod tests {
 
         let a = array![[4.0, 7.0], [2.0, 6.0]];
         let expected = array![[0.6, -0.7], [-0.2, 0.4]];
-        let got = matrix_inverse(a.clone());
+        let got = try_matrix_inverse(a.clone()).expect("matrix should not be singular");
         assert!(
             got.iter()
                 .zip_eq(expected.iter())
@@ -1057,7 +1033,7 @@ mod tests {
 
         let b = array![[0.0, 2.0], [3.0, 4.0]];
         let expected = array![[-2.0 / 3.0, 1.0 / 3.0], [0.5, 0.0]];
-        let got = matrix_inverse(b.clone());
+        let got = try_matrix_inverse(b.clone()).expect("matrix should not be singular");
         assert!(
             got.iter()
                 .zip_eq(expected.iter())
@@ -1066,13 +1042,16 @@ mod tests {
 
         let c = array![[1.0, 2.0, 3.0], [0.0, 1.0, 4.0], [5.0, 6.0, 0.0]];
         let reference = array![[-24.0, 18.0, 5.0], [20.0, -15.0, -4.0], [-5.0, 4.0, 1.0]];
-        let got = matrix_inverse(c);
+        let got = try_matrix_inverse(c).expect("matrix should not be singular");
         assert!(
             got.iter()
                 .zip_eq(reference.iter())
                 .all(|(x, y)| (x - y).abs() < 1e-10)
         );
 
+        let singular = array![[1.0, 2.0], [2.0, 4.0]];
+        assert!(try_matrix_inverse(singular).is_none());
+
         Ok(())
     }
 
@@ -1084,9 +1063,9 @@ mod tests {
             &[3999.5, 2999.5],
             &[-0.8, 0.156],
         );
-        let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
+        let transform = Transform::from_rotation(1.0, &[0.0, 0.0]);
-        let image_b = transform.interpolate_par::<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())?;
+        let image_c = transform.interpolate::<3, _, _>(image_a.view(), &Default::default())?;
 
         assert!(
             image_b
@@ -1105,8 +1084,8 @@ mod tests {
             &[3999.5, 2999.5],
             &[-0.8, 0.156],
         );
-        let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
+        let transform = Transform::from_rotation(1.0, &[0.0, 0.0]);
-        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(())
     }
@@ -1119,8 +1098,8 @@ mod tests {
             &[3999.5, 2999.5],
             &[-0.8, 0.156],
         );
-        let transform = Transform::from_rotation(1.0, &[3999.5, 2999.5]);
+        let transform = Transform::from_rotation(1.0, &[0.0, 0.0]);
-        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(())
     }
@@ -1139,8 +1118,7 @@ mod tests {
         .slice(s![.., 0])
         .mapv(|i| i as f64);
         let q = vec![0.85, 4.0];
-        let im_b =
+        let im_b = Transform::new(q).interpolate::<1, _, _>(&im_a, &Default::default())?;
-            Transform::new(q.clone(), vec![im_a.shape()[0]]).interpolate::<1, _, _>(&im_a)?;
 
         // Test with ALL points (like grad_check uses)
         let metric_all = MattesMetric::<ndarray::Ix1>::new(
@@ -1225,7 +1203,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)?;
+            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]);
@@ -1335,10 +1313,10 @@ mod tests {
 
         let q = vec![0.85, 2.0];
         let moving =
-            Transform::new(q.clone(), vec![fixed.shape()[0]]).interpolate::<3, _, _>(&fixed)?;
+            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(), 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.
@@ -1347,9 +1325,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 t_point = t.transform_point(&[point], &center);
-            let e_point = q_inv.transform_point(&[point]);
+            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;
@@ -1386,8 +1369,7 @@ mod tests {
         .slice(s![.., 0])
         .mapv(|i| i as f64);
         let q = vec![0.85, 2.0];
-        let im_b =
+        let im_b = Transform::new(q.clone()).interpolate::<1, _, _>(&im_a, &Default::default())?;
-            Transform::new(q.clone(), vec![im_a.shape()[0]]).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())?;
@@ -1401,7 +1383,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;
 
@@ -1506,9 +1488,8 @@ 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 =
+        let transform = Transform::<Ix2>::new_around_center(params.clone(), &center);
-            Transform::<Ix2>::new_with_center(params.clone(), center.to_vec(), shape.to_vec());
+        let moving = transform.interpolate::<3, _, _>(fixed.view(), &Default::default())?;
-        let moving = transform.interpolate::<3, _, _>(fixed.view())?;
 
         let q_inv = transform.inverse()?.parameters;
 
@@ -1540,28 +1521,26 @@ mod tests {
         use crate::metric::{MattesMetric, SamplingArg};
         use algos::ObjectiveFunction;
 
-        let shape = [200, 200];
-        let center = [99.5, 99.5];
         let im_a = julia_image(
-            &shape,
+            &[200, 200],
             &[1.0, 0.0, 0.0, 1.0, 0.0, 0.0],
-            &center,
+            &[99.5, 99.5],
             &[-0.8, 0.156],
         )
         .mapv(|i| i as f64);
-        let rotation = Transform::<Ix2>::from_rotation(f64::PI() / 4.0, &center);
+        let rotation = Transform::<Ix2>::from_rotation(f64::PI() / 4.0, &[0.0, 0.0]);
-        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 identity = vec![1.0, 0.0, 0.0, 1.0, 0.0, 0.0];
         let edge = 0.01;
 
-        // Test at full resolution with different sigma levels (matching elastix approach)
+        // Test at full resolution with different sigma levels (matching Elastix approach)
         // Elastix does NOT downsample — it smooths at full resolution with
         // sigma = 0.5 * factor where factor is from the pyramid schedule [8, 4, 2, 1]
         println!(
-            "=== Full resolution (no downsampling, matching elastix FixedSmoothingImagePyramid) ==="
+            "=== Full resolution (no downsampling, matching Elastix FixedSmoothingImagePyramid) ==="
         );
         for (level, sigma_val) in [4.0, 2.0, 1.0, 0.5].iter().enumerate() {
             let f = gaussian_smooth(im_a.view(), &[*sigma_val; 2])?;
@@ -1619,20 +1598,18 @@ mod tests {
         use crate::filter::gaussian_smooth;
         use crate::metric::{MattesMetric, SamplingArg};
         use algos::ObjectiveFunction;
-        let shape = [200, 200];
-        let center = [99.5, 99.5];
         let im_a = julia_image(
-            &shape,
+            &[200, 200],
             &[1.0, 0.0, 0.0, 1.0, 0.0, 0.0],
-            &center,
+            &[99.5, 99.5],
             &[-0.8, 0.156],
         )
         .mapv(|i| i as f64);
         let edge = 0.01;
         // 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 rotation = Transform::<Ix2>::from_rotation(angle_deg.to_radians(), &[0.0, 0.0]);
-            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])?;
@@ -1723,17 +1700,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 t_point = t.transform_point(&point, &center);
-                let e_point = expected_transform.transform_point(&point);
+                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;
@@ -1765,9 +1743,8 @@ 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())
+        let e = Transform::<Ix2>::new(vec![0.8, 0.0, 0.0, 1.0, 0.0, 0.0]).inverse()?;
-            .inverse()?;
+        let moving = e.interpolate::<3, _, _>(fixed.view(), &Default::default())?;
-        let moving = e.interpolate::<3, _, _>(fixed.view())?;
 
         let t = Transform::<Ix2>::register(
             fixed.view(),
@@ -1781,13 +1758,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 t_point = t.transform_point(&point, &center);
-                let e_point = e_inv.transform_point(&point);
+                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;