Demian Wassermann

Matt Toews

Marc Niethammer

William Wells

Probabilistic Diffeomorphic Registration: Representing Uncertainty

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 framework is demonstrated in the case of landmark-based image registration, including simulated data and annotated pre and intra-operative 3D images.

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 [START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF][START_REF] Joshi | Landmark matching via large deformation diffeomorphisms[END_REF]. An important body of literature focuses on deformations in the form of diffeomorphisms [START_REF] Ashburner | Diffeomorphic registration using geodesic shooting and Gauss-Newton optimisation[END_REF][START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Joshi | Landmark matching via large deformation diffeomorphisms[END_REF], 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 [START_REF] Simpson | A Bayesian Approach for Spatially Adaptive Regularisation in Non-rigid Registration[END_REF].

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 [START_REF] Simpson | A Bayesian Approach for Spatially Adaptive Regularisation in Non-rigid Registration[END_REF], 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 [START_REF] Markussen | Large deformation diffeomorphisms with application to optic flow[END_REF]Alternatively, the posterior may be investigated via sampling methods, e.g. Markov chain Monte Carlo (MCMC) [START_REF] Risholm | Bayesian characterization of uncertainty in intra-subject non-rigid registration[END_REF] or Hamiltonian Monte Carlo [START_REF] Zhang | Bayesian estimation of regularization and atlas building in diffeomorphic image registration[END_REF].

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 0, 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.

Methods

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 Ω M and Ω F respectively, and let φ : Ω M → Ω F be a mapping between the two. The registration problem seeks a posterior probability density over mappings φ conditioned on data (M, F ), which is expressed via Bayes theorem as p(φ|M, F ) = p(φ)p(M, F |φ)/p(M, F ).

(

) 1
In Eq. (1), p(φ) is a prior density over φ embodying geometrical constraints such as smoothness. p(M, F |φ) is the data attachment factor or likelihood of the map φ relating F and M . E.g. the probability that M deformed by φ, which we note φ • M , is similar to F . Finally, p(M, F ) is a normalizing constant. The direct calculation of the posterior density p(φ|M, F ) is a difficult problem. Hence, we use a variational method to estimate a distribution q(φ) (abbreviated as q) that is close to p(φ|M, F ) in the sense of the Kullback-Leibler divergence (KL [• •]). Specifically, we seek q minimizing: KL [q p(φ|M, F )] = KL [q p(φ)] -log p(M, F |φ)dq(φ) + log p(M, F ).

(

) 2
In the registration literature the data attachment factor p(M, F |φ) is typically modeled using a measure of similarity between the registered objects: m : Ω F × Ω F → R which is minimal when two objects are exactly the same and grows as they become different. Adopting the equality -log p(M, F |φ) = m(φ • M, F ), Eq. ( 2) may be rewritten as:

KL [q p(φ|M, F )] = KL [q p(φ)] + m(φ • M, F ) q + log p(M, F ), (3)
where m(φ • M, F ) 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 [START_REF] Ashburner | Diffeomorphic registration using geodesic shooting and Gauss-Newton optimisation[END_REF][START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF]. First, instead of seeking a single optimal deformation φ, e.g. the maximum a-posteriori (MAP) solution in the Bayesian formulation, we seek to obtain the full distribution q(φ). In this way, we obtain both the MAP deformation φ in addition to the uncertainty at any given point in space, which can be calculated from q(φ). Second, we obtain q(φ) by minimizing the data attachment term over a weighted combination of all possible deformation fields in the family of φ instead of only at a single deformation φ.

Probabilistic Diffeomorphic Deformations

The variational approximation to p(φ|M, F ) described in the previous section requires a parameterization for q(φ) 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 [START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF], common to many diffeomorphic registration approaches [START_REF] Joshi | Landmark matching via large deformation diffeomorphisms[END_REF]. 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 [START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF], many diffeomorphic deformation formulations seek an optimal registration solution φ, e.g. the MAP deformation in Eq. ( 1), by constraining the map φ to be the solution at t = 1 of the ordinary differential equation (ODE)

d dt φ t (x) = v t (φ t (x)), φ 0 (x) = x, t ∈ [0, 1]. (4)
and setting φ 1 in Eq. ( 4) to minimize

E(φ; M, F ) = E v (v) + E φ (φ; M, F ), E v (v) = 1 2 1 0 Ωm Lv t (x) 2 2 dxdt. ( 5
)
From the two terms of E, E v (v) regularizes the evolution of the time-varying velocity field and E φ (φ; M, F ) drives E such that the deformed moving object, φ 1 • 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 [START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF], Eq. ( 5) restricts φ t (x) defined as in Eq. ( 4) to the space of diffeomorphisms [START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF].

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 v p(v) = GP(0, Σ ts (x, y)), Σ ts (x, y) ∈ SPD d (6)
where the covariance function Σ 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 v induces a random process φ with density p(φ) on the deformation field of Eq. ( 4) that is a solution of the stochastic differential equation (SDE)

dφ t (x) = v t (φ t (x))dt + Σ t (φ t (x))dW t , (7)
W t ∼ GP(0, Θ ts (x, y)), Θ ts (x, y) = min(t, s) Id, (8)
where √ Σ is the square root matrix

√ Σ √ Σ = Σ; the GP W t ∈ R d is called
Brownian motion or a Wiener process [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]; v t (x) ∈ R d is a deterministic velocity field like in Eq. (4); and Σ t (x) Σ tt (x, x) ∈ SPD d5 , the covariance of the probabilistic prior defined in Prop. 1, is a consequence of Eq. (5). Proposition 3. For the stochastic process φ with density p(φ), defined in prop. 2, the mean φ and covariance Λ functions are solutions of the deterministic ODEs

d dt φt (x) = v t (φ t (x)) p (9a) Λ ts (x, y) = cov p [φ t (x), φ s (y)] = φ t (x)φ s (y) p -φt (x) φ s (y) (9b) d dt φ t (x)φ s (y) p = v t (φ t (x))φ s (y) p + φ t (x)v s (φ s (y)) p (9c) + Σ ts (φ t (x), φ s (y)) p .
Moreover, up to a first order approximation:

p(φ) = GP( φt (x), Λ ts (x, y)), (10)
The proofs for Propositions 1-3 are as follows: Proof of Prop. 1 This proposition has been proven by Joshi et. al. [START_REF] Joshi | Landmark matching via large deformation diffeomorphisms[END_REF], 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 v, p(v):

-log(p(v)) = E v (v) + const = 1 2 1 0 Ω M Lv t (x) 2 2 dxdt + const, (11)
To show that p(v) is a stochastic process with a particular distribution, we need to prove that any finite sample of the domain Ω M ∈ R d has the same parametric distribution [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. We take N samples X ∈ R N ×d in space and t ∈ [0, 1] N in time, and let

V ij = [v ti (X i )] j ∈ R N ×d .
Then we rewrite Eq. ( 11) as

-log(p(V)) = 1 2 (L vec V) L vec V + const = 1 2 vec V L L vec V + const, ( 12
)
where L is the matrix such that 12) is recognisable as the log probability of a centered multivariate Gaussian with covariance C = (L L) -1 and, therefore v(x) is a GP. The covariance function Σ ts (x, y) can be calculated as the matrix Green's function of the operator L [START_REF] Joshi | Landmark matching via large deformation diffeomorphisms[END_REF]; specifically, if x, y ∈ R d , Σ ts (x, y) ∈ SPD d where [Σ ts (x, y)] ij is the covariance between x i at time t and y j at time s.

[LV i• ] = [Lv ti (X i )] j . Eq. (
This shows that for a given velocity field 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 W t ∈ R d is white noise, which multiplied by

√
Σ t is a centered Gaussian random variable with covariance Σ t , a sample drawn from Eq. ( 6). We noted the stochastic velocity field in Eq. ( 7) v to distinguish it from the deterministic one v. Eq. ( 7) ceases to be an ODE as a sample path of 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 φ [9, Chap. 8].

Proof of Prop. [START_REF] Ashburner | Diffeomorphic registration using geodesic shooting and Gauss-Newton optimisation[END_REF] The ODE for the mean of the stochastic process φ, 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 W t as a zero-centered Wiener process in Eq. [START_REF] Zhang | Bayesian estimation of regularization and atlas building in diffeomorphic image registration[END_REF].

To obtain the ODE for the second moment of φ, shown in Eq. (9c), we use the Itō product rule [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] to obtain an expression for d(φ t (x)φ s (y)) and substitute it in Eq. ( 7) obtaining

d(φ t (x)φ s (y)) p = d φ t (x)φ s (y) p = v t (φ t (x))dtφ s (y) p + √ Σ t (φ t (x))dW t φ s (y) p + φ t (x)(v s (φ s (y))) ds p + φ t (x)( √ Σ t (φ s (y))dW s ) p + v t (φ t (y))dt + √ Σ t (φ t (y))dW t v s (φ s (y))ds + √ Σ s (φ s (y))dW s p , (13
) which, using the Itō identities for expected values of differentials [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] results in Eq. (9c).

Obtaining a parametric form of the density of φ, p(φ), satisfying the SDE [START_REF] Risholm | Bayesian characterization of uncertainty in intra-subject non-rigid registration[END_REF] 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 √ Σ are linear functions on their time and location parameters, and the initial condition φ t=0 is a GP, φ t (x) is known to be a GP [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]. With this purpose we define a locally linearized (LL) v and √ Σ centered at t 0 , x 0 [START_REF] Biscay | Local linearization method for the numerical solution of stochastic differential equations[END_REF]:

v t (x) ≈ v t0 (x 0 ) + ∂ t v t0 (x 0 )(t -t 0 ) + D x vt 0 (x 0 )(x -x 0 ) (14a) √ Σ t (x) ≈ √ Σ t0 (x 0 ) + ∂ t √ Σ t0 (x 0 )(t -t 0 ) + i ∂ xi √ Σ t0 (x 0 )(x -x 0 ) i (14b)
where D x vt is the Jacobian of v t (x) w.r.t. x and ∂ t v t its partial derivative w.r.t. t. Considering that L is assumed time-invariant in Prop. 1, the time derivative of √ Σ in Eq. (14b) is equal to 0. Then, using the LL equations Eqs. (14a) and (14b), we approximate Eq. ( 7) as

dφ t (x) ≈ (A t φ t (x) + a t )dt + i S i t φ t (x) i + R t dW t (15) A t D x vt 0 (x 0 ) a t -D x vt 0 (x 0 )x 0 + ∂ t v t0 (x 0 )(t -t 0 ) + v t0 (x 0 ) S i t ∂ xi √ Σ t0 (x 0 ) R t - i S i t • (x 0 ) i + √ Σ t0 (x 0 ). (16)
The LL approximations in Eqs. (14a) and (15) lead to an approximation of the mean function of φ, φ, by the solution of the ODE

d φt(x) dt ≈ v t ( φt (x)) where v t (x) ≈ (A t φt (x) + a t ), (17)
and its second moment φ t (x)φ s (y) when t = s by

d φ t (x)φ t (y) dt ≈ A t φ t (x)φ t (y) + φ t (x)φ t (y) A t + a t φ t (y) + φt (x)a t + ij S i t φ t (x)φ t (y) S j t + i S i t φt (x) i R t + R t i S i t φt (y) i + R t R t , (18)
where A t ; a t ; S i t ; and R 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 φ t=0 is GP, the linear approximation of φ t is a GP uniquely determined by φ and Λ [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]. Then, given a set of stochastic velocity fields v 0 . . . v t M -1 with t 0 = 0 and t M -1 = 1, the parameters of the stochastic process representing the transform φ are obtained integrating Eqs. ( 17) and (18) with the initial conditions φ t=0 ∼ GP( φt=0 , Λ t=0 ). Having characterized stochastic transformations representing a diffeomorphic deformation, we are in position to formulate our probabilistic diffeomorphic registration algorithm.

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 Σ; we use this model as q(φ 1 ), our variational distribution. In this section, we show how to compute the parameters of q(φ 1 ) minimizing Eq. ( 2) for a particular registration problem. Taking the approach of [START_REF] Ashburner | Diffeomorphic registration using geodesic shooting and Gauss-Newton optimisation[END_REF], 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 φ derived with Prop. 1 is the joint probability of the fields at each time t: p(φ 1 ) = 1 0 p(v t ) dt with p(v t ) ∼ GP(0, Σ 0 ). Then, we rewrite leftmost term of Eq. ( 2) as KL [q p(φ 1 )] = 1 0 KL [q(v t )|p(v t )] dt. We parameterize each stochastic velocity field v t by a N -point set represented as a matrix X t ∈ R N ×d rendering its mean equivalent to a spline model [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. This sets the distributions of the discretized velocity field prior to p(vec v t |X t ) = 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(vec v t |X t ) = G(µ 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, µ t (X) = vec v t (X) and [S t=0 (X)] di+k,dj+l = [Σ t=0 (X i , X j )] kl , i, j = 1 . . . N, k, l = 0 . . . d -1 [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. This leads to an objective which we minimize to obtain q(φ 1 ) representing the registration problem:

E(q(φ 1 )) = KL [q(φ 1 ) p(φ 1 )] + m(φ 1 • M, F ) q + log p(M, F ) (19a) KL [q(φ 1 ) p(φ 1 )] = 1 0 µ t S -1 t=0 µ t dt, s.t. d φt dt ≈ µ t ( φt ), φ0 = id . (19b)
Using the ideas of [START_REF] Ashburner | Diffeomorphic registration using geodesic shooting and Gauss-Newton optimisation[END_REF], Eq. ( 19) can be minimized through geodesic shooting [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF], i.e. it depends only on M , F and µ 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 µ. In fact, the problem formulation equations, shown in Eq. ( 4), stay the same as in [START_REF] Ashburner | Diffeomorphic registration using geodesic shooting and Gauss-Newton optimisation[END_REF], only the final condition E φ changes, which is then warped to t = 0 for a gradient descent with respect to the initial velocity µ t=0 leading to the objective function of q(φ 1 ) parameterized on µ t=0

argmin µt=0 E(q µt=0 (φ)) = 1 2 µ t=0 S -1 t=0 µ t=0 + m(φ 1 • M, F ) q -log p(M, F ). (20)
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 R N ×d ; m(M, F ) = M -F 2 2 ; and the random variable Φ t = vec(φ t • M ) vec φ t (M ). This allows us to rewrite

m(φ 1 • M, F ) q = m(Φ 1 , F ) q = tr Φ 1 Φ 1 q -2 Φ 1 F + tr F F . (21)
Replacing Eq. ( 21) in Eq. (20) leads to the gradient

∇ µt=0 E(q µt=0 (φ 1 )) = 1 2 S -1 t=0 µ t=0 -2 Φ1 -2F .
Having this gradient, we minimize E w.r.t. µ t=0 using a gradient descent algorithm.

Experiments

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

Σ ts (x, y) ij = δ(t -s) exp - xi-yj 2 2 2σ 2 ij ∈ R d×d (22)
where σ 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 φ 0 (x) ∼ GP(x, Σ t=0,s=0 (x, y)) and Σ t=0,s=0 (x, y) = 0.

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 φ [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]. 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 σ ∈ {.1, 2, 5}. We sampled from the SDE in Eq. ( 7) using the standard Euler-Maruyama method [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] 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 σ value, the mean arrival locations for both methods differed by .5±.02 for σ = .1; .1±.003 for σ = 2 and .012±.0003 for σ = 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 ± 1 for σ = .1; 3 ± .02 for σ = 2; and .5 ± .01 for σ = 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.

Synthetic Registration Experiment

In order to compare our diffeomorphic model with a stochastic short deformation model [START_REF] Simpson | A Bayesian Approach for Spatially Adaptive Regularisation in Non-rigid Registration[END_REF], 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 [START_REF] Ashburner | Diffeomorphic registration using geodesic shooting and Gauss-Newton optimisation[END_REF][START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF]. 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.

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 [START_REF] Mercier | Online database of clinical MR and ultrasound images of brain tumors[END_REF] which include 12 clinical cases of brain tumor re- section. 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 σ = 50mm. Priors with σ < 50mm, were not able to move the left-out landmark to the ideal position and had increased variance. Priors, σ > 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.

Discussion and Conclusion

In this paper we presented a probabilistic diffeomorphic registration methodology. By extending the usual diffeomorphic model of [START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF] 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

Fig. 1 -

 1 Fig.1-Comparison between small deformation and diffeomorphic registrations with equal parameter values, the uncertainty is represented by the Frobenius norm of the covariance. The small deformation model has a smaller variance in general at the expense of a possibly invalid deformation field away from the landmarks.

Fig. 2 -

 2 Fig. 2 -Registration of pre-operative and intra-operative images: (a) The pre-operative MRI linearly registered and projected onto the ultrasound space. (b) The intra-operative ultrasound image; and (c) the pre-operative MRI of (a) registered to (b) using our algorithm were we show the warping according to the average registration field. The crosshair indicator shows how correspondence between a-b is not as accurate as the one using deformable registration (c).

SP D d : set of symmetric positive definite matrices of dimension d

This work was supported by grants NIH P41EB015898, R01CA138419 and NSF EECS-1148870.

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.