A Model for Elastic Evolution on Foliated Shapes

TL;DR

This paper introduces a shape evolution model based on anisotropic elasticity, useful for modeling slow biological shape changes like atrophy, with a focus on layered shapes inspired by cortical structures.

Contribution

It presents a novel elastic evolution framework using optimal control for shape deformation, including a layered shape specialization inspired by cortical organization.

Findings

Preliminary experiments on synthetic layered shapes demonstrate elasticity effects.

The model effectively captures slow biological shape changes like atrophy.

Layered shape decomposition offers new insights into cortical volume deformation.

Abstract

We study a shape evolution framework in which the deformation of shapes from time t to t + dt is governed by a regularized anisotropic elasticity model. More precisely, we assume that at each time shapes are infinitesimally deformed from a stress-free state to an elastic equilibrium as a result of the application of a small force. The configuration of equilibrium then becomes the new resting state for subsequent evolution. The primary motivation of this work is the modeling of slow changes in biological shapes like atrophy, where a body force applied to the volume represents the location and impact of the disease. Our model uses an optimal control viewpoint with the time derivative of force interpreted as a control, deforming a shape gradually from its observed initial state to an observed final state. Furthermore, inspired by the layered organization of cortical volumes, we consider a special case of our model in which shapes can be decomposed into a family of layers (forming a "foliation"). Preliminary experiments on synthetic layered shapes in two and three dimensions are presented to demonstrate the effect of elasticity.

Tables

Table 1. Table 1 : Estimated force parameters compared to their ground truth.

$({\hat{c}}_{\mathrm{true}},h_{\mathrm{true}})$	(1.60, 0.50, 0.30)	(2.20, 0.50, 0.30)	(1.60, 0.65, 0.30)	(2.20, 0.65, 0.30)
$(\hat{c},h)$	(1.60, 0.50, 0.31)	(2.21, 0.49, 0.29)	(1.60, 0.65, 0.32)	(2.21, 0.65, 0.29)
$({\hat{c}}_{\mathrm{true}},h_{\mathrm{true}})$	(1.60, 0.50, 0.35)	(2.20, 0.50, 0.35)	(1.60, 0.65, 0.35)	(2.20, 0.65, 0.35)
$(\hat{c},h)$	(1.60, 0.49, 0.36)	(2.20, 0.50, 0.34)	(1.60, 0.65, 0.36)	(2.21, 0.65, 0.34)

Equations

$$\{\begin{array}{rl}& {\partial }_t\phi (t)=v(t,\phi (t)),\phi (0)=\mathrm{id}\\ & L_{\phi (t)}v(t)=j(t)dp\end{array}.$$

$$U_{\Lambda }(u)-{\int }_{M}F(p{)}^{\top }u(p)dp,$$

$$u=\underset{u^{\prime }\in V}{\mathrm{argmin}}\left(\frac{\delta }{2}\Vert u^{\prime }{\Vert }_V^2+U_{\Lambda }(u^{\prime })-{\int }_{M}F(p{)}^{\top }{u}^{\prime }(p)dp\right),$$

$$\epsilon (u)=\frac{1}{2}(Du+D{u}^{\top }).$$

$$U_{\Lambda }(u)={\int }_M{\Lambda }_p({\epsilon }_p(u))dp,$$

$$U_{\Lambda }(u)=\frac{1}{2}\left({\mathcal{A}}_{\Lambda }u|u\right),$$

$$\frac{1}{2}\left({\mathcal{L}}_{\Lambda }u|u\right)=\frac{\delta }{2}\Vert u{\Vert }_V^2+U_{\Lambda }(u).$$

$$\bar{\Lambda }(\xi )$$

$$+{\lambda }_{\mathrm{ang}}\left({\xi }_{13}^2+{\xi }_{23}^2\right).$$

$$C_{p}^{\prime}=\left[\begin{array}[]{ccc}\\[-5.0pt] T_{1}^{\prime}(p)&T_{2}^{\prime}(p)&\displaystyle\frac{S(p)}{|S(p)|}\\[-5.0pt] \phantom{a}\end{array}\right]=\left[\begin{array}[]{ccc}\\[-5.0pt] T_{1}(p)&T_{2}(p)&\displaystyle\frac{S(p)}{|S(p)|}\\[-5.0pt] \phantom{a}\end{array}\right]\left[\begin{array}[]{cc}\\[-5.0pt] \,G(p)&0\\[5.0pt] \,0&1\end{array}\right],$$

$$\bar{\Lambda }(\xi )$$

$$V_0=\{u\in V:u_{{|}_{{\mathcal{M}}_{\mathrm{bottom}}}}\equiv 0\},$$

$$u=\underset{u^{\prime }\in V_0}{\mathrm{argmin}}\left(\frac{\delta }{2}{\Vert u^{\prime }\Vert }_V^2+\frac{1}{2}\left({\mathcal{A}}_{\Lambda }{u}^{\prime }|u^{\prime }\right)-\left(Fdp|u^{\prime }\right)\right).$$

$$\tilde{\Phi }:[0,1]\times \phi ({\mathcal{M}}_{\mathrm{bottom}})\to \phi (M),(\nu ,\tilde{p})\mapsto \phi (\Phi (\nu ,{\phi }^{-1}(\tilde{p}))).$$

$$g(p;c,h,{\sigma }_{\tan },{\sigma }_{\mathrm{tsv}})=h\exp \left(-\frac{d_{\tan }^2(p,c)}{2{\sigma }_{\tan }^2}-\frac{d_{\mathrm{tsv}}^2(p,c)}{2{\sigma }_{\mathrm{tsv}}^2}\right),$$

$$d_{\tan }(p,c)=\frac{1}{2}\left(\mathrm{length}({\gamma }_{\nu })+\mathrm{length}({\gamma }_{{\nu }^{\prime }})\right)$$

$$d_{\mathrm{tsv}}(p,c)=\frac{1}{2}\left(\mathrm{length}({\alpha }_c)+\mathrm{length}({\alpha }_p)\right).$$

$$J(c,h)=d(\phi (1,M_0),M_1)\text{subject to}\{\begin{array}{rl}& {\partial }_t\phi (t)=v(t,\phi (t))\\ & {\mathcal{L}}_{{\Lambda }^{\phi (t)}}v(t)=j(t)dp\\ & j(t)=(D\phi (t)j_0)\circ \phi (t{)}^{-1}\\ & j_0=\nabla g\\ & c\in {\mathcal{M}}_{\frac{1}{2}},h>0\end{array}.$$

$$d(M,M^{\prime })=\mathrm{volume}(M△M^{\prime })={\int }_{\text{mathbbm}{R}^3}{\left(\text{mathbbm}{1}_M-\text{mathbbm}{1}_{M^{\prime }}\right)}^{2}dx,$$

Full text

1 Introduction

Understanding changes in anatomical shapes is an essential problem for the analysis of many diseases, and more specifically for their tracking over time. While a large variety of shape analysis methods have been successful in cross-sectional problems, in which one exhibits differences in some shape spaces between two populations, a longitudinal analysis of shape changes in a single subject requires a careful modeling of the process, including the structure of the deforming material and the process causing the changes.

Our interest is more specifically in the modeling of dynamical changes in biological tissues that typically occur over long periods of time (such as the progressive atrophy of an organ along the years of disease) thus resulting in sequences of slow alterations in the material properties. Unlike several past approaches for longitudinal analysis based on generic deformation analysis pipelines [16, 7], this work constitutes a step towards a model that could explain morphological variations through more physically interpretable features, while also enabling the estimation of the potential source location and severity of the pathology.

The setting we introduce in this paper is in part inspired from elastic models describing slow changes in biological shapes [18, 13, 6, 1, 19]. We propose a shape evolution framework in which the transition between time $t$ and $t+dt$ is governed by a regularized anisotropic elasticity model in which shapes at rest (or stress free) at time $t$ find a new elastic equilibrium at time $t+dt$ as a result of the application of a small force $dF$, reaching a configuration that becomes their new resting state for subsequent evolution.

The paper is organized as follows. Sect. 2 introduces a general mathematical setting for the shape evolution process, while Sect. 3 details the specific laminar elastic model used in our experiments. Sect. 4 provides some numerical illustrations of the method, focusing, so far, on synthetic data. We conclude the paper in Sect. 5.

2 Shape Evolution Paradigm

We first introduce a general shape evolution model which we shall make specific to layered elastic materials in the next section. Given an initial shape $M_{0}$ which we take as a compact domain of $\mathbbm{R}^{3}$, we consider deformed shapes $t\mapsto\varphi(t,M_{0})$ in which $\varphi$ is a time-dependent diffeomorphism obtained as the flow of a time-dependent velocity field $t\mapsto v(t,M_{0})$ satisfying the system:

[TABLE]

We shall assume that each $v(t)$ belongs to $V$, a reproducing kernel Hilbert space (RKHS) of vector fields on $\mathbbm{R}^{3}$, continuously embedded in the space $C^{2}_{0}(\mathbbm{R}^{3},\mathbbm{R}^{3})$ which denotes the Banach space of $C^{2}$ vector fields on $\mathbbm{R}^{3}$ such that $u$, $Du$, and $D^{2}u$ vanish at infinity, equipped with the norm defined by $\|u\|_{2,\infty}=\|u\|_{\infty}+\|Du\|_{\infty}+\|D^{2}u\|_{\infty}$. ($Du$ denotes the differential of $u$.) Also, in (1), $L_{\varphi}$ is a certain linear operator from $V$ to $V^{*}$, $j(t)$ is a vector field on $\mathbbm{R}^{3}$ supported by $\varphi(t,M_{0})$, and $j(t)\,dp$ is its associated element in $V^{*}$. In the rest of the paper, for any generic function $(t,x)\mapsto f(t,x)$, we adopt the convention of writing $f(t)$ instead of $f(t,\cdot)$.

Our objective consists in estimating $j(t)$ given an observed initial shape $M_{0}$ and target shape $M_{1}$, under the evolution model expressed in (1). In this general setting, such an inverse problem is typically formulated as an optimal control problem by adding either regularization penalties or constraints on $j$ for well-posedness.

In spite of the fact that $\varphi(t)$ is a diffeomorphism obtained as the flow a vector field, the important difference between this approach and standard registration of volumetric data in diffeomorphic frameworks such as [2] is that the evolution here is not controlled directly by the velocity $v$. Instead, $v$ is implicitly determined by $j$ through the operator $L_{\varphi}$. The definition of $L_{\varphi}$ thus considerably affects the nature and properties of optimal solutions. In the next section, we define $L_{\varphi}$ for layered structures based on anisotropic linear elasticity.

3 Elastic Evolution of Layered Shapes

We build the operator $L_{\varphi}$ step by step. First, for a fixed elastic shape with elastic tensor $\Lambda$, we define the operator $\mathcal{L}_{\Lambda}$ such that the solution of $\mathcal{L}_{\Lambda}u=F\,dp$ characterizes, up to a small regularization, the elastic displacement from resting state to equilibrium when a force density $F$ is applied. Next, we define an elastic tensor $\Lambda$ in the special case of layered shapes. We then specify how the elastic tensor is transformed to $\Lambda^{\varphi}$ when a diffeomorphism $\varphi$ acts on layered shapes, which is needed to keep track of the configuration as it gets progressively deformed.

By our assumption that the deformed shapes are at rest at each time, the infinitesimal displacement at each time is given by $\mathcal{L}_{\Lambda^{\varphi}}\delta u=\delta F\,dp$ for some infinitesimal force density $\delta F$. After normalization by $\delta t$, we finally obtain the elastic operator for layered shapes as $L_{\varphi}=\mathcal{L}_{\Lambda^{\varphi}}$ such that $\mathcal{L}_{\Lambda^{\varphi}}v=j\,dp$.

3.1 Linear Elastic Model for Small Deformation

Let $M\subset\mathbbm{R}^{3}$ be a compact domain, which is assumed to represent an elastic material at rest. When applying a small force $F$ to $M$, the displacement field $u$ of the small deformation $\psi=id+u$ can be obtained based on the minimum total potential energy principle, in which $u$ minimizes

[TABLE]

where $U_{\Lambda}(u)$ is the linear elastic energy resulting from the displacement $u$ on $M$, and $\Lambda$ is the elastic tensor describing elastic properties of the material. The exact expression of $U_{\Lambda}(u)$ is detailed below. Here, we assume that only a body force “inside the volume”, i.e., a force density, affects the material, but that no pressure acts on its boundary.

In order to ensure, eventually, a diffeomorphic evolution of the shape when passing to a time-dependent model, we add a regularization term to the total potential energy and look for a regularized response $u$ to the force $F$ given by

[TABLE]

in which $\delta>0$ is a small regularization parameter.

Following [10, 14, 4], we now specify the definition of the elastic energy and characterize the solution of (2). For a displacement field $u$ on $M$, we denote the linear strain tensor by

[TABLE]

Let $\varepsilon_{p}(u)$ denote the evaluation of $\varepsilon(u)$ at a point $p\in M$. The linear elastic energy corresponding to $u$ can then be written as

[TABLE]

where the elastic tensor $p\mapsto\Lambda_{p}$ is a mapping from $M$ to the set of non-negative quadratic forms over the space of $3$ by $3$ symmetric matrices. The function $\Lambda$ encodes the elastic properties of the material. For example, an isotropic linear elastic material has $\Lambda_{p}(\varepsilon)=\frac{\lambda}{2}\left(\sum_{i=1}^{3}\varepsilon_{ii}\right)^{2}+\mu\sum_{i,j=1}^{3}\varepsilon_{ij}^{2}$ for all $p\in M$, where $\lambda$ and $\mu$ are the Lamé parameters.

If we further assume that $p\mapsto\|\Lambda_{p}\|$ is integrable over $M$, then there exists a linear operator $\mathcal{A}_{\Lambda}:C^{1}_{0}(\mathbbm{R}^{3},\mathbbm{R}^{3})\to C^{1}_{0}(\mathbbm{R}^{3},\mathbbm{R}^{3})^{*}$ uniquely defined by

[TABLE]

where ${\left({\alpha}\,\left|\vphantom{\alpha}\,{u}\right.\right)}$ denotes the application of a linear form $\alpha$ to a vector $u$. Since $V\hookrightarrow C^{2}_{0}(\mathbbm{R}^{3},\mathbbm{R}^{3})$, $\mathcal{A}_{\Lambda}$ can also be seen as an operator from $V$ to $V^{*}$. Letting $L:V\to V^{*}$ denote the duality operator on $V$ such that ${\left\langle{v}\,,\,{w}\right\rangle}_{V}={\left({Lv}\,\left|\vphantom{Lv}\,{w}\right.\right)}$, one can define the operator $\mathcal{L}_{\Lambda}:V\to V^{*}$ by $\mathcal{L}_{\Lambda}=\delta L+\mathcal{A}_{\Lambda}$ and then write

[TABLE]

As a result, the solution of (2) is given by $u$ such that $\mathcal{L}_{\Lambda}u=F\,dp$.

3.2 Layered Shapes and Tangential Isotropic Elasticity

Besides being elastic, we assume that $M$ comes with a layered structure as illustrated in Fig. 1 and that the elastic behavior of $M$ is “isotropic tangentially to layers.” This is defined precisely below by specifying the particular form of $\Lambda$ in this configuration. One motivating example for such a construction is that of the cerebral cortex and its organization along cortical layers and columns.

For the compact domain $M\subset\mathbbm{R}^{3}$, we assume that there are two surfaces $\mathcal{M}_{\mathrm{bottom}}$, the bottom layer, $\mathcal{M}_{\mathrm{top}}$, the top layer, included in $\partial M$. Moreover, we assume that the top layer is obtained from a diffeomorphism $\Phi:[0,1]\times\mathcal{M}_{\mathrm{bottom}}\rightarrow M$, such that $\Phi(0)=\mathrm{id}$ and $\Phi(1,\mathcal{M}_{\mathrm{bottom}})=\mathcal{M}_{\mathrm{top}}$. We let $S$ denote the vector field defined by $\partial_{\nu}\Phi(\nu,p)=S(\Phi(\nu,p))$, and we assume that $S$ is continuously differentiable.

Define intermediate layers as $\mathcal{M}_{\nu}=\Phi(\nu,\mathcal{M}_{\mathrm{bottom}})$, so that $\{\mathcal{M}_{\nu}\}_{\nu\in[0,1]}$ forms a foliation of $M$ to which $S$ is transversal. Algorithms building such foliations of volumes have been introduced in brain imaging in order to estimate thickness, see e.g., [11, 5, 17]. Given such a transversal vector field $S$ and layers $\{\mathcal{M}_{\nu}\}_{\nu\in[0,1]}$, we now define an elastic tensor $\Lambda$ which is consistent with this layered structure and “isotropic on layers”. For $p\in M$, let $\{T_{1}(p),T_{2}(p)\}$ be any orthonormal basis of $T_{p}\mathcal{M}_{\nu}$ (the tangent plane to $\mathcal{M}_{\nu}$ at $p$). We define the skew linear change of coordinates at $p$ as the $3$ by $3$ matrix $C_{p}=\left[T_{1}(p)\ T_{2}(p)\ \frac{S(p)}{|S(p)|}\right]$. Now we consider $\Lambda$ of the form $\Lambda_{p}(\varepsilon)=\bar{\Lambda}(C_{p}^{\top}\varepsilon\,C_{p})$, where

[TABLE]

$\bar{\Lambda}$ is independent of the position $p$ if the parameters $\mu_{\mathrm{tan}}$, $\lambda_{\mathrm{tan}}$, $\lambda_{\mathrm{tsv}}$, and $\lambda_{\mathrm{ang}}$ are fixed, which is assumed in the following. Letting $\xi=C_{p}^{\top}\varepsilon_{p}\,C_{p}$, we note that $\xi_{11},\xi_{22}$ measure the stretches along tangential directions $T_{1}(p),T_{2}(p)$ respectively, and $\xi_{33}$ measures the stretch along the transversal direction, which explains the particular form we take for $\bar{\Lambda}$.

Let $\{T_{1}^{\prime}(p),T_{2}^{\prime}(p)\}$ be another orthonormal basis of $T_{p}\mathcal{M}_{\nu}$. Then we have

[TABLE]

where $G(p)$ is a $2$ by $2$ orthogonal matrix. It follows that $\mathrm{tr}(\xi)$, $\mathrm{tr}(\xi^{2})$, and $\xi_{11}^{2}+\xi_{22}^{2}+2\,\xi_{12}^{2}$ are invariant under this transformation, so $\Lambda_{p}$ does not depend on the choice of orthonormal basis of $T_{p}\mathcal{M}_{\nu}$. $\Lambda$ can thus be thought as “isotropic on $T_{p}\mathcal{M}_{\nu}$” for all $p\in M$.

Remark. One can make the same construction for a 2D domain $M$, with a foliation $\{\mathcal{M}_{\nu}\}_{\nu\in[0,1]}$ made of curves, instead of surfaces. One can then describe $T_{p}\mathcal{M}_{\nu}$ by a single unit vector $T(p)$ and define $C_{p}=\left[T(p)\ \frac{S(p)}{|S(p)|}\right]$, a 2 by 2 matrix. In this case, the elastic tensor simplifies to

[TABLE]

3.3 Boundary Condition

As the elastic energy is insensitive to the effect of rigid motions, one usually needs to specify certain boundary conditions in order to ensure that the displacement $u$ resulting from $F$ is well-defined. Here, we impose boundary conditions on $u$ such as, for example, $u=0$ on some subset of $M$ (typically the bottom layer $\mathcal{M}_{\mathrm{bottom}}$). Going back to the case of cortical volumes, this is consistent with the typical assumption that grey matter atrophy is balanced mainly by expansion of the cerebrospinal fluid (CSF) region [12], which in our case corresponds the top layer of $M$. The bottom layer represents the gray/white matter boundary and should remain relatively stable.

In the context of reproducing kernel Hilbert spaces, this boundary condition corresponds to replacing $V$ by its closed subspace

[TABLE]

which contains smooth vector fields that do not move $\mathcal{M}_{\mathrm{bottom}}$. Under the constraint, the displacement response $u$ to a force $F$ is now given by

[TABLE]

Because $V_{0}$ is a closed subspace of $V$, the induced norm makes it into an RKHS. Note that $\mathcal{A}_{\Lambda}$ can then be restricted as an operator from $V_{0}$ to $V_{0}^{*}$, so the same formal analysis applies, leading to an operator from $V_{0}$ to $V_{0}^{*}$ that we still denote by $\mathcal{L}_{\Lambda}$.

3.4 Action of Diffeomorphisms

If the volume $M$ is transformed by a diffeomorphism $\varphi$ to $\varphi(M)$, the deformed layered structure can be specified by

[TABLE]

In particular, the transversal vector field $S$ becomes $\widetilde{S}$ such that $\widetilde{S}(\varphi(p))=D\varphi(p)\,S(p)$. Let $\widetilde{\mathcal{M}}_{\nu}=\varphi(\mathcal{M}_{\nu})$ and $\{\widetilde{T}_{1}(\widetilde{p}),\widetilde{T}_{2}(\widetilde{p})\}$ be an orthonormal basis of $T_{\widetilde{p}}\widetilde{\mathcal{M}}_{\nu}$. One can then define a skew linear transformation $\widetilde{C}=\left[\widetilde{T}_{1},\widetilde{T}_{2},\frac{\widetilde{S}}{|\widetilde{S}|}\right]$ and the transformed elastic tensor $\Lambda^{\varphi}$ given by $\Lambda^{\varphi}_{\widetilde{p}}(\varepsilon)=\bar{\Lambda}(\widetilde{C}^{\top}_{\widetilde{p}}\varepsilon\,\widetilde{C}_{\widetilde{p}})$ with $\bar{\Lambda}$ unchanged. This directly provides operators $\mathcal{A}_{\Lambda^{\varphi}}$ and $\mathcal{L}_{\Lambda^{\varphi}}=L_{\varphi}$.

3.5 Evolution Model for Large Deformation

We now have all required elements in place to operate the large deformation model in (1) for the particular operator $L_{\varphi}$ defined in the previous sections. The connection between the small deformation model and the evolution equations can be made by considering that, at a given time $t$, the displacement $\delta u$ resulting from an infinitesimal force density $\delta F$ is given by $L_{\varphi}\hskip 1.0pt\delta u=\delta F\,dp$. The resulting velocity is $v=\frac{\delta u}{\delta t}$, and letting $j=\frac{\delta F}{\delta t}$ yields $L_{\varphi}\,v=j\,dp$.

The regularization in equation (2) ensures that $v(t)$ is smooth enough such that, under some conditions controlling the size of the control $j(t)$, it generates a flow of diffeomorphisms. This approach is inspired by methods in image and shape registration that use similar regularization, such as the large deformation diffeomorphic metric mapping algorithm [2, 20].

4 Numerical Results

In this section, we showcase several numerical experiments as applications of our model. In all these experiments, we consider a simplified case where $j$ is of the form $j(t,\varphi(t,x))=D\varphi(t,x)\,j_{0}(x)$, that is, $j$ is completely specified by its value $j_{0}$ at time zero. We describe in the next section a low dimensional model for $j_{0}$, specified by a small number of parameters. We then consider the inverse problem of estimating these parameters based on the observation of the original and deformed shapes.

4.1 Force Model

As a simple model of atrophy, we assume that $j_{0}$ is the gradient of a curved Gaussian which is defined by

[TABLE]

where

[TABLE]

and

[TABLE]

As shown in Fig. 2, $\gamma_{\nu}$ and $\gamma_{\nu^{\prime}}$ are geodesics on $\mathcal{M}_{\nu}$ and $\mathcal{M}_{\nu^{\prime}}$ joining $c$, $\widetilde{p}$ and $\widetilde{c}$, $p$ respectively, and $\alpha_{c}$ and $\alpha_{p}$ are integral curves joining $c$, $\widetilde{c}$ and $\widetilde{p}$, $p$ respectively.

4.2 Inverse Problem

We assume that $\sigma_{\mathrm{tan}}$ and $\sigma_{\mathrm{tsv}}$ are known parameters and that the center of Gaussian $c$ is constrained on the middle layer $\mathcal{M}_{\frac{1}{2}}$. This assumption is partially motivated by the case of Alzheimer’s disease in which histological observations suggest that tau protein accumulation preceding the onset occurs in internal cortical layers [3]. Under all these assumptions, the control $j$, and, as a consequence, the entire evolution, is completely determined by the location of $c$ and the height of peak $h$. Given an undeformed shape $M_{0}$ and a target shape $M_{1}$, the inverse problem can now be stated as minimizing

[TABLE]

We use in our experiments

[TABLE]

which computes the volume of non-overlapping region of two shapes. Here, we assume that one does not have access to the layered structure for $M_{1}$. Therefore, we do not include any information of layers in the discrepancy measure $d$.

Since the constraint $c\in\mathcal{M}_{\frac{1}{2}}$ burdens an optimization algorithm, we consider the coordinate representation of $J$ in $c$. In other words, let $(U,\psi)$ be a local chart of $c\in\mathcal{M}_{\frac{1}{2}}$. We define $\widehat{J}:\psi(U)\times\mathbbm{R}_{+}\rightarrow\mathbbm{R}$ by $\widehat{J}(\widehat{c},h)=J(\psi^{-1}(\widehat{c}),h)$. Instead of minimizing the objective function $J$, we minimize equivalently its coordinate representation $\widehat{J}$. By moving the domain of the objective function from the curved space $\mathcal{M}_{\frac{1}{2}}\times\mathbbm{R}_{+}$ to the Euclidean space, we can then utilize a derivative-free optimization algorithm with only box constraints.

Remark. Following the terminology introduced in [9], our model specifies a “deformation module,” that we call “elastic module” in the following. Such modules provide a deformation mechanism (here represented by $j_{0}$) that both drives the shape evolution and is advected by it (see [9] for more details). The free parameters for these modules are the control variables $(c,h)$ with additional geometric parameters $\theta_{f}=(\sigma_{\mathrm{tan}},\sigma_{\mathrm{tsv}})$ for the force and $\theta_{e}=(\delta,\mu_{\mathrm{tan}},\lambda_{\mathrm{tan}},\lambda_{\mathrm{tsv}},\lambda_{\mathrm{ang}})$ for the elastic properties of the volume. One can also relate our construction to that provided in [8], in which our choice for $j_{0}$ provides what is called a “sink” in the referenced paper.

4.3 Discretization

Recall the layered structure $\mathcal{M}_{\nu}=\Phi(\nu,\mathcal{M}_{\mathrm{bottom}})$ where $\Phi$ is a diffeomorphism. We use a discrete set of layers (hence letting $\nu$ be an integer) with consistent triangulations. More precisely, we assume that vertices in $\mathcal{M}_{\nu}$ are $(p_{1}^{\nu},\ldots,p_{N}^{\nu})$ (with $N$ independent of $\nu$) and faces (which are triples or integer in $\{1,\ldots,N\}$) are also independent of $\nu$. We also define transverse edges between $p_{k}^{\nu}$ and $p_{k}^{\nu+1}$ and let $S(p_{k}^{\nu})=p_{k}^{\nu+1}-p_{k}^{\nu}$ represent the discrete version of the transversal vector field (such a representation is provided, for example, by the algorithm introduced in [17]). This provides a decomposition of the volume $M$ into triangular prisms. Without adding vertices, we further decompose the prisms into a consistent tetrahedral mesh using a method similar to that introduced in [15].

For the evaluation of the objective function, we use Dijkstra’s algorithm to compute the length of geodesics needed to define the curved Gaussian map. Discretizing $j\,dp$ as a weighted sum of Dirac functions supported by vertices, one can show that the solution of $\mathcal{L}_{\Lambda^{\varphi}}v=j\,dp$ is uniquely specified by the values of $v$ at vertices. This results in a large linear system (including the constraint that $v=0$ on the bottom layer), which is solved using conjugate gradient.

4.4 Simulated Deformation

Figs. 3(b) to 3(d) display the deformed shapes $\varphi(1,M_{0})$ after applying the same $j_{0}$ to the undeformed shape $M_{0}$, which is shown in Fig. 5(a). The colors of landmarks from blue to red represent the values of curved Gaussian from low to high. The force parameters are $c=(0,0)$, $h=3$, $(\sigma_{\mathrm{tan}},\sigma_{\mathrm{tsv}})=(0.1,0.05)$. We used three different sets of elasticity parameters, that we call “tangent easy,” “transverse easy” and “angle easy.” All three cases use $V$ with a Matérn kernel of order three and width 0.01, $\delta=10^{-6}$, $(\lambda_{\mathrm{tan}},\lambda_{\mathrm{tsv}},\lambda_{\mathrm{ang}})=(1,3,3)$, $(3,1,3)$ and $(3,3,1)$ in the tangent-easy, transverse-easy, and angle-easy cases, respectively.

The different elastic parameters produce different responses to the same $j_{0}$: in Fig. 3(b), the landmarks mainly move tangentially; in Fig. 3(c), they mainly move transversely; while in Fig. 3(d), we see the least movement due to the constraint on the change of lengths between landmarks. Recall that in all cases the bottom layer is kept fixed by our model.

4.5 Simulated Inverse Problem

We now present 2D and 3D simulated inverse problems, in which we first generated deformed targets using curved Gaussian, and then we tried to retrieve the parameters $c$ and $h$. In our experiments, we used surrogateopt in Matlab to find a minimizer.

Two-dimensional results are shown in Figs. 5 and 7. We used a regular shape in Fig. 5, while Fig. 7 was created to resemble a lateral cut of cerebral cortex. The initial shapes are given in Figs. 5(a) and  7(a), and the simulated targets are shown in Figs. 5(b) and  7(b), in which the dashed curves indicate the positions of the initial shapes. We used a very small force in both cases for the purpose of demonstration. Recall that $c$ is restricted on the middle layer, and therefore only has one degree of freedom, which we use the arc length from the leftmost landmark in both cases as its coordinate representation $\widehat{c}$. The level curves of the objective function $\widehat{J}(\widehat{c},h)$ are presented in Figs. 5(c) and 7(c), which demonstrate that the global minimizers are very close to the force parameters generating the targets, even though the deformation is small. Fig. 7(c) also suggests that the optimization problem could be challenging for a highly curved shape. For example, if a gradient descent-like algorithm starts at $(\widehat{c},h)=(1,1)$, it would be difficult to locate the minimizer around $(2.2,0.3)$.

Our 3D experiments use the initial shape shown in Fig. 8, which is a synthetic cerebral cortex. Table 1 compares the estimated parameters $(\widehat{c},h)$ with the ground truth for a series of simulated deformations. Again, we see that the true solution for a simulated target can be accurately retrieved by our method.

5 Discussion

We have presented in this paper a new model describing shape evolution where the control can be interpreted as the derivative of a body force density in the deforming volume. We also have provided a preliminary set of experiments, using simulated data, and based on derivative-free optimization methods.

Current and future work include applying this approach to medical imaging data, and in particular to the study of atrophy due to neurodegenerative diseases. This will require a more general definition of the force field than the single sink model we have used, and an interesting problem will be to maintain parametric identifiability for such models, by injecting suitable biological priors. Using more sophisticated optimization methods will require being able to explicitly compute gradients, which, even though formally feasible, represent serious numerical challenges, probably involving a solution of the linear system $L_{\varphi}v=j\,dp$ on parallel hardware.

We have not discussed in this paper any theoretical result on the well-posedness of the considered problems, such as sufficient conditions on the existence of solutions of (1), or consistency of the discretization schemes. We however included what we believe to be suitable smoothness assumptions regarding, in particular, the RKHS $V$. Such results are under investigation and will be published in the near future.