374 lines
21 KiB
Markdown
374 lines
21 KiB
Markdown
# 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]([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<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:
|
||
|
||
```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::<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.
|
||
|
||
```rust
|
||
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):
|
||
|
||
```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::<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
|
||
|
||
```rust
|
||
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.
|
||
|
||
```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<Option<f64>>` 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::<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).
|
||
|
||
```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<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):
|
||
|
||
```rust
|
||
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](LICENSE-APACHE) or [MIT](LICENSE-MIT), at your option.
|