PDE-constrained optimization in medical image analysis

Abstract

PDE-constrained optimization problems find many applications in medical image analysis, for example, neuroimaging, cardiovascular imaging, and oncologic imaging. We review the related literature and give examples of the formulation, discretization, and numerical solution of PDE-constrained optimization problems for medical imaging. We discuss three examples. The first is image registration, the second is data assimilation for brain tumor patients, and the third is data assimilation in cardiovascular imaging. The image registration problem is a classical task in medical image analysis and seeks to find pointwise correspondences between two or more images. Data assimilation problems use a PDE-constrained formulation to link a biophysical model to patient-specific data obtained from medical images. The associated optimality systems turn out to be sets of nonlinear, multicomponent PDEs that are challenging to solve in an efficient way. The ultimate goal of our work is the design of inversion methods that integrate complementary data, and rigorously follow mathematical and physical principles, in an attempt to support clinical decision making. This requires reliable, high-fidelity algorithms with a short time-to-solution. This task is complicated by model and data uncertainties, and by the fact that PDE-constrained optimization problems are ill-posed in nature, and in general yield high-dimensional, severely ill-conditioned systems after discretization. These features make regularization, effective preconditioners, and iterative solvers that, in many cases, have to be implemented on distributed-memory architectures to be practical, a prerequisite. We showcase state-of-the-art techniques in scientific computing to tackle these challenges.

Appendices

Appendix 1: Computing $\det \nabla y$

We do not differentiate the deformation map y, but transport $\psi \mathrel {\mathop :}=\det \nabla \phi$ where $y^{-1}(x) = \phi (x,t=1)$. That is, we solve

$$\begin{aligned} \partial _t \psi + v\cdot \nabla \psi&= \psi \nabla \cdot v&\quad \text {in}\;\;{\varOmega }\times (0,1] \end{aligned}$$

(13a)

$$\begin{aligned} \psi&= 1&\quad \text {in}\;\;{\varOmega }\times \{0\} \end{aligned}$$

(13b)

with periodic boundary conditions on $\partial {\varOmega }$.

Remark 6

If we start from an identity map ${\text {id}}$ (zero velocity field) the determinant of the deformation gradient will have a value of one. Hence, we require that the determinant of the deformation gradient remains strictly positive for all $x \in {\varOmega }$ for the computed deformation map to be locally diffeomorphic. Strictly speaking, we also require that the map y is smooth, and has a smooth inverse for y to be a diffeomorphism. If we use (slightly more than) $H^2$-regularity for v we meet these requirements from a theoretical perspective (assuming that the images are smooth) (Beg et al. 2005; Vialard et al. 2012). However, we note that ensuring that the discrete deformation map y (which does not appear in our formulation) is a diffeomorphism, requires sophisticated discretization schemes (see, e.g., Burger et al. 2013). In our formulation, we transport the data. Since the deformation map y does not appear explicitly, we decided to compute the determinant of the deformation gradient by computing the solution of the transport problem in Appendix 1. We found that this is more stable than computing $\det \nabla y$ from y using our numerical scheme.

Appendix 2: Gauss–Newton Krylov algorithm

Here, we showcase the individual steps necessary for solving (1) in Algorithm 3 (Newton iterations; outer iterations) and Algorithm 4 (solution of the system in (12), inner iterations). We refer to Mang and Biros (2015, 2017), Gholami et al. (2016) for details on the particular implementation of our nonlinear solver. We rely on PETSc’s TAO package (Balay et al. 2016; Munson et al. 2017) for our distributed-memory solver (Mang et al. 2016; Gholami et al. 2017) (see Sect. 5.5 for details). The steps are the same as outlined in Algorithms 3 and 4, respectively.

Appendix 3: Newton step

Here, we describe the Newton step for the two problems described in Sects. 2.2 and 3.2. The Newton step is obtained by computing the second variation of (6). We will see that the structure of the PDE operators that appear in the reduced space Hessian is very similar to the first-order optimality systems in Sect. 5.1.

1.1 Diffeomorphic registration

The nonlinearity, ill-posedness, and the dimensionality of the search space makes the solution of this variational problem computational challenging. Most available implementations use first-order iterative schemes to compute a minimizer to the variational optimization problem (see, e.g., Avants et al. 2011; Beg et al. 2005; Chen and Lorenz 2011; Ha et al. 2010; Hart et al. 2009; Vialard et al. 2012). To increase the rate of convergence, they typically do not use $g_v$ in their iterative scheme, but the gradient in the Sobolev space induced by the regularization operator ${\mathcal {A}}$, i.e.,

$$\begin{aligned} \tilde{g}_v = v + (\beta {\mathcal {A}})^{-1}\int _0^1\lambda \nabla m{\mathrm {d}}t. \end{aligned}$$

First-order methods are usually inefficient; they require a lot of iterations to converge to a specific tolerance; we have analyzed this in Mang and Biros (2015). An alternative is to apply variants of Newton’s method to the KKT system. To derive the operators for applying Newton’s method we have to compute the second variations of (6). We arrive at the incremental state and adjoint equation

$$\begin{aligned} \partial _t \tilde{m} + \nabla \tilde{m} \cdot v + \nabla m \cdot \tilde{v}&= 0&{\mathrm{in}}\;\; {\varOmega } \times (0,1], \end{aligned}$$

(14a)

$$\begin{aligned} \tilde{m}&= 0&\quad {\mathrm{in}}\;\; {\varOmega } \times \{0\}, \end{aligned}$$

(14b)

$$\begin{aligned} -\partial _t \tilde{\lambda } - \nabla \cdot (v\tilde{\lambda } + \tilde{v}\lambda )&= 0&\quad \text {in}\;\; {\varOmega } \times [0,1), \end{aligned}$$

(14c)

$$\begin{aligned} \tilde{\lambda }&= - \tilde{m}&\quad {\mathrm{in}}\;\; {\varOmega } \times \{1\}, \end{aligned}$$

(14d)

and the expression for the reduced space Hessian matvec

$$\begin{aligned} {\mathcal {H}}(v)\tilde{v} \mathrel {\mathop :}=\beta {\mathcal {A}}[\tilde{v}] + {\mathcal {K}}[\int _0^1\tilde{\lambda }\nabla m + \lambda \nabla \tilde{m}{\mathrm {d}}t\;]&\quad\text {in}\,\, {\varOmega }, \end{aligned}$$

(14e)

with periodic boundary conditions on $\partial {\varOmega }$. The incremental control variable $\tilde{v}$ corresponds to the search direction of our scheme. The transport equations for the incremental state and adjoint fields $\tilde{m}$ and $\tilde{\lambda }$ are hidden in the integro-differential operator in (14e); we have to solve (14a) and (14c) every time we apply ${\mathcal {H}}$ to a new vector $\tilde{v}$. We use a Gauss–Newton approximation to the true Hessian to ensure that the operator is positive definite far away from the optimum. This corresponds to neglecting all terms that involve $\lambda$ in (14c) and (14e).

1.2 Data assimilation

The expressions for evaluating the Hessian operator are derived by computing the second variations for (9). The Hessian matvec is given by ${\mathcal {H}}\tilde{p} \mathrel {\mathop :}=\beta \tilde{p} - {\varPhi }^\mathsf {T}\tilde{\lambda }_0$, ${\mathcal {H}} : {\mathbf {R}}^{n_p} \rightarrow {\mathbf {R}}^{n_p}$. To be able to evaluate this expression we have to solve the incremental state and adjoint equations given by

$$\begin{aligned} \partial _t \tilde{m} - \nabla \cdot k \nabla \tilde{m} - \rho (1-2m)\tilde{m}&= 0&\quad {\mathrm{in}}\;{\varOmega }_B\times (0,1], \end{aligned}$$

(15a)

$$\begin{aligned} \tilde{m}&={\varPhi }\tilde{p}&\quad {\mathrm{in}}\;{\varOmega }_B\times \{0\}, \end{aligned}$$

(15b)

$$\begin{aligned} -\partial _t \tilde{\lambda } - \nabla \cdot k \nabla \tilde{\lambda } - \rho (1-2m - 2\lambda )\tilde{\lambda }&= 0&\quad {\mathrm{in}}\;{\varOmega }_B\times [0,1), \end{aligned}$$

(15c)

$$\begin{aligned} \tilde{\lambda }&= -{\mathcal {Q}}^\mathsf {T}{\mathcal {Q}} \tilde{m}&\quad {\mathrm{in}}\;{\varOmega }_B\times \{1\}, \end{aligned}$$

(15d)

with Neumann boundary conditions on $\partial {\varOmega }_B$. For the Gauss–Newton approximation to the true Hessian the terms that involve $\lambda$ in (15c) need to be dropped.