PDE-constrained LDDMM via geodesic shooting and inexact Gauss-Newton-Krylov optimization using the incremental adjoint Jacobi equations

The class of non-rigid registration methods proposed in the framework of PDE-constrained Large Deformation Diffeomorphic Metric Mapping is a particularly interesting family of physically meaningful diffeomorphic registration methods. Inexact Newton-Krylov optimization has shown an excellent numerical accuracy and an extraordinarily fast convergence rate in this framework. However, the Galerkin representation of the non-stationary velocity fields does not provide proper geodesic paths. In this work, we propose a method for PDE-constrained LDDMM parameterized in the space of initial velocity fields under the EPDiff equation. The derivation of the gradient and the Hessian-vector products are performed on the final velocity field and transported backward using the adjoint and the incremental adjoint Jacobi equations. This way, we avoid the complex dependence on the initial velocity field in the derivations and the computation of the adjoint equation and its incremental counterpart. The proposed method provides geodesics in the framework of PDE-constrained LDDMM, and it shows performance competitive to benchmark PDE-constrained LDDMM and EPDiff-LDDMM methods

Efficient algorithms for geodesic shooting in diffeomorphic image registration

Diffeomorphic image registration is a common problem in medical image analysis. Here, one searches for a diffeomorphic deformation that maps one image (the moving or template image) onto another image (the fixed or reference image). We can formulate the search for such a map as a PDE constrained optimization problem. These types of problems are computationally expensive. This gives rise to the need for efficient algorithms. After introducing the PDE constrained optimization problem, we derive the first and second order optimality conditions. We discretize the problem using a pseudo-spectral discretization in space and consider Heun's method and the semi-Lagrangian method for the time integration of the PDEs that appear in the optimality system. To solve this optimization problem, we consider an L-BFGS and an inexact Gauss-Newton-Krylov method. To reduce the cost of solving the linear system that arises in Newton-type methods, we investigate different preconditioners. They exploit the structure of the Hessian, and use algorithms to efficiently compute an approximation to its inverse. Further, we build the preconditioners on a coarse grid to further reduce computational costs. The different methods are evaluated for two-dimensional image data (real and synthetic). We study the spectrum of the different building blocks that appear in the Hessian. It is demonstrated that low rank preconditioners are able to significantly reduce the number of iterations needed to solve the linear system in Newton-type optimizers. We then compare different optimization methods based on their overall performance. This includes the accuracy and time-to-solution. L-BFGS turns out to be the best method, in terms of runtime, if we solve solving for large gradient tolerances. If we are interested in computing accurate solutions with a small gradient norm, an inexact Gauss-Newton-Krylov optimizer with the regularization term as preconditioner performs best

Combining the Band-Limited Parameterization and Semi-Lagrangian Runge–Kutta Integration for Efficient PDE-Constrained LDDMM

The family of PDE-constrained Large Deformation Diffeomorphic Metric Mapping (LDDMM) methods is emerging as a particularly interesting approach for physically meaningful diffeomorphic transformations. The original combination of Gauss–Newton–Krylov optimization and Runge–Kutta integration shows excellent numerical accuracy and fast convergence rate. However, its most significant limitation is the huge computational complexity, hindering its extensive use in Computational Anatomy applied studies. This limitation has been treated independently by the problem formulation in the space of band-limited vector fields and semi-Lagrangian integration. The purpose of this work is to combine both in three variants of band-limited PDE-constrained LDDMM for further increasing their computational efficiency. The accuracy of the resulting methods is evaluated extensively. For all the variants, the proposed combined approach shows a significant increment of the computational efficiency. In addition, the variant based on the deformation state equation is positioned consistently as the best performing method across all the evaluation frameworks in terms of accuracy and efficiency

CLAIRE: Scalable GPU-Accelerated Algorithms for Diffeomorphic Image Registration in 3D

We present our work on scalable, GPU-accelerated algorithms for diffeomorphic image registration. The associated software package is termed CLAIRE. Image registration is a non-linear inverse problem. It is about computing a spatial mapping from one image of the same object or scene to another. In diffeomorphic image registration, the set of admissible spatial transformations is restricted to maps that are smooth, one-to-one, and have a smooth inverse. We formulate diffeomorphic image registration as a variational problem governed by transport equations. We use an inexact, globalized (Gauss--)Newton--Krylov method for numerical optimization. We consider semi-Lagrangian methods for numerical time integration. Our solver features mixed-precision, hardware-accelerated computational kernels for optimal computational throughput. We use the message-passing interface for distributed-memory parallelism and deploy our code on modern high-performance computing architectures. Our solver allows us to solve clinically relevant problems in under four seconds on a single GPU. It can also be applied to large-scale 3D imaging applications with data that is discretized on meshes with billions of voxels. We demonstrate that our numerical framework yields high-fidelity results in only a few seconds, even if we search for an optimal regularization parameter

An inexact Newton-Krylov algorithm for constrained diffeomorphic image registration

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.Comment: 32 pages; 10 figures; 9 table

Partial Differential Equation-Constrained Diffeomorphic Registration from Sum of Squared Differences to Normalized Cross-Correlation, Normalized Gradient Fields, and Mutual Information: A Unifying Framework; 35632143

This work proposes a unifying framework for extending PDE-constrained Large Deformation Diffeomorphic Metric Mapping (PDE-LDDMM) with the sum of squared differences (SSD) to PDE-LDDMM with different image similarity metrics. We focused on the two best-performing variants of PDE-LDDMM with the spatial and band-limited parameterizations of diffeomorphisms. We derived the equations for gradient-descent and Gauss-Newton-Krylov (GNK) optimization with Normalized Cross-Correlation (NCC), its local version (lNCC), Normalized Gradient Fields (NGFs), and Mutual Information (MI). PDE-LDDMM with GNK was successfully implemented for NCC and lNCC, substantially improving the registration results of SSD. For these metrics, GNK optimization outperformed gradient-descent. However, for NGFs, GNK optimization was not able to overpass the performance of gradient-descent. For MI, GNK optimization involved the product of huge dense matrices, requesting an unaffordable memory load. The extensive evaluation reported the band-limited version of PDE-LDDMM based on the deformation state equation with NCC and lNCC image similarities among the best performing PDE-LDDMM methods. In comparison with benchmark deep learning-based methods, our proposal reached or surpassed the accuracy of the best-performing models. In NIREP16, several configurations of PDE-LDDMM outperformed ANTS-lNCC, the best benchmark method. Although NGFs and MI usually underperformed the other metrics in our evaluation, these metrics showed potentially competitive results in a multimodal deformable experiment. We believe that our proposed image similarity extension over PDE-LDDMM will promote the use of physically meaningful diffeomorphisms in a wide variety of clinical applications depending on deformable image registration