An Inexact Newton--Krylov Algorithm for Constrained Diffeomorphic Image Registration

Abstract

We propose numerical algorithms for solving large deformation diffeomorphic image registration problems. We formulate the nonrigid image registration problem as a problem of optimal control. This leads to an infinite-dimensional partial differential equation (PDE) constrained optimization problem. The PDE constraint consists, in its simplest form, of a hyperbolic transport equation for the evolution of the image intensity. The control variable is the velocity field. Tikhonov regularization on the control ensures well-posedness. We consider standard smoothness regularization based on H1- or H2-seminorms. We augment this regularization scheme with a constraint on the divergence of the velocity field (control variable) rendering the deformation incompressible (Stokes regularization scheme) and thus ensuring that the determinant of the deformation gradient is equal to one, up to the numerical error. We use a Fourier pseudospectral discretization in space and a Chebyshev pseudospectral discretization in time. The latter allows us to reduce the number of unknowns and enables the time-adaptive inversion for nonstationary velocity fields. We use a preconditioned, globalized, matrix-free, inexact Newton–Krylov method for numerical optimization. A parameter continuation is designed to estimate an optimal regularization parameter. Regularity is ensured by controlling the geometric properties of the deformation field. Overall, we arrive at a black-box solver that exploits computational tools that are precisely tailored for solving the optimality system. We study spectral properties of the Hessian, grid convergence, numerical accuracy, computational efficiency, and deformation regularity of our scheme. We compare the designed Newton–Krylov methods with a globalized Picard method (preconditioned gradient descent). We study the influence of a varying number of unknowns in time. The reported results demonstrate excellent numerical accuracy, guaranteed local deformation regularity, and computational efficiency with an optional control on local mass conservation. The Newton–Krylov methods clearly outperform the Picard method if high accuracy of the inversion is required. Our method provides equally good results for stationary and nonstationary velocity fields for two-image registration problems.

Understanding the Connections between Inexact Newton- Krylov PDE-Constrained LDDMM and LDDMM via Geodesic Shooting and Efficient Adjoint Calculation

“…The velocity v is referred as the control variable. The optimization of this problem has been approached using gradient-descent [3] and inexact Newton-Krylov methods [6]. The gradient-descent update equation is obtained from…”

A Semi-Lagrangian Two-Level Preconditioned Newton--Krylov Solver for Constrained Diffeomorphic Image Registration

“…Here, we follow up on our preceding work on constrained diffeomorphic image registration [58, 59]. In diffeomorphic image registration we require that the map y is a diffeomorphism , i.e., y is a bijection, continuously differentiable, and has a continuously differentiable inverse.…”

“…If v is sufficiently smooth it can be guaranteed that the resulting y is a diffeomorphism [10,31,70]. In our formulation, we augment this type of smoothness regularization by constraints on the divergence of v [58, 59]. For instance, for ∇ · v = 0 the flow becomes incompressible.…”

“…Velocity field formulations for diffeomorphic image registration can be distinguished between approaches that invert for a time dependent v [10, 31, 46, 58] and approaches that invert for a stationary v [48, 57, 59]. We invert for a stationary v .…”

“…We solve for a stationary velocity field $\mathrm{bold-italic\nu }\in \mathrm{scriptU}$ and a mass source as follows [59]: $$\underset{m,\mathrm{bold-italic\nu },w}{min}\frac{1}{2}{\Vert m_R-m_1\Vert }_{L^{2}false(\mathrm{normal\Omega }false)}^2+\frac{{\beta }_v}{2}{\Vert \mathrm{bold-italic\nu }\Vert }_{\mathcal{V}}^2+\frac{{\beta }_w}{2}{\Vert w\Vert }_{\mathcal{W}}^2$$subject to $${\partial }_{t}m+\nabla m\cdot \mathrm{bold-italic\nu }=0\text{in}\mathrm{normal\Omega }\times false(0,1false],$$ $$m=m_T\text{in}\mathrm{normal\Omega }\times \left\{0\right\},$$ $$\nabla \cdot \mathrm{bold-italic\nu }=w\text{in}\mathrm{normal\Omega },$$and periodic boundary conditions on ∂ Ω. In our formulation m 1 ( x ) := m ( x , t = 1)—-i.e., the solution of the hyperbolic transport equation (1b) with initial condition (1c)—is equivalent to m T ○ y ; the deformation map y can be computed from v in a post-processing step (see, e.g., [58, 59]). The weights β v > 0, and β w > 0 control the regularity of v .…”