A Lagrangian Gauss-Newton-Krylov Solver for Mass- and   Intensity-Preserving Diffeomorphic Image Registration

Andreas Mang; Lars Ruthotto

arXiv:1703.04446·cs.CV·November 2, 2017

A Lagrangian Gauss-Newton-Krylov Solver for Mass- and Intensity-Preserving Diffeomorphic Image Registration

Andreas Mang, Lars Ruthotto

PDF

TL;DR

This paper introduces an efficient Lagrangian Gauss-Newton-Krylov solver for diffeomorphic image registration that preserves mass and intensity, supporting both stationary and non-stationary velocity fields with scalable optimization techniques.

Contribution

It develops a novel PDE constraint elimination method using Lagrangian hyperbolic PDE solvers and integrates it into a scalable, efficient registration framework supporting various models.

Findings

01

Supports large-scale problems with up to 14.7 million degrees of freedom.

02

Achieves fast convergence with spectral preconditioners and Gauss-Newton methods.

03

Demonstrates effectiveness on synthetic and real-world datasets.

Abstract

We present an efficient solver for diffeomorphic image registration problems in the framework of Large Deformations Diffeomorphic Metric Mappings (LDDMM). We use an optimal control formulation, in which the velocity field of a hyperbolic PDE needs to be found such that the distance between the final state of the system (the transformed/transported template image) and the observation (the reference image) is minimized. Our solver supports both stationary and non-stationary (i.e., transient or time-dependent) velocity fields. As transformation models, we consider both the transport equation (assuming intensities are preserved during the deformation) and the continuity equation (assuming mass-preservation). We consider the reduced form of the optimal control problem and solve the resulting unconstrained optimization problem using a discretize-then-optimize approach. A key contribution is…

Tables2

Table 1. Table 1: Commonly used symbols and abbreviations.

DCT	discrete cosine transform
FAIR	Flexible Algorithms for Image Registration [54]
LDDMM	large deformation diffeomorphic metric mapping
MRI	magnetic resonance imaging
PDE	partial differential equation
PET	positron emission tomography
PIC	particle-in-cell (method)
SSD	sum-of-squared-differences
$ℛ (x)$	reference / fixed image
$𝒯 (x)$	template image (image to be registered)
$𝒞$	PDE constraint
$𝒟$	distance or similarity measure
$𝒮$	regularization model (smoother)
$α$	regularization weight
$x$	spatial coordinate; $x \in Ω$
$Ω$	spatial domain; $Ω \subset ℝ^{d}$
$v (x, t)$	velocity field
$u (x, t)$	transported image intensities
$y (x)$	transformation / mapping
$N$	number of time steps for computing the characteristic
$n$	number of unknowns (i.e., the dimension of the discretized velocity field)
$n_{t}$	number of cells in space-time grid
$I$	interpolation operator
$\nabla$	gradient operator
$\nabla \cdot$	divergence operator
$\partial_{t}$	time derivative

Table 2. Table 2: Registration quality. We report registration results as a function of the regularization weight α 𝛼 \alpha for the first two datasets of the NIREP repository. We compare the proposed method considering stationary and non-stationary velocity fields. We report results for the curvature ( H 2 superscript 𝐻 2 H^{2} ) regularization model (theoretically required to guarantee the existence of a diffeomorphic deformation map) and a diffusive ( H 1 superscript 𝐻 1 H^{1} ) regularization model. We compare the proposed method to a formulation for diffeomorphic image registration based on a hyperelastic regularization method [ 15 ] . We use the experimentally determined regularization weights α 1 = 100 subscript 𝛼 1 100 \alpha_{1}=100 (length regularizer), α 2 = 10 subscript 𝛼 2 10 \alpha_{2}=10 (surface regularizer), and α 3 = 100 subscript 𝛼 3 100 \alpha_{3}=100 (volume regularizer). We report values (from left to right) for the Dice coefficient after registration and the extremal values for the Jacobian ( min \min and max \max ). The initial value for the Dice coefficient for the considered datasets is 5.543 534 386 E − 01 5.543534386E-01 5.543\,534\,38610-01 .

method	run	$α$	dice	$\min (det \nabla y)$	$\max (det \nabla y)$	time (speedup)
hyperelastic	#1	100, 10, 100	$7.933 170 596 E - 01$	$4.513 155 512 E - 01$	$1.862 803 169 E + 00$	$3.706 054 800 E + 03$	( $1.000 000 000$ )
stationary LDDMM
curvature	#2	10	$7.764 123 715 E - 01$	$8.876 992 884 E - 02$	$8.793 200 978 E + 00$	$1.816 074 592 E + 03$	( $2.040 695 254$ )
	#3	25	$7.548 417 772 E - 01$	$1.193 601 635 E - 01$	$5.766 703 129 E + 00$	$1.855 497 767 E + 03$	( $1.997 337 246$ )
	#4	50	$7.336 621 571 E - 01$	$1.789 841 204 E - 01$	$4.735 466 508 E + 00$	$1.447 520 290 E + 03$	( $2.560 278 309$ )
	#5	100	$7.141 735 014 E - 01$	$2.369 603 205 E - 01$	$3.132 339 125 E + 00$	$1.474 501 181 E + 03$	( $2.513 429 523$ )
diffusion	#6	300	$7.860 163 795 E - 01$	$7.789 530 517 E - 02$	$6.380 539 417 E + 00$	$1.957 822 580 E + 03$	( $1.892 947 215$ )
	#7	400	$7.750 227 852 E - 01$	$1.017 359 434 E - 01$	$6.264 590 010 E + 00$	$1.941 301 868 E + 03$	( $1.909 056 423$ )
	#8	500	$7.659 294 175 E - 01$	$1.178 093 780 E - 01$	$6.166 939 898 E + 00$	$1.948 818 006 E + 03$	( $1.901 693 636$ )
non-stationary LDDMM
curvature	#9	5	$7.595 132 911 E - 01$	$9.545 388 132 E - 02$	$6.357 335 053 E + 00$	$7.640 378 580 E + 03$	( $0.485 061 671$ )
	#10	10	$7.595 373 122 E - 01$	$8.480 122 115 E - 02$	$5.726 130 378 E + 00$	$7.498 044 615 E + 03$	( $0.494 269 505$ )
	#11	25	$7.467 159 258 E - 01$	$1.310 324 124 E - 01$	$5.677 905 399 E + 00$	$6.178 465 357 E + 03$	( $0.599 834 196$ )
	#12	50	$7.350 425 329 E - 01$	$2.123 510 134 E - 01$	$3.680 210 984 E + 00$	$5.980 091 501 E + 03$	( $0.619 732 123$ )
diffusion	#13	300	$8.053 458 765 E - 01$	$7.705 591 708 E - 02$	$6.475 742 122 E + 00$	$3.293 675 673 E + 03$	( $1.125 203 319$ )
	#14	400	$7.925 344 925 E - 01$	$1.076 371 450 E - 01$	$5.395 472 375 E + 00$	$3.342 763 690 E + 03$	( $1.108 679 866$ )
	#15	500	$7.672 242 420 E - 01$	$1.528 603 048 E - 01$	$6.498 834 778 E + 00$	$2.588 212 381 E + 03$	( $1.431 897 485$ )

Equations46

v, u min {J (v, u) := D (u (\cdot, 1), R) + α S (v)} subject to C (v, u) = 0,

v, u min {J (v, u) := D (u (\cdot, 1), R) + α S (v)} subject to C (v, u) = 0,

\mathcal{C}(u,v)=\left\{\begin{array}[]{l}\partial_{t}u(x,t)+v(x,t)\cdot\nabla u(x,t)=0,\\ u(x,0)=\mathcal{T}(x)\end{array}\right.

\mathcal{C}(u,v)=\left\{\begin{array}[]{l}\partial_{t}u(x,t)+v(x,t)\cdot\nabla u(x,t)=0,\\ u(x,0)=\mathcal{T}(x)\end{array}\right.

\mathcal{C}(u,v)=\left\{\begin{array}[]{l}\partial_{t}u(x,t)+\nabla\cdot(u(x,t)v(x,t))=0,\\ u(x,0)=\mathcal{T}(x).\end{array}\right.

\mathcal{C}(u,v)=\left\{\begin{array}[]{l}\partial_{t}u(x,t)+\nabla\cdot(u(x,t)v(x,t))=0,\\ u(x,0)=\mathcal{T}(x).\end{array}\right.

\partial_{t} y (v, y_{0}, 0, t) = v (y (v, y_{0}, 0, t), t) and y (v, y_{0}, 0, 0) = y_{0} .

\partial_{t} y (v, y_{0}, 0, t) = v (y (v, y_{0}, 0, t), t) and y (v, y_{0}, 0, 0) = y_{0} .

S^{diff} (v) = \frac{1}{2} \int_{0}^{1} \int_{Ω} k = 1 \sum d ∣\nabla v^{k} (x, t) ∣^{2} d x d t

S^{diff} (v) = \frac{1}{2} \int_{0}^{1} \int_{Ω} k = 1 \sum d ∣\nabla v^{k} (x, t) ∣^{2} d x d t

S^{curv} (v) = \frac{1}{2} \int_{0}^{1} \int_{Ω} k = 1 \sum d ∣Δ v^{k} (x, t) ∣^{2} d x d t,

S^{curv} (v) = \frac{1}{2} \int_{0}^{1} \int_{Ω} k = 1 \sum d ∣Δ v^{k} (x, t) ∣^{2} d x d t,

y (v, x, 0, t_{k + 1}) = y (v, x, 0, t_{k}) + Δ t I (v, y (v, x, 0, t_{k}), t_{k}), \forall k = 0, 1, \dots, N - 1,

y (v, x, 0, t_{k + 1}) = y (v, x, 0, t_{k}) + Δ t I (v, y (v, x, 0, t_{k}), t_{k}), \forall k = 0, 1, \dots, N - 1,

u (x_{c}, 1) = T (y (v, x_{c}, 1, 0)) .

u (x_{c}, 1) = T (y (v, x_{c}, 1, 0)) .

u (x_{c}, 1) = F (y (v, x_{c}, 0, 1)) T (x_{c}) .

u (x_{c}, 1) = F (y (v, x_{c}, 0, 1)) T (x_{c}) .

b^{δ} (x, y_{j}) = ⎩ ⎨ ⎧ 1 - (y_{j} - x) / δ 1 + (x - y_{j}) / δ 0 for y_{j} - δ \leq x \leq y_{j}, for y_{j} < x \leq y_{j} + δ, else,

b^{δ} (x, y_{j}) = ⎩ ⎨ ⎧ 1 - (y_{j} - x) / δ 1 + (x - y_{j}) / δ 0 for y_{j} - δ \leq x \leq y_{j}, for y_{j} < x \leq y_{j} + δ, else,

u_{i} (v, 1) = j = 1 \sum n_{p} \int_{x_{i}}^{x_{i + 1}} u_{j} b^{δ} (x, y_{j}) d x = j = 1 \sum n_{p} u_{j} (B^{δ} (x_{i + 1}, y_{j}) - B^{δ} (x_{i}, y_{j})) = j = 1 \sum n_{p} u_{j} F_{ij},

u_{i} (v, 1) = j = 1 \sum n_{p} \int_{x_{i}}^{x_{i + 1}} u_{j} b^{δ} (x, y_{j}) d x = j = 1 \sum n_{p} u_{j} (B^{δ} (x_{i + 1}, y_{j}) - B^{δ} (x_{i}, y_{j})) = j = 1 \sum n_{p} u_{j} F_{ij},

B^{δ} (x, y_{j}) = ⎩ ⎨ ⎧ 0 x - \frac{1}{2 δ} (2 y_{j} x - x^{2}) x + \frac{1}{2 δ} (x^{2} - 2 y_{j} x) 1 for x < y_{j} - δ, for y_{j} - δ \leq x \leq y_{j}, for y_{j} < x \leq y_{j} + δ, for y_{j} + δ < x .

B^{δ} (x, y_{j}) = ⎩ ⎨ ⎧ 0 x - \frac{1}{2 δ} (2 y_{j} x - x^{2}) x + \frac{1}{2 δ} (x^{2} - 2 y_{j} x) 1 for x < y_{j} - δ, for y_{j} - δ \leq x \leq y_{j}, for y_{j} < x \leq y_{j} + δ, for y_{j} + δ < x .

v min {J (v) := D (u_{1} (v), R) + α S (v)},

v min {J (v) := D (u_{1} (v), R) + α S (v)},

D^{SSD} (u_{1} (v), R) = \frac{h ^{d}}{2} res (v)^{⊤} res (v), where res (v) := u_{1} (v) - R .

D^{SSD} (u_{1} (v), R) = \frac{h ^{d}}{2} res (v)^{⊤} res (v), where res (v) := u_{1} (v) - R .

d_{v} J (v) = d_{v} u_{1} (v) d_{u_{1}} D (u_{1} (v), R) + α d_{v} S,

d_{v} J (v) = d_{v} u_{1} (v) d_{u_{1}} D (u_{1} (v), R) + α d_{v} S,

d_{u_{1}} D^{SSD} (u_{1} (v), R) = h^{d} res (v) .

d_{u_{1}} D^{SSD} (u_{1} (v), R) = h^{d} res (v) .

d_{v} S (v) = Δ t h^{d} (B^{⊤} B) v .

d_{v} S (v) = Δ t h^{d} (B^{⊤} B) v .

H (v) \approx d_{2} J (v) = d_{v} u_{1} (v) d_{2} D (u_{1} (v), R) d_{v} u_{1} (v)^{⊤} + α Δ t h^{d} B^{⊤} B + γ I_{n},

H (v) \approx d_{2} J (v) = d_{v} u_{1} (v) d_{2} D (u_{1} (v), R) d_{v} u_{1} (v)^{⊤} + α Δ t h^{d} B^{⊤} B + γ I_{n},

d_{2} D^{SSD} = h^{3} I_{n} .

d_{2} D^{SSD} = h^{3} I_{n} .

d_{v} u_{1} (v) = d_{v} T (y (v, x_{c}, 1, 0)) = \nabla T (y (v, x_{c}, 1, 0)) d_{v} y (v, x_{c}, 1, 0) .

d_{v} u_{1} (v) = d_{v} T (y (v, x_{c}, 1, 0)) = \nabla T (y (v, x_{c}, 1, 0)) d_{v} y (v, x_{c}, 1, 0) .

d_{v} u_{1} (v) = d_{v} (F (y (v, x_{c}, 0, 1)) T (x_{c})) = d_{y} (F (y (v, x_{c}, 0, 1)) T (x_{c})) d_{v} y (v, x_{c}, 0, 1) .

d_{v} u_{1} (v) = d_{v} (F (y (v, x_{c}, 0, 1)) T (x_{c})) = d_{y} (F (y (v, x_{c}, 0, 1)) T (x_{c})) d_{v} y (v, x_{c}, 0, 1) .

d_{v} y (v, x_{c}, 1, t_{k + 1}) = d_{v} y (v, x_{c}, 1, t_{k}) + \dots + Δ t d_{v} I (v, y (v, x_{c}, 1, t_{k}), t_{k}) + \dots Δ t d_{y} I (v, y (v, x_{c}, 1, t_{k}), t_{k}) d_{v} y (v, x_{c}, 1, t_{k}),

d_{v} y (v, x_{c}, 1, t_{k + 1}) = d_{v} y (v, x_{c}, 1, t_{k}) + \dots + Δ t d_{v} I (v, y (v, x_{c}, 1, t_{k}), t_{k}) + \dots Δ t d_{y} I (v, y (v, x_{c}, 1, t_{k}), t_{k}) d_{v} y (v, x_{c}, 1, t_{k}),

H (v^{i}) δ v = - d_{v} J (v^{i}) .

H (v^{i}) δ v = - d_{v} J (v^{i}) .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

A Lagrangian Gauss–Newton–Krylov Solver for Mass- and Intensity-Preserving Diffeomorphic Image Registration

Andreas Mang Institute for Computational Engineering and Sciences, University of Texas at Austin, Austin, Texas, USA. ([email protected])

Lars Ruthotto Department of Mathematics and Computer Science, Emory University, Atlanta, Georgia, USA. ([email protected])

Abstract

We present an efficient solver for diffeomorphic image registration problems in the framework of Large Deformations Diffeomorphic Metric Mappings (LDDMM). We use an optimal control formulation, in which the velocity field of a hyperbolic PDE needs to be found such that the distance between the final state of the system (the transformed/transported template image) and the observation (the reference image) is minimized. Our solver supports both stationary and non-stationary (i.e., transient or time-dependent) velocity fields. As transformation models, we consider both the transport equation (assuming intensities are preserved during the deformation) and the continuity equation (assuming mass-preservation).

We consider the reduced form of the optimal control problem and solve the resulting unconstrained optimization problem using a discretize-then-optimize approach. A key contribution is the elimination of the PDE constraint using a Lagrangian hyperbolic PDE solver. Lagrangian methods rely on the concept of characteristic curves. We approximate these curves using a fourth-order Runge-Kutta method. We also present an efficient algorithm for computing the derivatives of the final state of the system with respect to the velocity field. This allows us to use fast Gauss-Newton based methods. We present quickly converging iterative linear solvers using spectral preconditioners that render the overall optimization efficient and scalable. Our method is embedded into the image registration framework FAIR and, thus, supports the most commonly used similarity measures and regularization functionals. We demonstrate the potential of our new approach using several synthetic and real world test problems with up to 14.7 million degrees of freedom.

keywords:

Diffeomorphic Image Registration, Large Deformation Diffeomorphic Metric Mapping, Optimal Control, PDE-Constrained Optimization, Lagrangian Methods

AMS:

68U10, 49J20, 35Q93, 65M32, 76D55, 65K10.

1 Introduction

In this paper, we present efficient numerical methods for diffeomorphic image registration in the framework of Large Deformation Diffeomorphic Metric Mapping (LDDMM) [65, 23, 9]. We use an optimal control formulation similar to the one in [11, 12, 10, 48, 50, 49]. Here, the task is to find a smooth velocity field $v$ such that the distance between two images (or densities), $\mathcal{T}$ (the template image) and $\mathcal{R}$ (the reference image), is minimized, subject to some regularization norm for $v$ and a transformation model, given by a hyperbolic partial differential equation (PDE) that models the deformation of $\mathcal{T}$ . We consider the transport equation (assuming intensities are related at corresponding points) and the continuity equation (assuming that mass is preserved) as constraints. The connection to traditional image registration formulations [53, 54] is that a sufficiently smooth velocity field $v$ gives rise to a diffeomorphism $y$ via the method of characteristics. Vice versa, representing diffeomorphisms through velocity fields has been used for efficient statistical analysis; see, e.g., [3].

Solving the variational problem associated with LDDMM is, in theory, known to yield a diffeomorphic transformation $y$ if $v$ is sufficiently smooth [23, 65, 9]. Although, the theory of diffeomorphic registration using LDDMM is well explored [52, 70, 71], efficient numerical optimization is not. Until recently [5, 39, 48, 50, 49] mostly first-order optimization methods were used; see, e.g., [9, 67, 17, 58]. A key component in LDDMM is the numerical method for solving the hyperbolic PDE. Hyperbolic PDE solvers can be roughly divided into Eulerian (in which the density is discretized at the same locations for each time point) and Lagrangian (in which the grid moves over time along the characteristic curves) solvers; see also [45, 25, 46]. Intermediates are Semi-Lagrangian (SL) methods, which follow the velocity for a short time step and then estimate the density at fixed points. In this work, we use a Lagrangian solver.

It is well known that explicit Eulerian methods for hyperbolic PDEs require the size of the time step to be sufficiently small to ensure numerical stability. The maximal admissible time step size depends on the accuracy of the spatial discretization and the magnitude of the velocities. In optimal control problems, like ours, the velocity field $v$ is not known a priori and, thus, it is difficult to come up with an efficient and stable choice of the time step. SL and Lagrangian solvers are explicit methods that, unlike explicit Eulerian methods, are stable without a restriction on the maximal admissible time step size. One drawback of SL methods is their memory requirements. For efficient derivative (or sensitivity) computations in Gauss–Newton type optimization schemes, we have to store intermediate images [51]. Further, SL methods require a repeated interpolation of the initial image and may therefore introduce severe dissipation if implemented naively.111We note, that interpolation errors and numerical diffusion can be minimized, by, e.g., applying high-order interpolation schemes and/or evaluating the interpolation on a finer grid; this is costly. As we will see, Lagrangian methods require only one image interpolation at the final time. Secondly, derivatives of the Lagrangian solver can be obtained efficiently, without storing intermediate variables. A feature of SL and Lagrangian methods is that they can be easily modified to solve intensity- and mass-preserving problems, since the characteristic curves coincide.

The key idea of our work is to use a discretize-then-optimize strategy based on a Lagrangian hyperbolic PDE solver to efficiently solve the reduced formulation of the PDE-constrained optimization problem arising in LDDMM. We show that Lagrangian methods lead to a finite dimensional optimization problem that can be solved using an inexact Gauss–Newton method. Our PDE solver requires numerical computation of the characteristic curves. We use a fourth-order Runge–Kutta (RK4) method to numerically approximate the transformation $y$ that is associated with a, in general non-stationary, velocity field $v$ . As we show, derivatives of the transformation with respect to $v$ can be derived analytically and computed efficiently. Due to the hyperbolic nature of the PDE, the derivatives can be represented as sparse matrices; the procedure can be paralellelized: characteristics starting at different points can be computed independently. Given these characteristics, the hyperbolic PDEs can be solved by a single interpolation step (for the advection equation) or the particle-in-cell method (for the continuity equation).

1.1 Contributions

•

We propose a discretize-optimize method for solving LDDMM using a Lagrangian hyperbolic PDE solver. Our scheme is based on an RK4 method to approximate the characteristic curves and we derive an efficient algorithm for computing the derivative of the solver with respect to the velocities.

•

The storage requirement of our method is independent on the number of time steps used in the numerical solver. Also, the Hessian of the objective function can be build explicitly at moderate costs, which is useful, e.g., to accelerate matrix-vector products. In this work, we use and numerically study the convergence of spectral preconditioners to iteratively solve the Gauss–Newton system.

•

We extend the LDDMM framework to mass-preserving registration, which has been proved to be an adequate model for many relevant biomedical applications involving with density images, e.g., in [16, 59, 21, 30, 61].

•

We derive a flexible framework supporting both stationary and non-stationary velocity fields. Our methods are embedded into the FAIR framework [54]. This allows us to consider different regularization norms and distance measures. Our implementation is freely available as an add-on to FAIR at:

https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM

•

We provide detailed numerical experiments on four different data sets that demonstrate the flexibility and effectiveness of our method. We show that our prototype implementation is competitive to state-of-the-art packages for diffeomorphic image registration [15, 60]. We study registration quality for synthetic benchmark problems and real-world applications leading to optimization problem with up to 14.7 million degrees of freedom.

1.2 Related Work

We limit this review to work closely related to ours. For a general insight into the area of image registration, its applications, and its formulation we refer to [53, 54, 63]. Our work builds upon the LDDMM framework described in [23, 65, 9], which is based on the pioneering work on velocity-based fluid registration described in [20]. We adopt an optimal control point of view; we also do not directly invert for the transformation $y$ but for its velocity $v$ . We arrive at a hyperbolic PDE-constrained optimization problem. We refer to [31, 43, 41, 13] for a general introduction into optimal control theory and developments in PDE-constrained optimization. Related optimal control formulations for diffeomorphic image registration can, e.g., be found in [11, 12, 38, 10, 17, 48, 50, 49, 51]. Other formulations for velocity-based diffeomorphic image registration are described in [4, 66, 5]. Our work also shares characteristics with optical flow formulations [11, 12, 17]. Our formulation for mass-preserving registration problems is related to the Monge–Kantorovich functional arising in optimal mass transport [10, 35].

Most work on large deformation diffeomorphic image registration still considers first-order methods for numerical optimization (see, e.g., [6, 7, 8, 11, 9, 38, 44, 17, 67]); the exceptions are [5, 48, 50, 49, 66]. First-order schemes for numerical optimization do in general require a larger number of iterations than Newton type optimization schemes.222We note that in LDDMM most implementations use a gradient descent scheme in the Sobolev space induced by the regularization operator (dual space). This leads to a significant speedup compared to standard gradient descent approaches. However, it has been demonstrated experimentally that Gauss–Newton–Krylov methods are superior [48]. The work in [5] uses geodesic shooting and estimates the initial value of a non-stationary velocity field that parameterizes the diffeomorphism $y$ . Other approaches that reduce the size of the optimization problem are based on stationary velocity fields; see, e.g., [40, 47, 50, 49, 51].

PDE-constrained optimization commonly requires a repeated solution of the forward problem. Thus, the design of an efficient forward solver is critical. The approaches described in [10, 17, 48, 50, 49, 38] are based on an Eulerian formulation. They employ explicit high-order schemes [11, 12, 48, 50, 38, 58], which suffer from a restriction on a maximally admissible time step, implicit schemes [10], or explicit SL schemes [49, 51, 17]. The latter were originally proposed in the context of weather prediction [64]. They are a hybrid between Lagrangian and Eulerian schemes, and unconditionally stable. SL schemes have been used in the context of Lagrangian formulations for diffeomorphic image registration (to compute the characteristics) [9, 40]. Conditionally stable schemes require small time steps, which can result in a significant amount of memory that needs to be allocated to store the time-space fields necessary to evaluate the gradient or Hessian. This makes a direct application of these type of methods to large-scale 3D problems challenging. One remedy is to turn to parallel architectures and use sophisticated checkpointing schemes in time to reduce the amount of memory that has to be allocated; see, e.g., [1]. Implicit schemes typically suffer from severe numerical diffusion. The same is true for straightforward implementations of SL schemes. Pure Lagrangian schemes for diffeomorphic image registration have been described in [6, 7, 8]. The time integration for computing the characteristics in [6, 7, 8] is based on a first-order explicit scheme.

What sets our work apart is the numerical solver and the generalization of our formulation to both the transport and the continuity equation. Most existing works on optimal control formulations for diffeomorphic image registration consider an optimize-then-discretize approach [17, 48, 50, 49, 51]. We use a discretize-then-optimize strategy instead.333Advantages and disadvantages of these two techniques are discussed, e.g., in [31]. Similar to [48], we describe a method that can handle stationary and non-stationary velocity fields. Our numerical scheme is, likewise to [6, 7, 8], based on a purely Lagrangian approach. We consider a reduced formulation of the PDE-constrained optimization problem arising in LDDMM, i.e., we eliminate the hyperbolic PDE constraint (state equation) from the variational problem. We use a Lagrangian solver to parameterize the final state in terms of the velocity. In doing so we avoid many of the complications we reviewed above: our method is unconditionally stable, limits numerical diffusion, and does not require the storage of multiple space-time fields for the evaluation of the gradient or Hessian operators. The work in [6, 7, 8] uses first-order information for numerical optimization and a first-order accurate explicit time integrator to compute the characteristic. We use a fourth-order, explicit RK scheme instead. We derive expressions for the exact derivative of the characteristics. We arrive at a Gauss–Newton–Krylov scheme that—combined with an efficient iterative linear solver—yields an approximate solution within a few iterations, and has an overall algorithmic complexity of $\mathcal{O}(n\log(n))$ , where $n$ is the dimension of the discretized velocity field (i.e., the number of unknowns).

2 Mathematical Formulation

We describe the variational optimal control formulation of the LDDMM problem next. We denote the image domain by $\Omega\subset\mathbb{R}^{d}$ , where $d\in\{2,3\}$ represents the spatial dimension. We assume that the template and the reference image, denoted by $\mathcal{T}:\Omega\to\mathbb{R}$ and $\mathcal{R}:\Omega\to\mathbb{R}$ , are compactly supported on $\Omega$ and continuously differentiable. Given these two images, the task of image registration is to find a plausible transformation $y:\Omega\to\mathbb{R}^{d}$ so that the transformed template image $\mathcal{T}\circ y$ becomes similar to the reference image $\mathcal{R}$ [54]. The definitions of plausibility, similarity, and the transformation model depend on the context; see [54, 53, 63] for examples. Many relevant applications, e.g., in medical imaging, require that plausible transformations are diffeomorphic, i.e., smooth mappings with a smooth inverse. One framework that contains the most commonly used definitions of these three terms within the medical imaging application domain is LDDMM [23, 65, 9]. The variational optimal control formulation of the LDDMM problem can, in general format, be stated as follows:

[TABLE]

where $\mathcal{D}$ is a distance (or similarity) measure, $\mathcal{S}$ is a regularizer (smoother), $v:\Omega\times[0,1]\to\mathbb{R}^{d}$ is the sought after velocity field, and $u:\Omega\times[0,1]\to\mathbb{R}$ is a time series of images. In an optimal control context, $v$ is commonly referred to as the control variable and $u$ as the state variable. Here, $\alpha>0$ is a regularization parameter that balances between minimizing the image distance and the smoothness of the velocity field (and consequently controls the properties of the resulting transformation). In the current work, we assume that $\alpha$ is chosen by the user and refer to [36, 34, 68] for some works on automatic selection of the parameters.444Examples for an automatic selection of the regularization parameter in the context of image registration can, e.g., be found in [33, 48].

In this work, we assume that the constraint $\mathcal{C}$ , which describes the transformation model, is given either by the advection (also called transport) equation

[TABLE]

or the continuity equation

[TABLE]

The former assumes intensity values are preserved during the transformation; the latter preserves the overall mass of the image. The choice of the transformation model depends on the application. Intensity preservation is commonly used, e.g., for registration of medical images acquired from different subjects [56]. Mass-preservation has been successfully used, e.g., for motion correction in position emission tomography (PET) [21, 30] or artifact correction of magnetic resonance imaging (MRI) [16, 61].

The models in (2) and (3) can be used to establish point-to-point correspondences between the template and the reference image. One way of showing this is the method of characteristics [46, 24]. To better illustrate this, we consider the advection equation (2), for which the intensity $u$ is constant along the characteristics. This means that, for all $y_{0}\in\Omega$ and $t\in[0,1]$ , it holds that $u(y(v,y_{0},0,t),t)=\mathcal{T}(y_{0})$ where the characteristic curve $t\mapsto y(v,y_{0},0,t)$ satisfies

[TABLE]

Similarly, the characteristics can be traced backwards in time to compute the state at some point $y_{1}\in\mathbb{R}^{d}$ at $t=1$ . The position of the point is given by $t\mapsto y(v,y_{1},1,t)$ , which satisfies (4) with final time condition $y(v,y_{1},1,1)=y_{1}$ . Clearly, both operations are inverse to one another. That is, for all $x\in\mathbb{R}^{n}$ it holds that $y(v,y(v,x,0,1),1,0)=x$ , i.e., a composition of both maps yields identity.

Note that (2) through (4) involve a non-stationary (i.e., time-dependent) velocity field $v$ . This is a key assumption in the original LDDMM formulation [9]. To make the problem computational tractable it is, however, often assumed that $v$ is stationary (stationary velocity field based registration) [3, 4, 40, 50, 51, 57].555Another strategy to reduce the computational burden is to invert for an initial momentum [67] that encodes the trajectory of the diffeomorphism. The numerical framework we propose in this paper can efficiently handle both stationary and non-stationary velocities. As demonstrated in our numerical experiments in Sec. 4.2, the stationary model is less flexible in that we can only invert for a subset of the deformation maps living on the manifold of diffeomorphisms. Our experiments also suggest that stationary velocity fields are adequate for registration problems involving two topologically similar images, yielding little to no difference in the recovered deformation map [4, 40, 48]. However, using non-stationary velocity fields may become critical in applications involving large and highly nonlinear transformations and/or the registration of time series of images with large motion between time frames (e.g., typically seen in tracking or optical flow problems).

2.1 Regularization functionals

Due to the ill-posedness of the image registration problem, the literature on regularization in image registration is rich; see, e.g., [53, 27, 54, 63] for extensive overviews. It is established that the existence and regularity of a diffeomorphic map $y$ depends on the smoothness of the velocity field $v$ as well as the smoothness of the images $\mathcal{R}$ and $\mathcal{T}$ [23, 65, 50, 17]. Modeling the images as functions of bounded variation, $H^{3}$ -regularity [17] is required (assuming that $v$ is divergence free). For continuous images we can relax the $H^{3}$ -regularity to an $H^{2}$ -regularity [9] or—under additional assumptions on the divergence of $v$ —even to an $H^{1}$ -regularity [17].666We use interpolation and padding of the discrete image data to obtain continuously differentiable and compactly supported functions (see Sec. 3). Most implementations for LDDMM consider an $H^{2}$ -norm for $\mathcal{S}$ in (1) (or an approximation based on a Gaussian kernel within a gradient descent scheme in the Sobolev space induced by the regularization norm); see, e.g., [5, 9, 40].

In our framework, we regard the regularizer as a modular component that can be replaced or extended. In practical applications we control the regularization parameter by monitoring the Jacobian in an attempt to generate transformations that are diffeomorphic (in a discrete setting) and yield high-fidelity (low mismatch) results. In our numerical examples, we consider $H^{1}$ - and $H^{2}$ -seminorms as regularization models. These type of regularization models are also referred to as diffusive or curvature regularizers in the context of traditional variational image registration formulations [26, 53, 54], respectively. Assuming that $v$ is non-stationary, we have

[TABLE]

and

[TABLE]

respectively. Regularization of stationary velocity fields is along the same lines, by simply dropping the time integration in the above equations. Our formulation can also be extended to enforce smoothness in time, as also done in [11].

3 Numerical Methods

In this section, we describe a discretize-then-optimize approach for solving the variational problem (1). We eliminate the hyperbolic PDE constraints (2) and (3) using Lagrangian methods. We then describe the discretization of the objective functional itself. Following [54], we consider a multilevel Gauss–Newton method that allows us to efficiently solve the discrete optimization problem. Finally, we give some details about our implementation as an extension to the FAIR toolbox [54].

Assume, for simplicity, that the domain $\Omega=(0,1)^{d}$ is divided into a regular mesh of $m$ cells of edge length $h=1/m$ along each coordinate direction. We use interpolation and padding of the discrete image data to obtain continuously differentiable and compactly supported functions. We approximate integrals in (1) by a midpoint quadrature rule, which requires evaluating the final state on the cell-centered points, $\mathbf{x}_{c}\in\mathbb{R}^{d\cdot m^{d}}$ , of the grid. Without loss of generality, we assume that the velocity field $v$ is discretized on the same mesh. However, the domain size and number of cells can be varied in practice. The discrete velocity field, denoted by $\mathbf{v}\in\mathbb{R}^{n}$ , is discretized in time at the nodes of a regular grid with $n_{t}$ cells and in space at cell-centered grid points. The total number of unknowns is $n=(n_{t}+1)\cdot d\cdot(m^{d})$ .

3.1 Lagrangian Methods for Hyperbolic PDEs

Lagrangian methods exploit the fact that solutions to hyperbolic PDEs evolve along characteristic curves [46, 24]. These methods are Lagrangian in the sense that the transport of the density is referred to in the moving coordinate system. Lagrangian methods typically consist of two steps: First, the characteristics are computed numerically. Then, in a second step, the final image (or density) is computed. While the characteristic curves are identical for the advection and the continuity equation, the computation of the second part is not.

Step 1 (Computing the Characteristics)

Here, we describe our implementation for computing the characteristic curve passing through a given point. The velocity field $v:\Omega\times[0,1]\to\mathbb{R}^{d}$ is known and represented by the coefficients $\mathbf{v}\in\mathbb{R}^{n}$ . We solve (4) using an RK4 method with $N$ equidistant time steps of size $\Delta t=\pm 1/N$ ; the sign of the time step depends on whether the characteristics are computed forward or backward in time. Note that the number of time steps $N$ for the RK4 method and the number of cells $n_{t}$ in the space-time grid do not necessarily have to be equal. The former is a parameter of the numerical solver and controls the accuracy of the characteristics. The latter is a modeling parameter and ultimately controls the search space for the transformation $y$ .

Our choice of an RK4 scheme is motivated by accuracy considerations for computing the characteristics. For simplicity, we illustrate the concept of integrating $v$ based on a first-order forward Euler scheme. The derivation of our RK4 scheme is along the same lines; it is outlined in Algorithm 1.

Let $\mathbf{x}\in\mathbb{R}^{d\cdot n_{p}}$ be the coordinates of the start (or end) points of the characteristics, e.g., the cell-centers of a regular mesh. Introducing the (time-dependent) transformation $y:\mathbb{R}^{n}\times\mathbb{R}^{d\cdot n_{p}}\times[0,1]^{2}\to\mathbb{R}^{d\cdot n_{p}}$ , and imposing the initial condition $y(\mathbf{v},\mathbf{x},0,0)=\mathbf{x}$ , we compute

[TABLE]

where $t_{k}=k\Delta t$ are the time points and $I$ interpolates the velocity field $\mathbf{v}$ at the transformed points $y$ . In our experiments, we use a bi- or trilinear interpolation model in space and a linear interpolation model in time, applied separately to each component of the velocity field $\mathbf{v}$ . We found by experimentation that using low-order interpolation schemes for the (smoothness regularized) velocity fields is sufficiently accurate for our numerical scheme to yield high-fidelity results in practical applications. We note that the interpolation is a modular component; it can be replaced by (computationally more expensive) higher-order methods. Notice, that the positions at previous time steps, $y(\mathbf{v},\mathbf{x},0,0),\ldots,y(\mathbf{v},\mathbf{x},0,t_{k-1})$ , do not enter the computation in (5); the memory requirements of our numerical scheme are independent of the number of time steps $N$ . The end points of the characteristics are then $y(\mathbf{v},\mathbf{x},0,1)$ , which, if $\mathbf{x}=\mathbf{x}_{c}$ , can also be interpreted and visualized as a deformed regular grid.

Step 2 (Solving the hyperbolic PDE)

While solutions to both the transport and the continuity equations evolve along the same characteristics, the steps for computing the transported quantity at $t=1$ vary. Thus, we discuss both cases separately. An illustration of both schemes can be found in Fig. 3.1.

Transport equation: Considering the transport equation (2), we compute the intensities of the advected image $\mathcal{T}$ on the deformed cell-centered grid $\mathbf{x}_{c}$ by following the characteristics backwards in time. This yields

[TABLE]

In general, $y(\mathbf{v},\mathbf{x}_{c},1,0)$ does not coincide with a grid point; the intensity has to be computed by interpolation. In our numerical experiments we use a bi- or tri-linear interpolation model and regularized cubic approximation methods provided in FAIR to obtain the intensity values of the deformed image; see [54] for implementation details and other common choices.

Continuity equation: To solve the continuity equation (3) we consider the Particle-In-Cell (PIC) method that pushes mass along the characteristics forward in time; see, e.g., [18]. For a given grid point $y_{0}\in\mathbb{R}^{d}$ , we introduce a particle with its mass given by the intensity value $\mathcal{T}(y_{0})$ . Then, we follow the trajectory of the particle along the characteristic to its final point $y_{1}:=y(\mathbf{v},y_{0},0,1)$ . In general, $y_{1}$ does not coincide with a grid point. We obtain the value of the final state in a given cell by integrating the mass of all particles whose support intersects the cell. Equivalently, we compute the final density by splitting the mass of each particle among the cells adjacent to its final location. Ideally, the particles are represented by Dirac delta functions. In practice, we consider bi- or tri-linear hat functions of a certain isotropic width $\delta>0$ as proposed in [18].

Using the push-forward matrix $\mathbf{F}$ , which has also been used in [29], this process can be written as

[TABLE]

In the following, we construct the push-forward matrix in a way that ensures mass-preservation at the discrete level for any choice of $\delta$ and $h$ . For ease of presentation, we describe the procedure for the one dimensional case ( $d=1$ ). The derivation extends to tensor meshes in higher dimensions in a straightforward way under the assumption that the basis functions are piecewise polynomials in the coordinate directions. Let us assume that for $j=1,2,\ldots,m$ particles are located at the cell-centered points $\mathbf{x}_{j}=(j-\frac{1}{2})h$ and their respective mass is given by $\mathbf{u}_{0,j}=\mathcal{T}(\mathbf{x}_{j})$ . The particles are advected to the points $\mathbf{y}_{j}=y(\mathbf{v},\mathbf{x}_{j},0,1)$ . For some $\delta>0$ the particles are represented by the shifted basis functions

[TABLE]

for each $j=1,2,\ldots,n_{p}$ . The mass of $u(\cdot,1)$ contained in the $i$ th interval $[\mathbf{x}_{i},\mathbf{x}_{i+1}]$ is given by

[TABLE]

where $B^{\delta}$ denotes an anti-derivative of $b^{\delta}$ and is given by

[TABLE]

Repeating the process outlined in (8) for all $i=1,2,\ldots,m$ yields the discrete transformed density, which is summarized in (7). Note that our scheme is mass-preserving at the discrete level by design since exact integration is performed. In other words: the columns in $\mathbf{F}$ sum to one regardless of the choices for $\delta$ and $h$ . Also note that $\mathbf{F}$ is sparse. The level of sparsity depends on the ration between $\delta$ and $h$ . If we choose $n_{p}=m$ , $\delta=h$ , and $\mathbf{x}_{j}=h(j+1/2)$ , it is easy to see that $\mathbf{F}$ is the transpose of the linear interpolation matrix [29]. Thus, the relation between (6) and (7) mirrors the adjoint relation between the continuity and the advection equation.

Bibliography71

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] V. Akcelik, G. Biros, and O. Ghattas. Parallel multiscale Gauss-Newton-Krylov methods for inverse wave propagation. In Proc ACM/IEEE Conference on Supercomputing , pages 1–15, 2002.
2[2] A. Alexanderian, N. Petra, G. Stadler, and O. Ghattas. A fast and scalable method for A-optimal design of experiments for infinite-dimensional Bayesian nonlinear inverse problems. SIAM Journal on Scientific Computing , 38(1):A 243–A 272, 2016.
3[3] V. Arsigny, O. Commowick, X. Pennec, and N. Ayache. A log-Euclidean framework for statistics on diffeomorphisms. 9(Pt 1):924–931, 2006.
4[4] J. Ashburner. A fast diffeomorphic image registration algorithm. Neuro Image , 38(1):95–113, Oct. 2007.
5[5] J. Ashburner and K. J. Friston. Diffeomorphic registration using geodesic shooting and Gauss–Newton optimisation. Neuro Image , 55(3):954–967, Apr. 2011.
6[6] B. Avants, P. T. Schoenemann, and J. C. Gee. Lagrangian frame diffeomorphic image registration: Morphometric comparison of human and chimpanzee cortex. Medical Image Analysis , 10:397–412, 2006.
7[7] B. B. Avants, C. L. Epstein, M. Brossman, and J. C. Gee. Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain. Medical Image Analysis , 12(1):26–41, 2008.
8[8] B. B. Avants, N. J. Tustison, G. Song, P. A. Cook, A. Klein, and J. C. Gee. A reproducible evaluation of AN Ts similarity metric performance in brain image registration. Neuro Image , 54:2033–2044, 2011.