Probabilistic Diffeomorphic Registration: Representing Uncertainty

TL;DR

This paper introduces a Bayesian variational framework for representing uncertainty in large deformation diffeomorphic image registration, enabling estimation of the full posterior distribution of deformations.

Contribution

It proposes a novel SDE-based Gaussian process prior for deformations, allowing direct approximation of the full posterior distribution instead of MAP or Monte Carlo methods.

Findings

Effective uncertainty quantification demonstrated on landmark-based registration.

Applicable to simulated and real 3D medical images.

Provides a probabilistic interpretation of deformation fields.

Abstract

This paper presents a novel mathematical framework for representing uncertainty in large deformation diffeomorphic image registration. The Bayesian posterior distribution over the deformations aligning a moving and a fixed image is approximated via a variational formulation. A stochastic differential equation (SDE) modeling the deformations as the evolution of a time-varying velocity field leads to a prior density over deformations in the form of a Gaussian process. This permits estimating the full posterior distribution in order to represent uncertainty, in contrast to methods in which the posterior is approximated via Monte Carlo sampling or maximized in maximum a-posteriori (MAP) estimation. The frame-work is demonstrated in the case of landmark-based image registration, including simulated data and annotated pre and intra-operative 3D images.

Figures

Equations

$$p(\varphi |M,F)=p(\varphi )p(M,F|\varphi )/p(M,F).$$

$$\text{KL}\left[q\Vert p(\varphi |M,F)\right]$$

$$\text{KL}\left[q\Vert p(\varphi |M,F)\right]=\text{KL}\left[q\Vert p(\varphi )\right]+{⟨m(\varphi \circ M,F)⟩}_q+\log p(M,F),$$

$$\frac{d}{dt}{\varphi }_t(x)=v_t({\varphi }_t(x)),{\varphi }_0(x)=x,t\in [0,1].$$

$$E(\varphi ;M,F)=E_v(v)+E_{\varphi }(\varphi ;M,F),E_v(v)=\frac{1}{2}{\int }_0^1{\int }_{{\Omega }_m}\Vert L{v}_t(x){\Vert }_2^{2}dxdt.$$

$$p(\mathbf{v})=\mathcal{G}\mathcal{P}(0,{\Sigma }_{ts}(x,y)),{\Sigma }_{ts}(x,y)\in {\mathrm{SPD}}^d$$

$$d{\varphi }_t(x)$$

$${\mathbf{W}}_t$$

$$\frac{d}{dt}{\bar{\varphi }}_t(x)$$

$${\Lambda }_{ts}(x,y)$$

$$\frac{d}{dt}{⟨{\varphi }_t(x){\varphi }_s^{⊺}(y)⟩}_p$$

$$+{⟨{\Sigma }_{ts}({\varphi }_t(x),{\varphi }_s(y))⟩}_p.$$

$$p(\varphi )=\mathcal{G}\mathcal{P}({\bar{\varphi }}_t(x),{\Lambda }_{ts}(x,y)),$$

$$-\log (p(\mathbf{v}))=E_v(v)+const=\frac{1}{2}{\int }_0^1{\int }_{{\Omega }_M}\Vert L{\mathbf{v}}_t(x){\Vert }_2^{2}dxdt+const,$$

$$-\log (p(\mathbf{V}))=\frac{1}{2}{\left(L\mathrm{vec}\mathbf{V}\right)}^{⊺}L\mathrm{vec}\mathbf{V}+\text{const}=\frac{1}{2}\mathrm{vec}{\mathbf{V}}^{⊺}{L}^{⊺}L\mathrm{vec}\mathbf{V}+\text{const},$$

$${⟨d({\varphi }_t(x){\varphi }_s^{⊺}(y))⟩}_p=d{⟨{\varphi }_t(x){\varphi }_s^{⊺}(y)⟩}_p={⟨v_t({\varphi }_t(x))dt{\varphi }_s^{⊺}(y)⟩}_p+{⟨{\sqrt{\Sigma }}_t({\varphi }_t(x))d{\mathbf{W}}_t{\varphi }_s^{⊺}(y)⟩}_p$$

$$+{⟨{\varphi }_t(x)(v_s({\varphi }_s(y)){)}^{⊺}ds⟩}_p+{⟨{\varphi }_t(x)({\sqrt{\Sigma }}_t({\varphi }_s(y))d{\mathbf{W}}_s{)}^{⊺}⟩}_p$$

$$\displaystyle\quad+\left\langle\bigg{(}v_{t}(\boldsymbol{\mathbf{\phi}}_{t}(y))dt+\sqrt{\Sigma}_{t}(\boldsymbol{\mathbf{\phi}}_{t}(y))d\boldsymbol{\mathbf{W}}_{t}\bigg{)}\bigg{(}v_{s}(\boldsymbol{\mathbf{\phi}}_{s}(y))ds+\sqrt{\Sigma}_{s}(\boldsymbol{\mathbf{\phi}}_{s}(y))d\boldsymbol{\mathbf{W}}_{s}\bigg{)}^{\mathstrut\intercal}\right\rangle_{p},$$

$$v_t(x)$$

$${\sqrt{\Sigma }}_t(x)$$

$$d{\varphi }_t(x)\approx (A_t{\varphi }_t(x)+a_t)dt+\left(\sum _i{S}_t^i{\varphi }_t(x{)}_i+R_t\right)d{\mathbf{W}}_t$$

$$\begin{array}{rlrl}{A}_t& \triangleq {\mathsf{D}}_{v_{t_0}}^x(x_0)& a_t& \triangleq -{\mathsf{D}}_{v_{t_0}}^x(x_0)x_0+{\partial }_t{v}_{t_0}(x_0)(t-t_0)+v_{t_0}(x_0)\\ S_t^i& \triangleq {\partial }_{x_i}{\sqrt{\Sigma }}_{t_0}(x_0)& R_t& \triangleq -\sum _i{S}_t^i\cdot (x_0{)}_i+{\sqrt{\Sigma }}_{t_0}(x_0).\end{array}$$

$$\frac{d{\bar{\varphi }}_t(x)}{dt}\approx v_t({\bar{\varphi }}_t(x))\text{ where }{v}_t(x)\approx (A_t{\bar{\varphi }}_t(x)+a_t),$$

$$\begin{array}{r}\frac{d⟨{\varphi }_t(x){\varphi }_t^{⊺}(y)⟩}{dt}\approx A_t⟨{\varphi }_t(x){\varphi }_t^{⊺}(y)⟩+⟨{\varphi }_t(x){\varphi }_t^{⊺}(y)⟩A_t^{\prime ⊺}+a_t{\bar{\varphi }}_t^{⊺}(y)+{\bar{\varphi }}_t(x)a_t^{\prime ⊺}\\ +\sum _{ij}{S}_t^i⟨{\varphi }_t(x){\varphi }_t^{⊺}(y)⟩{\left(S_t^{\prime j}\right)}^{⊺}+\left(\sum _i{S}_t^i{\bar{\varphi }}_t(x{)}_i\right)R_t^{\prime ⊺}+R_t{\left(\sum _i{S}_t^{\prime i}{\bar{\varphi }}_t(y{)}_i\right)}^{⊺}+R_t{R}_t^{\prime ⊺},\end{array}$$

$$\mathcal{E}(q({\varphi }_1))=\text{KL}\left[q({\varphi }_1)\Vert p({\varphi }_1)\right]+⟨m({\varphi }_1\circ M,F){⟩}_q+\log p(M,F)$$

$$\text{KL}\left[q({\varphi }_1)\Vert p({\varphi }_1)\right]={\int }_0^1{\mu }_t^{⊺}{S}_{t=0}^{-1}{\mu }_{t}dt,\text{ s.t. }\frac{d{\bar{\varphi }}_t}{dt}\approx {\mu }_t({\bar{\varphi }}_t),{\bar{\varphi }}_0=\mathrm{id}.$$

$$\begin{array}{c}\underset{{\mu }_{t=0}}{\mathrm{argmin}}\mathcal{E}(q_{{\mu }_{t=0}}(\varphi ))=\frac{1}{2}{\mu }_{t=0}^{⊺}{S}_{t=0}^{-1}{\mu }_{t=0}^{⊺}+⟨m({\varphi }_1\circ M,F){⟩}_q-\log p(M,F).\end{array}$$

$${⟨m({\varphi }_1\circ M,F)⟩}_q={⟨m({\mathbf{\Phi }}_1,F)⟩}_q=\mathrm{tr}{⟨{\mathbf{\Phi }}_1{\mathbf{\Phi }}_1^{⊺}⟩}_q-2{\bar{\Phi }}_1^{⊺}F+\mathrm{tr}F{F}^{⊺}.$$

$${\nabla }_{{\mu }_{t=0}}\mathcal{E}(q_{{\mu }_{t=0}}({\varphi }_1))=\frac{1}{2}{S}_{t=0}^{-1}{\mu }_{t=0}-\left(2{\bar{\Phi }}_1-2F\right).$$

$${\Sigma }_{ts}(\vec{x},\vec{y}{)}_{ij}=\delta (t-s){\left[\exp \left(-\frac{\Vert {\vec{x}}_i-{\vec{y}}_j{\Vert }_2^2}{2{\sigma }^2}\right)\right]}_{ij}\in {\mathbb{R}}^{d\times d}$$

Full text

1 Introduction

Deformable image registration seeks to identify a deformation field that aligns two images, and is a key component of image analysis applications such as computational anatomy [1, 2]. An important body of literature focuses on deformations in the form of diffeomorphisms [3, 4, 2], one-to-one mappings between image coordinate systems that are smooth and invertible. These properties help in ensuring biologically plausible deformations, and avoiding phenomena such as folding or tearing that may occur in non-diffeomorphic registration approaches [5].

While a good deal of literature has focused on identifying optimal diffeomorphic registration solutions, it would be useful to quantify the inherent uncertainty in these solutions when interpreting the results of registration. Quantification of deformable registration uncertainty, particularly at point locations throughout the image, remains an open problem. The Bayesian approach quantifies probabilistic uncertainty via a posterior distribution over deformations conditioned on image data. Estimating the full posterior in the case of large deformation diffeomorphisms is desirable but computationally challenging, and has typically been avoided. Simpson et al. propose a Bayesian variational framework based on small deformation kinematics [5], however this does not address the general case of large deformations. Markussen proposed a stochastic differential equation (SDE) model for large deformations, however only the computation of the maximum a posteriori deformation is provided lacking the estimation of a distribution on the deformations [6]Alternatively, the posterior may be investigated via sampling methods, e.g. Markov chain Monte Carlo (MCMC) [7] or Hamiltonian Monte Carlo [8].

This paper introduces a novel mathematical framework that allows representing and computing of the full Bayesian posterior in the case of large deformation diffeomorphisms. Our framework considers a SDE modeling the deformation field as the evolution of a time-varying velocity field, with additive noise in the form of a Wiener process. A Gaussian process (GP) density results from a locally linear approximation of the SDE and taking the initial deformation field to be Gaussian process distributed. Deformation field uncertainty is quantified by the point-wise covariance of the deformation field throughout the image, and can be summarized, e.g., via the Frobenius norm of the covariance (FC). This can be pictured through the following example: if the FC at a point approaches [math], the marginal density of the transform approaches an impulse function denoting the existence of a single probable solution. On the other hand, when FC is large, the density becomes “broader” denoting a larger set of solutions with high probability at that point. Hence the point-wise FC is a model of uncertainty. Experiments demonstrate our framework in the context of landmark correspondences, where a heteroscedastic model accounts for variable uncertainty in landmark localization. This is particularly useful when estimates of landmark localization variability are available.

2 Methods

2.1 Variational Approximation to Registration

We start by posing the registration problem in a probabilistic framework. Let $M$ and $F$ be moving and fixed objects with domains in $\Omega_{M}$ and $\Omega_{F}$ respectively, and let $\phi:\Omega_{M}\mapsto\Omega_{F}$ be a mapping between the two. The registration problem seeks a posterior probability density over mappings $\phi$ conditioned on data $(M,F)$, which is expressed via Bayes theorem as

[TABLE]

In Eq. 1, $p(\boldsymbol{\mathbf{\phi}})$ is a prior density over $\boldsymbol{\mathbf{\phi}}$ embodying geometrical constraints such as smoothness. $p(M,F|\boldsymbol{\mathbf{\phi}})$ is the data attachment factor or likelihood of the map $\phi$ relating $F$ and $M$. E.g. the probability that $M$ deformed by $\phi$, which we note $\phi\operatorname{\circ}M$, is similar to $F$. Finally, $p(M,F)$ is a normalizing constant.

The direct calculation of the posterior density $p(\boldsymbol{\mathbf{\phi}}|M,F)$ is a difficult problem. Hence, we use a variational method to estimate a distribution $q(\boldsymbol{\mathbf{\phi}})$ (abbreviated as $q$) that is close to $p(\boldsymbol{\mathbf{\phi}}|M,F)$ in the sense of the Kullback-Leibler divergence ($\textrm{KL}\left[\cdot\|\cdot\right]$). Specifically, we seek $q$ minimizing:

[TABLE]

In the registration literature the data attachment factor $p(M,F|\boldsymbol{\mathbf{\phi}})$ is typically modeled using a measure of similarity between the registered objects: $m:\Omega_{F}\times\Omega_{F}\mapsto\mathbb{R}$ which is minimal when two objects are exactly the same and grows as they become different. Adopting the equality $-\log p(M,F|\boldsymbol{\mathbf{\phi}})=m(\phi\operatorname{\circ}M,F)$, Eq. 2 may be rewritten as:

[TABLE]

where $\left\langle m(\boldsymbol{\mathbf{\phi}}\operatorname{\circ}M,F)\right\rangle_{q}$ is the expected value of $m$ with respect to the density $q$. There are two main differences of this formulation with respect to common diffeomorphic registration approaches [3, 4]. First, instead of seeking a single optimal deformation $\phi$, e.g. the maximum a-posteriori (MAP) solution in the Bayesian formulation, we seek to obtain the full distribution $q(\boldsymbol{\mathbf{\phi}})$. In this way, we obtain both the MAP deformation $\phi$ in addition to the uncertainty at any given point in space, which can be calculated from $q(\boldsymbol{\mathbf{\phi}})$. Second, we obtain $q(\boldsymbol{\mathbf{\phi}})$ by minimizing the data attachment term over a weighted combination of all possible deformation fields in the family of $\phi$ instead of only at a single deformation $\phi$.

2.2 Probabilistic Diffeomorphic Deformations

The variational approximation to $p(\boldsymbol{\mathbf{\phi}}|M,F)$ described in the previous section requires a parameterization for $q(\boldsymbol{\mathbf{\phi}})$ over which Eq. 3 can be minimized. In this section we derive a novel parameterization in the form of a Gaussian Process (GP). The theoretical basis for our derivation lies in a stochastic interpretation of the work of [1], common to many diffeomorphic registration approaches [2]. Here, we begin by outlining the relevant elements of this work, then we present our derivation in three propositions and their proofs, with our primary contributions being in Propositions 2 and 3.

Following the work of [1], many diffeomorphic deformation formulations seek an optimal registration solution $\phi$, e.g. the MAP deformation in Eq. 1, by constraining the map $\phi$ to be the solution at $t=1$ of the ordinary differential equation (ODE)

[TABLE]

and setting $\phi_{1}$ in Eq. 4 to minimize

[TABLE]

From the two terms of $E$, $E_{v}(v)$ regularizes the evolution of the time-varying velocity field and $E_{\phi}(\phi;M,F)$ drives $E$ such that the deformed moving object, $\phi_{1}\operatorname{\circ}M$, becomes as similar to $F$ as possible. The regularization term $E_{v}(v)$, is driven by $L$, a linear differential operator. The key insight here is that given suitable $L$ according to [1], Eq. 5 restricts $\phi_{t}(x)$ defined as in Eq. 4 to the space of diffeomorphisms [1]. From this, we make the following propositions:

Proposition 1

Under a probabilistic interpretation, the regularization term $E_{v}(v)$ in Eq. 5 corresponds to the negative logarithm of the GP prior on the stochastic velocity field $\boldsymbol{\mathbf{v}}$

[TABLE]

where the covariance function $\Sigma_{{ts}}(x,y)$, representing the relation between the point $x$ at time $t$ and the point $y$ at time $s$, is determined by the operator $L$ in Eq. 5.

Proposition 2

Interpreting the energy $E_{v}(v)$ in Eq. 5 as the negative logarithm of the density of a stochastic process $\boldsymbol{\mathbf{v}}$ induces a random process $\boldsymbol{\mathbf{\phi}}$ with density $p(\boldsymbol{\mathbf{\phi}})$ on the deformation field of Eq. 4 that is a solution of the stochastic differential equation (SDE)

[TABLE]

where $\sqrt{\Sigma}$ is the square root matrix $\sqrt{\Sigma}\sqrt{\Sigma}^{\mathstrut\intercal}=\Sigma$; the GP $\boldsymbol{\mathbf{W}}_{t}\in\mathbb{R}^{d}$ is called Brownian motion or a Wiener process [9]; $v_{t}(x)\in\mathbb{R}^{d}$ is a deterministic velocity field like in Eq. 4; and $\Sigma_{t}(x)\triangleq\Sigma_{tt}(x,x)\in\operatorname{SPD}^{d}$ 111$SPD^{d}$: set of symmetric positive definite matrices of dimension $d$, the covariance of the probabilistic prior defined in Prop. 1, is a consequence of Eq. 5.

Proposition 3

For the stochastic process $\phi$ with density $p(\phi)$, defined in prop. 2, the mean $\bar{\phi}$ and covariance $\Lambda$ functions are solutions of the deterministic ODEs

[TABLE]

Moreover, up to a first order approximation:

[TABLE]

The proofs for Propositions 1-3 are as follows:

Proof of Prop. 1 This proposition has been proven by Joshi et. al.[2], here we provide a sketch of the relevant points. We start by relating $E_{v}(v)$ in Eq. 5 to a probability density on velocity fields as stochastic processes $\boldsymbol{\mathbf{v}}$, $p(\boldsymbol{\mathbf{v}})$:

[TABLE]

To show that $p(\boldsymbol{\mathbf{v}})$ is a stochastic process with a particular distribution, we need to prove that any finite sample of the domain $\Omega_{M}\in\mathbb{R}^{d}$ has the same parametric distribution [10]. We take $N$ samples $X\in\mathbb{R}^{N\times d}$ in space and $t\in[0,1]^{N}$ in time, and let $\boldsymbol{\mathbf{V}}_{ij}=\left[\boldsymbol{\mathbf{v}}_{t_{i}}\left(X_{i}\right)\right]_{j}\in\mathbb{R}^{N\times d}$. Then we rewrite Eq. 11 as

[TABLE]

where $L$ is the matrix such that $[L\boldsymbol{\mathbf{V}}_{i\cdot}]=[L\boldsymbol{\mathbf{v}}_{t_{i}}(X_{i})]_{j}$. Eq. 12 is recognisable as the log probability of a centered multivariate Gaussian with covariance $C=(L^{\mathstrut\intercal}L)^{-1}$ and, therefore $\boldsymbol{\mathbf{v}}(x)$ is a GP. The covariance function $\Sigma_{{ts}}(x,y)$ can be calculated as the matrix Green’s function of the operator $L$ [2]; specifically, if $x,y\in\mathbb{R}^{d}$, $\Sigma_{{ts}}(x,y)\in\operatorname{SPD}^{d}$ where $[\Sigma_{{ts}}(x,y)]_{ij}$ is the covariance between $x_{i}$ at time $t$ and $y_{j}$ at time $s$.

This shows that for a given velocity field ${\boldsymbol{\mathbf{v}}}(x)$, random perturbations according to the regularization term $E_{v}(v)$ in Eq. 5 or prior Eq. 11 follow a GP, therefore the velocity fields according to $E_{v}(v)$ in Eq. 5 have the density specified in Eq. 6 proving Prop. 1.

Proof of Prop. 2 A formal proof of Prop. 2 is beyond the scope of this paper. Instead, using Prop. 1 we argue its validity and provide appropriate references. In Prop. 1 we characterized the density of random perturbations of velocity fields according to Eq. 5. Adding such random perturbations to Eq. 4 leads to Eq. 7.

The second term in Eq. 7 comes from considering the velocity fields $v_{t}$ of Eq. 4 as a stochastic process according to Prop. 1. We achieve this by perturbing the right hand side of Eq. 4 with noise. The factor $\boldsymbol{\mathbf{W}}_{t}\in\mathbb{R}^{d}$ is white noise, which multiplied by $\sqrt{\Sigma}_{t}$ is a centered Gaussian random variable with covariance $\Sigma_{t}$, a sample drawn from Eq. 6. We noted the stochastic velocity field in Eq. 7 ${\boldsymbol{\mathbf{v}}}$ to distinguish it from the deterministic one $v$. Eq. 7 ceases to be an ODE as a sample path of $\boldsymbol{\mathbf{W}}_{t}$ is almost surely not differentiable. Alternatively, using the Itō interpretation of Eq. 7 leads to the SDE in Eq. 7, whose solution is the density on $\boldsymbol{\mathbf{\phi}}$ [9, Chap. 8].

Proof of Prop. 3 The ODE for the mean of the stochastic process $\boldsymbol{\mathbf{\phi}}$, Eq. 9a, is obtained by calculating the expectation on both sides of Eq. 7. It is a consequence of the linearity of the expected value and the derivative operator and the definition of $\boldsymbol{\mathbf{W}}_{t}$ as a zero-centered Wiener process in Eq. 8.

To obtain the ODE for the second moment of $\boldsymbol{\mathbf{\phi}}$, shown in Eq. 9c, we use the Itō product rule [9] to obtain an expression for $d(\boldsymbol{\mathbf{\phi}}_{t}(x)\boldsymbol{\mathbf{\phi}}_{s}^{\mathstrut\intercal}(y))$ and substitute it in Eq. 7 obtaining

[TABLE]

which, using the Itō identities for expected values of differentials [9] results in Eq. 9c.

Obtaining a parametric form of the density of $\phi$, $p(\phi)$, satisfying the SDE (7) in the general case is an open problem and a wide field of study. However, in the case where the drift $v$ and diffusion coefficient $\sqrt{\Sigma}$ are linear functions on their time and location parameters, and the initial condition $\phi_{t=0}$ is a GP, $\phi_{t}(x)$ is known to be a GP [9]. With this purpose we define a locally linearized (LL) $v$ and $\sqrt{\Sigma}$ centered at $t_{0},x_{0}$ [11]:

[TABLE]

where $\mathsf{D}_{v_{t}}^{x}\!$ is the Jacobian of $v_{t}(x)$ w.r.t. $x$ and $\partial_{t}{v_{t}}\!$ its partial derivative w.r.t. $t$. Considering that $L$ is assumed time-invariant in Prop. 1, the time derivative of $\sqrt{\Sigma}$ in Eq. 14b is equal to [math]. Then, using the LL equations Eqs. 14a and 14b, we approximate Eq. 7 as

[TABLE]

The LL approximations in Eqs. 15 and 14a lead to an approximation of the mean function of $\boldsymbol{\mathbf{\phi}}$, $\bar{\phi}$, by the solution of the ODE

[TABLE]

and its second moment $\left\langle\boldsymbol{\mathbf{\phi}}_{t}(x)\boldsymbol{\mathbf{\phi}}_{s}(y)^{\mathstrut\intercal}\right\rangle$ when $t=s$ by

[TABLE]

where $A^{\prime}_{t}$; $a^{\prime}_{t}$; $S^{\prime i}_{t}$; and $R^{\prime}_{t}$ are the same as $A_{t}$; $a_{t}$; $S^{i}_{t}$; and $R_{t}$ in Eq. 16 substituting $y$ and $y_{0}$ for $x$ and $x_{0}$.

As long as the initial condition $\boldsymbol{\mathbf{\phi}}_{t=0}$ is a GP, the linear approximation of $\boldsymbol{\mathbf{\phi}}_{t}$ is a GP uniquely determined by $\bar{\phi}$ and $\Lambda$ [9]. Then, given a set of stochastic velocity fields $v_{0}\ldots v_{t_{M-1}}$ with $t_{0}=0$ and $t_{M-1}=1$, the parameters of the stochastic process representing the transform $\boldsymbol{\mathbf{\phi}}$ are obtained integrating Eqs. 17 and 18 with the initial conditions $\boldsymbol{\mathbf{\phi}}_{t=0}\sim\mathcal{GP}(\bar{\phi}_{t=0},\Lambda_{t=0})$. Having characterized stochastic transformations representing a diffeomorphic deformation, we are in position to formulate our probabilistic diffeomorphic registration algorithm.

2.3 Probabilistic Diffeomorphic Registration

The stochastic diffeomorphic deformation model of Section 2.2 leads to a GP approximation on deformation fields, whose parameters are determined by $v$ and $\Sigma$; we use this model as $q(\boldsymbol{\mathbf{\phi}}_{1})$, our variational distribution. In this section, we show how to compute the parameters of $q(\boldsymbol{\mathbf{\phi}}_{1})$ minimizing Eq. 2 for a particular registration problem. Taking the approach of [3], we focus on operators $L$ for the energy Eq. 5 regularizing in space but not in time. Due to the time-independent regularization, the prior on velocity fields of $\boldsymbol{\mathbf{\phi}}$ derived with Prop. 1 is the joint probability of the fields at each time $t$: $p(\boldsymbol{\mathbf{\phi}}_{1})=\prod_{0}^{1}p(\boldsymbol{\mathbf{v}}_{t})^{dt}$ with $p(\boldsymbol{\mathbf{v}}_{t})\sim\mathcal{GP}(0,\Sigma_{0})$. Then, we rewrite leftmost term of Eq. 2 as $\textrm{KL}\left[q\|p(\boldsymbol{\mathbf{\phi}}_{1})\right]=\int_{0}^{1}\textrm{KL}\left[q(\boldsymbol{\mathbf{v}}_{t})|p(\boldsymbol{\mathbf{v}}_{t})\right]dt$. We parameterize each stochastic velocity field $\boldsymbol{\mathbf{v}}_{t}$ by a $N$-point set represented as a matrix $X_{t}\in\mathbb{R}^{N\times d}$ rendering its mean equivalent to a spline model [10]. This sets the distributions of the discretized velocity field prior to $p(\operatorname{vec}\boldsymbol{\mathbf{v}}_{t}|X_{t})=\mathcal{G}(0,S_{t=0}|X_{t})$. As in usual LDDMM approaches, we keep the $L$ operator, hence the covariance $S$, fixed. Hence, the parameterized form of variational approximation to the posterior of the velocity fields becomes $q(\operatorname{vec}\boldsymbol{\mathbf{v}}_{t}|X_{t})=\mathcal{G}(\mu_{t},S_{t=0}|X_{t})$. Due to GP properties given the mean and covariance functions for the GP, we can characterize the mean and covariance for the discretized velocity field as, $\mu_{t}(X)=\operatorname{vec}v_{t}(X)$ and $[S_{t=0}(X)]_{di+k,dj+l}=[\Sigma_{t=0}(X_{i},X_{j})]_{kl},i,j=1\ldots N,k,l=0\ldots d-1$ [10]. This leads to an objective which we minimize to obtain $q(\boldsymbol{\mathbf{\phi}}_{1})$ representing the registration problem:

[TABLE]

Using the ideas of [3], Eq. 19 can be minimized through geodesic shooting [4], i.e. it depends only on $M$, $F$ and $\mu_{t=0}$. The shooting equations for the proposed probabilistic diffeomorphic registration can be derived using Eq. 19 in combination with the evolution equation based on the most probable velocity field $\mu$. In fact, the problem formulation equations, shown in Eq. 4, stay the same as in [3], only the final condition $E_{\phi}$ changes, which is then “warped” to $t=0$ for a gradient descent with respect to the initial velocity $\mu_{t=0}$ leading to the objective function of $q(\boldsymbol{\mathbf{\phi}}_{1})$ parameterized on $\mu_{t=0}$

[TABLE]

Up to this point the framework we presented is general for cases where $M$ and $F$ are images or landmarks. Henceforth, we specialize the treatment of the registration problems for the landmark case where $M$ and $F$ are matrices in $\mathbb{R}^{N\times d}$; $m(M,F)=\|M-F\|^{2}_{2}$; and the random variable $\boldsymbol{\mathbf{\Phi}}_{t}=\operatorname{vec}(\boldsymbol{\mathbf{\phi}}_{t}\operatorname{\circ}M)\triangleq\operatorname{vec}\boldsymbol{\mathbf{\phi}}_{t}(M)$. This allows us to rewrite

[TABLE]

Replacing Eq. 21 in Eq. 20 leads to the gradient

[TABLE]

Having this gradient, we minimize $\mathcal{E}$ w.r.t. $\mu_{t=0}$ using a gradient descent algorithm.

3 Experiments

We are now in position to perform experiments using our probabilistic diffeomorphic registration algorithm. For all our experiments, we chose the covariance function

[TABLE]

where $\sigma^{2}$ is the model parameter. To conclude the specification of the model, we know with certainty that the starting point of the registration algorithm is the identity transform, hence $\phi_{0}(\vec{x})\sim\mathcal{GP}(\vec{x},\Sigma_{t=0,s=0}(\vec{x},\vec{y}))$ and $\Sigma_{t=0,s=0}(\vec{x},\vec{y})=\vec{0}$.

3.1 Validity of the Locally Linear Approximation

To test the validity of our GP model for diffeomorphic deformations we compared the GP through the LL method with one of the standard numerical solver for SDEs which does not assume a parametric density on $\phi$ [9]. We generated two sets of landmarks, as shown in Fig. 1, a circle and one resembling a flower, both of radius $10mm$. Then we generated random initial velocity fields with the covariance function in Eq. 22 with a range of $\sigma\in\{.1,2,5\}$. We sampled from the SDE in Eq. 7 using the standard Euler-Maruyama method [9] and then calculated the mean and covariance of the samples at the end time of the simulation. On the other side we calculated the mean and covariance at the same end time using the ODEs in Eqs. 17 and 18. After generating $100$ experiments per landmark set and $\sigma$ value, the mean arrival locations for both methods differed by $.5\pm.02$ for $\sigma=.1$; $.1\pm.003$ for $\sigma=2$ and $.012\pm.0003$ for $\sigma=5$ all at least two orders of magnitude smaller than the radius of the datasets; the Frobenius norm of the difference between covariances was $11\pm 1$ for $\sigma=.1$; $3\pm.02$ for $\sigma=2$; and $.5\pm.01$ for $\sigma=5$ which is small in comparison with the original variance of the points $74$. This shows good agreement between the LL and the Euler-Maruyama methods.

3.2 Synthetic Registration Experiment

In order to compare our diffeomorphic model with a stochastic short deformation model [5], we implemented our model and then registered the landmarks in the circle shown in green in Fig. 1 to those of the “flower” shown in blue. The results for the short deformation model are illustrated in Fig. 1(a) and those of the diffeomorphic in Fig. 1(b). It is noticeable that in the short deformation model the domain has been warped into a non-invertible deformation which is not possible in the diffeomorphic case [3, 1]. We also show the uncertainty in the transform as modeled by the of the deformation field at each point. In Figs. 1(a) and 1(b) it is noticeable how, as expected, the uncertainty is lower close to the landmarks and it grows as we move far away from them. Moreover, in both models the FC values are comparable, showing that the increased complexity of the diffeomorphic model has not increased the uncertainty in the results.

3.3 Registration of Pre-operative and Intra-operative Images

We illustrate the strength of our method in the case of multi-modal registration. We use publicly available images [12] which include $12$ clinical cases of brain tumor resection. For these cases T1-MRI images have been acquired pre-surgically, manually annotated with between 20 and 37 anatomical landmarks and a tumor delineation and then intra-operative 3D ultrasound (US) reconstructions were acquired for the same subjects before tumor resection. The same experts annotated the US images with the same landmarks as the MRI.

We tested the accuracy of our registration algorithm on areas were there is no explicit information. For this, we used a leave-one-out (LOO) validation. For each subject we took one of the landmarks out, registered all others and then measured the distance of the landmark that we left out with the solution that were obtained by including it in the registration. We show the results in Fig. 3(a). The results are over $12$ subjects with between $20$ and $37$ landmarks per subject. We obtained the best results with $\sigma=50mm$. Priors with $\sigma<50mm$, were not able to move the left-out landmark to the ideal position and had increased variance. Priors, $\sigma>50mm$, had a closer distance to the ideal location but an increased number of outliers. Finally, we register these subjects using all the available landmarks and, through visual inspection, we are able the see that the deformable registration improves the image matching as shown in Fig. 2. Moreover, we also show how a prior enforcing a stronger smoothness constraint increases the certainty in of the registration in the whole image. We illustrate this in Fig. 3(b) where the increase of the low uncertainty (blue) area of the image correlates with the increase of the smoothness parameter.

4 Discussion and Conclusion

In this paper we presented a probabilistic diffeomorphic registration methodology. By extending the usual diffeomorphic model of [1] from a deterministic ODE formulation to a stochastic one, we were able to include in our model the registration error, or uncertainty. To the best of our knowledge, this is the first algorithm proposing a probabilistic diffeomorphic registration using a parametric density of the diffeomorphic deformations including a numerical method to calculate the parameters. Having presented our model, we devised an algorithm to implement it through a locally linear approximation to a parametric density. We successfully tested this approximation against usual methods for SDEs where a parametric density is not available. Then, we analyzed the performance of our algorithm in synthetic and human data. Our experiments showed that our algorithm produces good results. We measured this quantitatively through a LOO experiment as well as qualitatively by visual assessment of 12 registrations between MRI and US modalities.