Medical Image Registration

Medical image registration finds the spatial correspondence between medical images: a CT and an MRI of the same head, two lung CT scans taken at inspiration and expiration, a follow-up scan and the baseline it is compared with, or a patient’s brain and a labelled atlas. The result is a transform that maps every point of one image onto the point of the other that shows the same anatomy, so that the images can be fused, compared voxel by voxel, or used to carry labels from one to the other. The task is also called alignment, coregistration, matching or, when the target is a standard template, spatial normalization; registration that bends the image locally is called deformable, nonrigid or nonlinear registration.

The general method (transform models, similarity measures, search, and feature-based versus intensity-based alignment) is described in image registration. Medical images add volumes in physical coordinates, soft-tissue and inter-patient deformation, intensities that differ between CT and MRI contrasts, the requirement that a deformation stay physically plausible, and the difficulty of knowing whether a registration is right.

Problem Definition

One image is held fixed, $F$, and the other, the moving image $M$, is resampled onto it. Registration estimates a transform $\varphi$ that maps each point $\mathbf{x}$ of the fixed image domain to a point $\varphi(\mathbf{x})$ of the moving image, so that the warped image $M \circ \varphi$ resembles $F$. Intensity-based methods write this as an energy, in the form VoxelMorph’s authors use (Balakrishnan et al., 2019):

$$\hat{\varphi} = \arg\min_{\varphi} \; \mathcal{L}_{\text{sim}}(F,\, M \circ \varphi) + \lambda\, \mathcal{L}_{\text{smooth}}(\varphi),$$

where $\mathcal{L}_{\text{sim}}$ penalizes differences in appearance, $\mathcal{L}_{\text{smooth}}$ penalizes implausible deformations, and $\lambda$ balances the two. For a deformable transform, $\varphi(\mathbf{x}) = \mathbf{x} + \mathbf{u}(\mathbf{x})$ with a displacement $\mathbf{u}$ at every voxel. The direction matters: SimpleITK’s documentation, for example, states that the estimated transform maps points from the fixed image coordinate system to the moving one, which is what resampling the moving image onto the fixed grid needs (image warping). Labels, landmarks and contours move with the same transform, labels with nearest-neighbour interpolation, as in Klein et al.’s (2009) evaluation.

Coordinates are physical. Scanners record where each voxel lies in millimetres and how the voxel grid is oriented (the DICOM image position, orientation and spacing, or the NIfTI affine; see MRI), and voxels are often anisotropic. A transform estimated in voxel indices of one scan is meaningless for another unless both grids are related through these headers.

Medical registration usually asks for more than a good image match. A deformation of anatomy should not tear or fold tissue, so many methods restrict $\varphi$ to diffeomorphisms, smooth invertible maps whose inverses are also smooth; such transforms preserve topology and maintain invertibility (Balakrishnan et al., 2019). Locally, the Jacobian matrix $J_\varphi(\mathbf{x}) = \nabla \varphi(\mathbf{x})$ describes how a small neighbourhood is stretched, sheared and rotated, and its determinant $\det J_\varphi$ is the local volume change. A map is locally invertible and orientation-preserving only where $\det J_\varphi > 0$; voxels with $\det J_\varphi \le 0$ mark folding (Dalca et al., 2019).

Input and Output

• Input: two 3D (sometimes 2D or 4D) images with their spatial headers, a transform model, often masks (the brain, the lungs, a region of interest), an initial alignment, and for some methods landmarks or segmentations of the training or input images.
• Output: a transform (a $4 \times 4$ matrix for rigid and affine models, a grid of control-point displacements, a dense displacement or velocity field), the warped moving image, warped labels, and ideally the inverse transform. Thirion’s (1998) bijective demons, for example, return the forward and an approximate inverse transform together.

Variants

• Intra-patient and inter-patient. Images of one patient (follow-up scans, two contrasts of one session, inspiration and expiration) share the same anatomy, and many differences are motion or deformation. Images of different patients differ in anatomy itself, and correspondence may be only partial.
• Monomodal and multimodal. Same-contrast pairs can be compared with squared differences or correlation; CT against MRI, or MRI against ultrasound, needs a measure that does not assume equal intensities, such as mutual information.
• Rigid, affine and deformable. A head scanned twice moves rigidly (6 parameters in 3D), and affine registration (12 parameters) absorbs global size differences; Klein et al. preceded all nonlinear registrations with linear registration by FSL’s FLIRT. Soft tissue, breathing and anatomy between patients need deformable transforms.
• Pairwise, atlas and groupwise. A scan can be registered to another scan, to an atlas or template, or a whole set can be registered jointly to build a template.
• Image-guided intervention. Pre-operative images are registered to images acquired during surgery or radiotherapy; one Learn2Reg task registers pre-operative MRI to intra-operative ultrasound to correct tissue shift during brain tumour resection (Hering et al., 2023).

Challenges

• No ground truth. For real patients the true correspondence is unknown; Klein et al. (2009) had to treat manual labels as “silver standards”.
• Missing and changed anatomy. A tumour may grow, tissue may be resected, and an organ may be outside one field of view; Learn2Reg lists missing correspondences, partial visibility and differing fields of view among the difficulties of its tasks (Hering et al., 2023).
• Large deformations of small structures. In its lung task, Hering et al. note that breathing motion can be larger than the vessels and airways that must be aligned.
• Intensity differences. Between modalities, but also within one: MRI intensities have no absolute scale and vary with protocol and coil.
• Plausibility. An energy driven mainly by image match can fold the field; physical constraints (incompressible tissue, rigid bone) are often wanted but rarely known exactly (Sotiras et al., 2013).
• Cost. A dense 3D field has three unknowns per voxel; Klein et al. estimated about 28 million degrees of freedom for SyN on their data.

Approaches

Similarity measures

Squared differences suit monomodal pairs with equal intensities. Local (normalized) cross-correlation, computed in a small window around each voxel, tolerates smooth intensity variations; SyN in Klein et al.’s study used it with a $5 \times 5 \times 5$ window, and VoxelMorph with a $9^3$ window (Balakrishnan et al., 2019). Mutual information and normalized mutual information handle different modalities: their history, from Viola and Wells and the Leuven group in 1995–1997, is told in Alignment by Maximization of Mutual Information, and their estimation in image registration. Learn2Reg entries also used normalized gradient fields and MIND descriptors (Hering et al., 2023). Several Learn2Reg methods drove registration by segmentation overlap or keypoints instead of intensities; its authors read this as a partial return from intensity-based to surface-based registration.

Regularization and physical models

Early deformable methods modelled the moving image as a physical material. Thirion (1998) reported results qualitatively similar to the deformable-grid methods of Bajcsy and Kovačič and of Christensen and colleagues, the latter modelling a viscous fluid exactly and, in its 1996 form, run on a massively parallel computer; he cites Bro-Nielsen and Gramkow for a fast convolution implementation of fluid registration and for showing that the Gaussian smoothing used in demons approximates linear elasticity. Sotiras et al. (2013) organize deformation models into those derived from physical models (elastic bodies, viscous fluids, diffusion, curvature, flows of diffeomorphisms), those derived from interpolation theory (radial basis functions, free-form deformations and others), knowledge-based models, and task-specific constraints such as topology, volume or rigidity preservation.

Demons

Thirion (1998) described matching as diffusion: the object boundaries of the fixed image act as semi-permeable membranes, and “demons” on them push the moving image, treated as a deformable grid, through. His “demons 1” variant puts a demon at every voxel with a non-zero gradient, computes a displacement like an optical flow estimate at each (the formula is in image registration), and smooths the whole displacement field with a Gaussian at each iteration, inside a multi-scale pyramid. He used $\sigma = 1$, four pyramid levels and four iterations at the finest level for all experiments; a 2D sinusoidal deformation took 30 s of CPU time, and matching two 3D brain MR images of $256 \times 256 \times 128$ voxels, with the bijective scheme described next, took about 30 minutes on a DEC Alpha workstation. Instead of enforcing a positive Jacobian, which he judged too expensive, he computed forward and backward fields together and removed half of their composition’s residual from each at every iteration, keeping them close to inverses.

Later analysis gave demons a variational reading. As Sotiras et al. (2013) summarize, Pennec and colleagues showed that Thirion’s force corresponds to a second-order gradient descent on the sum of squared differences, and that Gaussian smoothing of the field corresponds to one gradient step on a regularizer. Because the force is driven by intensity differences, the original demons are monomodal; multimodal variants alternate registration with intensity correction or replace the measure (Sotiras et al., 2013).

Free-form deformations

Free-form deformations place a coarse grid of control points over the image and interpolate their displacements with cubic B-splines, as in the tensor-product form

$$\mathbf{u}(\mathbf{x}) = \sum_{l=0}^{3} \sum_{m=0}^{3} \sum_{n=0}^{3} B_l(\mu_x)\, B_m(\mu_y)\, B_n(\mu_z)\, \mathbf{d}_{i+l,\, j+m,\, k+n},$$

where $\mathbf{d}$ are the control-point displacements, $B_l$ the cubic B-spline basis functions, and $\mu$ the position of $\mathbf{x}$ within its grid cell (Sotiras et al., 2013). Sotiras et al. credit the combination with cubic B-splines, as in Rueckert et al.’s (1999) registration of breast MR images, for the wide adoption of FFDs in medical imaging. The field is smooth and has few parameters, but topology is not guaranteed; Sotiras et al. report that Rueckert and colleagues later obtained diffeomorphic FFDs by limiting control-point displacements to 0.4 of the grid spacing. In Klein et al.’s (2009) study, the IRTK package (Rueckert and colleagues) combined an affine and a B-spline FFD transform with normalized mutual information, and was run without regularization.

Diffeomorphic registration

Flows of diffeomorphisms. Instead of a displacement, the transform is generated by a velocity field $v_t$ integrated over time, $\partial \varphi_t / \partial t = v_t(\varphi_t)$ with $\varphi_0$ the identity, and the regularizer penalizes the velocity, $\int_0^1 \lVert v_t \rVert_V^2 \, dt$, in a norm that enforces smoothness (Sotiras et al., 2013). In this Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework (for image volumes, Sotiras et al. cite Beg et al., 2005, among others), the optimal path is a geodesic and its length defines a distance between anatomies, used to study anatomical variability. Sotiras et al. note the price: integrating the velocity over time costs computation and memory, and the usual gradient descent converges slowly.

Stationary velocity fields. A cheaper option fixes the velocity over time, so that $\varphi = \exp(v)$, the flow at time one. The exponential is computed by scaling and squaring: start from $\varphi^{(1/2^T)} = \mathrm{Id} + v/2^T$ and compose the map with itself $T$ times, $\varphi^{(1/2^{t-1})} = \varphi^{(1/2^t)} \circ \varphi^{(1/2^t)}$ (Dalca et al., 2019). Ashburner’s fast diffeomorphic algorithm of 2007 (Sotiras et al., 2013), distributed as SPM’s DARTEL toolbox (Klein et al., 2009), and diffeomorphic demons take this route. Vercauteren et al. (2009) recast Thirion’s algorithm as an optimization over displacement fields, derived from that analysis a theoretical advantage for the symmetric-forces variant and confirmed in experiments that it converges faster, and made demons diffeomorphic by replacing the addition of each update with a composition with its exponential, at the cost of a few compositions per iteration; they released an implementation in the Insight Toolkit (ITK).

Symmetric normalization. Avants et al. (2008) introduced SyN, which Sotiras et al. describe as a symmetric LDDMM driven by cross-correlation. As described in Klein et al.’s (2009) supplementary material, SyN parameterizes the shortest diffeomorphic path between two images symmetrically, optimizes both halves by bidirectional gradient descent so that the result does not depend on which image is fixed, smooths the velocity field with a Gaussian, and returns forward and inverse maps. SyN is implemented in ANTs.

Atlases and atlas-based segmentation

An atlas is a labelled reference image. Registering it to a new scan and carrying its labels along gives a segmentation: Thirion (1998) already applied his 3D inter-patient field to a labelled brain and noted that the inferred labels were mostly correct only where the two subjects had no local morphological differences. Iglesias and Sabuncu (2015) describe the history: single-atlas, registration-based segmentation was limited because one atlas cannot capture anatomical variation; probabilistic atlases, with label statistics precomputed in a single atlas space, dominated early atlas-guided segmentation; and in 2003–2004 Rohlfing and colleagues, followed by Klein et al. (2005) and Heckemann et al. (2006), registered every atlas of a set to the new image and fused the propagated labels, the approach called multi-atlas segmentation. Fusion ranges from majority voting, which picks the most frequent label at each voxel and ignores intensities, to weighted voting with global or local weights from image similarity, and probabilistic fusion such as STAPLE, originally a model of manual segmentations as noisy observations of a hidden true segmentation. The price is one deformable registration per atlas (Iglesias & Sabuncu, 2015).

Deep networks now segment many structures directly (medical image segmentation), but atlas registration still supplies a common coordinate frame for population studies, and the two can be combined: for its inter-patient tasks, ConvexAdam, the overall winner of Learn2Reg, trained a segmentation network and used its output in the cost function of an optimization-based registration (Hering et al., 2023).

Learned registration

Learned methods replace the per-pair optimization with a network that predicts the field. Early networks were trained on ground-truth fields from simulation or from classical registration (Balakrishnan et al., 2019). VoxelMorph trains a convolutional network without ground-truth fields: it takes the affinely aligned fixed and moving volumes, predicts a displacement field, warps the moving volume with a differentiable spatial transformer, and minimizes the energy above, with mean squared error or local cross-correlation and a diffusion regularizer on the gradients of the displacement. Its authors call this amortization: the parameters are optimized once over a dataset of pairs, and a new pair costs one forward pass. Optionally, segmentations available only at training time add a Dice term. The authors’ comparison with ANTs SyN and NiftyReg on 3D brain MRI is given in image registration; in the same table, all methods produced non-positive Jacobian determinants at fewer than 1% of the 5.2 million brain voxels (0.14% for ANTs SyN, 0.37% for VoxelMorph with cross-correlation, 0.60% for NiftyReg). Training with the segmentation term alone, without image and smoothness terms, produced 9.9–13.1% folding voxels (their Table IV, arXiv v3; the text rounds this to 10–13%).

A diffeomorphic variant predicts a stationary velocity field and integrates it with scaling-and-squaring layers inside the network (Dalca et al., 2019). On the same brain data its authors reported an average Dice of 0.754 against 0.749 for ANTs SyN, and fewer than one folding voxel per test scan on average, against thousands for the baselines (their Table 1, arXiv v2). The two papers report the same Dice and runtimes for the three baselines but different mean folding counts per scan: 9,662 and 7,523 voxels for ANTs SyN, 41,251 and 33,838 for NiftyReg, and 19,077 and 19,715 for VoxelMorph with cross-correlation. Neither paper explains the difference.

Learn2Reg found that hybrids did well: methods that decouple learned feature extraction from optimization, or refine a network’s prediction with per-pair instance optimization (Hering et al., 2023).

Software

Learn2Reg’s authors name SyN, elastix, NiftyReg, deeds and VoxelMorph as public methods commonly used as baselines (Hering et al., 2023). As of October 2026, three open-source toolkits, all under the Apache-2.0 licence, cover classical registration:

• ANTs (Advanced Normalization Tools) is a C++ library used through the command line, with R and Python interfaces (ANTsR, ANTsPy) and examples for template construction. Its latest release is v2.6.5, January 2026.
• elastix is an ITK-based toolbox for rigid and nonrigid registration, with a modular design and a command-line interface for batch processing, and Python access through ITKElastix; its lead developers are Stefan Klein and Marius Staring. Its latest release is 5.3.1, March 2026.
• SimpleITK is a simplified layer over ITK for Python and other languages. Its ImageRegistrationMethod is a multi-resolution framework with metrics including mean squares, correlation, ANTs neighbourhood correlation and two mutual information estimators, and it implements several demons variants separately. Its latest release is v2.5.7, 30 October 2026.

VoxelMorph’s code is public in the authors’ repository, also under Apache-2.0.

Example: Multimodal Rigid and Deformable Registration

The first script builds a smoothed 2D “head” from the Shepp–Logan phantom bundled with scikit-image, at $128 \times 128$ pixels. The moving image is the same anatomy after a known rigid motion (12°, 5.5 and −3.2 pixels), passed through a non-monotonic intensity map so that the skull, bright in the fixed image, is dark in the moving one: a stand-in for a second modality, not a model of any scanner. An exhaustive grid search over angle and translation, refined with Powell’s method, is run once with mutual information (32-bin joint histogram) and once with the mean squared difference. TRE is the mean target registration error at five points inside the head. Run with NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2.

import numpy as np
from scipy import ndimage, optimize
from skimage.data import shepp_logan_phantom

rng = np.random.default_rng(0)
N, s = 128, 400 / 128                          # image size; high-res pixels per pixel
H = ndimage.gaussian_filter(shepp_logan_phantom(), s / 2)   # band-limited "anatomy"
yy, xx = np.mgrid[0:N, 0:N].astype(float)
c = (N - 1) / 2

def anatomy(Y, X):                             # sample the anatomy at pixel coordinates
    return ndimage.map_coordinates(H, [(Y + 0.5) * s - 0.5, (X + 0.5) * s - 0.5],
                                   order=1, mode="constant", cval=0.0)

def second_contrast(v):                        # a non-monotonic intensity mapping
    return np.interp(v, [0.0, 0.1, 0.2, 0.3, 0.4, 1.0],
                        [0.05, 0.5, 0.7, 0.3, 0.9, 0.15])

def rigid(p, Y=yy, X=xx):                      # T(x) = R(x - c) + c + t, p = (deg, tx, ty)
    th = np.radians(p[0])
    return (np.sin(th) * (X - c) + np.cos(th) * (Y - c) + c + p[2],
            np.cos(th) * (X - c) - np.sin(th) * (Y - c) + c + p[1])

truth = np.array([12.0, 5.5, -3.2])            # moving(T(x)) shows the anatomy at x
th = np.radians(truth[0])                      # build the moving image with T^-1
dX, dY = xx - c - truth[1], yy - c - truth[2]
fixed = anatomy(yy, xx) + 0.01 * rng.standard_normal((N, N))
moving = second_contrast(anatomy(np.cos(th) * dY - np.sin(th) * dX + c,
                                 np.cos(th) * dX + np.sin(th) * dY + c))
moving += 0.01 * rng.standard_normal((N, N))

def pairs(p, step=1):                          # fixed and resampled moving over the overlap
    Y, X = rigid(p, yy[::step, ::step], xx[::step, ::step])
    w = ndimage.map_coordinates(moving, [Y, X], order=1, mode="constant", cval=np.nan)
    ok = ~np.isnan(w)
    return fixed[::step, ::step][ok], w[ok]

def ssd(a, b):
    return np.mean((a - b) ** 2)

def neg_mi(a, b, bins=32):
    p = np.histogram2d(a, b, bins=bins)[0]
    p /= p.sum()
    pa, pb = p.sum(1, keepdims=True), p.sum(0, keepdims=True)
    nz = p > 0
    return -(p[nz] * np.log(p[nz] / (pa @ pb)[nz])).sum()

landmarks = np.array([[64, 40], [64, 100], [45, 64], [83, 64], [64, 64]], float)  # (x, y)
def tre(p):                                    # mean landmark error in pixels
    a = np.array(rigid(p, landmarks[:, 1], landmarks[:, 0]))
    b = np.array(rigid(truth, landmarks[:, 1], landmarks[:, 0]))
    return np.linalg.norm(a - b, axis=0).mean()

grid = [(a, tx, ty) for a in range(-30, 31, 2)
        for tx in range(-8, 9) for ty in range(-8, 9)]
print("        rot     tx     ty   measure at: estimate   truth  identity   TRE (px)")
for name, cost, sign in [("MI", neg_mi, -1), ("SSD", ssd, 1)]:
    f = lambda p, step=1: cost(*pairs(p, step))
    start = min(grid, key=lambda p: f(p, step=2))          # exhaustive coarse search
    p = optimize.minimize(f, start, method="Powell",       # local refinement
                          options={"xtol": 1e-3, "ftol": 1e-7}).x
    v = [sign * f(q) for q in (p, truth, [0, 0, 0])]
    print(f"{name:4s} {p[0]:7.2f} {p[1]:6.2f} {p[2]:6.2f}   {v[0]:20.4f} {v[1]:7.4f}"
          f" {v[2]:8.4f}   {tre(p):8.2f}")
print(f"truth{truth[0]:7.2f} {truth[1]:6.2f} {truth[2]:6.2f}")

Output:

        rot     tx     ty   measure at: estimate   truth  identity   TRE (px)
MI     12.15   5.45  -3.17                 1.0577  1.0520   0.3762       0.08
SSD     6.12   4.14  -2.01                 0.0982  0.1042   0.1139       2.64
truth  12.00   5.50  -3.20

The “measure” columns give MI in nats (higher is better) and the mean squared difference (lower is better) at the estimate, the true transform and the identity. Mutual information recovers the motion to a tenth of a pixel at the landmarks. The squared difference is lower at the truth than at the identity, but lower still at a pose rotated by about half the true angle: with remapped intensities its minimum is not where the anatomy aligns, so a better optimizer would not help. This is the lesson of the affine example in image registration, here with an intensity map that is not monotonic.

The second script tests a deformable step on a known deformation. A “follow-up” image is the same anatomy in the same contrast, warped by a sinusoidal displacement of amplitude 3 pixels and period 32 pixels, like Thirion’s (1998) synthetic test; the exact registration field is computed by inverting the warp with a fixed-point iteration. Registration uses Thirion’s demon force with Gaussian smoothing of the field, either adding each update, as classical demons do, or composing the field with the exponential of the update (four scaling-and-squaring steps), in the manner of diffeomorphic demons (Vercauteren et al., 2009), at two smoothing widths. Each runs on the smooth phantom and on a version with fine random “tissue texture” inside the head. The report gives the RMS intensity residual in the head, the mean field error over the head and on structure boundaries, the Dice coefficient of two ellipses carried from the moving to the fixed image, the minimum Jacobian determinant and the number of folded pixels ($\det J \le 0$).

import sys
import numpy as np
from scipy import ndimage
from skimage.data import shepp_logan_phantom

rng = np.random.default_rng(int(sys.argv[1]) if len(sys.argv) > 1 else 1)
N, s = 128, 400 / 128
phantom = shepp_logan_phantom()
flat = ndimage.gaussian_filter(phantom, s / 2)
tex = ndimage.gaussian_filter(rng.standard_normal(phantom.shape), 2 * s)
textured = flat + 0.05 * (flat > 0.05) * tex / tex.std()   # add fine "tissue texture"
yy, xx = np.mgrid[0:N, 0:N].astype(float)

def sample(img, Y, X, order=1):                # image at pixel coordinates (Y, X)
    return ndimage.map_coordinates(img, [(Y + 0.5) * s - 0.5, (X + 0.5) * s - 0.5],
                                   order=order, mode="constant", cval=0.0)

def interp(f, uy, ux):                         # f(x + u(x)), bilinear, edge-clamped
    return ndimage.map_coordinates(f, [yy + uy, xx + ux], order=1, mode="nearest")

def compose(ay, ax, by, bx):                   # displacement of (x + a) o (x + b)
    return by + interp(ay, by, bx), bx + interp(ax, by, bx)

def jac_det(uy, ux):                           # det of the Jacobian of x + u(x)
    dyy, dyx = np.gradient(uy)
    dxy, dxx = np.gradient(ux)
    return (1 + dxx) * (1 + dyy) - dxy * dyx

# Known deformation: follow-up(x) = baseline(x + w(x)), sinusoidal, period 32 pixels
A, P = 3.0, 32.0
wy, wx = A * np.sin(2 * np.pi * xx / P), A * np.sin(2 * np.pi * yy / P)
# The true registration field u satisfies x + u(x) = (x + w)^-1(x): fixed-point iteration
uy_t, ux_t = np.zeros((N, N)), np.zeros((N, N))
for _ in range(100):
    uy_t, ux_t = -interp(wy, uy_t, ux_t), -interp(wx, uy_t, ux_t)

lab_img = ndimage.median_filter(phantom, 3)
labels = np.round(sample(lab_img, yy, xx, 0), 1) == 0.3              # two ellipses
lab_moving = np.round(sample(lab_img, yy + wy, xx + wx, 0), 1) == 0.3
head = sample(flat, yy, xx) > 0.05
gy, gx = np.gradient(sample(flat, yy, xx))
boundaries = head & (np.hypot(gy, gx) > 0.02)                        # structure edges

def demons(fixed, moving, sigma, diffeomorphic, iters=300):
    gy, gx = np.gradient(fixed)
    uy, ux = np.zeros((N, N)), np.zeros((N, N))
    for _ in range(iters):
        diff = interp(moving, uy, ux) - fixed
        den = gx**2 + gy**2 + diff**2
        den[den < 1e-9] = np.inf
        dy, dx = -diff * gy / den, -diff * gx / den        # Thirion's demon force
        if diffeomorphic:                                   # u <- u o exp(d)
            dy, dx = dy / 2**4, dx / 2**4                   # scaling and squaring
            for _ in range(4):
                dy, dx = compose(dy, dx, dy, dx)
            uy, ux = compose(uy, ux, dy, dx)
        else:                                               # u <- u + d
            uy, ux = uy + dy, ux + dx
        uy, ux = ndimage.gaussian_filter(uy, sigma), ndimage.gaussian_filter(ux, sigma)
    return uy, ux

def report(name, fixed, moving, uy, ux):
    err = np.hypot(uy - uy_t, ux - ux_t)
    warped_lab = interp(lab_moving.astype(float), uy, ux) > 0.5
    dice = 2 * (warped_lab & labels).sum() / (warped_lab.sum() + labels.sum())
    J = jac_det(uy, ux)
    rms = np.sqrt(np.mean((interp(moving, uy, ux) - fixed)[head] ** 2))
    print(f"  {name:26s} {rms:6.4f}  {err[head].mean():5.2f} {err[boundaries].mean():6.2f}"
          f"  {dice:5.3f}  {J[head].min():5.2f}  {(J[head] <= 0).sum():4d}")

print("                                RMS   field error   Dice    min  fold")
print("                             resid.   head  edges         det J    px")
for name, img in [("flat phantom", flat), ("textured phantom", textured)]:
    fixed = sample(img, yy, xx) + 0.01 * rng.standard_normal((N, N))
    moving = sample(img, yy + wy, xx + wx) + 0.01 * rng.standard_normal((N, N))
    print(name)
    report("no registration", fixed, moving, 0 * uy_t, 0 * ux_t)
    report("true field", fixed, moving, uy_t, ux_t)
    for sigma in (1.0, 0.5):
        for diffeo in (False, True):
            method = ("diffeomorphic" if diffeo else "additive") + f", sigma = {sigma}"
            report(method, fixed, moving, *demons(fixed, moving, sigma, diffeo))

Output:

                                RMS   field error   Dice    min  fold
                             resid.   head  edges         det J    px
flat phantom
  no registration            0.2138   2.89   3.13  0.811   1.00     0
  true field                 0.0267   0.00   0.00  0.981   0.74     0
  additive, sigma = 1.0      0.0365   2.38   2.24  0.951   0.32     0
  diffeomorphic, sigma = 1.0 0.0383   2.39   2.27  0.937   0.48     0
  additive, sigma = 0.5      0.0359   2.60   2.30  0.971  -1.33   422
  diffeomorphic, sigma = 0.5 0.0253   2.47   2.25  0.968   0.10     0
textured phantom
  no registration            0.2168   2.89   3.13  0.811   1.00     0
  true field                 0.0269   0.00   0.00  0.981   0.74     0
  additive, sigma = 1.0      0.0342   0.83   1.17  0.963   0.47     0
  diffeomorphic, sigma = 1.0 0.0354   0.87   1.24  0.956   0.45     0
  additive, sigma = 0.5      0.0349   1.11   1.42  0.945  -3.93   164
  diffeomorphic, sigma = 0.5 0.0226   0.87   1.24  0.962   0.11     0

On the smooth phantom every method brings the intensity residual close to that of the true field and raises the Dice of the carried ellipses from 0.81 to 0.94–0.97, yet the field is still wrong by more than 2 pixels on average, even on boundaries. The images agree because the phantom’s flat regions constrain nothing and its edges constrain only the displacement across them (our reading: the aperture problem of optical flow, which the demon force inherits); many fields make the images match equally well. With texture the same algorithms recover the field to about a pixel or better in the head. Image similarity and label overlap rose almost as much in both cases, so neither told the two situations apart; only the known field did.

The weakly smoothed additive scheme ($\sigma = 0.5$) folds hundreds of pixels on the smooth phantom and over a hundred on the textured one, while matching the images and labels about as well as the others; on the smooth phantom it even has the highest Dice of the four. The compositional update with the same smoothing keeps the determinant positive everywhere, and so do both schemes with $\sigma = 1$. The script takes the noise and texture seed as an argument; with seeds 2, 3 and 4 the pattern is the same, folding appearing only for the additive scheme at $\sigma = 0.5$ (83 to 395 pixels across the two images and three seeds). Strictly, Gaussian smoothing of a composed field is not guaranteed to keep it invertible; here it did.

Datasets and Benchmarks

Brain MRI evaluations. Learn2Reg’s authors trace benchmark papers back to RIRE, which compared rigid alignment of head MR, PET and CT. Klein et al. (2009) ran what they describe as the largest evaluation of nonlinear brain registration ever conducted: 14 algorithms (AIR, ANIMAL, ART, diffeomorphic demons, FNIRT, IRTK, JRD-fluid, ROMEO, SICLE, SyN and four SPM5 methods, including DARTEL), more than 45,000 registrations between 80 manually labelled brains from four label sets (LPBA40, IBSR18, CUMC12 and MGH10), scored by 8 error measures. The relative performance of the methods changed little with subject population, labelling protocol and overlap measure; ART, SyN, IRTK and DARTEL gave the best results, and ART and SyN were the most consistent. The algorithms’ developers were co-authors, but the study cautions that the images were older and of lower quality than some developers were used to, that each brain was labelled only once, and that overlap measures are blind to misregistration inside a region and to folding.

Challenges. Hering et al. (2023) name earlier single-task challenges: EMPIRE10 (lung CT, scored with landmarks, fissure segmentations and Jacobian determinants), CuRIOUS (intra-operative ultrasound and MRI) and ANHIR (histology), and the Continuous Registration Challenge of 2018. DIR-Lab’s lung CT data, with public landmarks, are a common benchmark, which Hering et al. note makes them prone to overfitting.

Learn2Reg. Held with MICCAI 2020 and 2021, Learn2Reg combined six 3D tasks: CuRIOUS (MRI to intra-operative ultrasound, landmarks), Hippocampus MR (394 scans, two labels), Abdomen CT-CT (inter-patient, 13 organs), Abdomen MR-CT (intra-patient, with 4 organs labelled for training and 9 for testing), OASIS (416 training brain MRIs, 35 FreeSurfer structures) and Lung CT (expiration to inspiration, 100 landmarks per test case). The organizers provided pre-processed, pre-aligned data and evaluation code, received more than 65 entries from more than 20 teams, and ranked them per metric by pairwise Wilcoxon signed-rank tests, combining metric ranks by a geometric mean, in the manner of the Medical Segmentation Decathlon (Hering et al., 2023). Their findings:

• No single approach was best on all tasks. ConvexAdam, a hybrid of learned features and GPU-accelerated optimization, ranked first overall; the network LapIRN second; MEVIS, mainly optimization-based, third; and corrField, which uses no labels, fourth and best for intra-patient registration.
• Training with anatomical labels helped inter-patient registration clearly and intra-patient registration hardly at all. With ConvexAdam on Abdomen CT-CT, Dice over eight test organs was 69% with all 13 labels in training, 55% with 4 and 45% with none; labels used in training also improved nearby unlabelled organs.
• Evaluating on structures that were not used in training exposed label bias: on the MR-CT task, LapIRN, PIMed and MEVIS lost accuracy on the nine organs not labelled for training, while ConvexAdam and corrField did not.
• Optimization on GPUs took a few seconds per 3D pair, so the authors found virtually no speed difference between learning- and optimization-based methods.

Learn2Reg has continued with new tasks each year. As of October 2026, its website announces the 2026 edition at MICCAI 2026 in Strasbourg, with tasks on longitudinal whole-body PSMA PET/CT and inspiratory–expiratory chest CT, hosted on CodaBench. Data licences differ by task and are listed on its datasets page (for example, CC BY 4.0 for its NLST lung task).

Evaluation Metrics

• Landmark error (TRE). The distance in millimetres between corresponding anatomical landmarks after registration (image registration defines it); Learn2Reg also reports TRE30, the 30th percentile of the largest landmark distances, as a robustness score.
• Overlap of propagated labels. Dice between the fixed labels and the moving labels warped with the transform; Learn2Reg adds DSC30 (the 30th percentile of Dice over structures and cases) and the 95th-percentile Hausdorff distance. Klein et al. (2009) used target overlap, the intersection over the target region’s volume, and mean overlap, which is Dice. Overlap is most informative on labels the method did not see in training, which Learn2Reg’s DSC9 for its MR-CT task measures.
• Plausibility. The number or fraction of voxels with $\det J_\varphi \le 0$ (Balakrishnan et al., 2019; Dalca et al., 2019), and the standard deviation of $\log \det J_\varphi$ (SDlogJ), computed by central differences, which Learn2Reg uses as its smoothness score. Klein et al. also name inverse-consistency error and transitivity error as measures of a transform’s intrinsic properties.
• Runtime, measured on fixed hardware for Learn2Reg’s Docker submissions.

As the example shows, a registration can match intensities and propagated labels while its field is wrong inside structures or folds, which is why Learn2Reg combines accuracy, robustness, plausibility and runtime.

Applications

• Multimodal fusion of CT, MRI and PET of the same patient, and pre-operative with intra-operative images in image-guided surgery, as in the variants above.
• Longitudinal comparison. Follow-up scans are aligned to baselines; Learn2Reg’s 2026 PSMA PET/CT task targets therapy planning, dosimetry and response assessment, according to its website.
• Motion estimation. Thirion (1998) tracked heart motion in image sequences; Learn2Reg’s lung task serves monitoring disease, estimating motion in radiotherapy planning and assessing lung ventilation (Hering et al., 2023).
• Population studies, atlases and segmentation. Comparative morphometry and group analysis of brain images need correspondences across brains (Klein et al., 2009); abdominal registration across patients can model organ variation and provide a canonical atlas space (Hering et al., 2023); and atlases carry expert labels to new scans, as described above.

Open Problems

• Correspondence itself. Klein et al. (2009) note that one cannot expect every brain to map one-to-one onto every other, and that image correspondence is easily mistaken for anatomical correspondence.
• Validation without ground truth. Landmarks are sparse, labels are coarse and biased toward what was annotated, and surrogate scores can reward implausible fields, as VoxelMorph’s segmentation-only training shows.
• Generalization. Learned methods might lose accuracy on data other than their training set, a concern Learn2Reg tested; it found top methods robust on multi-site brain MRI and on a second lung dataset, but its authors caution that conclusions apply mainly to the selected tasks (Hering et al., 2023). Iglesias and Sabuncu (2015) also name robustness to changes in MRI contrast as an open need for multi-atlas methods.
• Pre-alignment and fields of view. Learn2Reg supplied pre-aligned images and left out affine or rigid pre-alignment, differing fields of view and voxel resolutions, and very high-resolution scans; its authors write that pre-alignment remains an active problem, particularly for deep-learning methods.
• Pathology and change. When a lesion appears, grows or is removed, the energy above has, by our reasoning, no correct answer at those voxels, and a registration flexible enough to match them erases the change it was meant to measure.