image-registration/README.md

image-registration

Find an ND affine transformation that registers an ND test image to a fixed image. Then apply this transformation to other images.

image-registration is a pure-Rust, dependency-light implementation of intensity-based image registration in the spirit of SimpleElastix. It finds an affine (rigid + scale + shear + translation) transform that aligns a moving image to a fixed image by maximizing the Mattes mutual-information 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 arrays, works in arbitrary dimensions (1D, 2D, 3D, …), and is parallelized with 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:

use image_registration::error::Error;
use image_registration::transform::Transform;
use ndarray::Array2;

fn register(fixed: &Array2<f64>, moving: &Array2<f64>) -> Result<Transform, 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(
        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:

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:

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

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::<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]);

    // 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.

use image_registration::error::Error;
use image_registration::transform::Transform;

fn main() -> Result<(), Error> {
    // `image` is an Array2<f64> (200×150, say)
    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
    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):

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], &center, &[-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::<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())?;

    // Register (all 6 parameters free, default steps)
    let t = Transform::<Ix2>::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

use image_registration::error::Error;
use image_registration::transform::Transform;

fn main() -> Result<(), Error> {
    let (t, steps) = Transform::<ndarray::Ix2>::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.

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<Option<f64>> as fixed_mu: None optimizes a parameter, Some(v) keeps it fixed at v. For example, register for translation only:

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::<ndarray::Ix2>::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).

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<f64>
    let fixed = ndarray::Array2::<f64>::zeros((64, 64));
    let moving = ndarray::Array2::<f64>::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):

use image_registration::error::Error;
use image_registration::transform::Transform;
use std::path::PathBuf;

fn main() -> Result<(), Error> {
    let t = Transform::<ndarray::Ix2>::from_rotation(0.2, &[99.5, 99.5]);
    t.to_file(PathBuf::from("transform.yaml"))?;
    let back = Transform::<ndarray::Ix2>::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<N, D> 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<u8>).

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 or MIT, at your option.