CLAIRE: A Distributed-Memory Solver for Constrained Large Deformation Diffeomorphic Image Registration

Abstract

We introduce CLAIRE, a distributed-memory algorithm and software for solving constrained large deformation diffeomorphic image registration problems in three dimensions. We invert for a stationary velocity field that parameterizes the deformation map. Our solver is based on a globalized, preconditioned, inexact reduced space Gauss-Newton-Krylov scheme.We exploit state-of-the-art techniques in scientific computing to develop an effective solver that scales to thousand of distributed memory nodes on high-end clusters. Our improved, parallel implementation features parameter-, scale-, and grid-continuation schemes to speedup the computations and reduce the likelihood to get trapped in local minima. We also implement an improved preconditioner for the reduced space Hessian to speedup the convergence.We test registration performance on synthetic and real data. We demonstrate registration accuracy on 16 neuroimaging datasets. We compare the performance of our scheme against different flavors of the DEMONS algorithm for diffeomorphic image registration. We study convergence of our preconditioner and our overall algorithm. We report scalability results on state-of-the-art supercomputing platforms. We demonstrate that we can solve registration problems for clinically relevant data sizes in two to four minutes on a standard compute node with 20 cores, attaining excellent data fidelity. With the present work we achieve a speedup of (on average) 5× with a peak performance of up to 17× compared to our former work. Introduction.Deformable registration is a key technology in the medical imaging. It is about computing a map y that establishes a meaningful spatial correspondence between two (or more) images m R (the reference (fixed) image) and m T (the template (deformable or moving) image; image to be registered) of the same scene [43,96]. Numerous approaches for formulating and solving image registration problems have appeared in the past; we refer to [43,60,96,97,116] for lucid overviews. Image registration is typically formulated as a variational optimization problem that consists of a data fidelity term and a Tikhonov regularization functional to over-come ill-posedness [40,43]. A key concern is that y is a diffeomorphism, i.e., the map y is differentiable, a bijection, and has a differentiable inverse. We require that the determinant of the deformation gradient det ∇y does not vanish or change sign. An intuitive approach to safeguard against nondiffeomorphic maps y is to add hard and/or soft constraints on det ∇y to the variational optimization problem [29,58,103,108]. An alternative strategy to ensure regularity is to introduce a pseudo-time variable t and invert for a smooth velocity field v that encodes the map y [16,37,94,126]; existence of a diffeomorphism y can be guaranteed if v is adequately smooth [16,30,37,121]. This model is, as opposed to many traditional formulations that directly invert for y [58,88,95,108], adequate for recovering large, highly nonlinear deformations. Our approach falls into this category. T...

Incompressible Image Registration Using Divergence-Conforming B-Splines

“…There are two approaches to discretise the constraint: 1) relaxing the constraint into an additive soft constraint [1,3,5,9]. 2) using a specific parameterisation for the transformation [2,6,7]. Soft constraints introduce hyperparameters that are difficult to tune reliably and cannot guarantee an incompressible transformation.…”

“…Yet, this approach only provides a first order approximation of an incompressible deformation, still requiring a soft constraint. Recently [6,7] proposed to parameterise the transformation by a divergence-free stationary velocity field (SVF). The transformation is obtained via the Lie exponential mapping exp : v → Φ that maps any divergence-free SVF v into an incompressible transformation Φ [6,7].…”

“…Recently [6,7] proposed to parameterise the transformation by a divergence-free stationary velocity field (SVF). The transformation is obtained via the Lie exponential mapping exp : v → Φ that maps any divergence-free SVF v into an incompressible transformation Φ [6,7]. This reduces the non-convex constraint in (1) to a linear constraint in the continuous domain.…”

“…This reduces the non-convex constraint in (1) to a linear constraint in the continuous domain. However, in [6,7] the SVF is parameterised as linear B-splines, and the linear constraint is discretised by imposing the constraint only on the points of the deformation grid. As a result, guarantees to obtain a continuous incompressible transformation do not hold anymore.…”

“…As a result, guarantees to obtain a continuous incompressible transformation do not hold anymore. To mitigate this issue, [6] proposed to work with images of higher resolution with a finer grid for the linear B-spline leading to the need for distributed super-computing. Moreover, [6,7] methods are limited to using sum of squared differences (SSD) as image similarity metric, which limits their applicability to images with similar intensity distributions.…”