From 687c37360d0b10bc5dc51792af080c0fb75dd572 Mon Sep 17 00:00:00 2001 From: "w.pomp" Date: Mon, 3 Aug 2026 13:59:09 +0200 Subject: [PATCH] - readme --- README.md | 374 +++++++++++++++++++++++++++++++++++++++++++++++++- src/metric.rs | 4 - 2 files changed, 372 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index 7b8be61..ff74007 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,373 @@ -# Image registration +# image-registration + Find an ND affine transformation that registers an ND test image to a fixed image. Then apply this transformation to -other images. \ No newline at end of file +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 +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. + +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). + +## Features + +- **Affine registration** in any number of dimensions — rotation, scale, shear, and translation parameters + (`Transform::register`). +- **B-spline interpolation** of order up to 5 (linear, cubic, …) for resampling/warping images under a transform + (`interpolate`, `interpolate_par`). +- **Mattes mutual-information metric** with B-spline joint-histogram Parzen estimation, analytic gradients, and + configurable sampling (fixed, fixed-at-points, or random) and binning. +- **Multi-resolution registration** with a fixed-smoothing image pyramid with Gaussian smoothing. +- **Two optimizers**: L-BFGS (default, for smooth/consistent gradients) and ASGD (adaptive stochastic gradient descent, + robust to noisy gradients). +- **Flexible parameter control**: fix any subset of transform parameters during registration (`FixedMu`). +- **Transform algebra**: composition, inverse, matrix forms, coordinate transforms, serialization to/from YAML. +- **Parallel** interpolation and metric evaluation via rayon. +- Pure Rust — no BLAS/LAPACK or system libraries required (matrix inversion is implemented internally with LU + factorization). + +## Quick start + +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; + +fn register(fixed: &Array2, moving: &Array2) -> Result { + // 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::::register( + fixed.view(), + moving.view(), + vec![None; 6], + None, + None, + ) +} +``` + +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 +let transform = register(&fixed, &moving)?; +let registered = transform.interpolate_par::<3, _, _>(moving.view())?; // cubic B-spline, parallel +``` + +## Concepts + +### Parameterization of an affine transform + +An ND affine transform is stored as a flattened `N×N` linear part plus `N` translation parameters, for `N*N + N` values +in total: + +- 1D: `[scale, translation]` +- 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```: + +```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` +lets you override it. + +### Registration convention + +The registration minimizes the negative Mattes mutual information between the fixed image and the warped moving image. +For a moving image that was created by applying transform `q` to a fixed image, registration recovers `q⁻¹` (the +transform that maps fixed coordinates back onto the moving image). + +## Examples + +### 1. Build transforms and transform coordinates + +```rust +use image_registration::error::Error; +use image_registration::transform::Transform; +use ndarray::Ix2; + +fn main() -> Result<(), Error> { + // Rotation of 45° about the image center of a 200×200 image + let rot = Transform::::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]); + + // Transform a single point (the center is subtracted and re-added) + let p = t.transform_point(&[200, 300]); + + // 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())?; + + // 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)); + Ok(()) +} +``` + +### 2. Warp an image (resampling) + +Interpolate an image under a transform, producing an output image of the same shape. Pixels whose transformed +coordinates fall outside the input image are set to zero. + +```rust +use image_registration::error::Error; +use image_registration::transform::Transform; + +fn main() -> Result<(), Error> { + // `image` is an Array2 (200×150, say) + let image = ndarray::Array2::::zeros((200, 150)); + let transform = Transform::::from_rotation(0.3, &[99.5, 74.5]); + + // Order-1 (linear) interpolation, single-threaded + let warped_linear = transform.interpolate::<1, _, _>(image.view())?; + + // Order-3 (cubic) B-spline interpolation, parallel + let warped_cubic = transform.interpolate_par::<3, _, _>(image.view())?; + Ok(()) +} +``` + +### 3. Full affine registration with a known ground truth + +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 +use image_registration::error::Error; +use image_registration::julia_image; +use image_registration::transform::Transform; +use ndarray::Ix2; + +fn main() -> Result<(), Error> { + // Fixed image (a Julia fractal, included for testing) + let shape = [200, 200]; + let center = [99.5, 99.5]; + let fixed = julia_image(&shape, &[1.0, 0.0, 0.0, 1.0, 0.0, 0.0], ¢er, &[-0.8, 0.156]) + .mapv(|i| i as f64); + + // Ground-truth affine: 10° rotation, 5% scale change, (12, -8) px translation + 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 = + Transform::::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())?; + + // Register (all 6 parameters free, default steps) + let t = Transform::::register(fixed.view(), moving.view(), vec![None; 6], None, None)?; + + // The registered transform should equal the inverse of the ground truth + let q_inv = ground_truth.inverse()?.parameters; + let sse: f64 = t + .parameters + .iter() + .zip(q_inv.iter()) + .map(|(a, b)| (a - b).powi(2)) + .sum(); + assert!(sse < 1.0); + Ok(()) +} +``` + +### 4. Registration with inspection of every optimization step + +```rust +use image_registration::error::Error; +use image_registration::transform::Transform; + +fn main() -> Result<(), Error> { + let (t, steps) = Transform::::register_debug( + fixed.view(), moving.view(), vec![None; 6], None, None, + )?; + + for step in steps { + println!( + "sigma={:?} n_bins={} optimizer_iters={} converged={} -> {:?}", + step.sigma_fixed, step.n_bins, step.iterations, step.converged, step.optimal_point + ); + } + Ok(()) +} +``` + +### 5. Custom registration steps (multi-resolution pyramid) + +`RegistrationStep` gives you full control over the smoothing, sampling, binning, and optimizer at each level. The +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 +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( + Sigma::Absolute(vec![4.0, 4.0]), // Gaussian sigma per dimension + SamplingArg::Fixed(3000), // 3000 sampled points (cached) + 32, // histogram bins + 1e-4, // tolerance + 0.05, // edge fraction + 512, // max iterations + 1.0, // learning rate + ); + coarse.optimizer = Optimizer::ASGD; // robust to noisy gradients + + let mut fine = RegistrationStep::new( + Sigma::Absolute(vec![0.5, 0.5]), + SamplingArg::Fixed(3000), + 32, + 1e-8, + 0.05, + 2048, + 1.0, + ); // optimizer stays LBFGS (the default) + + let mut reg = Registration::new(vec![None; 6]); + reg.set_steps(vec![coarse, fine]); + let t = reg.register(fixed.view(), moving.view())?; + Ok(()) +} +``` + +> **Tip:** `RegistrationStep::default_steps(ndim, n)` returns a sensible default pyramid: a 5-level schedule that starts +> with ASGD on heavily smoothed images and finishes with L-BFGS on the unsmoothed image. + +### 6. Fix (pin) some parameters during registration + +Pass a `Vec>` as `fixed_mu`: `None` optimizes a parameter, `Some(v)` keeps it fixed at `v`. For example, +register for translation only: + +```rust +use image_registration::error::Error; +use image_registration::transform::Transform; + +fn main() -> Result<(), Error> { + // 2D: fix the linear part to identity, optimize tx and ty + let t = Transform::::register( + fixed.view(), + moving.view(), + vec![Some(1.0), Some(0.0), Some(0.0), Some(1.0), None, None], + None, + None, + )?; + Ok(()) +} +``` + +For the common "translation only" case there is also the convenience method +`Transform::register_translation(fixed, moving)`. + +### 7. Evaluate the Mattes metric directly + +Build a `MattesMetric` and evaluate the (negative) mutual information and its gradient at arbitrary parameters. +`SamplingArg` controls which points are used: + +- `Random(n)` — `n` random continuous positions (drawn from a thread-local RNG, cached per metric). +- `Fixed(n)` — `n` random positions, frozen at construction. +- `FixedAt(points)` — an explicit list of sample positions (fully deterministic; recommended for reproducible results). + +```rust +use algos::ObjectiveFunction; +use image_registration::bspline::BSpline; +use image_registration::error::Error; +use image_registration::metric::{MattesMetric, SamplingArg, Sigma}; + +fn main() -> Result<(), Error> { + // `fixed` and `moving` are Array2 + let fixed = ndarray::Array2::::zeros((64, 64)); + let moving = ndarray::Array2::::zeros((64, 64)); + + let f = Sigma::Absolute(vec![2.0, 2.0]).smooth(fixed.view())?; + let m = Sigma::Absolute(vec![2.0, 2.0]).smooth(moving.view())?; + + let metric = MattesMetric::new( + BSpline::<0, _>::new(f.view()), // fixed image B-spline (nearest-neighbor) + BSpline::<3, _>::new(m.view()), // moving image B-spline (cubic) + SamplingArg::FixedAt((0..100).map(|i| vec![i as f64, i as f64]).collect()), + 32, + 0.05, + )? + .with_fixed_mu(vec![None; 6]); // which parameters are variable + + let id = [1.0, 0.0, 0.0, 1.0, 0.0, 0.0]; + let val = metric.evaluate(&id); + let grad = metric.gradient(&id).unwrap(); + println!("mi({:?}) = {:.6}, grad = {:?}", id, val, grad); + Ok(()) +} +``` + +### 8. Save and load transforms + +Transforms serialize to YAML (via serde): + +```rust +use image_registration::error::Error; +use image_registration::transform::Transform; +use std::path::PathBuf; + +fn main() -> Result<(), Error> { + let t = Transform::::from_rotation(0.2, &[99.5, 99.5]); + t.to_file(PathBuf::from("transform.yaml"))?; + let back = Transform::::from_file(PathBuf::from("transform.yaml"))?; + assert_eq!(t, back); + Ok(()) +} +``` + +## Module reference + +| Module | Contents | +|----------------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------| +| `image_registration::transform` | `Transform`, the affine transform type: construction (`new`, `new_with_center`, `from_translation`, `from_scaling`, `from_rotation`), algebra (`inverse`, `with_translation`, `with_scaling`, `with_rotation`, `with_rotation_around`, `Mul`), coordinate transforms (`transform_point`, `transform_points`, `matrix`, `dmatrix`), resampling (`interpolate`, `interpolate_par`), registration (`register`, `register_debug`, `register_affine`, `register_translation`), serialization (`to_file`, `from_file`), plus the free function `transform_point` and a blas-free `matrix_inverse`. | +| `image_registration::register` | `Registration`, `RegistrationStep`, `RegistrationResult`, `Optimizer` (LBFGS/ASGD). | +| `image_registration::metric` | `MattesMetric` (Mattes mutual information with analytic gradient), `SamplingArg` (Fixed/FixedAt/Random), `Sigma` (None/Absolute/Relative), `FixedMu`. | +| `image_registration::bspline` | `BSpline` B-spline interpolation of arbitrary order, `BSplineTrait`, `BSplineMem`, Parzen kernel helpers. | +| `image_registration::filter` | Gaussian smoothing and FFT helpers: `gaussian_smooth`, `gaussian_kernel`, `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. | +| `image_registration::julia_image` | Test-image generator (Julia fractal on an `Array2`). | + +## How it works + +1. **Image pyramid.** Each registration level smooths the fixed and moving images with a Gaussian of the configured + `Sigma` (and optionally downsamples them). Coarse levels use large sigma so the metric is smooth and far-sighted; + fine levels use small sigma for precision. +2. **Metric.** `MattesMetric` samples points in the fixed image, warps them into the moving image through the current + transform, and builds a joint intensity histogram using cubic B-spline Parzen windows. It minimizes (the negative of) + the mutual information with an analytic gradient computed via the chain rule (image Jacobian × transform Jacobian). +3. **Optimization.** ASGD is used at coarse levels (robust to noisy gradients, with automatic parameter estimation), + L-BFGS at fine levels (quadratic convergence on smooth objectives). +4. **Result.** The final transform maps fixed-image coordinates into moving-image coordinates; its inverse maps moving + coordinates back into the fixed frame. + +## Known limitations + +- With `SamplingArg::Random` / `SamplingArg::Fixed`, sample positions are drawn from a thread-local RNG, so different + runs can converge to different (local) optima — especially on small or nearly uniform images where the MI signal is + weak. For reproducible results, prefer `SamplingArg::FixedAt` with explicit sample points, and use enough samples to + cover the structure of interest. +- The metric is defined only where the warped sample falls inside the moving image; points near the edge are smoothly + weighted out of the histogram, so registration can be unreliable if the two images overlap only slightly. + +## License + +Licensed under either of [Apache-2.0](LICENSE-APACHE) or [MIT](LICENSE-MIT), at your option. diff --git a/src/metric.rs b/src/metric.rs index 0d1ef6e..45ec974 100644 --- a/src/metric.rs +++ b/src/metric.rs @@ -214,8 +214,6 @@ where /// memory for (point, value, gradient) combination, refcell gives interior mutability metric: RefCell, edge: f64, - // mem_f: RefCell, - // mem_m: RefCell, } impl MattesMetric @@ -231,8 +229,6 @@ where ) -> Result { let shape = fixed.shape().iter().map(|&i| i as f64).collect::>(); let center = shape.iter().map(|i| (i - 1.0) / 2.0).collect(); - // let mem_f = fixed.get_mem(); - // let mem_m = moving.get_mem(); let (min, max) = fixed .iter() .chain(moving.iter())