  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 $H^1$- or
$H^2$-seminorms. We augment this regularization scheme with a constraint on the
divergence of the velocity field rendering the deformation incompressible 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. 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. 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 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.
