Loïc Le Folgoc

Hervé Delingette

Antonio Criminisi

Nicholas Ayache

Sparse Bayesian registration of medical images for self-tuning of parameters and spatially adaptive parametrization of displacements

Keywords: Non-rigid registration, Bayesian modelling, Sparse structured prior, Variational Bayes, ARD, Cardiac Imaging

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Non-rigid image registration is the ill-posed task of inferring a deformation Ψ from a pair of observed (albeit noisy), related images I and J. Classical approaches propose to minimize a functional which weighs an image similarity criterion D against a regularizing (penalty) term R:

arg min Ψ E(Ψ) = D(I, J, Ψ) + λ • R(Ψ) (1)
Prior knowledge to precisely model the space of plausible deformations or the regularizing energy is generally unavailable.

The optimal trade-off between the image similarity term and the regularization prior is itself difficult to find. Typically the user would manually adjust the trade-off until a qualitatively good fit is achieved, which is time consuming and calls for some degree of expertise. Alternatively if quantitative benchmarks are available on a similar set of images, they can serve as a metric of reference on which to optimize parameters, under the assumption that the value that achieves optimality is constant across the dataset. Unfortunately, this assumption generally does not hold. Probabilistic interpretations of registration recently emerged as a way to automate the process [START_REF] Richard | A SAEM algorithm for the estimation of template and deformation parameters in medical image sequences[END_REF][START_REF] Simpson | Probabilistic inference of regularisation in non-rigid registration[END_REF][START_REF] Risholm | Bayesian characterization of uncertainty in intra-subject non-rigid registration[END_REF]. [START_REF] Gee | Elastic matching: Continuum mechanical and probabilistic analysis[END_REF] first noted that, in a Bayesian paradigm, the two terms in Eq. ( 1) relate respectively to a likelihood and prior on the latent transformation Ψ. In fact the trade-off parameter itself can be treated as a hidden random variable, equipped with a broad prior distribution, and jointly inferred with Ψ or integrated out. In practice, analytical inference is precluded and various strategies are devised for approximate inference. [START_REF] Risholm | Bayesian characterization of uncertainty in intra-subject non-rigid registration[END_REF] characterize the distributions of interest from MCMC samples. This is a most principled and accurate approach provided that enough samples can be drawn within the available computational bud- get. Aside from monitoring the progress of the scheme, two difficulties arise: crafting an efficient proposal distribution over Ψ and computing the acceptance probability of the proposed sample. To circumvent this latter issue, the authors sample from an approximate posterior distribution derived in a variational free-energy framework. Alternatively, the full inference can be tackled in a variational Bayes framework [START_REF] Simpson | Probabilistic inference of regularisation in non-rigid registration[END_REF][START_REF] Simpson | Probabilistic non-linear registration with spatially adaptive regularisation[END_REF]. This offers an appealing compromise between the computational burden and the quality of the estimates, depending on the chosen family of variational (approximate) posterior distributions. In this article, we propose to extend the Bayesian framework of registration to automatically select the optimal location and scale of bases parametrizing the transformation. Spatial refinement of the parametrization was previously handled heuristically [START_REF] Rohde | The adaptive bases algorithm for intensity-based nonrigid image registration[END_REF], or led to alternative formulations of registration via spatially anisotropic filtering [START_REF] Stefanescu | Grid powered nonlinear image registration with locally adaptive regularization[END_REF]. Dynamic refinements of the displacement space have also been proposed by [START_REF] Glocker | Optical flow estimation with uncertainties through dynamic mrfs[END_REF]; [START_REF] Parisot | Concurrent tumor segmentation and registration with uncertainty-based sparse non-uniform graphs[END_REF] for MRF-based discrete registration. The displacement quantization is refined using local, min-marginal based estimates of uncertainty. Dynamic quantization (see also [START_REF] Tang | Random walks with efficient search and contextually adapted image similarity for deformable registration[END_REF] and [START_REF] Heinrich | Deformable image registration by combining uncertainty estimates from supervoxel belief propagation[END_REF]) does not affect the registration energy however: its purpose is simply to accelerate convergence towards the optimum. In our work, the registration cost function forces the complexity of the mapping to adapt to the underlying dataset: finer bases are introduced only in the presence of coherent image information and motion, while coarser bases ensure better extrapolation of the motion to textureless, uninformative regions. In that spirit of model selection, [START_REF] Stewart | The dual-bootstrap iterative closest point algorithm with application to retinal image registration[END_REF] use an information criterion to choose regionally among a limited pool of deformation models (e.g. similarity, affine, quadratic), but do not address combinatorial issues arising in fully non-rigid registration from the number of possible parametrizations. [START_REF] Shi | Temporal sparse free-form deformations[END_REF] couple sparse optimization with a multiscale free-form representation of deformations, demonstrating gains in registration accuracy.

Here and to our knowledge, for the first time, basis selection in registration is approached on principled grounds within a probabilistic framework.

We propose a Bayesian model of registration that allows to automatically infer from the data the optimal parametrization of displacements, along with all model parameters. The inference scheme is efficient and tractable for real scale non-rigid registration tasks. The model and inference strategy are based on the Relevance Vector Machine [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF][START_REF] Tipping | Fast marginal likelihood maximisation for sparse bayesian models[END_REF], a generic approach to sparse regression and classification. To make it suitable for registration, where smooth solutions are looked for, we extend it to richer Gaussian priors with arbitrarily structured covariance, at no cost in algorithmic complexity. We also generalize the approach to multivalued regression (regression of vector fields), so as to preserve the natural invariance of the problem to changes of coordinate frames.

This article expands on earlier work of the authors (Le [START_REF] Folgoc | Sparse bayesian registration[END_REF] in several ways. Inference is fully presented within a variational Bayes framework. We propose a different approximation of the likelihood term, effectively removing a computational bottleneck: the voxelwise, local optimization of the image similarity via dense block-matching. It is replaced by a step where the registration energy is optimized w.r.t. the reduced parametrization. Finally we introduce a flexible noise model that is more robust to acquisition noise and artefacts, adapting over a range of image modalities.

Section 2 describes the statistical model of pairwise registration. The inference strategy is exposed in section 3. Section 4 reports experiments on tasks of motion tracking on real cardiac data, specifically time sequences of 3D cine or tagged MR images and echocardiographic images.

Statistical Model of Registration

Image registration assumes images to relate via some transformation of space such as, in a medical context, when imaging the motion of organs throughout a sequence of time frames. Registration then aims at recovering the unknown transformation of space from the observed data, which is formally an inference problem. We now specify a generative model of the data given the transformation Ψ, along with a sparse structured prior over the admissible set of transformations. Fig. 2 provides a graphical depiction of the model.

The sparse structured prior model was previously proposed by [START_REF] Sabuncu | The relevance voxel machine (rvoxm): A self-tuning bayesian model for informative image-based prediction[END_REF] for image-based tasks of classification and regression. Inference followed the guidelines of [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF], and was later accelerated for specific priors exploiting sparsely connected graphs [START_REF] Ganz | An improved optimization method for the relevance voxel machine[END_REF]. In Section 3 we develop alternative inference schemes that are applicable with no restriction: the structured part of the prior, irrelevant to the algorithmic complexity, may be arbitrarily defined. The gain in algorithmic complexity reflects the way in which the later work of [START_REF] Tipping | Fast marginal likelihood maximisation for sparse bayesian models[END_REF] accelerates the original Relevance Vector Machine [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF]. This ef- fectively renders the approach applicable to non-rigid registration.

Data Likelihood

A good transformation Ψ should adequately map the datasets, up to some misalignment and residual error attributable to the data formation process. The knowledge of this process is captured in a likelihood model, which assigns a probability p(D|Ψ; P) for the data D to be observed under some transformation Ψ (often conditioned on a set of hyperparameters P). The likelihood typically assumes the form of a Boltzmann distribution:

p(D|Ψ; P) ∝ exp -D(D, Ψ; P) , (2)
which explicitely bridges the gap with the classic optimization framework of Eq. (1). For pairwise registration of images, the simplest and most common image similarity term is the sum of squared difference (SSD) of voxel intensities, which can be improved upon by modeling spatially varying noise levels [START_REF] Simpson | A bayesian approach for spatially adaptive regularisation in non-rigid registration[END_REF] and artefacts [START_REF] Hachama | Bayesian technique for image classifying registration[END_REF], or by relaxing assumptions over the intensity mapping between imagese.g. to a piecewise constant mapping [START_REF] Richard | A SAEM algorithm for the estimation of template and deformation parameters in medical image sequences[END_REF], to a locally affine mapping [START_REF] Cachier | Iconic feature based nonrigid registration: the pasha algorithm[END_REF] or to a more complex, non-linear (Parzen-window type) intensity mapping [START_REF] Janoos | Bayesian characterization of uncertainty in multi-modal image registration[END_REF]. Mutual information is another popular image similarity, especially in the context of registering images of different modalities [START_REF] Wells | Multimodal volume registration by maximization of mutual information[END_REF], and has been successfully applied to the registration of cardiac images [START_REF] Chandrashekara | Analysis of 3-d myocardial motion in tagged mr images using nonrigid image registration[END_REF]. SSD is a simple yet efficient image similarity term for registration of monomodal cardiac images. It naturally lends itself to a probabilistic interpretation and eases mathematical derivations. The target image J is modeled as the warped source image I • Ψ -1 further corrupted by additive, independent identically distributed (i.i.d.) noise e i ∼ N(0, β) at each voxel

i = 1 • • • N: J = I • Ψ -1 + e (3)
where e ∼ N(0, βI), I the N × N identity matrix. β is a global scaling parameter: it stands for the inverse variance (precision) of the noise across the image. The SSD model can be described in a more familiar manner by the energy of Eq. ( 4), where

{v i } N i=1
is the list of voxel centers in the fixed image and

V i = Ψ -1 (v i )
are the paired coordinates in the moving image.

D β (J; I, Ψ) = β 2 N i=1 J v i -I V i 2 (4)
Since the SSD is quadratic w.r.t to intensity differences of paired voxels, both the penalty for intensity discrepancies and the rate at which it grows can become arbitrarily high. As seen in Fig. 3, this renders registration vulnerable to strong local intensity biases, introduced for instance by acquisition artefacts or by topology changes in the imaged objects. In addition residual misalignments between structures of interest tend to yield higher intensity residuals than those observed at background voxels (see for instance Fig. 4a). Sources of model bias and acquisition noise cannot be captured together in a plausible manner with a single, spatially uniform noise level. In other words, the SSD noise model is neither robust nor flexible enough.

To address this limitation we propose to model the noise

e i ∼ 1≤l≤L π l N(0, β l ) at each voxel i = 1 • • • N with a mix- ture of L Gaussian distributions. Implicitely at each voxel, the residue e i = J[v i ] -I[V i ]
is independently assigned to one of the L components, with π l the probability of being assigned to the lth component N(0, β l ). Introducing a set of binary assignment variables {z i1 • • • z iL } for the ith voxel, such that z il = 1 if assigned to the lth component and z il = 0 otherwise, the above can be summarized as: with

p(J|I, Ψ, β, z) ∝ N i=1 exp - 1 2 ( L l=1 z il β l ) J v i -I V i 2 (5) p(z|π) ∝ N i=1 L l=1 π z il l (6)
β = {β 1 • • • β l }, π = {π 1 • • • π l } and z = {z il } i=1•••N,l=1•••L .
This yields a spatially varying model of noise that is better suited to render the complexity of noise patterns in medical images. Unlike in previous work [START_REF] Simpson | A bayesian approach for spatially adaptive regularisation in non-rigid registration[END_REF][START_REF] Folgoc | Sparse bayesian registration[END_REF] the noise here is not assumed to vary smoothly across the image, as patterns arising from misalignment and imaging artefacts are local in nature. Integrating over assignment variables, we explicitely retrieve the mixture-of-Gaussian structure:

p(J|I, Ψ, β, π) = N i=1 L l=1 π l Z l exp - β l 2 J v i -I V i 2 (7)
where Z l = 2π/β l stands for the normalizing constant for the Gaussian probability distribution function. The corresponding data matching energy is given in Eq. ( 8):

D β,π,L (J; I, Ψ) = - N i=1 log L l=1 π l Z l exp - β l 2 J v i -I V i 2 (8)
Fig. 4b shows the histogram of intensity residuals for the example registration of Fig. 3, along with the learned Gaussian mixture (jointly fit during registration). The profiles of the standard SSD loss and the Gaussian mixture (GMM) loss are displayed in Fig 4c . The characteristic inflexion of the GMM loss, with a reduced growth rate as the intensity residual becomes higher, is responsible for its robustness towards intensity artefacts compared to the standard SSD quadratic loss. Mixtures of Gaussian (or Student-t) distributions have long been used as building blocks for robust autoregressive models [START_REF] Roberts | Variational bayes for generalized autoregressive models[END_REF][START_REF] Tipping | Variational inference for student-t models: Robust bayesian interpolation and generalised component analysis[END_REF], including in medical imaging [START_REF] Penny | Robust bayesian general linear models[END_REF] but remain uncommon for registration [START_REF] Leventon | Multi-modal volume registration using joint intensity distributions[END_REF].

A limitation of SSD shared by all aforementioned variants is to assume that voxelwise intensity residuals are independent. This assumption does not hold [START_REF] Simpson | Probabilistic inference of regularisation in non-rigid registration[END_REF]. In practice, the residual between the warped image I • Ψ -1 and its counterpart J exhibits local spatial correlations, either intrinsic to the image acquisition and pre-processing (e.g. image presmoothing, image upsampling) or introduced as a consequence of registration misalignments. Ignoring local correlations in the noise pattern leads to an artificial increase in the number of independent observations and induces over-confidence in the data term. On the other hand, modeling precisely the noise structure would come at a significant computational cost. Here, we follow [START_REF] Simpson | Probabilistic inference of regularisation in non-rigid registration[END_REF] in artificially downweighting the data term by a factor α that captures redundancies in voxelwise observations, based on a virtual decimation procedure suggested by [START_REF] Groves | Linked independent component analysis for multimodal data fusion[END_REF].

Representation of displacements

We proceed in a small deformation framework, Ψ -1 = Id + u, with a parameterized representation of the displacement field u : x ∈ R d → u(x) ∈ R d . The displacement field is expressed over a dictionary {φ k } M k=1 of Gaussian radial basis functions, φ k (x) = K S k (x k , x) I where I is the d × d identity matrix and

K S (x, y) = exp - 1 2 (x -y) S -1 (x -y) . (9)
In other words, the displacement field u is parametrized by a set of weights w k ∈ R d associated to each basis φ k :

u w (x) = 1≤k≤M φ k (x)w k = φ(x)w . (10) φ(x) = φ 1 (x) • • • φ M (x) and w = w 1 • • • w M are respectively the concatenation, for k = 1 • • • M, of φ k (x) and w k .
The basis centers x k span a predefined regular grid of points, typically the whole range of voxel centers. The kernel width S k is also allowed to vary and spans a user-predefined set of values S 1 , S 2 , • • • , S q . This yields a redundant, multiscale representation of displacements. Larger kernels make the representation more compact, whereas smaller kernels allow to capture finer local details. The genericity of the approach w.r.t. the choice of dictionary is discussed in AppendixH.

Transformation Prior

In non-rigid registration, the displacement u w is insufficiently constrained by the data and some regularizing prior has to be imposed over its parameters. This prior distribution encapsulates our knowledge of the deformation and our modeling assumptions (see for instance [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF] for an exhaustive review of deformation priors). We will consider Gaussian priors of the form

p(w|λ, {A k }) ∝ exp - 1 2        λw Rw + M k=1 w k A k w k        (11)
where λ and {A k } k=1•••M are model parameters. The motivation for such a prior is two-fold.

Regularity control. Gaussian priors in the form of Eq. ( 12) let us penalize physically implausible deformations. They have been commonly used in the literature starting with [START_REF] Broit | Optimal registration of deformed images[END_REF], both because of their natural interpretability and soundness in mechanical terms, and their convenience from an algorithmic and computational standpoint.

q(w|λ) ∝ exp - 1 2 λw Rw (12)
The structure of the precision matrix R can be adjusted to penalize the magnitude Du 2 of the first derivative of the displacement field [START_REF] Gee | Elastic matching: Continuum mechanical and probabilistic analysis[END_REF] or higher order derivatives [START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast mr images[END_REF][START_REF] Ashburner | A fast diffeomorphic image registration algorithm[END_REF][START_REF] Ashburner | Symmetric diffeomorphic modelling of longitudinal structural MRI[END_REF], effectively encoding a wide range of priors. We recall in AppendixH how to compute R efficiently using Fourier analysis for general families of basis functions instead of relying on costly numerical integration. With the parametrization of displacements given in section 2.2, classical energies are in fact implemented in closed-form. In this work, we specifically rely on a bending (thin-plate) energy (D = ∇ 2 ).

Basis selection. The second factor in our prior, recalled in Eq. ( 13), induces the desired basis selection mechanism. The additional term w k A k w k for each basis φ k lets us penalize independently the recourse to this basis to capture the displacement, by penalizing high magnitudes of its associated weight w k . Each A k is an arbitrary d × d symmetric positive matrix, so that w k can be penalized in a different manner along different orientations. The improper limit case of infinite A k actually constrains w k to be null and thus forbids the use of φ k to represent the signal. In section 3 we determine optimal values of the set {A k } k=1•••M in a principled manner, from which most of them turn out to be infinite: we thus obtain a sparse representation of the displacement from the initial, over-complete dictionary. Introducing A diag(A 1 • • • A M ) as the block diagonal matrix whose kth d × d diagonal block is A k , the full prior takes the more compact form p(w|λ, A) ∝ exp -1 2 w (A + λR)w.

q(w|{A k }) ∝ m k=1 exp - 1 2 w k A k w k (13)

Hyperpriors

The value of model parameters β, π, λ and {A k } k=1•••M is unknown. We regard them as additional model variables and endow them with prior distributions. When possible, the choice of conjugate priors facilitates inference. The noise levels β = {β 1 • • • β L } for each component of the Gaussian mixture are assigned independent Gamma priors Γ(

β l |c 0 , d 0 ) ∝ β c 0 -1 l e -d 0 β l . The noise mixture proportions π = {π 1 • • • π L } are equipped with a Dirichlet prior Dir(π|η 0 ) ∝ l π η 0 -1 l .
λ is endowed with a Gamma prior Γ(λ|a 0 , b 0 ) ∝ λ a 0 -1 e -b 0 λ . In absence of strong prior knowledge, broad uninformative priors can be chosen

(a 0 = b 0 = c 0 = d 0 = η 0 → 0).
An improper uniform prior is taken over basis penalties A k , with the added benefit of making inference invariant to rescaling of basis functions. Moreover given the inference strategy of section 3, AppendixB and AppendixE, optimality conditions state that

A -1 k = α -1 k n k n k is at most rank-one, with α k ∈ R + ∪ {+∞}.
We found advantageous to further restrain α k to be either 0 or +∞. In other words, along any given direction, A k either constrains w k to be null or does not constrain it whatsoever. This prevents direct competition between the two regularization mechanisms of Eq. ( 12) and (13). Moreover, the Gamma prior over λ then becomes conjugate to p(w|A, λ).

Related work: Sparse Coding & Registration

Sparsity-inducing priors have a two-fold motivation. The first benefit is in terms of algorithmic complexity. Unless resorting to low parametric models, the size of the parametrization makes direct optimization cumbersome without the recourse to sophisticated solvers. The computation of exact covariance matrices that are typically involved in probabilistic approaches also becomes unfeasible, while diagonal approximations used in their stead discard significant interactions induced by the data and priors. Secondly, basis selection mechanisms adaptively constrain the space of deformations, automatically tuning the degrees of freedom to the smallest set sufficient to capture the observed displacement. Coupled with a multiscale set of basis functions, this yields a data-driven, automatic spatial refinement of the granularity of the displacement field that complements the otherwise scale-blind L 2 regularization. Adaptive, multiscale regularization was shown to yield stateof-art results e.g. in denoising natural scenes [START_REF] Fanello | Filter forests for learning data-dependent convolutional kernels[END_REF], but also in medical image registration [START_REF] Shi | Temporal sparse free-form deformations[END_REF]. Fig. 6 gives a naive insight into the key impact of scale when limited data is available.

L 1 priors have been widely used in all areas of sparse coding, including for registration [START_REF] Shi | Temporal sparse free-form deformations[END_REF]. Other sparsityinducing norms such as k-support norms and variants [START_REF] Argyriou | Sparse prediction with the k-support norm[END_REF]Belilovsky et al., 2015a), that improve over the performance of the L 1 norm w.r.t. the degree of sparsity in presence of strongly correlated explanatory variables, have recently been proposed. They were shown to be attractive on tasks of functional MR imaging [START_REF] Jenatton | Multi-scale Mining of fMRI data with Hierarchical Structured Sparsity[END_REF]Belilovsky et al., 2015b). Here, we turn instead towards sparse Bayesian learning, with the prospect of joint estimation of model parameters and that of uncertainty quantification. For an extensive review of sparse methods, we refer the reader to the work of [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF], and to that of [START_REF] Mohamed | Bayesian and l1 approaches for sparse unsupervised learning[END_REF] for a benchmark of L 1 and bayesian sparse learning methods. The prior of Eq. ( 13) was first introduced by [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF] for regression and classification tasks with the so called Relevance Vector Machine. The authors demonstrated its relevance for sparse coding when used in conjunction with the framework of Automatic Relevance Determination [START_REF] Mackay | Bayesian interpolation[END_REF]. [START_REF] Bishop | Variational relevance vector machines[END_REF] offer an alternative sparse Bayesian learning (SBL) view on the Relevance Vector Machine, where they opt for a Variational Bayes treatment. [START_REF] Wipf | A new view of automatic relevance de-termination[END_REF] further investigate links between the SBL and ARD frameworks and resulting schemes. Alternatively, Eq. ( 11) can be interpreted as a generalized spike-and-slab prior [START_REF] Mitchell | Bayesian variable selection in linear regression[END_REF] despite using a different parametrization, provided that each A k is constrained to a binary state -either null or infinite.

Model Inference

Bayesian inference summarizes both prior knowledge p(θ) and data-driven information p(D|θ) on model parameters and hyperparameters θ = {w, π, β, z, λ, A} within a posterior distribution

p(θ|D) = p(D|θ) • p(θ) p(D) . ( 14
)
The goal of Bayesian inference is thus to characterize the joint posterior p(θ|D) or to characterize marginals of interest, such as the marginal posterior distribution p(w|D) of transformation parameters for the purpose of registration. Exact inference is precluded and we proceed in the framework of variational Bayes (VB) inference [START_REF] Bishop | Pattern recognition and machine learning[END_REF].

Variational Bayes inference

VB inference approximates the true posterior p(θ|D) among a restricted family of variational posterior distributions q(θ) that benefits analytical and computational derivations. The objective under VB inference is to minimize the Kullbach-Leibler divergence KL(q(θ) p(θ|D)) or equivalently to maximize a lower bound L(q) of the log-evidence. This equivalence follows from the following identity, where on the left-hand side the evidence p(D) for the model is constant w.r.t. q:

log p(D) = KL q(θ) p(θ|D) - θ q(θ) log q(θ) p(θ, D) dθ L(q) (15)
When the true posterior lies within the variational family, the (non negative) Kullbach-Leibler divergence is minimized (zero) when q * (θ) = p(θ|D). In practice however, the choice of variational family q reflects a trade-off between the accuracy of the approximation q * (θ) and its actual tractability.

The mean-field approximation assumes q to factorize over subsets of model variables. In our case, we consider variational distributions for which the transformation parameters and individual penalties, the regularization level, the noise levels, the mixture proportions and the voxel assignments factorize: q(θ) = q w,A (w, A)q λ (λ)q β (β)q π (π)q z (z) .

(16)

Let q I be any one of the individual factors and q -I the product of remaining factors, e.g. q λ and q w,A (w, A)q β (β)q π (π)q z (z) respectively. Let θ I and θ -I denote corresponding subsets of variables within θ. Exploiting the factorization, each factor q I (θ I ) ≈ p(θ I |D)1 can be seen to give an approximation of a given marginal of interest. The optimum q * among variational posteriors compatible with this factorization is known from calculus of variations to satisfy, for each individual factor q I :

log q * I = log p(θ, D) q * -I + const . (17)
Here, • q -I denotes expectation w.r.t. q -I (θ -I ). Eq. ( 17) naturally suggests inference schemes that update each factor q I in turn until convergence, guaranteeing decrease of the objective L(q) at each iteration. Moreover when the Bayesian model uses conjugate exponential distributions, mean-field VB updates are considerably simplified. Each factor q * I (θ I ) lies in the same exponential family as the corresponding prior p(θ I |θ -I ) so that VB inference resolves into much more practical updates of the exponential distribution parameters. For instance, the optimal variational posterior for λ is again Gamma distributed, q * λ (λ) = Γ(λ|a, b), with closed form expressions for a, b. In the proposed model, most conditional probabilities do belong to conjugate exponential families. One exception is the likelihood p(D|w, β, z). To enable analytical derivations of q β , q z and q w,A , a Gaussian approximation of the likelihood is used (section 3.4). In addition the family of variational posteriors q w,A (w, A) is further restricted (section 3.2).

Constraining q w,A (w, A) for fast sparse Bayes inference

From the factorization of Eq. ( 16) and optimality conditions of Eq. ( 17), q * w,A (w, A) unfortunately does not simplify into a convenient distribution. Without loss of generality yet, q w,A (w, A) = q(w|A)q A (A), where q A (A) w q w,A (w, A) dw is a variational approximation to p(A|D). We propose to constrain the variational posterior q A (A) to take the form of a Dirac distribution q A (A) δ Â(A) with all its mass assigned at the value A k = Âk for k = 1 . . . M, so that q w,A (w, A) = q(w| Â)δ Â(A). Under this assumption, the optimum

q * w,A = arg max Â, q(w| Â) L(q) q * (w|A * )δ A * (A) (18)
can be derived by calculus of variations in two steps. Given any Â, q * (w| Â) satisfies optimality conditions similar to Eq. ( 17):

log q * (w| Â) = log p(θ -w,-Â, w, Â, D) q -w,-Â + const . ( 19
)
Reinjecting this expression into Eq. ( 18) turns it into a maximization w.r.t. Â only, so that the maximizer A * of L(q) can be shown to maximize the following quantity:

A * = arg max  w exp log p(θ -Â, Â, D) q -Â,-w dw (20) = arg max  p(  | D, λ q λ , β q β , z q z ) , (21)
= arg max  p(D | Â, λ q λ , β q β , z q z ) . (22)
Eq. ( 21) uses the fact that p(w|A, λ) and p(D|w, β, z) belong to exponential families. Eq. ( 22) follows from Bayes' rule with the improper prior p(A) ∝ 1, and shows that maximizing L(q) w.r.t. Â is the same as maximizing the conditional evidence p(D | Â, λ q λ , β q β , z q z ). In fact to increase the value of the objective L(q) w.r.t. q w,A we merely need to increase (not necessarily maximize) the conditional evidence w.r.t. Â, then update q(w| Â) according to the optimality condition of Eq. ( 19).

Based on this remark, we derive an active set method that greedily improves on the objective functional L(q). The active set refers to the subset S of basis functions φ k for which Âk is finite along at least a direction, as opposed to the inactive set of basis functions for which Âk is infinite and constrains w k = 0. The scheme starts with an arbitrary active set (typically S = ∅) and proceeds by updating one Âk at a time, maximizing the quantity of Eq. ( 22)) w.r.t. this basis function only. This results in adding a new basis to the active set if Âk is made finite along at least a direction, or removing a previously active basis if Âk becomes infinite. The hyperparameter Âk that is selected for an update is the one, among all indices k = 1 . . . M, that provides the highest gain w.r.t. the objective. AppendixB shows that, in the case of a Gaussian likelihood, all necessary updates can be performed efficiently using rank-1 linear algebra identities.

Related work

In a simplified setting (λ = 0), [START_REF] Bishop | Variational relevance vector machines[END_REF] use the factorization q w,A (w, A) = q w (w)q A (A) to derive closedform updates for all factors. Unfortunately, the smoothnessinducing prior destroys model conjugacies on which the authors rely. In addition, the resulting updates have a complexity O(M 3 ) that does not scale favorably w.r.t the number of basis functions. Finally basis functions that are numerically pruned from the model cannot be reintroduced at a later stage. As an alternative, the evidence-maximization criterion of Eq. ( 21) was also proposed by [START_REF] Tipping | Fast marginal likelihood maximisation for sparse bayesian models[END_REF] on the grounds of type-II maximum likelihood inference. Our active set method generalizes their fast marginal likelihood maximization procedure in presence of a smoothness-inducing prior.

Gaussian approximation of the likelihood

Although voxel intensities in the warped and fixed images are related by assumption via Gaussian noise (or a mixture thereof), the transformation Ψ -1 w acts non-linearly on intensity profiles and the resulting likelihood w.r.t. w does not belong to a standard family. To retrieve the required conjugacies during updates of q β , q z and q w,A , a Gaussian approximation of the data likelihood is used. It is derived from an efficient second-order Taylor expansion of the log-likelihood (AppendixA). The Taylor expansion is local: it depends on the point w around which it is computed. For updates of q β and q z , the approximation is used around the known mode of q(w|A). For updates of q w,A however, the approximation is taken at the mode of the true posterior p(w|D, A, β q β , π q π , λ q λ ), which must first be computed. This is done by quasi Newton optimization (L-BFGS) w.r.t. the subset of active variables (AppendixB).

Algorithm overview

The scheme proceeds according to Algorithm 1. We start with no active bases, S = ∅. We cycle between updates of the noise mixture parameters, of the transformation parametrization and parameters, and of the regularization parameter. Prior to updating q w,A (w, A), we update the approximation of the data likelihood. The global objective L(q) provides an always increasing lower-bound to the evidence p(D) and can be used to monitor convergence. Alternatively, the scheme can simply stop after a certain number of updates to the set of active bases has been performed.

Algorithm 1 Sparse Bayesian registration algorithm 1: Initialize Âk = ∞ for all k (S = ∅) and q λ 2: repeat 3:

for T iterations do Update q z to arg max q z L(q) following AppendixD.

5:

Update q β to arg max q β L(q) following AppendixD.

6:

Update q π to arg max q π L(q) following AppendixD. Update the likelihood approximation (AppendixA). 9: Update A (active set method) to greedily increase L(q) then set q(w|A) to arg max q(w|A) L(q) (AppendixB, AppendixE, AppendixF, AppendixG) 10:

Update q λ to arg max q λ L(q) according to AppendixC. 11: until no significant increase in L(q) or maximum number of iterations reached.

Algorithmic complexity

Updates of the mixture parameters take O(L • N) per pass on the image. The cost is dominated by the computation, for each of the N voxels, of soft-assignments to the L mixture components. Several O(L • N) passes are typically performed. The regularization level λ is updated in O(|S| 2 ), where |S| is the number of active bases. Updates to the parametrization A occur one basis at a time (a single A k is changed): each update takes O(|S| 2 + M|S| + N log N) to maintain necessary statistics, exploiting rank-1 linear algebra identities. As a byproduct, an update of q(w|A) is obtained. An overhead of O(|S| 3 + M|S| 2 + |S| • N log N) adds up to this, since statistics must be recomputed once in full after the mixture parameters and the likelihood approximation are updated. The likelihood approximation itself involves minimizing a registration energy w.r.t. the subset of active basis parameters (|S| M) to find the posterior mode, then a O(N) cost to compute the Gaussian approximation around this mode.

As a point of comparison, a single gradient descent step when optimizing the classical registration energy of Eq. ( 1) w.r.t. the full set of variables costs O(M 2 + N log N), where the left-hand term stems from the gradient of the regularization energy and the right-hand term from the gradient of the dataenergy. Exact Hessian computation in absence of sparsity costs O(M 2 + MN log N) and Hessian inversion is O(M 3 ).

Experiments & Results

Material

We experiment with the proposed framework on tasks of cardiac motion tracking. The goal is to recover the motion of the cardiac muscle over the course of the cardiac cycle from a time series of 3D images. The first experiment gives insight into the empirical behaviour of the proposed algorithm on a simple example of 2D pairwise registration. Other experiments involve full 3D + t motion tracking on various imaging modalities.

The first dataset consists of synthetic sequences of 3D ultrasound data provided as part of the registration challenge organized for the 2012 MICCAI workshop on Statistical Atlases and Computational Models of the Heart (STACOM). Details on the challenge methodology can be found in De Craene et al. (2013). These synthetic images count approximately 10 million voxels each, at a very fine isotropic resolution of 0.33mm. To avoid further optimization of our code in terms of RAM management, we downsampled them by a factor of 2. We thus worked at a resolution of 0.66mm at the finest level. The second and third datasets are hosted by the Cardiac Atlas Project. They were made available following the cardiac motion analysis challenge [START_REF] Tobon-Gomez | Benchmarking framework for myocardial tracking and deformation algorithms: An open access database[END_REF] organized for the 2011 MICCAI workshop on Statistical Atlases and Computational Models of the Heart (STACOM). The data includes a set of 15 sequences of real 3D tagged MR images at roughly 1mm × 1mm × 1mm resolution (1 million voxels), and a set of 15 sequences of real cine MR images at about 1.25mm × 1.25mm × 8mm resolution. The tagged sequences contain 20 to 30 frames each, the cine MR sequences 30 frames each. Fig. 8 displays example 2D slices from frames of each modality.

Details of the experimental setting

The experimental setup is identical across all modalities. The multiscale parametrization of the displacement field consists of isotropic Gaussian kernels of respective variance S 1 = 20 2 mm 2 and S 2 = 10 2 mm 2 , plus an anisotropic Gaussian kernel of variance 10 2 mm 2 in the short axis plane and 20 2 mm 2 along the long axis. The proposed framework imposes no restriction on the parametrization of the displacement field and we expect the anisotropy to be of potential relevance given the ventricle anatomy. All hyperparameters are set to uninformative values (a 0 = b 0 = c 0 = d 0 = η 0 → 0). The registration is accelerated with a multiresolution pyramidal scheme, starting with downsampled (smoothed) versions of images I and J and progressively moving through the pyramid of images to the images I and J at full resolution. Three resolution levels are introduced, downsampling by a factor of 2 at each level. Note that we do not make use of pre-segmentations of regions of interest. Computations were run on an Intel Xeon processor X5660 (@2.80GHz, 6 cores, 12 threads) and took 15 -30min per image pair for cine MRI, 30 -45min for tagged MRI and ∼ 90min for 3D US.

Self-tuning registration algorithm: an analysis

We use the example 2D registration of Fig. 3 and Fig. 7 to give some insight into the registration algorithm.

Basis selection & regularization. Fig. 9 demonstrates how basis selection mechanisms empirically combine with parameter re-estimation throughout iterations. A heuristic provides a large initial value for the regularization level λ. This initially discourages the addition of finer dictionary bases, whose impact on the signal regularity is too high at this stage. Coarse bases are added instead to capture the global trends in the observed displacement. The regularization level is consequently refined to reflect the regularity of the inferred displacement. As λ decreases towards a more sensible value, finer bases are incorporated in the active set to capture finer local details of the visible motion, or to ensure that these finer details of the inferred motion blend smoothly with the rest of the displacement field. In case of significant overlap between a subset of fine bases and a coarse basis, the basis at the coarsest scale may be deemed no longer to contribute towards a better explanation of the data and removed. Towards the last iterations, most actions consist in updating the orientation of active bases rather than in additions or deletions from the active set, and λ reaches a plateau as well. Fig. 7 further illustrates this mechanism of basis selection. The location of active bases is shown at two points in time: in the initial steps of the algorithm and at the end. The noise model learned on downsampled, smoothed images has a higher probability of low noise (higher peak around 0) but also of high noise (higher tails) due to the increased misalignment at the beginning of registration.

Noise model estimation. The noise model is jointly estimated

over the course of the algorithm. In all experiments it displayed a fast convergence. The Gaussian mixture also adapts quickly to changes in the distribution of intensity residuals that arise from the multiresolution pyramidal scheme, when hopping from a smoothed downsampled image to the next level in the pyramid of images, as seen from Fig 10.

Robustness w.r.t. initialization. Fig. 11 provides evidence towards the empirical robustness of the estimated level of regularity λ w.r.t. its initial value. The initial value of λ spans 4 orders of magnitude, whereas its final estimate varies by at most a factor of 4 across runs. Empirically, we observe that the regularity level λ decreases monotonically from its starting value towards a reasonable local optimum. As a limitation, it follows that the scheme will typically not recover the expected regularity level if initialized from too low a value of λ.

Synthetic 3D Ultrasound Cardiac Dataset

The appeal of this benchmark is to offer a dense ground truth in terms of motion and strain inside the cardiac muscle. The workflow of image synthesis uses the output of a cardiac electromechanical model to prescribe displacements in the myocardium. For each sequence of images, the ground truth consists of a sequence of meshes of the left and right ventricles end-diastole frame end-systole frame deformed over the cardiac cycle. The data extracted from such ground truth meshes can be compared to that obtained by deforming the mesh at a reference time point (namely, end diastole) throughout the cardiac cycle with the transformation output by the proposed registration approach. The visual and qualitative behaviour of the proposed approach was found to be satisfactory, even in terms of extrapolation: the inferred motion remained consistent in areas of the right ventricle that fall outside of the field of view (Fig. 12). This hints at an effective regularization mechanism, despite being automatically tuned.

We evaluate the accuracy of the proposed approach on a first subset of sequences that image the same motion at various Signal-to-Noise Ratios (SNRs). Because the proposed approach infers a consistent motion both inside and outside of the field of view, we find natural to assess its accuracy from statistics based on the whole mesh. This slightly departs from the methodology of De Craene et al. ( 2013) where part of the left ventricle only is considered. Fig. 13 reports the median pointto-point error in the inferred displacement for each time frame, where the median statistics is computed from every node in the mesh. At the best SNR, the highest error is observed around end systole with a median of 0.83mm, although the spread of error values becomes wider in the last frames. This falls in the same range as that reported for challenge participants by De Craene et al. (2013) -although slightly higher than the most accurate methodology. Of course part of the error is likely to be attributable to the use of downsampled, smoothed images with a resolution of 0.66mm as opposed to 0.33mm. Besides as the signal to noise ratio degrades, we observe as expected a global trend of increased error magnitude. As seen from Fig. 14, the increased SNR impacts the noise model (with a higher prevalence of small intensity residuals at high SNR) learned by the proposed approach, which in turn becomes more conservative in its estimates of displacements.

Fig. 15a reports (Green Lagrangian) strain measures at end systole averaged over AHA segments. This provides indirect evidence of the relevance of the automatically tuned regularity level λ and of the displacement parametrization. Ground truth values of strain obtained from the corresponding ground truth mesh are compared to those estimated from the output of registration. Variations in the strain across segments are generally well captured, even more so for its longitudinal and circumferential components. Similarly to most methodologies however, the radial strain -which captures the thickening of the muscle during the contraction -appears to be globally somewhat underestimated in the left ventricle. This might indicate a slight under-estimation of the endo-and/or epicardium displacement, due to a coarse parametrization or over-regularized transformation. The following table provides statistics on the number of bases of each scale used for the parametrization of the displacement field, for the normal case at highest SNR. The number of active bases on these sequences is typically smaller than that used in our experiments on cine and tagged data, with a lesser reliance on fine-scale bases. It may evidence increased conservatism in the estimated displacements, as well as indicate greater regularity of the synthetic ground truth motion. The benchmark also provides datasets that aim at reproducing pathological cardiac function, including a case where certain AHA segments become quasi akinetic due to ischemy. Fig. 15b summarizes estimated regional strains for this case, with qualitative retrieval of the ischemic segments (bolded contours), as emphasized by the comparison with the normal case. The accuracy on the ischemic case is similar to that of the normal case at identical SNR, with a median error at end systole of 0.80mm.

STACOM 2011 tagged MRI benchmark

On this image modality the grid-like tags allow to follow the motion of keypoints on the boundary of, or inside the cardiac Fig. 16 summarizes challengers' results along with ours at End-Systole (ES). The proposed approach achieves stateof-art results on this benchmark with a median accuracy of 1.46mm. As a point of comparison, the variability in the landmark tracking was estimated as part of the challenge methodology at 0.84mm. We perform two simple statistical tests to quantify the statistical significance of the increase in accuracy of our methodology compared to the challengers: a pairwise Student-t test and a pairwise Kolmogorov-Smirnov test. The tests are run for each pair of samples involving the proposed approach against a challenge participant's. The Student-t test aims at detecting significant differences in the true mean error of our method versus a challenger's, whereas the Kolmogorov-Smirnov test more generally aims at detecting whether the underlying distribution of errors differ. Figures are reported in Table 2 and provide some evidence towards a significant improvement from at least 3 of the 4 methodologies. The proposed formulation appears to achieve in a quasi automatic manner results qualitatively and quantitatively on par with the state of the art. In particular we insist that most parameters involved in the proposed formulation -the noise model and regularity level λ, the active parametrization of the displacement field -were automatically determined during registration. The strain maps and mesh deformations produced by the proposed scheme, as illustrated for instance in Fig. 17, also appear to be qualitatively on par with the best challenge results in that respect, and superior to that of the closest competing methodology accuracy wise (please refer to Tobon-Gomez et al. ( 2013) for a direct counterpart to Fig. 17). This again hints at the practical viability of the automatically adjusted trade-off between data and regularization energies. We report in Table 3 the number of bases of each scale used for the parametrization of the displacement field.

Cine MRI dataset: qualitative results and uncertainty

Original images had a low inter-slice resolution of 8mm compared to the in-plane resolution of 1.25mm. We upsampled them (typically by a factor of 5) prior to the registration process to prevent a degradation of numerical accuracy. To obtain a ground truth by direct manual tracking of landmarks over time was deemed difficult for this image modality. Instead the accuracy of the proposed algorithm was evaluated by crosscomparison with direct 3D+t segmentation results. Specifically, the endocardium was delineated over time on 2D slices using the freely available software Segment2 [START_REF] Heiberg | Time resolved three-dimensional automated segmentation of the left ventricle[END_REF], yielding a 3D point set of discretized contours. A 3D surface was then reconstructed as the zero level set of a signed distance map computed by radial basis interpolation, after estimating the normal to the surface at every point in the set from a local neighborhood3 . We then assessed the discrepancy between the reference end diastole segmentation transported over time via the output of registration, and the surface estimated by direct segmentation of the endocardium at each time step.

Fig. 18 summarizes the distribution of errors over time, pooled over all 15 sequences and all contour points, displaying the evolution of key quantile-based statistics. The median error reaches a satisfactory maximum of 1.82mm for frame 10, which roughly coincides with the end systole time for all volunteers. As a point of comparison, the volumes under consideration have a spacing of 1.25mm in the short-axis plane (i.e. within slices) and 8mm along the long-axis (i.e. inter-slice). The wide spread of error values partly reflects the challenge in obtaining a 3D segmentation of the endocardium that remains consistent over time (e.g. due to the variable appearance of papillary muscles). Misalignment of short-axis slices in 3D volumes, which may arise from the (slice by slice) image acquisition process, also accounts for some of the largest discrepancies. We observed no evident spatial pattern in the distribution of errors, although the segmentation rarely reached the very tip of the apical region.

The mixture-of-Gaussians noise model captured variations in the level of noise of an order of magnitude between distinct components (a factor of 10 between the standard deviations of extreme components). The visual aspect of the cardiac muscle changes drastically over time, so that these regions tend to be assigned higher noise levels than the baseline acquisition noise Figure 19: Example registration for the cine SSFP dataset of volunteer V5. We propagate the segmentation from the reference frame to the rest of the time-series with the output of the registration. The resulting mesh is overlayed on a 2D slice and visualized at three representative timesteps. The 3D mesh attests to the regularity of the underlying transform, and to its coherence over the cardiac cycle. level. Voxels in basal slices, with visible outflow tracts and apparent topology changes, also tend to fall in the noisiest components. Fig. 20 attests to the high variability (several orders of magnitude) of the optimal model parameter λ for varying sequences and time steps, which would render its manual estimation via a trial-and-error or cross-validation approach cumbersome. The apparent bimodality of the histogram might reflect the fact that cardiac phases with significant contraction or relaxation, around end systole, alternate with phases of lesser motion around end diastole.

Discussion

Adequacy of modelling assumptions

Despite using generic RBFs, the inferred parametrizations of 3D displacements were highly sparse, typically involving no more than a hundred degrees of freedom. This shows that the proposed sparsity-inducing mechanism is potent, and benefits both algorithmic complexity and memory usage. Finer bases were used more often in experiments with tagged MRI; in addition to the higher resolution of these volumes compared to cine MRI data, tags may have been regarded as reliable, informative structures along all directions of motion. While good accuracy was achieved on synthetic echocardiographic time series with a reduced number of bases, the synthetic motion from which the sequences were reconstructed is likely to have enjoyed greater regularity as well.

We did not explicitely make use of temporal regularization (as in e.g. De Craene et al. ( 2012)), but the temporal (and spatial) consistency of deformations remained satisfactory (Fig. 1(a) and Fig. 19). Incorporating temporal regularization and moving towards a large deformation framework [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Arsigny | A log-euclidean framework for statistics on diffeomorphisms[END_REF] with geodesic-by-part trajectories may still be advantageous, although technicalities should be addressed to maintain a reasonable computational complexity.

The proposed smoothness-inducing prior is widely used, and realizes a pragmatic trade-off between quality and complexity: it is stationary and involves a single hyperparameter λ. Spatially varying levels of regularization have been used in the recent literature [START_REF] Simpson | Probabilistic non-linear registration with spatially adaptive regularisation[END_REF][START_REF] Gerig | Spatially varying registration using gaussian processes[END_REF] to account for imaging artefacts or heterogeneous image content (e.g. tissue, blood), at the cost of additional technicalities and approximations in the variational inference. Here, spatially varying regularization occurs indirectly via a spatially varying noise level, and directly thanks to the adaptive parametrization of displacements. The latter results in a non-stationary prior on transformations (with coarser bases naturally enforcing higher regularization).

The sparsity-inducing prior over individual basis functions is parametrized by a positive symmetric precision matrix A k . This allows to activate a given basis along a single direction while constraining w k to be null in the orthogonal plane. This type of (at most) rank one A -1 k = α -1 k n k n k is in fact optimal under variational Bayes inference. To simplify derivations, an alternative all or nothing parametrization A k = α k I could have been chosen, where a basis function is artificially constrained to be either fully in use or fully pruned along all directions, at the cost of artificially tripling the number of active degrees of freedom.

For mono-modal registration the assumption that intensities between source and target images coincide up to spatially varying noise mostly holds. In this context, the mixture-of-Gaussian model of residuals is flexible yet simple enough to be efficiently and robustly fit jointly during the registration of images of interest, as opposed to beforehand on training data [START_REF] Leventon | Multi-modal volume registration using joint intensity distributions[END_REF][START_REF] Zhou | Boostmotion: Boosting a discriminative similarity function for motion estimation[END_REF][START_REF] Lee | Learning similarity measure for multi-modal 3d image registration[END_REF][START_REF] Tang | Locally-adaptive similarity metric for deformable medical image registration[END_REF]. For higher interpretability of the mixture, the variational Bayes approach allows to select the optimal number of components [START_REF] Penny | Robust bayesian general linear models[END_REF][START_REF] Archambeau | Robust bayesian clustering[END_REF], although this was not pursued here. For multi-modal registration, the mapping between source and target image intensities can also be regressed [START_REF] Guimond | Three-dimensional multimodal brain warping using the demons algorithm and adaptive intensity corrections[END_REF][START_REF] Janoos | Bayesian characterization of uncertainty in multi-modal image registration[END_REF].

The assumed independence of voxelwise intensity residuals is not realistic. To avoid placing too much confidence on data, the virtual decimation procedure downweights the data term (by up to two factors of magnitude in 3D experiments) in a mostly ad-hoc but empirically viable manner. Designing spatially varying, correlation-aware models of image discrepancies would address the matter more elegantly, but is outside of the scope of this paper.

Modelling uncertainty in the interpolation of discrete intensity profiles (cf. AppendixA) proved to be appropriate. This is examplified by cine MR data, due to the lower long-axis resolution. If not accounting for it the regularity of the inferred transform in the long axis direction was systematically found, upon visual inspection, to be of lesser quality. This behaviour is expected if the scheme is unaware of its increased reliance on interpolation to match intensity values between images. Image upsampling prior to registration also involves interpolation and was accounted for in an identical manner.

Implementation & computational load

Running times for the proposed approach are on the same order of magnitude as the state-of-the-art in cardiac motion tracking [START_REF] De Craene | 3d strain assessment in ultrasound (STRAUS): A synthetic comparison of five tracking methodologies[END_REF][START_REF] Tobon-Gomez | Benchmarking framework for myocardial tracking and deformation algorithms: An open access database[END_REF]. Although our implementation is CPU based (with partial multithreading), most computations involve statistics at the level of individual voxels or basis functions and are highly parallelizable. The active set method used to update the active parametrization is technical but computationally inexpensive. Significant improvements in running time may be obtained instead by basic optimization of the energy minimization problem solved to compute the posterior mode for the Laplace approx-imation4 , by partial parallelization of the VB updates for the noise mixture (AppendixD) and by optimizing convolution filters that we rely on (AppendixH). Limitations of the proposed approach include the technicality of its implementation and its memory consumption.

Upon inspection, the expressions derived for λ q λ , β q β and other hyperparameter expectations relate to intuitive quantities, such as averages of voxelwise square intensity residuals, or the energy in the estimated displacement field. Simplified versions of these updates can be used independently of the specifics of the registration scheme as ad-hoc recipes for automated parameter tuning.

Atlas of motion

In the present work, the reduced parametrization simply benefits the algorithmic complexity of the schemes and the quality of interpolation. If an anatomically relevant parametrization were desired however, we note that the active set method naturally extends to a multi-subject setting: it may be used to yield an optimal, joint parametrization of displacements. This could constitute a basis to learn a parametric atlas of motion from a small dataset of 3D+t images in the spirit of e.g. [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF]; [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF]; [START_REF] Gori | Bayesian atlas estimation for the variability analysis of shape complexes[END_REF].

Uncertainty quantification

The proposed method returns a Gaussian variational approximation q(w|A) N(µ, Σ) of the posterior distribution p(w|I, J) of transformation parameters. The |S| × |S| covariance matrix Σ is of small size and readily computed over the course of the algorithm. The covariance on transformation parameters can be turned into directional estimates of uncertainty at any point in space by simple linear algebra, or can be sampled from at a marginal cost to efficiently explore the joint variability of the full transformation. Sampling the transformation itself, unlike sampling displacements independently at each point in space, preserves correlations in the displacement of close-by points. This allows to derive empirical uncertainty estimates on integral geometrical quantities.

Fig. 1(b) reports estimates of uncertainty in the volume enclosed over time by the endocardium surface, as segmented on the reference frame (at time 0), for a cine MRI sequence (volunteer 5). For the same volunteer, Fig. 1(c) summarizes, in the form of a tensor map, the uncertainty in the inferred displacement field at end-systole, accounting for uncertainty in the output of each frame-to-frame registration between end-diastole and end-systole. Tensors are rasterized at the voxel centers of the end-systole frame. Each tensor encodes (the square root of) the 3 × 3 covariance matrix of the pointwise displacement and is elongated in directions of higher uncertainty. Due to voxel anisotropy, the direction of higher uncertainty is, consistently across space, aligned with the long-axis. The color scheme thus encodes the second principal direction of highest uncertainty. Steep intensity gradients in the underlying image typically translate into directions where tensors are least elongated.

Tensor magnitude and principal directions vary smoothly across space, as estimates of uncertainty incorporate information of a local (and in fact, global) nature. The yellow dashed line gives a visual cue as to the position of the left ventricle endocardium boundary. The agreement between exact and approximate posteriors should be explored in future work.

Conclusions

In this paper we proposed an approach to data-driven, spatially adaptative, multiscale parametrizations of deformations for registration. It uses larger kernels in regions of high uncertainty due to e.g. lack of image gradients or incoherent information in registered images, and uses smaller kernels where a finer motion can be estimated with confidence from local cues in paired images. This is achieved in a Bayesian framework so that the approach retains natural advantages of probabilistic formulations. It is self-tuning, with hyperparameters being jointly inferred during registration. Inference is tractable on real-scale data thanks to an efficient Variational Bayes method.

The core methodological contribution is a procedure for fast marginal likelihood maximisation in sparse Bayesian models, that relaxes the assumptions made by [START_REF] Tipping | Fast marginal likelihood maximisation for sparse bayesian models[END_REF] for the fast Relevance Vector Machine. The prior is allowed to encode correlations between explanatory variables so as to favor smooth solutions. The proposed structured sparse Bayesian model itself, and variants thereof, are relevant to a variety of generalized regression problems, including image-based classification and regression tasks.

Acknowledgments. The first author was partly funded by the Microsoft Research -Inria Joint Centre, and part of this work was funded by the European Research Council through the ERC Advanced Grant MedYMA (2011-291080) on Biophysical Modeling and Analysis of Dynamic Medical Images.

AppendixA. Second order Taylor expansion of the model likelihood

Taking a Gaussian approximation of the likelihood facilitates the variational Bayes inference, as the variational posterior q(w|A) ≈ N(µ, Σ) is then approximately Gaussian (Ap-pendixB). We seek such an approximation of the likelihoodor equivalently we look for a quadratic approximation of the log-likelihood. In this appendix we derive one from the secondorder Taylor expansion of D β,z (J; I, w) =log p(D|w, β, z). Around the point w * , it takes the form:

D β,z (J; I, w) ≈ D β,z (J; I, w * ) + ∇[D β,z ] w * (w -w * ) + 1 2 (w -w * ) H[D β,z ] w * (w -w * ) (A.1)
where ∇[D β,z ] w * stands for the gradient at w * and H[D β,z ] w * for the Hessian at w * . For the gradient we obtain:

∇[D β,z ] w * = N i=1 βi (I[V i ] -J[v i ]) φ(v i ) T ∇I(V i ) , (A.2)
where βi = 1≤l≤L z il β l . {v i } N i=1 is the list of voxel centers in the fixed image and V i Ψ -1 w * (v i ) = v i + φ(v i )w * are the paired coordinates in the moving image. For the Hessian:

H[D β,z ] w * = N i=1 βi φ(v i ) ∇I(V i )∇I(V i ) φ(v i ) + N i=1 βi (I[V i ] -J[v i ]) φ(v i ) H[I](V i )φ(v i ) . (A.3)
After dropping the term involving the Hessian of the image5 , we arrive at the following approximation of D β,z (J; I, w): .4) t i ∈ R d and H i ∈ M d×d only depend on the point w * around which the approximation is taken, and are respectively given by Eq. (A.5) and Eq. (A.6), noting u * (v i ) φ(v i )w * .

D β,z (J; I, w) ≈ 1 2 N i=1 (t i -φ(v i )w) βi H i (t i -φ(v i )w) . (A
t i = u * (v i ) - I(V i ) -J(v i ) ∇I(V i ) 2 ∇I(V i ) , (A.5) H i = ∇I(V i )∇I(V i ) . (A.6)
Eq. (A.4) demonstrates that up to a second order local approximation, registration can be recast as a regression task with a set of (virtual) observations t i and associated (heteroscedastic, anisotropic) confidence H i . This is known as a generalized linear model. In block form, the data likelihood can be expressed as p(D|w, β, z) ≈ p(t|w, β) with:

p(t|w, β) ∝ exp - 1 2 (t -Φw) βH(t -Φw) . (A.7)
t denotes the concatenation of all t i . βH is a block diagonal matrix with the ith diagonal block equal to ( βH) i βi H i .

In our particular instance, the virtual pairings relate to the optical flow: if we dropped the confidence tensors H i , Eq. (A.4) would yield an approximation of Eq. (B.6) much in the spirit of the demons algorithm [START_REF] Thirion | Image matching as a diffusion process: an analogy with maxwell's demons[END_REF][START_REF] Cachier | Isotropic energies, filters and splines for vector field regularization[END_REF]. The tensors H i vary sharply across the image however, e.g. as edges or boundaries are crossed. They assign anisotropic, spatially varying confidence in voxelwise pairings and account for how informative and structured the image is at the point of interest. The local approximation of Eq. (A.4) transforms an image-based criterion into a landmark-based one, and the proximity to formulations in the related literature [START_REF] Rohr | Spline-based elastic image registration: integration of landmark errors and orientation attributes[END_REF] is indeed striking in this form.

The confidence βi H i in the virtual voxelwise pairing t i + v i can grow arbitrarily high for arbitrarily high intensity gradients. These expressions result from linearizing the intensity profile I around the current pairing V i , and are blind to interpolation uncertainty in evaluating I(V i ) and ∇I(V i ). To address this shortcoming, we propose to replace βi

H i = βi ∇I(V i )∇I(V i ) by ( βi H i ) -1 + D int -1 = 1 1 + tr[ βi H i D int ]
• βi H i (A.8) which implements a soft upper threshold on the precision, as a heuristic for interpolation uncertainty. D int acts as a minimum covariance: it is a diagonal matrix set to the square of -say-half the voxel spacing to prevent unreasonable subvoxel confidence.

AppendixB. VB inference -optimization of q w,A

We recall that the variational posterior q w,A is constrained to lie in the family q w,A (w, A) = q(w| Â)δ Â(A). AppendixB.1 restates exact optimality conditions and exhibits where the likelihood is involved. Under a Gaussian approximation of the likelihood, (approximate) optimality conditions yield tractable updates of q(w| Â) and suggest an efficient active set scheme to update  (AppendixB.2). Technicalities for the Gaussian approximation are clarified in AppendixB.3.

AppendixB.1. Exact optimality conditions

From optimality conditions given in section 3.2 by Eq. ( 19) and keeping explicit only the terms that depend on w, we obtain:

q * (w| Â) = p D|w, β q β , z q z • p w| Â, λ q λ p D| Â, λ q λ , β q β , z q z , (B.1)
where we exploited the linearity of the logarithm of distributions of interest w.r.t. β, z and λ to take the expectations of those variables w.r.t. the associated variational factors inside those distributions. The numerator of Eq. (B.1) is simply the product of the likelihood times the prior on transformation parameters, with λ, β and z fixed at their expected values with respect to their respective variational posteriors. Again from Eq. ( 22), the optimal A * maximizes the denominator of Eq. (B.1):

p D | Â, λ q λ , β q β , z q z = w p D|w, β q β , z q z • p w| Â, λ q λ dw (B.2)
and q * w,A (w, A) = q * (w|A * )δ A * (A). Maximizing Eq. (B.2) w.r.t. Â and evaluating Eq. (B.1) however is hard since the likelihood p(D|w, β q β , z q z ) is non trivial as a function of w. For convenience we drop the expectations from notations in what follows (and the hat over A), as it does not affect derivations.

AppendixB.2. Approximate optimality conditions

The tractability of the inference relies on approximating the data likelihood following AppendixA as a Gaussian distribution p(D|w, β, z) ≈ N(t|Φw, βH), with virtual data t interpretable as the concatenation of noisy voxelwise displacements. The prior p(w|A, λ) = N(0, A + λR) is also Gaussian, yielding in turn a Gaussian variational posterior q * (w|A) ≈ N(µ, Σ) with:

µ = ΣΦ βHt Σ = (Φ βHΦ + λR + A) -1 . (B.3)
Moreover the marginal likelihood p(D|A, λ, β, z) ≈ p(t| β, A, λ) is Gaussian as well:

p(t| β, A, λ) = | 2π C | -1/2 • exp - 1 2 t C -1 t (B.4)
where by identification C -1 = βH -( βH)ΦΣΦ ( βH). In other words, the marginal distribution of the virtual data t conditioned on the hyperparameters {A, λ, β} is Gaussian N(0, C). Furthermore, it follows from the Woodbury matrix identity that

C = ( βH) -1 + Φ(A + λR) -1 Φ . (B.5)
The objective is to increase Eq. (B.4) w.r.t. A. The two factors in Eq. (B.4) have antagonistic effects: while the left hand term penalizes covariance matrices C that waste mass (via |C|), the right hand term gives incentive to spend mass to better explain the data t. This compromise mechanically leads to sparsity. Indeed looking at the form of C in Eq. (B.5), we see that part of the data is explained for free by the contribution ( βH) -1 of the noise to C regardless of the right hand term; thus only a few degrees of freedom need be active (A k < ∞) to fully explain the data.

Algorithm 2 Optimization of A 1: if the likelihood approximation was updated then 2:

Recompute µ, Σ in full from Eq. (B.3).

3:

Recompute statistics q k , s k and κ k in full for all k. 4: end if 5: for p iterations do 6: ∀k ∈ S (resp. k S), compute the gain max A k ∆l(A k ) in log-evidence obtained by updating or deleting k from S (resp. adding k to S) from AppendixE.

7:

Select the most favorable action l such that:

max A l ∆l(A l ) ≥ max A k ∆l(A k ) 8:
Set A * l = arg max A l l(A l ) and update S.

9:

Update µ, Σ via rank-one identities (AppendixF).

10:

Update q k , s k and κ k for all k using rank-one updates (AppendixG). 11: end for

In AppendixE we single out the contribution ∆l(A k ) of each basis φ k to the log marginal likelihood given the state of other bases. AppendixF and AppendixG show that the statistics q k , s k and κ k required to express this contribution can be updated efficiently using rank-one linear algebra identities. This is the basis for the proposed active set method, which optimizes A by updating one A k at a time according to Algorithm 2. AppendixB.3. Around which point is the likelihood approximated? We have not yet specified around which point the Gaussian approximation is taken. The mode µ MP of the (true) posterior p(w|D, A, β q β , π q π , λ q λ ) is chosen. This departs from [START_REF] Simpson | Probabilistic inference of regularisation in non-rigid registration[END_REF], who use a Taylor expansion around the mode of the (Gaussian) variational posterior q w , and from Le [START_REF] Folgoc | Sparse bayesian registration[END_REF] who use a different strategy to derive a quadratic approximation. Here, as in energy-based registration we numerically solve for the minimizer µ MP of Eq. (B.6), using the current estimate of hyperparameters A, λ , β , π :

E(w) = D β , π (J; I, Φw) + 1 2 w (A + λ R)w (B.6)
with the notable difference, in terms of computational complexity, that Eq. (B.6) only involves the sparse subset S of active bases. Other weights are effectively constrained to zero due to an infinite penalty A k . Note that we marginalize over soft assignments z of voxels to components of the noise mixture. This comes at little computational cost and can reasonably be expected to accelerate convergence. The mode is found by quasi Newton (BFGS) optimization, using closed form expressions for the energy and its gradient. Expressions for the energy follow Eq. (B.6) and Eq. ( 8). The gradient of the data term is given by:

∇ w [D β , π ] = N i=1 βi (I[V i ] -J[v i ]) φ(v i ) T ∇I(V i ) . (B.7)
where βi = 1≤l≤L ρ il β l can be seen as an effective noise level for the ith voxel, and

ρ il ∝ π l Z l exp -β l 2 (J[v i ] -I[V i ]
) 2 (such that l ρ il = 1) can be seen as a soft assignment of the ith voxel to the lth component.

AppendixC. VB inference -optimization of q λ

Taking terms that do not depend on λ in the constant, the optimality condition of Eq. ( 17) rewrites as:

log q λ (λ) = log p(w|A, λ)p(λ) q w,A q π q z q β + const (C.1) = log p(w|A, λ)p(λ) q w,A + const (C.2) = log p(w| Â, λ) q(w| Â) + log p(λ) + const. (C.3)
Recall that • q(θ) = θ • q(θ)dθ denotes expectation w.r.t. q(θ).

The second equation uses the fact that the integrand does not depend on either π, z nor β. The third uses the fact that q w,A (w, A) = q(w|A)δ Â(A). The prior on λ is Gamma distributed, p(λ) = Γ(λ|a 0 , b 0 ), so that:

log p(λ) = -b 0 λ + (a 0 -1) log λ + const(λ) . (C.4) p(w| Â, λ) = N(0| Â + λR) is a degenerate Gaussian.
For all inactive bases k S, Â-1 k = 0, and the Gaussian degenerates to a Dirac at 0 along the d corresponding dimensions. For all active bases k ∈ S, the Gaussian degenerates along d -1 directions since Â-1 k = α -1 k n k n k is rank one (section 2.4, AppendixE). By assumption α k = 0 so that, noting w S = (w k n k ) k∈S and w -S = w\w S :

p(w| Â, λ) = δ 0 (w -S ) • | λ 2π R S | 1/2 exp - λ 2 w S R -1 S w S , (C.5)
where R S is an |S|×|S| matrix whose k, lth coefficient is n k Rn l . Finally, since q(w| Â) = N(µ, Σ):

log q λ (λ) = |S| 2 log λ - 1 2 µ R S µ + tr(ΣR S ) λ -b 0 λ + (a 0 -1) log λ + const(λ) . (C.6)
In other words, q λ (λ) = Γ(λ|a, b) is Gamma distributed with hyperparameters given by Eq. (C.7) and (C.8), and

λ q λ = a/b. a = a 0 + |S|/2 (C.7) b = b 0 + µ R S µ + tr(ΣR S ) /2 . (C.8)
AppendixD. VB inference -optimization of q π , q β , q z VB updates are stated without proof. Derivations follow the same strategy as previously, starting from the optimality condition of Eq. ( 17). They follow the same outline as those of [START_REF] Archambeau | Robust bayesian clustering[END_REF]. The variational posterior for assignments z of voxels to a component of the Gaussian mixture is a product of categorical distributions q z = 1≤i≤N C(z i |ρ i1 • • • ρ iL ) with parameters ρ il interpretable as softassignments, summing to 1 (∀i, l ρ il = 1) and given by Eq. (D.1). The variational posterior q π for the noise mixture proportions

π = {π 1 • • • π L } follows a Dirichlet distribution Dir(π|η), with η = (η 1 • • • η L )
given by Eq. (D.2). The variational posterior for noise levels q β = 1≤l≤L Γ(β l |c l , d l ) is a product of Gamma distributions over each individual mixture component, with parameter updates given by Eq. (D.3) and (D.4).

ρ il ∝ πl β1/2 l exp -1 2 β l q β e 2 i (D.1) η l = N πl + η 0 (D.2) c l = 1 2 N πl + c 0 (D.3) d l = 1 2 1≤i≤N ρ il e 2 i + d 0 (D.4)
The quantities required for these updates are as follows:

log πl log π l q π = ψ(η l ) -ψ( 1≤l≤L η l ) (D.5) log βl log β l q β = ψ(c l ) -log d l (D.6) πl 1 N 1≤i≤N ρ il (D.7) e 2 i (I[Ψ -1 w (v i )] -J[v i ]) 2 q(w| Â) ≈ I[Ψ -1 µ (v i )] -J[v i ] 2 + tr φ(v i ) Σφ(v i ) H i (D.8) β l q β = c l /d l (D.9)
where ψ( • ) stands for the digamma function and tr[ • ] for the trace. Eq. (D.8) follows from the quadratic approximation of AppendixA, noting in addition that q(w| Â) = N(µ, Σ) (cf. Ap-pendixB and Eq. (B.3) for the updates of µ, Σ). Each factor q z , q π , q β is updated in turn for T iterations; each iteration essentially computes then aggregates voxelwise statistics in a single pass over the image of squared intensity residuals e 2 i .

AppendixE. Contribution of a basis to the log marginal likelihood

From Eq. (B.4) we see that the log marginal likelihood L = log p(t| β, A, λ) is given up to additive constant by

L = - 1 2 log |C| + t C -1 t (E.1)
with C = ( βH) -1 + ΦLΦ , where we define

L (A + λR) -1 . (E.2)
Noting that C exclusively depends on the basis k via the kth (d× d block-) diagonal coefficient of A and the kth column (of d × d elements) of Φ, we would like to single out the contribution l(A k ) of any such basis φ k to the log marginal likelihood in the form:

L = l(A k ) + L -k (E.3)
where L -k does not depend on the basis k. If we denote by L -k the inverse of the matrix obtained by removing the kth column6 from L -1 = A + λR (or equivalently by setting A k = +∞7 in L), we see from the Woodbury rank one matrix identity 6 that L

= L -k + U k L kk U k , with U k = (λL -k R k ) I and L kk = (A k + κ k ) -1
, where the d × d matrix κ k is defined as:

κ k λR kk -(λR k ) L -k (λR k ) (E.4)
By injecting this latter decomposition of L into the expression of C, we derive a decomposition of C into the sum of a term that does not depend on the kth basis and of a rank one term:

C = C -k + (ΦU k )(A k + κ k ) -1 (ΦU k ) (E.5) Letting C -1 -k (C -k ) -1
, a second application of rank one update identities for the determinant and the inverse gives the two following expressions E.6 and E.7 for the two terms in the righthand side of the log marginal likelihood expression E.1:

|C| = |C -k | • |A k + κ k | -1 • |A k + κ k + s k | (E.6) t C -1 t = t C -1 -k t -q k (A k + κ k + s k ) -1 q k (E.7)
We introduced the statistics s k ∈ M d,d and q k ∈ R d respectively defined as:

s k (ΦU k ) C -1 -k (ΦU k ) (E.8) q k (ΦU k ) C -1 -k t (E.9)
Eq. (G.1) and Eq. (G.2) in AppendixG provide alternative, more easily interpretable expressions for these quantities. This yields the following expression for l(A k ):

l(A k ) = log |A k + κ k | -log |A k + κ k + s k | + q k {A k + κ k + s k } -1 q k . (E.10)
It is of practical significance to the algorithmic complexity of our schemes that the quantities involved (L, µ, Σ) do not actually depend on bases that are not in the active set S (i.e. all bases s.t. A k = +∞). Similarly s k , κ k and q k only involve the set of active bases S augmented with the kth basis, due of the form of U k .

The maximization of Eq. (E.10) under constraint that A k is a symmetric positive semidefinite d × d matrix involves the gradient of the (unconstrained) function l(A k ):

∇l(A k ) = -σ k q k q k -s k -s k (A k + κ k ) -1 s k σ k (E.11)
where σ k is shorthand for (A k + κ k + s k ) -1 , and ∇l(A k ) is a d ×d matrix. Since σ k is symmetric positive definite and q k q k is of rank one, ∇l(A k ) has at most one negative eigenvalue. More precisely, if q k q k -s k is negative then ∇l(A k ) is positive definite for all A k and the improper maximizer of l(A k ) lies at infinity A k → +∞. Otherwise there is exactly one negative eigenvalue and we look for maximizers of the form A -1 k = α -1 k n k n k . This is consistent with the intuitive comment that A -1 k ∈ M d,d cannot be fully determined from a single "observation" q k ∈ R d and should be degenerate. Rewriting Eq. (E.10) as a function of α, n leads to maximizing E.12 under constraint that α is positive (dropping the index k for convenience). Note also that E.12 is invariant under reparametrization n → νn, α → α/ν 2 . l(α, n) =log 1 + n sn α + n κn + (q n) 2 α + n (κ + s)n (E.12)

At a maximizer α * , n * = arg max α,n l(α, n) the constraint is either active (α * = 0) or inactive (α * > 0). If inactive, the solution actually maximizes the unconstrained function E.12 and is given by α * = ᾱ( n), n * = n where ᾱ(n) = (n sn) 2 (q n) 2 -n sn -n κn , (E.13) n = s -1 q . (E.14)

In this case l(α * , n * ) is simply equal to l(n * ), where l(n) is defined by E.15 with ξ(n) (q n) 2 /n sn.

l(n) -log ξ(n) + ξ(n) -1 (E.15)
In addition l( n) can be shown to always provide an upper bound to the maximum contribution of a basis to the evidence, max α,n l(α, n). If ᾱ( n) < 0, the upper bound is not reachable and the constraint is active, α * = 0. In that case, we numerically optimize over the unit sphere in R d to find n * . The case occurs when the l 2 -norm regularization is by itself sufficient along the direction n * , and no additional shrinkage is deemed necessary. It is also the default case if we purposedly restrict α to be either null or infinite. To save on unnecessary computations, we first check that the upper bound l( n) to the maximum contribution of the basis k to the evidence is superior to the current best contribution among bases already handled.

AppendixF. Update of µ, Σ, L

Updates of the moments of the posterior distribution µ, Σ = (Φ ( βH)Φ + A + λR) -1 and of L = (A + λR) -1 upon deletion from the model, update or addition to the model of a basis l are done similarly to [START_REF] Tipping | Fast marginal likelihood maximisation for sparse bayesian models[END_REF] and follow from Woodbury identities. Denoting updated quantities with a tilde, we get in the case of deletion8 : In particular, μl = Σll q l . The case of the update of a basis l is treated as a deletion followed by an addition, updating s l and κ l in-between these actions as they are needed in Eq. (F.7).

Σ = Σ -Σ l Σ -1 ll Σ l , (
AppendixG. Update of s k , κ k and q k

From the resolvent identity, we note that ΣΦ ( βH)ΦL = L-Σ. Using such relationships after developing the factors in E.8 and E.9, we derive alternative expressions for s k and q k :

s k = φ k ( βH)φ k -Π k Σ -k Π k + (λR k ) L -k (λR k ) (G.1) q k = φ k ( βH)t -Π k Σ -k Φ ( βH)t (G.2)
where Π k is a column vector of d × d matrices defined by Π k = Φ ( βH)φ k + λR k . Π k can be interpreted as the inner product of basis k with all the active bases w.r.t an appropriate metric, in the sense that its jth d × d coefficient is given by: Π jk = φ j ( βH)φ k + λ Dφ j |Dφ k . q k is the projection on the basis φ k of the optimal residual t -Φµ -k for the active set S\{k}, with a correction to account for regularization 9 : it is an unnormalized indicator of the relevance of basis φ k to better explain the data. s k + κ k = φ k ( βH)φ k + λR kk -Π k Σ -k Π k measures how much confidence basis φ k aggregates, taking into account its overlap with bases already in the active set (the confidence already captured by overlapping active bases is withdrawn). In the specific case where λ = 0, we retrieve the quantities and expressions derived by [START_REF] Tipping | Fast marginal likelihood maximisation for sparse bayesian models[END_REF] for the RVM. We 9 Indeed, q k = φ k ( βH)(t -Φµ -k ) -λR k µ -k .

found useful to introduce surrogate quantities t k and r k respectively defined according to G.3 and G.4:

t k φ k ( βH)φ k -Π k ΣΠ k + (λR k ) L(λR k ) (G.3) r k φ k ( βH)t -Π k ΣΦ ( βH)t (G.4)
These quantities merely differ from s k and q k in that the index -k was dropped from Σ -k and L -k . Our underlying motivation is to update simpler quantities t k and r k that still retain a straightforward link to the statistics s k and q k of interest for the computation of l(A k ). Indeed, for a basis φ k that does not lie in the model, Σ -k = Σ and L -k = L. Therefore, the quantities under consideration coincide: s k = t k and q k = r k . For a basis k that lies in the model and noting that Σ -k = Σ -Σ k Σ -1 kk Σ k , we obtain the statistics of interest efficiently as:

s k = t k + Π k Σ k Σ -1 kk Π k Σ k -[(λR k ) L k ] L -1 kk [(λR k ) L k ]
(G.5) where rk and tk denote updated quantities, as opposed to quantities prior to the update r k and t k . The quantities indexed by l are computed (once for all bases) following AppendixF. Similarly we maintain λR kk -(λR k ) L(λR k ) instead of κ k (index -k dropped from L -k ).

q k = r k + Π k Σ k Σ -1 kk µ k (G.

AppendixH. Basis functions: Computational complexity with translation invariant kernels

The section discusses the generality of the proposed approach w.r.t. the basis functions φ k parametrizing the displacement field. Smooth radial basis functions as well as many types of splines can be used without affecting the computational complexity. We show how to efficiently compute elements R kl of the quadratic regularizer and inner products φ k ( βH)Φ, which are required in the scheme (see for instance AppendixF and AppendixG).

Regularization and Fourier analysis. Let K(x, y) = κ(xy) be a translation invariant positive definite kernel such that κ is integrable on R d as well as its Fourier transform κ. We consider the space H of integrable d-vector fields f : R d → R d such that f 2 = f | f < +∞, where

f |g = 1 (2π) d ξ f (ξ) † g(ξ) κ-1 (ξ)dξ . (H.1) f (ξ) F [ f ](ξ)
x exp{-ιx ξ} f (x)dx is the Fourier transform of f . If κ is sufficiently smooth, the successive partial derivatives ∂ x i 1 •••x ip f of elements f ∈ H exist and all lie in H [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF]. As such, we can consider families of regularizers of the form R D ( f ) = D f 2 , where D is a linear differential operator. κ enables additional filtering in the frequency domain and can be set to κ = 1 if no such penalty is desired.

Figure 1 :

 1 Figure 1: (a) Trajectories of points on the endocardium, following the registration of a time series of cardiac MR images by the proposed approach. (b) LV volume over time and 99.7% confidence interval. (c) Tensor visualization of directional uncertainty at end-systole, rasterized at voxel centers of a 2D slice.

Figure 2 :

 2 Figure 2: Graphical model of registration. The generative model of data D involves a transformation Ψ of space, and noise governed by a set of underlying parameters P. Hyperpriors (with hyperparameters H P ) are in turn imposed over the noise parameters. The transformation is parametrized as a linear combination of predefined basis functions {φ k , k = 1 • • • M} with associated weights w k . Priors on the transformation smoothness and on the relevance of individual bases introduce additional parameters λ and A k . Random variables are circled, hyperparameters are squared. Arrows capture conditional dependencies. Shaded nodes are observed variables or fixed hyperparameters. The transformation Ψ is fully determined by its parent nodes (the φ k and w k ), hence the doubly circled node. The content of plates is replicated (M times).

Figure 3 :

 3 Figure 3: We illustrate the appeal of a robust variant of the SSD image loss based on a mixture-of-Gaussians model (GMM). Images (c,d) display the output warped images obtained after registering images (a) and (b), using respectively the SSD-based likelihood or the GMM-based likelihood (section 2.1). The arrows point towards a specific region that highlights the limitations of the SSD: the subset of hypo-intense voxels bordering the myocardium in the fixed image has no evident counterpart in the moving image. The SSD still drives the motion towards the best matches intensity-wise, which induces implausible tangential stretch of the myocardium. The GMM, on the other hand, incorporates a natural mechanism to downweight regions that cannot be reliably paired from image to image based on intensity values. The inferred motion is qualitatively closer to our expectations.

Figure 4 :

 4 Figure 4: (Left) Example residual image following registration. Artefacts and structures that changed appearance from one image to the other stand out much unlike ambient noise. Note that the intensity of the cardiac muscle itself differed in the pair of images. (Middle) Histogram of intensity residuals, with SSD and GMM fits overlayed respectively in red and green. (Right) Energy profiles for the SSD (grey dashes) and GMM (black line). The voxelwise penalty is plotted as a function of the intensity residual. GMM achieves robustness thanks to concave inflexions that result in a soft threshold on the penalty incurred for large intensity residuals.

Figure 5 :

 5 Figure 5: Impact of the regularization model. Displacements are parameterized by isotropic Gaussian kernels of set width σ = 0.25. From left to right, the regularizer varies. The data consists of 4 points regularly sampled on the unit circle, forming an axis aligned square, pulled twice as far away from the origin as they initially were. The warped grid obtained by regression is displayed along with the ground truth displacement.

Figure 6 :

 6 Figure 6: Impact of the basis scale on the inferred transform. From left to right, the displacement field is parameterized by isotropic Gaussian kernels of increasing width. The data consists of 8 points regularly sampled on the unit circle. The underlying motion is a rotation of π/4 radian. Only four of these eight rotated samples are displayed for readibility. The scale of the bases used to represent the transform affects the area of influence of the data points and the scale at which the regressed transform resembles a global rotation.

Figure 7 :

 7 Figure 7: Basis selection mechanism displayed on an example 2D registration between slices of cardiac MR images (cf. sec 4.3), respectively at ES (a) and ED (b). (d,e) Bases selected in the initial steps of the algorithm vs. at the end. The locations and scales of the Gaussian RBFs are indicated by circles (isocontour at 1 std). (c) Inverse displacement field output by the algorithm (scale factor: 2), smoothly varying across the whole image.

Figure 8 :

 8 Figure 8: Example slices for the cardiac imaging modalities that we experiment on, with artefacts and patterns peculiar to each modality. (Left) 3D tagged MR image. (Middle) 3D echocardiographic image. (Right) 3D cine MR image.

Figure 9 :

 9 Figure 9: Basis selection mechanism and its coupling with the jointly estimated regularization level, across iterations. (Top) Addition, update or deletion of dictionary bases in the active parametrization of the displacement field across iterations. Three distinct scales are used in the representation of displacements (1 curve per scale). (Bottom) Regularization parameter λ, updated every few iterations, plotted against the number of iterations run since the beginning of the registration.

Figure 10 :

 10 Figure 10: Inferred noise model at the beginning of registration and at the end.The noise model learned on downsampled, smoothed images has a higher probability of low noise (higher peak around 0) but also of high noise (higher tails) due to the increased misalignment at the beginning of registration.

Figure 11 :

 11 Figure 11: Robustness of the inferred regularity level w.r.t. its initial estimation. The 2D registration is run 3 times, and initialized each time with a differing level of regularity λ (respectively 10 4 , 10 6 , 10 8 ). Each curve shows the evolution of λ over the course of the associated run.

Figure 12 :

 12 Figure 12: Ground truth mesh (green transparent surface) vs. reference mesh transported via registration (overlayed black wireframe). The extrapolated motion out of the field of view (where the arrow points) remains close to the ground truth. The maximum error does not exceed 4mm. Best seen by zooming in.

Figure 13 :

 13 Figure 13: Accuracy benchmark on the 3D US STACOM 2012 normal dataset, reporting the median tracking error over time for varying SNRs (blue, green, red curves). For the reference SNR (in blue), quartiles are overlayed (boxplots) to picture the dispersion of error values.

Figure 14 :

 14 Figure 14: Evolution of the inferred noise model for increasing Signal-to-Noise Ratios.

Figure 15 :

 15 Figure 15: Bull's eye plots of the radial, longitudinal and circumferential strain components at end-systole, averaged over AHA segments: estimated (top) and ground truth (bottom). A healthy case (left) and an ischemic case (right) are reported.

Figure 16 :

 16 Figure 16: Accuracy benchmark on the 3D tag STACOM 2011 dataset, reporting box-plots of tracking errors on all methodologies. The doted black line represents the average inter-observer variability.

Figure 17 :

 17 Figure 17: Strain at ES, computed from the 3D tag data of volunteer V9.

Figure 18 :

 18 Figure 18: Accuracy benchmark on the cine MR dataset, reporting median error over time along with quartiles. Surfaces reconstructed from slice-by-slice 3D+t segmentation serve as ground truth. Points on the discrete contours at time 0 are transported over time using the registration output. For each time step, pointto-surface distances over all 15 sequences and all contour points are pooled. Errors at time 0 are induced by the surface reconstruction.

Figure 20 :

 20 Figure 20: Histogram of inferred values for the regularity hyperparameter λ, pooled over all 15 sequences and 30 frames per sequence.

F

  -Σ l (Σ -1 ll µ l ) .(F.3)These rank one updates carefully avoid matrix-matrix products and have a O(|S| 2 ) complexity. In the case of the addition of a basis, we first compute the new column of Σ (resp. L) before updating its full body as:s l + κ l + A l ) -1 , Lll = (κ l + A l ) -1 . (F.7)Π l is the column vector of d × d matrices defined by Π l = Φ ( βH)φ l + λR l . The counterpart of Eq. (F.3) for addition is given by Eq. (F.8): μ = µ + Σl q l . (F.8)

  6)Thus, we always maintain the quantities t k and r k (for every basis) and recompute s k and q k either at no cost for inactive bases or, for bases in the active set S, in O(|S| • d). Updates of t k and r k upon deletion from the model, update or addition to the model of a basis l are done similarly to[START_REF] Tipping | Fast marginal likelihood maximisation for sparse bayesian models[END_REF], in O(|S| • d) per basis. For instance, in the addition case, it follows from Woodbury identities that tk = t k -Π k Σl Σk -Π k Σl r l (G.8)

Table 1 :

 1

Number of bases at each scale in the active parametrization of the displacement field (pooled over all frames and all sequences). Median, first and third quartiles are reported.

Table 3 :

 3 Number of bases at each scale in the active parametrization (pooled over all cases). Median, first and third quartiles are reported.

Table 4 :

 4

Number of bases in the active parametrization (pooled over all cases), at each scale. Median, first and third quartiles are reported.

q I (θ I ) = q I (θ I ) θ -I q -I (θ -I ) dθ -I = θ -I q(θ) dθ -I ≈ θ -I p(θ|D) dθ -I

http://segment.heiberg.se

http://hdl.handle.net/10380/3149

The minimization is w.r.t. the active, reduced parametrization and should be significantly faster and easier than its classic counterpart.

This outer product approximation conveniently guarantees positivity of the Hessian and is discussed in e.g.[START_REF] Bishop | Pattern recognition and machine learning[END_REF].

This is an abuse of speech for the sake of convenience, when we in fact manipulate blocks of d × d elements, lines and columns of such d × d elements.

We abuse notations here again for convenience. We mean to say thatA -1 k = 0 or that ∀n ∈ R d , n A k n → ∞.

Whenever necessary we identify square matrices (resp. vectors) of differing dimensionalities if one can be obtained from the other by padding with zeros, e.g. in (F.1) (F.2) (F.3), the rank-one correction leaves the lth column and line (resp. coefficient) of the right-hand side of the equation equal to zero.

Representation on a finite dictionary. The displacement field u(x) = 1≤k≤M φ k (x)w k is assumed to be expressed as a linear combination of bases φ k ∈ H (with associated weight w k ∈ R d ). By linearity of the representation, R D ( f ) = w Rw where the k, lth block element R kl is a d × d matrix, whose i, jth coefficient is D(φ k e i ) | D(φ l e j ) , with e 1 • • • e d the canonical frame. We wish to compute these coefficients efficiently.

Computation of R kl . If the basis functions are regularly translated versions of a single kernel φ, i.e. φ k (x) = φ(xx k ), then we have from the properties of the Fourier transform and Eq. (H.1):

where F -1 is the inverse Fourier transform and δx = x kx l .

Values for all offsets can be computed at once using the fast Fourier transform. If the basis offset is the voxel size, this runs in O(N log N). The values can be stored in an image, so that retrieving a coefficient R kl at a later stage is done in O(1), after recomputing the corresponding offset. If the dictionary of basis functions is generated from a few (say, K) such translation invariant kernels, the same process is repeated K 2 times, once for each couple of generating kernels.

Case of the Gaussian kernel. For Gaussian basis functions and standard differential operators D, closed form expressions can be obtained instead of relying on the fast Fourier transform. We illustrate this on the bending energy (D = ∆). Recall that for any p-uplet of integers

The Fourier transform of a Gaussian basis

ξ}, which is again Gaussian. Using Eq. (H.3), regrouping the exponential factors and using Eq. (H.4), we finally obtain for two Gaussian bases φ k and φ l (with respective variance S k and S l and centered at x k and x l ):

The derivative (-∆) 2 [K S k,l ] of the Gaussian kernel at any point is known in closed form.

Inner products φ k ( βH)Φ. For a block-diagonal matrix βH (with N blocks) and an arbitrary set of M basis functions, computing the inner product of φ k with all other bases functions is O(MN). If the dictionary of bases is generated from K translation invariant kernels, these inner products can be computed in O(KN log N). φ k ( βH) is computed in O(N), then K convolutions (one with each generating kernel, in O(N log N) each) yield the desired result. Similarly all of the φ k ( βH)φ k , k = 1 . . . M, can be computed by convolution of the image whose voxels are the d × d tensors βi H i with the square of the K generating kernels.