Selective metamorphosis for growth modelling with applications to landmarks

TL;DR

This paper introduces a shape matching framework in computational anatomy that allows localized control of deformations, enabling growth modeling in specific image regions with applications to landmarks.

Contribution

It develops a mathematical model for localized template deformation in shape matching, incorporating control over growth regions, and analyzes its well-posedness and geodesic equations.

Findings

Preliminary numerical results demonstrate the framework's potential.

The model allows specifying or learning growth locations in images.

Analysis confirms the mathematical well-posedness of the approach.

Abstract

We present a framework for shape matching in computational anatomy allowing users control of the degree to which the matching is diffeomorphic. This control is given as a function defined over the image and parameterises the template deformation. By modelling localised template deformation we have a mathematical description of growth only in specified parts of an image. The location can either be specified from prior knowledge of the growth location or learned from data. For simplicity, we consider landmark matching and infer the distribution of a finite dimensional parameterisation of the control via Markov chain Monte Carlo. Preliminary numerical results are shown and future paths of investigation are laid out. Well-posedness of this new problem is studied together with an analysis of the associated geodesic equations.

Equations

$$\dot{{\phi }_t}=u_t\circ {\phi }_t,{\phi }_0=\text{id}.$$

$$S(u)={\int }_0^1\frac{1}{2}\Vert u_t{\Vert }_V^2\text{d}t+\frac{1}{2{\lambda }^2}F(I_0\circ {\phi }_1^{-1},I_1)⟶\text{min. subject to eqrefdiffeo},$$

$$\frac{\delta l}{\delta u}=\mathbf{m}(x)=\sum _{i=1}^M{p}_t^i\delta (x-q_t^i),$$

$$u_t(x)=\sum _{i=1}^M{p}_t^{i}K(x-q_t^i),$$

$${\mathbf{q}}_t={\phi }_t{\mathbf{\eta }}_t.$$

$$\mathbf{z}={\phi }_t\dot{\mathbf{\eta }}$$

$$\begin{array}{rl}{S}_m({\mathbf{q}}_t,{\mathbf{p}}_t,{\mathbf{z}}_t)=& {\int }_0^1\frac{1}{2}\left(\Vert u_t{\Vert }_V^2+\frac{1}{{\sigma }^2}\sum _{i=1}^M|z_t^i{|}^2\right)\text{d}t,\end{array}$$

$$\dot{{\mathbf{q}}_t}=u_t({\mathbf{q}}_t)+{\mathbf{z}}_t,$$

$$\mathbf{m}(x)=\frac{1}{{\sigma }^2}\sum _{i=1}^M{z}_t^i\delta (x-q_i)\Rightarrow \mathbf{z}={\sigma }^2{\mathbf{p}}_i,$$

$$\begin{array}{rl}\dot{{\mathbf{p}}_t}& =-\nabla u_t({\mathbf{q}}_t{)}^T{\mathbf{p}}_t\\ \dot{{\mathbf{q}}_t}& =u_t({\mathbf{q}}_t)+{\sigma }^2{\mathbf{p}}_t,\end{array}$$

$$\begin{array}{rl}{S}_{sm}^{\nu }({\mathbf{q}}_t,u_t,{\mathbf{z}}_t)=& {\int }_0^1\frac{1}{2}\left(\Vert u_t{\Vert }_V^2+\sum _{i=1}^M\frac{1}{\nu (q_t^i)}|z_t^i{|}^2\right)\text{d}t,\end{array}$$

$$\mathbf{m}(x)=\sum _{i=1}^M\frac{1}{\nu (q_t^i)}{z}_t^i\delta (x-q_i)\Rightarrow z_t^i=\nu (q_t^i)p_t^i\forall i,$$

$$\begin{array}{rl}{S}_{sm}^{\nu }({\mathbf{q}}_t,u_t,{\mathbf{p}}_t)=& {\int }_0^1\frac{1}{2}\left(\Vert u_t{\Vert }_V^2+\sum _{i=1}^M\nu (q_t^i)|p_t^i{|}^2\right)\text{d}t.\end{array}$$

$$\begin{array}{rl}{\dot{\mathbf{p}}}_t& =-\nabla u_t({\mathbf{q}}_t{)}^T{\mathbf{p}}_t-\frac{1}{2}\nabla \nu ({\mathbf{q}}_t)|{\mathbf{p}}_t{|}^2\\ {\dot{\mathbf{q}}}_t& =u_t({\mathbf{q}}_t)+\nu ({\mathbf{q}}_t){\mathbf{p}}_t,\end{array}$$

$$h(\mathbf{q},\mathbf{p})=h_l({\mathbf{q}}_t,{\mathbf{p}}_t)+\frac{1}{2}\sum _i\nu (q_i)|p_i{|}^2,$$

$$h_l(\mathbf{q},\mathbf{p})=\frac{1}{2}\sum _{i,j=1}^{M}K(q_i-q_j)p_i\cdot p_j.$$

$$\underset{\begin{array}{c}u\in L^2([0,1],V)\\ \mathbf{q}\in H^1([0,1],{\mathbb{R}}^{d\times M})\\ \mathbf{w}\in L^2([0,1],{\mathbb{R}}^{d\times M})\end{array}}{\inf }$$

$$\dot{{\mathbf{q}}_t^i}=u_t({\mathbf{q}}_t)+\sqrt{\nu ({\mathbf{q}}_t)}{\mathbf{w}}_t$$

$${\mathbf{q}}_0,{\mathbf{q}}_1\text{ fixed}$$

$$\displaystyle\Big{(}\sum_{i=1}^{M}|w_{t}^{i}|^{2}\Big{)}^{\frac{1}{2}}$$

$$⟨\sqrt{\nu ({\mathbf{q}}_t)}{\mathbf{w}}_t-\sqrt{\nu ({\mathbf{q}}_t^n)}{\mathbf{w}}_t^n,\varphi ⟩\lesssim {\nu }_{\text{inf}}⟨{\mathbf{w}}_t-{\mathbf{w}}_t^n,\varphi ⟩\to 0,\forall \varphi \in L^2([0,1],{\mathbb{R}}^{d\times M}).$$

$$⟨u_t({\mathbf{q}}_t)-u_t^n({\mathbf{q}}_t^n),\varphi ⟩$$

$$⟨u_t^n({\mathbf{q}}_t)-u_t^n({\mathbf{q}}_t^n),\varphi ⟩\le \text{Lip}(u_t^n)⟨{\mathbf{q}}_t-{\mathbf{q}}_t^n,\varphi ⟩\to 0$$

$$\Vert \nabla v(x)-\nabla w(y){\Vert }_2\lesssim \Vert v{\Vert }_V\Vert x-y{\Vert }_2+\Vert v-w{\Vert }_V\Vert y{\Vert }_2$$

$$\begin{array}{rl}& (q,p)\mapsto \nu (q)p\\ & (q,p)\mapsto \nabla \nu (q)|p{|}^2\end{array}$$

$${\nu }_h(x)=\sum _{k=1}^K{e}^{-{\sigma }_{\nu }^{-2}\Vert h_k-x{\Vert }_{{\mathbb{R}}^d}^2}.$$

Full text

1 Introduction

In computational anatomy [11, 12] one of the most fundamental problems is to continuously deform an image or shape into another and thereby obtain a natural notion of distance between them as the energy required for such a deformation. Common methods to compute image deformations are based on diffeomorphic deformations which assume that the images are continuously deformed into one another with the additional property that the inverse deformation is also continuous. This is a strong requirement for images which implies that the ’mass’ of any part of the image is conserved: we cannot create or close ’holes’. This is also a crucial property in fluid mechanics and in fact the theory of diffeomorphic matching carrying the moniker Large Deformation Diffeomorphic Metric Mapping (LDDMM) [24, 5] has been inspired by fluid mechanics. Indeed, Arnold [4] made the central observation that the geodesic equations for the diffeomorphism group induced by divergence-free vector fields corresponded to that of incompressible flows. If a strictly diffeomorphic matching is not possible or necessary, an extension of LDDMM called metamorphosis [26, 15] is available which introduces a parameter $\sigma^{2}$ parameterising the deviation from diffeomorphic matching allowing for topological variations e.g. growth via image intensity. In particular, if $\sigma^{2}=0$ the deformation is purely diffeomorphic as in LDDMM. See [23, 25, 19] for technical details pertaining to the construction of the metamorphosis problem. While diffeomorphic paths always exist for landmark problems [13] this theory allows one to match images of shapes with different topological features, which is ill-conditioned for standard LDDMM. Indeed, even inexact matching in LDDMM for such problems yields large energies and spurious geodesics that do not contribute to an intuitive matching, see figure 1. As observed here, introducing $\sigma^{2}>0$ regularises the problem and qualitative improves the matching.

In this work, we modify metamorphosis to include a spatially dependent control parameter $x\mapsto\nu(x)$ in order to selectively allow non-diffeomorphic (metamorphic) matching in parts of the image. For $\nu(\cdot)=\sigma^{2}$ our theory recovers the standard metamorphosis model. However, with a localised control (e.g. a Gaussian centred at a point in $\mathbb{R}^{d}$), we can selectively introduce metamorphosis in an image and model local topological effects such as growth phenomena. The difficulty of this problem is to infer the function $\nu(\cdot)$ without prior knowledge of the location of the topological effects. This problem is similar to the one treated in [3], where such functions were parameterising the randomness in LDDMM matching of shapes. We will use a Markov chain Monte Carlo (MCMC) approach to infer appropriate functions $\nu(\cdot)$, such that the topological effects are well described and a large part of the matching remains diffeomorphic. In this paper, we focus on landmark matching but aim to extend the theory of selective metamorphosis to data structures amenable to classical metamorphosis or LDDMM theory can handle.

Structure This paper is organised as follows. We review the theory of classical metamorphosis in section 2 and extend it to selective metamorphosis in section 3. We then introduce a Bayesian framework for inferring the metamorphic control parameter $\nu(\cdot)$ in section 4 and apply this theory to a few landmark examples in section 5. Section contains concluding remarks.

2 Metamorphosis for Landmarks

In this paper we are concerned with diffeomorphometric approaches to image and shape matching. To this end, we use time-dependent velocity fields $u_{t}$ occupying some Hilbert space $u_{t}\in V$, where $V$ is continuously embedded in $\textsf{C}_{0}^{k}(\mathbb{R}^{d})$, $k\geq 1$ inducing a curve $\varphi_{t}$ on a subgroup $\text{Diff}_{V}(\mathbb{R}^{d})$ of diffeomorphisms [4, 27] via the following ordinary differential equation

[TABLE]

This is often used in a minimisation problem where the objective is to match two images $I_{0}$ and $I_{1}$:

[TABLE]

where $F$ denotes a similarity measure between the deformed initial image $I_{0}\circ\varphi_{1}$ and the target image $I_{1}$ to allow inexact matching parameterised by $\lambda^{2}$. The LDDMM approach takes $F$ as an $L^{2}$ norm of the difference between its arguments. In order to simplify the exposition, we will consider singular solutions of this problem, which are given by the Ansatz

[TABLE]

for $M$ landmarks with position $q_{t}^{i}\in\mathbb{R}^{d}$ and momenta $p_{t}^{i}\in\mathbb{R}^{d}$. The vector field is thus

[TABLE]

where $K:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ is the kernel associated to the norm $\|\cdot\|_{V}$. This parameterisation holds throughout this paper and we set $d=2$. For metamorphosis, we introduce a discrete template variable $\boldsymbol{\eta}_{t}$ such that the deformation of a set of landmarks is written as the composition of the template position and deformation as

[TABLE]

Then, we can define the template velocity as

[TABLE]

and extend the action functional (2) to

[TABLE]

where now the reconstruction relation is

[TABLE]

obtained from (6) and (5) together with $u=\dot{\varphi}_{t}\varphi^{-1}$, see [15] for more details.

By taking variations carefully, see again [15], we directly obtain a relationship between the momentum variable $\mathbf{p}$ and the template variable $\mathbf{z}$ as

[TABLE]

and the equation of motions are

[TABLE]

where $u_{t}$ is defined in (4).

3 Selective Metamorphosis for Landmarks

We can now extend the metamorphosis setting to be able to locally control the amount of non-diffeomorphic evolution. For this, we introduce a function $\nu:\mathbb{R}^{2}\to\mathbb{R}$ replacing the parameter $\sigma^{2}$ such that $\nu(x)=\sigma^{2}$ corresponds to the classic landmark metamorphosis. The action for selective metamorphosis thus becomes

[TABLE]

which we minimise subject to the reconstruction equation (8) and the boundary conditions $\mathbf{q}_{0}$ and $\mathbf{q}_{1}$ at time $t=0,\,1$. In the case of landmarks we have as before that

[TABLE]

so we can eliminate the template variable $\mathbf{z}_{t}$ and write

[TABLE]

The problem defined by (13) yields the following equations for selective metamorphosis for landmarks:

[TABLE]

with $\mathbf{q}_{0},\,\mathbf{q}_{1}\text{ fixed}$. Again, the velocity is fully described by $\mathbf{p}$ and $\mathbf{q}$ via (4). The landmark dynamics follow standard LDDMM trajectories as $\nu(x)$ vanishes in parts of the domain. This can also be seen in the relation (12), where $\nu(x)=0\Rightarrow\mathbf{z}_{t}=0$ implies that the template variable remains fixed. Notice that these equations are Hamilton’s equations for the Hamiltonian $\mathbf{p}_{t}$ and $\mathbf{q}_{t}$

[TABLE]

where we have used (4) for the standard landmark Hamiltonian

[TABLE]

A practical procedure for solving (14) with the landmark Hamiltonian above is called shooting, where we replace the end-point condition $\mathbf{q}_{1}$ with a guess for $\mathbf{p}_{0}$, and iteratively update $\mathbf{p}_{0}$ using automatically computed adjoint (or backward) equations until $\mathbf{q}_{1}$ compares to $\mathbf{q}(1)$ below a certain tolerance. We will perform this procedure directly with an automatic differentiation package Theano [22], see [18, 17] for more details on the implementation. This section concludes with some theoretical results.

Theorem 3.1

Let $\nu$ be bounded from below away from zero by $\nu_{\text{inf}}\in\mathbb{R}$ and from above by $0<\sigma^{2}\in\mathbb{R}$. Then there exists a minimiser of (13) admissible to (8).

Proof

The functional in (13) is not convex so we work with a reformulation to ensure the required lower semi-continuity. Define a variable $w^{i}_{t}=\sqrt{\nu(q_{t}^{i})}p^{i}_{t}$ in the problem:

[TABLE]

First, note that owing to the constraint effectively being a boundary value problem, we cannot always find a $\mathbf{q}$ for arbitrary pairs of $(u,\,\mathbf{w})$. We define a bounded operator $(\mathbf{q},\,u_{t})\mapsto\frac{\dot{\mathbf{q}_{t}}-u_{t}(\mathbf{q}_{t})}{\sqrt{\nu(\mathbf{q}_{t})}}\triangleq\mathbf{w}$:

[TABLE]

From this we generate a minimising sequence $(\mathbf{q}^{n},u^{n},\mathbf{w}^{n})_{n\geq 0}$ admissible to (3). The rest of the proof is standard, see e.g. [27]. We show the constraint equation is continuous with respect to the weak topology on $X\triangleq H^{1}([0,1],\,\mathbb{R}^{d\times M})\times L^{2}([0,1],\,V)\times L^{2}([0,1],\,\mathbb{R}^{d\times M})$ i.e. $e(\mathbf{q}_{t}^{n},\,\mathbf{w}_{t}^{n},\,u_{t}^{n})\rightharpoonup e(\mathbf{q}_{t},\,\mathbf{w}_{t},\,u_{t})$ where $e(q,\,w,\,u)\triangleq\dot{q}-u(q)-\sqrt{\nu(q)}w$. Then,

[TABLE]

Further, for $\phi\in L^{2}([0,1],\,V)$,

[TABLE]

The first term vanishes trivially, while for the second we see

[TABLE]

Since linear operators are naturally compatible with the weak topology the required continuity follows. Passing to subsequences where necessary we can by classic results extract bounded subsequences converging to weak limits where necessary to obtain a minimiser. Convexity of $S$ implies weak lower semi-continuity concluding the proof. $\square$

Theorem 3.2

Assume $\nu\in W^{2,\infty}(\mathbb{R}^{d})$ and $V$ is embedded in $\textsf{C}_{0}^{k}(\mathbb{R}^{d})$, $k\geq 1$ (continuous functions with continuous derivatives to order $k$ vanishing at infinity). Then, given $\mathbf{p}_{0},\,\mathbf{q}_{0},\in\mathbb{R}^{d\times M}$, (14) with (4) are integrable for all time.

Proof

Establishing appropriate Lipschitz conditions implies integrability of the system akin to [6, Theorem 5]. We note that the kernel in (4) is Lipschitz in $(p_{t},q_{t})$ by assumption, so the composition $(p,\,q)\mapsto u\circ q$ is also Lipschitz. $u(q)\mapsto\nabla u(q)^{T}$ consider $v,\,w\in V$ and $x,\,y\in\mathbb{R}^{d}$:

[TABLE]

so the mapping is Lipschitz in both the position and velocity. Given the conditions on $\nu$ the mappings

[TABLE]

are locally Lipschitz. Consequently we verify that for any $(\mathbf{p}_{0},\,\mathbf{q}_{0})\in B(0,r)\subset\mathbb{R}^{d\times M}\times\mathbb{R}^{d\times M}$, the system (14) is locally Lipschitz with constant $L_{r,t_{0}}$ for some $t_{0}>0$. By the conservation of the Hamiltonian we can extend the existence of solutions to arbitary $t>t_{0}$. $\square$

4 Bayesian Framework

We now place a stochastic model on $\nu$ inspired by the approach taken in [6], see also [21, 1, 2] for similar Bayesian approaches in computational anatomy. The goal is to develop an algorithm to infer $\nu$, passing via the deterministic problem seen above. First, we present some preliminaries on the Bayesian approach to inverse problems in section 4.1, essentially quoting concepts and results from [9], or [7] for an exposition of algorithmic aspects of function space MCMC. Section 4.2 then formally describes how we apply this stochastic approach to inverse problems to $\nu$ by a finite-dimensional family of parameterisations.

4.1 General Framework

The general framework is based on the idea that we can cast optimisation problems in a probabilistic framework where minimisers, roughly speaking, correspond to modes of a certain distribution over function space. In the context of optimisation we define a likelihood $\Phi:X\rightarrow Y$, mapping some control variable in $X$ to an observable in $Y$. A maximum a posteriori (MAP) estimator $f^{*}$ satisfies $f^{*}=\arg\max_{f\in X}\Phi(f)p(f)$ where $p$ is a density over $X$. Equipped with a norm $\|\cdot\|_{X}$ we can then define the Gaussian density by $p(\cdot)\propto e^{-\|\cdot\|_{X}^{2}}$. Supposing further that the likelihood is on the form $\log\Phi(f)=-\|f-\lambda\|_{Y}^{2}$, for some desiderata $\lambda\in Y$ then, at least formally, the MAP estimator minimises $J=\|f^{*}\|_{X}^{2}+\|f^{*}-\lambda\|_{Y}^{2}$.

For general inverse problems on function space, several key properties documented in [9] must be verified before the inverse problem is well-posed. Beyond showing the existence of the MAP estimator minimising $J$ above, the infinite-dimensional version of Bayes’ rule must also be checked i.e. the Radon-Nikodym derivative of the prior with respect to the posterior must exist and be absolutely continuous. Finally, we request continuity of the posterior distribution w.r.t the initial data corresponding in a sense to Hadamard well-posedness in a probabilistic framework. Rigorously treating this Bayesian inverse problem is subject to further study in forthcoming works.

4.2 Finite-dimensional Parameterisation

We now introduce the main problem of this paper in the setting above. We consider $\nu$ as a random variable. If the growth location is not known a priori, then this is an appropriate framework as it allows for a qualitative evaluation of selective metamorphosis. Moreover, it can account for some observation error by its probabilistic nature. As a simple example, we consider the case where $\nu$ is given by a sum of Gaussians

[TABLE]

Here the finite family $h_{k}$ of centroids in $\mathbb{R}^{d}$ together with the uniform length-scale $\sigma_{\nu}$ fully determine $\nu_{h}$, thus greatly reducing the complexity of sampling. We defer sampling from function space to future work. Defining a density $p_{sm}\propto e^{-S_{sm}^{\nu}}$ over the space of triples $(\nu,\,\mathbf{q}_{\nu},\,\mathbf{p}_{\nu})$ leads to the preconditioned Crank-Nicholson algorithm 1, see e.g. [14].

Here, $N$ denotes the desired number of samples and $K$ the number of terms in (19). randomUnit() denotes a randomly generated number in $[0,\,1]$. The coefficient $\beta$ scales between the previous sample and the new step $\xi$. This needs to be calibrated as a too low value may increase the acceptance rate by taking shorter steps at the cost of slow exploration of the state space. Conversely, a too high value of $\beta$ results in lower acceptance and thus convergence. The next section shows the algorithm above in practice.

5 Numerical examples

This section displays some numerical results for our method. We apply algorithm 1 to infer a distribution for the growth location using the landmark boundary conditions seen in figure 1. The parameters and results for the first configuration is shown in figure 2. These preliminary results demonstrate that even for a small number of samples the density of accepted samples corresponds at least heuristically to the analytical density histogram obtained by computing the value of the metamorphosis functional in (11).

We arrive at the same conclusion for the second example, for which the results are shown in figure (3). Moreover, we note that the geodesic equations for $p$ and $q$ are time-reversible meaning that the configuration in figure 3 corresponds to both particle collapse as well as hole creation.

It is numerically relatively simple to control the behaviour of $\nu$ by simple scaling or by adding regularisation terms to (11) to e.g. penalise having $\nu$’s far away from the support of the images. Such cost can easily be added to the MCMC algorithm, depending on the prior information one can have on the shape matching problem.

6 Conclusion

We have presented a preliminary approach for selectively allowing photometric variation in a diffeomorphic image matching. We analysed the selective metamorphosis problem, the associated geodesic equations and demonstrated a proof of concept MCMC algorithm inferring a simple parameterisation of $\nu$. This generalises LDDMM and metamorphosis and could provide a first-order exploratory tool for physicians to see if the development of a biological feature stems from a few violations of diffeomorphic evolution. This paper paves the way towards surgically investigating growth phenomena between topologically different images.

Future work is manifold. Firstly, we aim to extend the equations of section 3 to images e.g. using the kernel framework in [20] or developing a space-time method. We also aim to find an explicit solution to the geodesic equations for $\mathbf{p}$ and $\mathbf{q}$ and with the additional terms involving $\nu$ à la [25] to eliminate the need $\nu$ to be bounded from below. Further, as outlined in 4 there are many aspects of the probabilistic framework that need rigorous treatment. Beyond the references therein, see also the work in [8]. Natural extensions of our probabilistic approach also include fully treating $\nu$ as a function and interpreting the resulting inverse problem through the appropriate measure-theoretical lens. Adding a time-dependency to $\nu$ can also be explored. Determining a truncated Fourier series of $\nu$ could lead to efficient numerical methods. More generally, we hope to reconcile our attempts to model growth here with the mathematically elegant approach described in [16] and with the more general mathematics of growth [10].

Finally, it is our hope that we can extend the probabilistic approach developed here to encompass classic metamorphosis as well; that is to say, to develop the necessary theory in order to place a stochastic model on the state space consisting of velocities and source functions and sample from function space. This provides a derivative-free method of solving classic metamorphosis (or other problems in shape analysis) at the expense of interpreting the results probabilistically. In this work, we only used a simple MCMC algorithm, but a Metropolis-adjusted Langevin algorithm or Hamiltonian Monte-Carlo algorithm may be more appropriate to solve this problem.

Acknowledgements

AA acknowledges EPSRC funding through award EP/N014529/1 via the EPSRC Centre for Mathematics of Precision Healthcare.