ElasticMM: Mass-Preserving Image Registration

• ElasticMM is an image registration framework that models source and target images as mass densities and enforces exact mass preservation via the Monge mass transport law.
• It minimizes linear elastic strain energy subject to the hard constraint ρ₀(x)=ρ₁(φ(x))·det(∇φ(x)), ensuring smooth and physically realistic deformations.
• The method utilizes a fully staggered finite-difference grid and an inexact SQP strategy to achieve high-accuracy registration for imaging modalities where intensities represent conserved physical quantities.

Searching arXiv for the primary paper and closely related registration work to ground the article. ElasticMM, short for Elastic Mass-preserving Matching, is a PDE-constrained image registration framework in which the source and target images are interpreted as mass densities and the deformation is required to satisfy the exact mass-preservation law from optimal transport while minimizing linear elastic strain energy. In "Elastic image registration with exact mass preservation" (Wlazło et al., 2016), the method is formulated for bounded Lipschitz domains $\Omega \subset \mathbb{R}^d$, with $d=2$ or $3$, and seeks a diffeomorphic map $\phi = id + u$ that matches $\rho_0$ to $\rho_1$ under the constraint

$\rho_0(x)=\rho_1(\phi(x))\,\det(\nabla \phi(x)).$

The framework is intended for modalities in which intensities represent conserved mass or quantity rather than brightness alone.

1. Registration model and mass-preservation principle

ElasticMM treats the source image $I_0$ and target image $I_1$ as mass densities $\rho_0$ and $d=2$0. The deformation is a map $d=2$1 with displacement $d=2$2 and Dirichlet boundary condition $d=2$3 on $d=2$4. The central modeling assumption is that the images reflect conserved mass, so the Jacobian determinant of the deformation must account exactly for local compression and expansion (Wlazło et al., 2016).

The exact mass-preservation law is

$d=2$5

or, in the displacement variable,

$d=2$6

A numerically convenient equivalent form is

$d=2$7

This constraint is the Monge mass conservation law from optimal transport. Classical $d=2$8 optimal transport would minimize $d=2$9 under this constraint. ElasticMM instead minimizes linear elastic strain energy, thereby enforcing spatial smoothness and physical plausibility beyond mere displacement magnitude. The analysis assumes positivity and continuity of densities, specifically a $3$0 such that $3$1 and $3$2 on $3$3, with $3$4 continuous, $3$5 and continuous up to the boundary, together with total mass equality

$3$6

The solution space is $3$7. For $3$8, the Jacobian determinant has favorable weak continuity properties; for $3$9, additional regularity or $\phi = id + u$0 with $\phi = id + u$1 is needed.

2. Variational formulation and optimality system

The objective functional is the linear elastic energy for a small-deformation, isotropic, homogeneous material,

$\phi = id + u$2

with Lamé parameters $\phi = id + u$3. ElasticMM solves

$\phi = id + u$4

The continuous Lagrangian is

$\phi = id + u$5

where the multiplier $\phi = id + u$6 enforces exact mass preservation (Wlazło et al., 2016). An augmented formulation with a penalty $\phi = id + u$7 can be written, but ElasticMM keeps mass as a hard constraint.

First-order optimality yields a saddle-point system:

• stationarity with respect to $\phi = id + u$8:

$\phi = id + u$9

• feasibility:

$\rho_0$0

Here $\rho_0$1 is the linear elasticity operator and $\rho_0$2 is the Jacobian of the constraint. The resulting structure is a classical KKT saddle point system, and the nonlinearity enters entirely through the mass-preservation constraint.

A common misconception is to treat mass-preserving registration as a brightness-based registration problem with an elastic regularizer added afterward. ElasticMM does not do this. The Jacobian determinant is part of the matching model itself, so intensity changes induced by local expansion and compression are explained geometrically rather than absorbed into a generic dissimilarity term.

3. Discretization on a fully staggered grid

The discretization uses a fully staggered finite-difference layout modeled on MAC grids. Scalar quantities such as $\rho_0$3, $\rho_0$4, and the multiplier $\rho_0$5 are stored at cell centers, while displacement components are stored on edge-centered staggered grids: $\rho_0$6 on horizontal edges and $\rho_0$7 on vertical edges in two dimensions. This arrangement stabilizes divergence-like constraints and prevents grid-scale checkerboarding (Wlazło et al., 2016).

Discrete divergence and curl are built with mimetic properties,

$\rho_0$8

and second-order finite differences are combined with bilinear projection operators between sub-grids. Rather than collocating $\rho_0$9 at cell centers, ElasticMM computes a discrete local cell-volume change

$\rho_1$0

which gives a second-order approximation to $\rho_1$1. The discrete constraint becomes

$\rho_1$2

For differentiability and accuracy, $\rho_1$3 is represented by cubic B-splines,

$\rho_1$4

so that analytic gradients are available in the linearization of $\rho_1$5. The elastic operator is discretized mimetically through

$\rho_1$6

with discrete energy

$\rho_1$7

This preserves the continuous Helmholtz-type decomposition and improves multigrid performance, especially near incompressibility.

4. Inexact SQP, saddle-point solves, and preconditioning

The nonlinear constrained problem is solved by an inexact SQP method embedded in a multiresolution strategy. At iterate $\rho_1$8, ElasticMM linearizes the constraint and solves the Newton system

$\rho_1$9

In practice, the method sets $\rho_0(x)=\rho_1(\phi(x))\,\det(\nabla \phi(x)).$0, omitting second derivatives of $\rho_0(x)=\rho_1(\phi(x))\,\det(\nabla \phi(x)).$1 to avoid prohibitive cost (Wlazło et al., 2016).

Each SQP iteration assembles the KKT matrix, approximately solves the linear system with GMRES and a block triangular preconditioner, performs an Armijo line search on the merit function

$\rho_0(x)=\rho_1(\phi(x))\,\det(\nabla \phi(x)).$2

and applies a diffeomorphic safeguard by shrinking the step length until

$\rho_0(x)=\rho_1(\phi(x))\,\det(\nabla \phi(x)).$3

Stopping is based on either the KKT residual or the constraint norm.

The preconditioner has the approximate form

$\rho_0(x)=\rho_1(\phi(x))\,\det(\nabla \phi(x)).$4

The inverse of $\rho_0(x)=\rho_1(\phi(x))\,\det(\nabla \phi(x)).$5 is approximated by a multigrid V-cycle with Gauss-Seidel smoothing on an augmented system involving

$\rho_0(x)=\rho_1(\phi(x))\,\det(\nabla \phi(x)).$6

while the Schur complement inverse is approximated by a commutator-based formula,

$\rho_0(x)=\rho_1(\phi(x))\,\det(\nabla \phi(x)).$7

with a direct LU applied to $\rho_0(x)=\rho_1(\phi(x))\,\det(\nabla \phi(x)).$8. The intent is to avoid inverting the full KKT system while retaining mesh- and parameter-independent multigrid performance.

5. Workflow, parameters, and empirical behavior

ElasticMM uses coarse-to-fine multiresolution. Images are normalized to densities, smoothed, clamped to $\rho_0(x)=\rho_1(\phi(x))\,\det(\nabla \phi(x)).$9, and scaled to equal total mass. The description gives $I_0$0 as an example for avoiding near-zero densities that destabilize the mass constraint (Wlazło et al., 2016).

At each resolution level, densities are restricted, $I_0$1 are initialized either at zero or by prolongation from the coarser level, and the SQP loop iterates until either $I_0$2 or the residual is below tolerance. The final output consists of the displacement $I_0$3, the deformation $I_0$4, the warped image $I_0$5 for visualization, and the registered density $I_0$6 for density consistency.

The reported experiments include a tailored example on $I_0$7 with $I_0$8, a uniform $I_0$9 grid, and parameters $I_1$0, $I_1$1. Metrics are

• $I_1$2,
• $I_1$3,
• $I_1$4.

The mass-preserving solution matches the ground-truth mass constraint exactly, and the elastic energy is bounded by the energy of the exact deformation. A real-world example uses brain cortical tissue slices on a $I_1$5 grid, pre-smoothed, rescaled to $I_1$6 with $I_1$7, and solved over levels $I_1$8 with $I_1$9, $\rho_0$0. About 36 SQP iterations total were reported, with deformation grids showing near-perfect alignment while preserving mass with high accuracy. A regularized baseline with $\rho_0$1 did not reach the same DMP tolerance; smaller $\rho_0$2 was needed but increased ill-conditioning.

The ratio $\rho_0$3 controls the preference for divergence-free components. Large $\rho_0$4 allows more compressible solutions, whereas $\rho_0$5 biases the displacement toward divergence-free behavior.

6. Relation to standard registration and neighboring elastic methods

ElasticMM differs fundamentally from standard elastic registration of the form

$\rho_0$6

because it does not rely on SSD, cross-correlation, or mutual information to explain intensity changes. Instead, it enforces the exact conservation law

$\rho_0$7

and is therefore parameter-free in the strict constrained form with respect to the usual regularization weight $\rho_0$8 (Wlazło et al., 2016).

Its closest conceptual relative is optimal transport, with which it shares the exact mass constraint. The difference is that ElasticMM replaces the classical $\rho_0$9 transport cost with elastic strain energy, thereby promoting smooth, physically plausible, diffeomorphic mappings. In the special case $d=2$00, the constraint reduces to

$d=2$01

and the KKT structure resembles a Stokes-type system.

Within elasticity-based matching more broadly, other approaches minimize the symmetric difference of shapes rather than imposing exact mass preservation. "Elasticity-based Matching by Minimizing the Symmetric Difference of Shapes" (Simon et al., 2015), for example, formulates 2D matching through linearized elasticity and a symmetric-difference objective on contours. That comparison makes clear that ElasticMM is not a generic elastic deformation model; its defining feature is the hard mass-preservation law.

The principal limitations follow directly from that choice. ElasticMM requires positive, calibrated density images with equal mass; ensuring $d=2$02 throughout the domain can be challenging; the problem is nonconvex and therefore susceptible to local minima; and three-dimensional use requires additional regularity and substantially heavier numerics. These constraints also define its proper domain of application: settings in which image intensity is a conserved physical quantity and local compression or expansion should be modeled explicitly rather than treated as appearance variation.