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.

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 (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, Ix2};

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::<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)?;
// 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

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 c that is supplied where the transform is used:

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_point/transform_points take the center as an argument, and resampling (interpolate, interpolate_par) 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

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

    // Matrix form (3×3 homogeneous) and exact inverse
    let m = t.matrix();
    let inv = t.inverse()?;
    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(())
}

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 (out-of-range pixels → 0.0)
    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(), &Default::default())?;
    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_around_center(params.clone(), &center);

    // Moving image = ground truth applied to the fixed image
    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)?;

    // 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 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
        Optimizer::ASGD,                 // robust to noisy gradients
    );

    let fine = RegistrationStep::new(
        Sigma::Absolute(vec![0.5, 0.5]),
        SamplingArg::Fixed(3000),
        32,
        1e-8,
        0.05,
        2048,
        1.0,
        Optimizer::LBFGS, // the default optimizer
    );

    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; // direct dependency of image-registration; add it to your crate to name the trait
use image_registration::bspline::{BSpline, BSplineTrait};
use image_registration::error::Error;
use image_registration::metric::{MattesMetric, SamplingArg, Sigma};
use image_registration::transform::Transform;

fn main() -> Result<(), Error> {
    // `fixed` and `moving` are Array2<f64>; moving is a rotated + shifted copy of fixed
    let center = [31.5f64, 31.5];
    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 sigma = Sigma::Absolute(vec![2.0, 2.0]);
    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)
        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_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 (new_affine, new_translation, with_steps, with_optimizer, …), RegistrationStep, RegistrationResult, Optimizer (LBFGS/ASGD).
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, DefaultValue (fill value for out-of-range pixels), Parzen kernel helpers.
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.
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.