Diffeomorphic Surface Registration with Atrophy Constraints

1. Introduction.

Over the last couple of decades, multiple studies have provided evidence of anatomical differences between control groups and cognitively impaired groups at the population level for a collection of diseases, including schizophrenia, depression, Huntington’s, and dementia [45, 19, 11, 46, 66, 53, 1, 36, 52, 64, 41]. In the particular case of neurodegenerative diseases, a repeated objective has been to design anatomical biomarkers, measurable from imaging data, that would allow for individualized detection and prediction. This goal has become even more relevant with the recent emergence of longitudinal studies involving patients at early stages, or “converters,” who showed that, when the effect is measured at the population level, anatomical changes caused by diseases like Alzheimer’s or Huntington’s were happening several years before cognitive impairment could be detected on individual subjects.

Shape analyses from medical imaging data rely most of the time on a registration step that places the subjects’ anatomies in a common coordinate system, often associated with a template, or average shape, to which all images are coregistered [15, 56, 58, 8, 5, 6, 9, 10, 38, 39, 59, 60]. A large number of registration methods have been proposed in the literature, with an extensive survey proposed in [51], including more than 400 references for image registration alone. Surface registration has also been extensively developed [27, 26, 57, 18, 34, 37, 4, 61, 68, 7, 13, 35, 67, 12].

We focus in this paper on surface registration using the large deformation diffeomorphic metric mapping (LDDMM) algorithm, which has been used extensively to analyze shape variation in regions of interest (ROIs) represented by triangulated surfaces. While some of its main advantages are its flexibility and its ability to render smooth, diffeomorphic, free-form shape changes, there are situations in which prior information is available and should be used to enhance the results of the shape analysis and make it more robust to noise. In this paper, we focus on situations in which no tissue growth is expected to occur. In such contexts, it is natural to ensure that shape analysis should only detect atrophy, even when noise and inaccuracy in the ROI segmentation process may lead in the other direction. (Here, we mean “atrophy” in the general sense of local volume loss.) In this paper, we introduce an atrophy-constrained registration algorithm, which includes some of the ideas introduced in [2], while extending them to inequality constraints associated to the problem we consider. The primary application of the proposed method is for longitudinal data (time series), for which our approach can be seen as a generalization of monotonic nonlinear regression. This adds a new tool to the growing set of methods developed for regression in infinite-dimensional shape spaces, which include, among others, kernel regression, geodesic regression, polynomial regression, and shape spline regression [16, 43, 55, 20, 21, 30, 31, 29, 32, 40].

The atrophy-constrained registration algorithm will be described in section 2, with our numerical approach discussed in section 3. Some theoretical results on existence of solutions and consistency of discrete approximations are provided in section 4. An extension of the algorithm to include affine alignment is provided in section 5. The application of this algorithm to time series is discussed in section 6. Finally, experimental results are provided in section 7.

2. Atrophy-constrained LDDMM.

Although time series, discussed in section 6, will be the primary application of the atrophy-constrained framework, we will start by discussing the approach with the standard two-shape case (“template to target”), since it includes all the new developments in a simpler context.

3. Numerical method.

Problem 3 is solved using augmented Lagrangian methods, introducing, as described in [44], slack variables to complete inequality constraints. More precisely, let

and let C(q) = (C_kl(q), k, l = 1, … , n) be the associated n × 3n matrix. Let K(q) be the 3n × 3n matrix formed from blocks ($K(q_k,q_l){\text{Id}}_{{\mathbb{R}}^3},k,l=1,\dots ,n$). For a vector u, let u⁺ denote the vector formed with the positive parts of each of the coordinates of u. Let ∣N∣ denote the n-dimensional vector with kth coordinates equal to ∣N_k∣.

The augmented Lagrangian is defined by

Here, the parameter μ is a small positive number, $\lambda \in {\mathbb{R}}^n$ is a Lagrange multiplier, and q is considered as a function of α via the state equation ∂_tq = K(q)α. The augmented Lagrangian optimization procedure iterates the following steps (starting with initial values (α₀, λ₀)):

1. Minimize α ↦ F(α, λ_n) to obtain a new value α_n+1 (and q_n+1 via the state evolution equation).

2. Update λ with λ_n+1 = (μλ_n − C(q_n+1)α_n+1 + ε∣N∣)⁺ /μ.

3. If needed (e.g., if ∣(C(q_n+1)α_n+1 − ε∣N∣)⁺∣ is larger than a threshold δ_n), replace μ by a smaller number, ρμ, with ρ < 1.

The most expensive step is, of course, the first one, which requires solving an optimal control problem equivalent in complexity to the unconstrained problem. The computation of ∇_αF(α, λ) (the gradient of F with respect to α) uses a back-propagation algorithm, also called the adjoint method. Since similar computations are described in multiple places [57, 65, 14, 4], we briefly summarize it here.

1. Given α, compute the associated state q via ∂_tq = K(q)α and evaluate the matrix C(q).

2. Introduce a costate $\mathit{p}(t)=(p_1,\dots ,p_n)\in ({\mathbb{R}}^3{)}^n$, t ∈ [0, 1], defined by

   $$\begin{gathered} \mathit{p}(1)={\nabla }_q{D}_{\delta }((q,F_0),({\mathit{q}}_1,F_1){)}_{\mid q=\mathit{q}(1)}\text{and} \\ {\partial }_t\mathit{p}=-{\nabla }_{\mathit{q}}\left({\mathit{p}}^{T}K(\mathit{q})\alpha -\frac{1}{2}{\alpha }^{T}K(\mathit{q})\alpha -\frac{1}{2\mu }\mid (C(\mathit{q})\alpha -\epsilon \mid \mathit{N}\mid -\mu \lambda {)}^{+}{\mid }^2\right). \end{gathered}$$

   (Note that the square of the positive part of a number is differentiable.)

3. Define ${\nabla }_{\alpha }F=K(\mathit{q})(\alpha -\mathit{p})+\frac{1}{\mu }{\left((C(\mathit{q})\alpha -\epsilon \mid \mathit{N}\mid -\mu \lambda {)}^{+}\right)}^{T}C(\mathit{q})$.

Remark. Ensuring that the total volume of the surface decreases (instead of enforcing inward motion at every point) can be done very similarly to the full atrophy constraint, using the single constraint ${\sum }_{k=1}^{n}v(t,q_k(t){)}^T{N}_k(t)\le \epsilon$ or, after reduction,

where N_k is given by (8). It is important here to use the area-weighted normal to discretize the continuous constraint ${\int }_{S(t)}v(t,\cdot {)}^{T}N(t,\cdot )d{\mathcal{A}}_{S(t)}\le \epsilon$, which provides the derivative of the total volume with respect to time. This constraint can be rewritten as $1_n^{T}C(\mathit{q}(t))\alpha (t)\le \epsilon$, where 1_n is the n-dimensional vector with all coordinates equal to 1. The reformulation of the augmented Lagrangian method for this scalar constraint is straightforward and left to the reader.

4. Existence of solutions and convergence.

The following theorems address two important properties of the proposed method. The first one shows that solutions of the constrained problem within our space of interest exist, and the second one considers the issue of consistency of the discretization approach, showing that under suitable assumptions, the solutions obtained with discrete surfaces converge to solutions of the continuous problem. These results are proved in the appendices.

Theorem 1. Assume that V is continuously embedded in $C_0^p({\mathbb{R}}^3,{\mathbb{R}}^3)$ for p ≥ 2 and that the data attachment term D is such that φ ↦ D(φ(S₀), S₁), defined over all C^p-diffeomorphisms φ such that φ − $id\in C_0^p$, is continuous for the uniform C^p convergence over compact sets. Then Problems 1, 2, and 3 always have an optimal solution v ∈ L²(0, 1; V).

The assumption on D applies to the one defined in (3) as soon as p ≥ 1.

We now introduce some assumptions to address the consistency of discrete approximations. We say that a sequence of triangulations ((qⁿ, Fⁿ)) converges nicely to a surface S of class C² if there exists an n ≥ n₀ such that the following conditions are satisfied.

1. The vertices of Sⁿ belong to $S:{\mathit{q}}^n=(q_1^n,\dots ,q_{m_n}^n)$ with $q_k^n\in S$ for k = 1, … , m_n.

2. The triangulations are nested:

   $$\{q_1^{n_1},\dots ,q_{m_{n_1}}^{n_1}\}\subset \{q_1^{n_2},\dots ,q_{m_{n_2}}^{n_2}\}.$$

   In particular, we can order the collection of all vertices of the triangulations in a countable sequence (q₁, q₂, …), and for each q_j there exists an integer n_j such that q_j is a vertex in qⁿ for all n ≥ n_j.

3. The maximum edge length in (qⁿ, Fⁿ) goes to 0.

4. If η(qⁿ, f) is the length of the largest edge of a face f, the ratio

   $$\frac{\text{area}(S({\mathit{q}}^n,f))}{\eta ({\mathit{q}}^n,f{)}^2},$$

   f ∈ Fⁿ, is bounded away from 0 (uniformly in n).

5. The triangulations are close enough to S to ensure that every point x ∈ S(qⁿ, Fⁿ) has a unique closest point ξ(x) on S, such that ξ is onto.

These hypotheses imply the convergence of normals as follows [42].

1. In addition

   $$\underset{f\in F^n,p\in S({\mathit{q}}^n,f)}{\max }\mid \frac{N({\mathit{q}}^n,f)}{\mid N({\mathit{q}}^n,f)\mid }-N(\xi (p))\mid \to 0$$

   when n goes to infinity, where N is the unit normal to S.

We have the following result.

Proposition 1. Assume that (qⁿ, Fⁿ) converges nicely to S and (${\mathit{q}}_1^n$, $F_1^n$) converges nicely to S₁. Then

when n goes to infinity, where D and D_δ are defined in (3), (2), and (9).

We can now state our consistency theorem.

Theorem 2. Fix ε^∞ > 0, and assume that V is continuously embedded in $C_0^p({\mathbb{R}}^3,{\mathbb{R}}^3)$ for p ≥ 2. Let S₀ and S₁ be two C² surfaces, and let ($S_0^n$), ($S_1^n$) be two sequences of triangulations that converge nicely to S₀ and S₁, respectively. Then there exists a decreasing sequence of real numbers εⁿ > 0 with εⁿ → ε^∞, such that if vⁿ ∈ L²(0, 1; V) solves Problem 2 with ε = εⁿ, initial surface $S_0^n$, and target $S_1^n$, then the sequence $(v^n{)}_{n\in \mathbb{N}}$ is bounded in L²(0, 1; V), and any weak limit point v* of vⁿ is a solution to Problem 1 with ε = ε^∞.

5. Affine alignment.

Because the RKHS V is embedded in $C_0^p({\mathbb{R}}^3,{\mathbb{R}}^3)$, vector fields v ∈ V vanish at infinity. This implies, in particular, that affine transformations do not belong to this Hilbert space, and that the diffeomorphic registration does not incorporate any rigid alignment. If such an alignment is needed, one can include it in the optimal control framework by completing the control with the corresponding vector fields.

Let the registration be computed over $G⋉{\mathbb{R}}^3$, where G is a subgroup of $G{L}_3(\mathbb{R})$, ⋉ referring to the semidirect product extending G with translations to obtain affine transformations. Let g be the Lie algebra of G, with basis E₁, … , E_h. Instead of v ∈ V, we use a control given by (v, β₁, … , β_h, τ) with ${\beta }_1,\dots ,{\beta }_h\in \mathbb{R}$ and $\tau \in {\mathbb{R}}^3$, and the state equation

with associated cost $\frac{1}{2}{\int }_0^1(\Vert v(t){\Vert }_V^2+{\sum }_{k=1}^h{c}_k{\beta }_k(t{)}^2+c_0\mid \tau (t){\mid }^2)dt$ for some nonnegative numbers c₀, c₁, … , c_h. The extension of the numerical procedure to this setting is rather straightforward and not detailed here. It is also natural to include non-Euclidean affine components in the atrophy constraint v · N ≤ ε∣N∣.

The Euclidean case (rigid correction) is, however, the most interesting for application. This corresponds to G = SO₃, the rotation group (so that h = 3). We will assume that Euclidean transformations act as isometries on V, which means that, for all v ∈ V, R ∈ SO₃, $b\in {\mathbb{R}}^3$, the vector field $\tilde{v}:x\mapsto R^{T}v(Rx+b)$ also belongs to V and $\Vert \tilde{v}{\Vert }_V=\Vert v{\Vert }_V$. Euclidean-invariant RKHSs of vector fields are extensively described in [22], to which we refer the reader for more details. In the case of scalar kernels $\mathbb{K}(x,y)=K(x,y){\text{Id}}_{{\mathbb{R}}^3}$, Euclidean invariance implies that K is a radial kernel, i.e., that K(x, y) = γ(∣x − y∣²) for some function γ (additional conditions on γ are needed to ensure that K is a positive kernel; see [22]). Assume finally that the endpoint cost D is invariant by the Euclidean transformation: D(S, S′) = D(RS + b, RS′ + b). In this case, the optimal control problem (using c₀ = ⋯ = c₃ = 0) that minimizes $\frac{1}{2}{\int }_0^1\Vert v(t){\Vert }_V^{2}dt+D(S(1),S_1)$ in v, β, τ, subject to

is equivalent to minimizing $\frac{1}{2}{\int }_0^1\Vert \tilde{v}(t){\Vert }_V^{2}dt+D(\tilde{S}(1),R_1{S}_1+b_1)$ in $\tilde{v}$, R₁, b₁, subject to

via the change of variable $q(t)=R(t)\tilde{q}(t)+b(t)$, $v(t,x)=R(t{)}^{-1}\tilde{v}(t,Rx+b)$, $N(t)=R(t)\tilde{N}(t)$, with ${\partial }_{t}R={\sum }_{l=1}^3{\beta }_l(t)E_{l}R$, ${\partial }_{t}b={\sum }_{l=1}^3{\beta }_l(t)E_{l}b(t)+\tau (t)$, R₁ = R(1)⁻¹, and b₁ = −R(1)⁻¹b(1). In other words, Euclidean registration via optimal control using (16) is equivalent to the original atrophy-constrained LDDMM optimizing its target over an orbit under the action of rigid transformations.

Note that c₀ = ⋯ = c_h = 0 should not be used with general affine transformations when noncompact components of the affine group are added to rotations. Intuitively, this would allow small deformations to be scaled up to larger ones at zero cost, and one can conjecture that the associated optimal control problem has no solutions and admits minimizing sequences of controls with vanishing cost at the limit. The equivalence with a problem in which the target is allowed to vary over its orbit via affine transformations is not true either for groups larger than the Euclidean one, essentially because invariant kernels do not exist in such cases.

Some attention should be paid to the discretization in time t, in particular in the rigid case. In our implementation, we use a simple Euler scheme to discretize the equation ∂_tq = v(t, q); i.e., we take q(t + δt) = q(t) + δt v(t, q(t)). When using rigid registration, however, we take, with $A(t)={\sum }_{l=1}^3{\beta }_l(t)E_l$ a skew-symmetric matrix,

to discretize ∂_tq = v(t, q) + Aq + τ, with the explicit expression $e^U=\text{Id}+\frac{\sin c_U}{c_U}U+\frac{1-\cos c_U}{c_U^2}{U}^2$, $c_u=\sqrt{-\text{tr}(U^2)}$ for a 3 × 3 skew-symmetric matrix U. This ensures that the rigid registration remains a displacement, even for large values of the coefficients β_l, which are made possible by the absence of cost on this transformation.

6. Application to time series.

Longitudinal studies generally include more than one followup scan, and we now extend our algorithm to the case of multiple targets observed at successive times. We still have a baseline surface S₀ and now assume p followup surfaces S₁, … , S_p. A direct generalization of the previous approach replaces (5) in Problem 1 by

where 0 < t₁ ≤ ⋯ ≤ t_n = 1 are attachment times, which can be fixed a priori (e.g., t_k = k/n) or also optimized. The expression in (18) should be minimized, starting from a baseline S₀, with the same constraints as in the single-target case.

As pointed out in [48, 49, 33, 3], this formulation of longitudinal registration (with or without atrophy) gives a special role to the baseline image and creates a bias in the way the information is obtained by breaking the symmetry within the acquisition and segmentation protocols that provided the images or surfaces. To address this, we add an extra error term, D(S(t₀), S₀) with t₀ = 0, making the sum start at k = 0 in (18). The trajectory baseline, S(0), becomes a new variable to estimate. To ensure that the topology of the estimated trajectory remains consistent with the data (as was originally enforced by S(0) = S₀), we assume that S(0) is a diffeomorphic transformation of another surface (a template) of known topology. Letting ${\overline{S}}_0$ denote the template, this leads to the minimization of

where λ is a regularization parameter and $\overline{v}$ controls a smooth deformation between the template ${\overline{S}}_0$ and the initial point of the trajectory, S(0). More precisely, we let $S(0)=\overline{S}(1)$, where $\overline{S}(t)$ is a template-to-baseline trajectory satisfying $\overline{S}(0)={\overline{S}}_0$ and ${\partial }_t\overline{S}=\overline{v}(t,\overline{S}(t))$. The cost function is then minimized in the pair of chained controls ($\overline{v}$, v). While atrophy constraints remain applicable to v, no such constraint is enforced on $\overline{v}$. This algorithm has been implemented in Python, and examples of trajectories are provided in the next section.

7. Experimental results.

We now present experiments that compare time series, as described in section 6. We start with simulated sequences, which we use to analyze the performance of the method in handling noise, as well as possible bias that could occur from the asymmetric constraint (atrophy vs. expansion). We will then describe real data results using the BIOCARD dataset.

8. Discussion.

In this paper, we have introduced a new approach targeting atrophy in the registration of longitudinal data, in view of applications to the analysis of neurodegenerative diseases. We have formulated the problem via constrained optimal control and proposed an algorithm that combined the adjoint method in optimal control with augmented Lagrangian methods. We also provided theoretical results on the existence of solutions for the problems, and on the consistency of our numerical scheme for approaching the solution.