Pietro Gori
email: pietro.gori@inria.fr

Olivier Colliot

Linda Marrakchi-Kacem

Yulia Worbe

Cyril Poupon

Andreas Hartmann

Nicholas Ayache

Stanley Durrleman

A Bayesian Framework for Joint Morphometry of Surface and Curve meshes in Multi-Object Complexes

Keywords: shape, Bayesian, varifolds, fiber bundle, morphometry, complex, multi-object, atlas

We present a Bayesian framework for atlas construction of multi-object shape complexes comprised of both surface and curve meshes. It is general and can be applied to any parametric deformation framework and to all shape models with which it is possible to define probability density functions (PDF). Here, both curve and surface meshes are modelled as Gaussian random varifolds, using a finite-dimensional approximation space on which PDFs can be defined. Using this framework, we can automatically estimate the parameters balancing data-terms and deformation regularity, which previously required user tuning. Moreover, it is also possible to estimate a well-conditioned covariance matrix of the deformation parameters. We also extend the proposed framework to data-sets with multiple group labels. Groups share the same template and their deformation parameters are modelled with different distributions. We can statistically compare the groups' distributions since they are defined on the same space. We test our algorithm on 20 Gilles de la Tourette patients and 20 control subjects, using three sub-cortical regions and their incident white matter fiber bundles. We compare their morphological characteristics and variations using a single diffeomorphism in the ambient space. The proposed method will be integrated with the Deformetrica software package, publicly available at www.deformetrica.org.

Introduction

Morphological analysis of human organs based on Magnetic Resonance Imaging (MRI) or Computed Tomography (CT) scans is an important field in medical imaging. An example of clinical application is the identification and quantification of the effect of a pathology on the anatomy of the brain (i.e., hippocampal atrophy in the Alzheimer's disease [START_REF] Gerardin | Multidimensional classification of hippocampal shape features discriminates Alzheimer's disease and mild cognitive impairment from normal aging[END_REF]) or of the heart (i.e., Tetralogy of Fallot [START_REF] Mansi | A Statistical Model for Quantification and Prediction of Cardiac Remodelling: Application to Tetralogy of Fallot[END_REF], stroke or myocardial ischemia). This provides precious insights into the pathological development giving the possibility to identify biomarkers of disease progression, quantitative rules for disease evaluation, therapy planning and monitoring.

The analysis can be carried out either directly on images or on selected structures which are previously segmented as 3D meshes. The second strategy has the drawback that it depends on an accurate segmentation, which can be challenging in some cases, but it has also some advantages. First, it is possible to focus only on certain structures 10 belonging to a particular area of an organ. Second, it allows visualisation and direct analysis of the entire 3D geometry of a structure without scrolling through a stack of images. And third, it becomes possible to combine different imaging modalities since one can analyse together meshes obtained from structural MRI images, CT scans and diffusion MRI images. Furthermore, the two strategies can be combined together in an iconic-geometric setting [START_REF] Siless | Medical Image Computing and Computer-Assisted Intervention -MICCAI[END_REF][START_REF] Graciano Fouquier | Iconic-Geometric Nonlinear Registration of a Basal Ganglia Atlas for Deep Brain Stimulation Planning[END_REF].

The structures of an organ can be modelled as 3D surfaces or 3D curves. In the brain, surfaces are used to model sub-cortical nuclei or cortical surface [START_REF] Durrleman | Statistical models of sets of curves and surfaces based on currents[END_REF][START_REF] Auzias | Diffeomorphic Brain Registration Under Exhaustive Sulcal Constraints[END_REF]. In the heart they are used as models of the left, right ventricle or of the entire myocardium [START_REF] Mansi | A Statistical Model for Quantification and Prediction of Cardiac Remodelling: Application to Tetralogy of Fallot[END_REF][START_REF] Palit | Computational modelling of left-ventricular diastolic mechanics: Effect of fibre orientation and right-ventricle topology[END_REF][START_REF] Lombaert | Joint Statistics on Cardiac Shape and Fiber Architecture[END_REF]. Lungs and liver can also be modelled as 3D surfaces [START_REF] Gorbunova | Lung CT registration combining intensity, curves and surfaces[END_REF]. Curves are instead employed in the brain to model trajectories of ensembles of neural tracts, commonly called fiber bundles (O'Donnell et al., 2009; Durrleman et al., 2011a). In the heart they can be used as models of cardiac myofibers [START_REF] Palit | Computational modelling of left-ventricular diastolic mechanics: Effect of fibre orientation and right-ventricle topology[END_REF] whereas in the lungs or liver they may be employed to model blood vessels [START_REF] Gorbunova | Lung CT registration combining intensity, curves and surfaces[END_REF]. The geometrical representation of an entire organ may thus combine surface and curve meshes into a single multi-object complex, which we call shape complex.

There are several examples of mesh based morphometry in the literature. Most of them use a single-object approach since they select and analyse only one particular structure of an organ (Golland et [START_REF] Cury | Analysis of anatomical variability using diffeomorphic iterative centroid in patients with Alzheimer's disease[END_REF]. This strategy limits the extent of the clinical conclusion to the chosen object, thus neglecting the information given by the surrounding structures. It does not seem appropriate to study intricate and complex systems of inter-connected objects such as the human organs. There is a growing awareness that multi-object studies could be more befitting. In this paper, we focus on a general morphometric approach called atlas construction. It consists of estimating an average shape complex of a population called template complex and the deformations of the embedding space which warp the template complex to the shape complexes of every subject. Deformations put into correspondence the template with every shape complex. The template captures the common morphological characteristics of the population and the deformations capture the variations in shape across the population. Every deformation is a single 3D diffeomorphism (a smooth deformation with smooth inverse) which transforms the entire 3D space. All meshes are therefore deformed together whatever their number or kind. This kind of deformation preserves the spatial organization of the structures since no intersection, folding or shearing may occur. It is thus possible to study a human organ as a whole, and not as a set of independent components, analysing not only the shape of each structure, but also their relative position.

The deformations align each structure of the template complex with the corresponding ones of the subject shape complex. To this end, one needs to define a similarity measure for every structure. Different metrics have been proposed for surfaces and curves. Some of them assume that it is possible to find point-correspondences between homologous structures, called landmarks. This is usually hard (or even impossible) especially when working with groups of curves. In this paper, we have opted for a correspondence-free metric based on the framework of varifolds [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF] which can be used with both surfaces and curves. It requires neither point-correspondences nor curve-correspondences. Moreover, surfaces and curves are treated as instances of the same mathematical object. This simplifies the atlas construction and the following statistical analysis.

A standard procedure to estimate a multi-object atlas (Durrleman et al., 2014; Avants and Gee, 2004; Ma et al., 2010) involves an optimization scheme where all template structures are optimised together. Every structure is weighted by a scalar which balances the importance of the structure with respect to the other structures and with respect to a regularity term on the deformations. These weights are fixed by the user and the results are rather sensitive to them. Optimal values should be determined by cross-validation which is time consuming and computationally intensive. Moreover, the computational load increases with the number of structures and subjects under study, making sometimes infeasible the automatic estimate of the weights via cross-validation.

Another strategy consists in embedding the atlas procedure into a statistical setting. Using a generative model, every shape can be interpreted as the sum between a deformation of the template and noise [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF]. Moreover, considering the deformation parameters, template and noise as random variables, it is possible to use Maximum Likelihood (ML) or Maximum A Posteriori (MAP) estimations to infer the atlas. If the noise follows a Gaussian distribution, it turns out that the weight of a structure can be interpreted as the variance of the noise and that it can be automatically estimated with ML or MAP techniques [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF][START_REF] Zhang | Bayesian Estimation of Regularization and Atlas Building in Diffeomorphic Image Registration[END_REF][START_REF] Simpson | Probabilistic inference of regularisation in non-rigid registration[END_REF][START_REF] Folgoc | Sparse Bayesian Registration[END_REF]. Other authors also proposed to infer the level of regularisation from the data by modelling it as the variance of the deformation parameters or deformation field [START_REF] Simpson | Probabilistic inference of regularisation in non-rigid registration[END_REF][START_REF] Folgoc | Sparse Bayesian Registration[END_REF][START_REF] Risholm | Bayesian characterization of uncertainty in intra-subject non-rigid registration[END_REF][START_REF] Wassermann | Probabilistic Diffeomorphic Registration: Representing Uncertainty[END_REF]. However, most of these works are about registration and atlas construction of images. Few of them focus on shapes. Moreover, they are usually based on single-object complexes and the variance of the data noise as well as the regularity level are usually fixed by the user. Three examples are [START_REF] Ma | A Bayesian Generative Model for Surface Template Estimation[END_REF] where the authors introduced a surface single-object diffeomorphic atlas procedure using a hypertemplate and Durrleman et al. (2011a[START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF] where the authors employed a generative model adapted for complexes of only curves or only surfaces respectively.

Here, we propose to embed a diffeomorphic multi-object atlas construction into a Bayesian framework. Shape complexes can be composed of curves and surfaces and they are both modelled as varifolds. The noise of each structure is assumed to be a random Gaussian varifold which belongs to an infinite dimensional space. We define a finite-dimensional approximation space where the noise is projected to and where it is possible to define probability distributions. Furthermore, we use a prior distribution on the variance of the noise as in [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF] adapted for shape complexes modelled as varifolds. This allows us to automatically estimate all balancing weights, one for each structure, in one single atlas construction. Moreover, we assume that also the deformation parameters follow a Gaussian distribution. We use a prior distribution on their covariance matrix similarly to [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF] which allows us to estimate it automatically even with a small sample size. Furthermore, the estimated covariance matrix is always well-conditioned whatever the number of subjects and it can be directly employed to compute statistics about the population under study.

This paper extends our previous work [START_REF] Gori | Bayesian Atlas Estimation for the Variability Analysis of Shape Complexes[END_REF] by providing a more detailed methodological description and more comprehensive results. Furthermore, we take advantage of the proposed Bayesian framework to define a multi-population atlas construction. We use the subjects' clinical diagnosis to subdivide the data-set into different groups. Their deformation parameters are then modelled as different Gaussian distributions. Since all sets of deformation parameters are defined on the same space, we can quantitatively compare them and the moments of their distributions such as the covariance matrices. Thanks to the proposed Bayesian framework, this multi-population analysis can be performed at no additional cost, namely without tuning any new balancing weight.

The paper is organised as follows. In Section 2, we first formulate the proposed Bayesian framework in a rather general way, showing that it could be employed with different shape models and deformation settings. We then present how to compute a multi-population atlas construction. Afterwards, we describe how to model multi-object shape complexes with varifolds and how to define a varifold random variable. Then, we show the diffeomorphic framework used in the experiments and how to integrate it in our statistical setting. In the last part of Section 2, we present the gradient-descent scheme used to optimise the cost function and how to initialise the template of both surfaces and curves. Eventually, we test the robustness of our algorithm with respect to the hyperparameters of the prior distributions and we also show how the proposed method can be employed to assess morphological differences between two groups of subjects.

Methods

Bayesian atlas construction

We define the mesh of structure j belonging to subject i as S i j . The total number of subjects is N and the number of analysed structures is M. Every subject shape complex S i is defined, for the moment, in a generic way as the ensemble of meshes S i j of subject i. We model it as a deformed template complex φ i (T) plus a residual i . Both T and i are also defined as the ensembles of templates T j and residuals i j . This formulation is known as forward model where we assume that all elements belong to an algebraic structure where addition is defined (Durrleman et al., 2011b;[START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF][START_REF] Ma | Bayesian template estimation in computational anatomy[END_REF] :

S i = φ i (T) + i (1)
The deformation φ i , proper to subject i, can belong to any diffeomorphic framework present in the literature. The only requirement, for the scope of this paper, is that it deforms the entire ambient space, namely all structures of the template complex simultaneously. The goal of the atlas construction is to estimate the template complex and its morphological variations within the population of shape complexes. The variations are described by the ensemble of deformations {φ i } and each one of them is parametrised by a set of parameters α i ∈ R q . We assume that these parameters follow a Gaussian distribution with zero mean and covariance matrix Γ α :

α i ∼ N(0, Γ α ) p(α i |Γ α ) ∝ 1 |Γ α | 1/2 exp - 1 2 α T i Γ -1 α α i (2)
From Eq.2, we can notice that the distribution of α i is completely described by Γ α . We can thus rephrase the goal of the atlas construction as estimating the template complex T and the covariance matrix Γ α knowing the shape complexes of the population {S i } or, more formally, maximizing their joint posterior distribution:

{T * , Γ * α } = argmax T,Γ α p(T, Γ α |{S i }) (3)
This maximization is constrained by the fact that the deformed template φ i (T) should resemble to the shape complex S i i.e., the residual i should be small. This is a common problem in statistical learning and it is usually tackled by assuming that the residual follows a Gaussian distribution centred at zero. Maximizing its likelihood is equivalent to minimizing its squared norm. We will start by modelling the structures of the shape complexes S i and template complex T with landmarks. Every structure j has a number of landmarks equal to Λ j and the norm of the difference between two meshes of the same structure is defined as the square root of the sum of squared differences between pairs of landmarks (L 2 -norm, || • || 2 ). The likelihood of the residuals modelled with landmarks is thus defined as:

i j ∼ N(0, σ 2 j Id Λ j ) p( i j |σ 2 j ) ∝ 1 |σ 2 j | Λ j /2 exp        - 1 2σ 2 j ||S i j -φ i (T j )|| 2 2        (4)
In section 2.3, we will make clear how to adapt this equation to varifolds by changing only the norm || • || and the definition of the parameter Λ j . The variance of the noise of structure j depends only on the parameter σ 2 j which is structure-dependent and modelled as a random variable.

It is important to notice that in Eq.1 a value of α i equal to zero means no deformation. The mean of α i is set to zero because the template is supposed to represent the average of the shape complexes (i.e. S i j ∼ N(φ i (T j ), σ 2 j Id Λ j )). In other words, using the terminology of mixed models, the deformation parameters α i are the random effects associated to the fixed effect T which is the expected value of the shape complexes.

Assuming independence between all random variables and considering α i as an unobserved nuisance variable, we can rewrite Eq.3 as:

{T * , Γ * α , σ 2 * j } = arg max T,Γ α ,σ 2 j         N i M j p(T j , Γ α , σ 2 j , α i , S i j )dα i         (5)
The maximization of Eq.5 is not tractable analytically. A possible solution is to employ an iterative algorithm like the EM (Expectation Maximization) and to approximate the conditional distribution of the E step with a Dirac distribution at its mode. Let Θ = {Γ α , T, {σ 2 j }} be the parameters of interest and Z = {{α i }} the unobserved nuisance variables, it results:

{Θ * , Z * } = argmin Θ,Z -log p({S i }/Z, Θ) -log p(Z/Θ) -log p(Θ) (6)
This algorithm is an approximation of the MAP estimator which does not converge when applied to a dataset with low signal-to-noise ratio (SNR), as demonstrated in [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF]. In order to ensure the convergence of the iterative estimation scheme one may use sampling algorithms like MCMC (Markov Chain Monte Carlo) [START_REF] Allassonnière | Construction of Bayesian deformable models via stochastic approximation algorithm: a convergence study[END_REF]. Unfortunately, these methods require a great computational load and execution time.

Not using priors for σ 2 j and Γ α is equivalent to a ML estimation and can produce degenerate estimates. Instead, as demonstrated in [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF], the introduction of inverse Wishart distributions makes it possible to obtain good estimates even when the number of subjects is small. They regularize the estimates of σ 2 j and Γ α . As in [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF], the probability density functions employed in this paper are adapted versions of the standard inverse Wishart distributions:

σ 2 j ∼ W -1 (P j , w j ) p(σ 2 j ; P j , w j ) ∝ (σ 2 j ) -w j 2 exp        - 1 2 w j P j σ 2 j        (7)
Γ α ∼ W -1 (P α , w α ) p(Γ α ; P α , w α ) ∝ |Γ α | -wα 2 exp - 1 2 w α T r(P T α Γ -1 α ) (8)
The scalars w j , P j and w α are strictly positive and P α is a positive symmetric matrix. The entire probabilistic model is graphically described in Fig. 1. Using these priors for σ 2 j and Γ α , and assuming that the template T has a non-informative prior distribution, Eq.6 becomes:

M j=1 N i=1 1 2σ 2 j ||S i j -φ i (T j )|| 2 2 + P j w j N + M j=1 1 2 (w j + Λ j N) log(σ 2 j )+ 1 2 N i=1 (α i ) T (Γ α ) -1 α i + 1 2 (w α + N) log(|Γ α |) + w α 2 tr((Γ α ) -1 P α ) (9)
The framed terms refer respectively to the data-terms and to the regularity of the deformations. If we fix both σ 2 j and Γ α and we set w j = 0, their sum is the cost function minimised in standard atlas construction procedures (Avants and Gee, 2004; Durrleman et al., 2014) where a statistical framework is not employed. This minimisation is often referred to as the computation of the Fréchet mean. It can be noticed that σ 2 j weights the contribution of structure j as mentioned in the Introduction. The shape S i j is generated by a deformation φ i of the template T j parametrised by α i . The noise i j is parametrised by a scalar value σ 2 j and a constant matrix Σ j which depends on the computational model employed for S i j .

Multi-population atlas construction

In the previous subsection we assumed that all subjects belonged to the same group. In a clinical application, a subject can be characterised by a label which can be for instance his/her clinical status. Depending on this label, subjects can be divided into different groups. Here we propose to take advantage of this information computing a multi-population atlas construction.

Let L be the set of labels that can be given to a subject shape complex S i . We assume that L has only two elements for simplicity: L = {c, p}, they can be interpreted as "control" and "patient". The extension to more than two labels is straightforward and it will not be considered here. With this notation, we can rewrite Eq.1 for both groups as:

S c i = φ c i (T) + c i S p i = φ p i (T) + p i ( 10
)
where we assume that the template complex T is shared between the two populations. Instead, the deformation parameters α i are divided into two classes which can be modelled as distinct Gaussian distributions:

α c i ∼ N(µ c , Γ c α ) α p i ∼ N(µ p , Γ p α ) (11)
We can quantitatively compare α c i with α p i since we assume that they have the same dimensionality. To avoid identifiability issues, we force the sum of the means µ c and µ p to equal zero: µ c =µ and µ p =-µ. This constraint forces the estimated template complex to lie halfway between the two groups' averages.

We can then assume either that the two groups share the same variability (Γ c α =Γ p α =Γ α ) or that their within-group variations are different (Γ c α Γ p α ). The cost functions for both cases can be found in Appendix A. Using the second assumption, we could test the hypothesis of equality between the two covariance matrices or quantitatively compare their respective modes. On the contrary, if we ignored the labels of the subjects computing a single-population atlas, we would implicitly assume that the two different groups belong to the same population and that they share the same variability. Furthermore, we could also adopt another strategy by constructing an atlas for each group and obtaining therefore two different template complexes. However, in this case, we could neither quantitatively compare the deformation parameters of the two groups nor directly compare the resulting template complexes.

The multi-population extension is computationally feasible only when using the proposed Bayesian framework. As previously shown, we use the log-likelihood of the deformation parameters as regularity term of the deformations. Dividing the deformation parameters into two different groups entails the addition of another regularity term into the cost function (see Appendix A). This new term should be balanced by a trade-off parameter tuned using crossvalidation. If we increased the number of groups, we would also augment the number of regularity terms and balancing values to tune. With the proposed statistical setting, there is no need to tune these weights since they are automatically taken into account in the estimates of the covariance matrices.

The proposed statistical framework, for both single-and multi-population, is general and it could be applied to any shape model and parametric diffeomorphic setting, provided that one defines a probability density function for the noise. Modelling shapes with landmarks simplifies the computations and it also eases the definition of random variables. Unfortunately, the choice of corresponding points between structures is not always an easy task and in some cases, like for fiber bundles in the brain, it is almost impossible. It is better to opt for a correspondence-free metric like the one based on the framework of varifolds [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF]).

Varifolds

The framework of varifolds is an extension of the one of currents [START_REF] Vaillant | Surface Matching via Currents[END_REF] and it can be used for both surfaces and curves. It does not require point-correspondence between structures or curve-correspondence between ensemble of tracts. This means that the same structure may have a different number of vertices among subjects. It is robust to mesh imperfections such as holes, spikes or tract interruptions. Moreover, it does not require a consistent orientation of the normal or tangent vectors among the population since the distance is invariant to a change of orientation of some normals or tangent vectors. Another important characteristic is that varifolds prevent what is known as the "currents cancelling effect". This happens when two surface/curve elements with opposite orientation cancel reciprocally their effect on the estimation of the template. For a more detailed discussion about varifolds the user is referred to Durrleman et al. (2014); [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF]. Let X and Y be two curves (oriented or unoriented rectifiable curves), they can be modelled as varifolds based on:

V X (ω) = X ω(x, ← -→ β(x))|β(x)| 2 dx V Y (ω) = Y ω(y, ← -→ γ(y))|γ(y)| 2 dy ( 12
)
where V X (ω) and V Y (ω) are scalars, β(p) and γ(y) are the tangents of X and Y at the point x and y respectively and ← → β p is the unoriented unit vector with the same direction of β p (respectively ← → γ q for γ q ). The set of unoriented unit vectors ← → β ∈ RP 2 can be formally defined as the quotient of the unit sphere in R 3 by the two elements group {±Id 3 }. Instead, the test field ω ∈ W is a function of both the points in the ambient space R 3 and the unoriented unit tangent vectors in RP 2 .

The curves X and Y can be seen as polygonal lines of P and Q segments respectively. Every segment of X is completely described by its center point x p ∈ R 3 and tangent vector β p ∈ R 3 centred at x p (respectively y q ∈ R 3 and γ q ∈ R 3 for Y). Let c p be the length of the tangent vector β p (respectively d q for γ q ), one can approximate Eq.12 with:

V X (ω) ≈ P p=1 ω(x p , ← → β p )c p and V Y (ω) ≈ Q q=1 ω(y q , ← → γ q )d q .
It is important to notice that, even if this approximation seems reasonable, it is still an open question how to ensure its convergence when the number of segments tends to infinite, i.e., with a more accurate sampling of the curve [START_REF] Charlier | The fshape framework for the variability analysis of functional shapes[END_REF]. Furthermore, Eq.12 shows that the space of varifolds W * is a linear functional from W to R. This makes W * a vector space and therefore the union of meshes is equal to a sum in the space of varifolds. But, contrary to currents, inverting the orientation of a mesh does not change the varifold representation. One chooses W to be a separable Reproducible Kernel Hilbert Space (RKHS) on R 3 xRP 2 in order to have an explicit definition of distance between two meshes modelled as varifolds. More precisely, since W is a product space, one defines its kernel K W as the tensor product between two kernels k x , Gaussian, and k β , Cauchy-Binet, defined on R 3 and on RP 2 respectively. Thus, the inner product between X and Y modelled as varifolds is defined as:

V X , V Y W * = P p=1 Q q=1 k x (x p , y q )k β ( ← → β p , ← → γ q )c p d q = P p=1 Q q=1 exp       -||x p -y q || 2 λ 2 W              β T p γ q c p d q        2 c p d q (13)
The distance between V X and V Y is therefore:

||V X -V Y || 2 W * = V X , V X W * + V Y , V Y W * -2 V X , V Y W * .
This metric is completely parametrized by the standard deviation λ W of the Gaussian kernel k X which is fixed by the user. This framework can be easily extended also to surface meshes by replacing x p (y q ) and β p (γ q ) with the centres and normals of the faces of the mesh, as shown in [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF]. The computations of the gradient with respect to the points of the mesh can be found in [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF]; [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF].

Varifolds random variable. In Eq.4, we defined the residual as a multivariate Gaussian random variable on the finitedimensional Euclidean space of landmarks. Here we extend this definition to the framework of varifolds following the same line of reasoning presented in [START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF] for the framework of currents. We define a random Gaussian varifold as a linear map between every test field ω ∈ W to a real random Gaussian variable G(ω) such that, given two test fields ω 1 and ω 2 , E[G(ω 1 )] = 0 and E[G(ω 1 ), G(ω 2 )] = ω 1 , ω 2 W . This shows that the kernel K W of the space W completely defines the covariance matrix of the Gaussian varifold. However, since it is infinite-dimensional, it has no probability density function and therefore we can not simulate instances of it. To tackle this problem, we define for each structure j a finite-dimensional space W * Λ j on which we project the template T j and all the shapes {S i j } i=1...N modelled as varifolds. This finite dimensional space W * Λ j is defined as the span of a set of delta Dirac varifolds: Span{δ (x u , ← → β k ) } where both the points {x u } and the unoriented unit vectors { ← → β k } are constrained to belong to two predefined grids, respectively Υ x and Υ β . The first one is a linearly spaced grid in the ambient space and Υ β is a regular sampling of the half unit sphere in R 3 . The number of points of the two grids is respectively Λ j

x and Λ β and Λ j is their product where the index j refers to the j-th structure.

For every finite-dimensional space W * Λ j , it is possible to define a block matrix K W,Λ j whose blocks are the RKHS kernel

K W,Λ j (x u , ← → β k ), (•, •) u=1,...,Λ j x k=1,...,Λ β between every possible combination of the couples {x u , ← → β k }. The matrix K W,Λ j has a dimension of [Λ j xΛ j ].
To project a delta Dirac varifold δ (y, ← → α ) |α| onto W * Λ j , using for instance the closest neighbour projection, we look for the closest point to y of the grid Υ x and for the closest direction to ← → α among the ones given by Υ β and we assign the scalar |α| to that particular couple of grid points. The projection is therefore completely parametrised by the vector {|α|} of size Λ j . Using this scheme, the squared norm of a projected varifold

|| L l=1 δ (y l , ← → α l |α l |)|| 2 W * Λ j is equal to u∈Υ x p∈Υ x k∈Υ β q∈Υ β c uk K W,Λ j (x u , ← → β k ), (x p , ← → β q ) c
pq where the scalars c uk and c pq refer to the values obtained at the end of the projection of the varifold in {x u , ← → β k } and {x p , ← → β q } respectively. Scaling the matrix K W,Λ j with a scalar σ 2 j , we can define the likelihood of the residuals modelled as varifolds in the following way:

p( i j |σ 2 j ) ∝ 1 |σ 2 j | Λ j /2 exp        - 1 2σ 2 j ||S i j -φ i (T j )|| 2 W * Λ j        ( 14
)
where it is important to notice that K W,Λ j is a constant matrix that is computed only once at the beginning of the atlas construction. This equation has the same formulation of the likelihood of the residuals modelled as landmarks in Eq.4. The only two differences are the definitions of the norm and of the parameter Λ j . These are also the only changes to apply to the cost function in Eq.9 when using the framework of varifolds instead of landmarks.

In practice, due to the heavy computational load and time required to calculate || • || 2

W * Λ j
, we choose to compute the exact expression of the norm ||(•)|| 2 W * using a fast GPU (CUDA) implementation. Moreover, we define Υ j x as a regular linearly-spaced grid containing both shapes and template of structure j. The distance between every couple of points is given by λ W . The resulting number of points Λ j

x is used to compute Λ j =Λ j x . In this way, Λ j depends only on λ W which is fixed by the user and it is related to the distance at which two shapes can be compared and at the same time to the degree of detail of the anatomical features one wants to detect. In Sec.3.3, we show that the number of grid points, if above a reasonable threshold, does not affect the results and the proposed choice seems a good one.

Diffeomorphic Transformations

We define here how to compute the diffeomorphic deformations of the template complex. Our approach relies on the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework based on the control point formulation presented in Durrleman et al. (2011b). Template transformations are built by integrating a time-varying vector field v t (x) over t ∈ [0, 1] where v t (x) represents the instantaneous velocity of every point x belonging to the ambient space at time t. Calling φ t (x) the position of a point at time t which was located in x at time t = 0, its evolution is given by:

∂φ t (x) ∂t = v t (φ t (x)) with φ 0 (x) = x, which produces a flow of deformations {φ t } t∈[0,1]
. Furthermore, v t belongs to a RKHS with Gaussian kernel K V . In [START_REF] Glaunès | Transport par diffomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numrique[END_REF] it is shown that if such a vector field is square integrable in the interval [0,1], then every deformation of the flow {φ t } t∈[0,1] is a diffeomorphism. Here, we define v t using a dynamical system of C p control points c={c k ∈ R 3 } scattered in the ambient space and a set of time-varying vectors called momenta α={α k ∈ R 3 } linked to each control point. The velocity of every point in the space is given by:

ẋ(t) = v t (φ t (x)) = C p p=1 K V (x(t), c p (t))α p (t) φ 0 (x) = x(0) = x ( 15
)
where φ t (x)=x(t). With this definition, every φ t is a diffeomorphism provided that all α(t) are square integrable for every t. The evolution of a point depends only on its initial position and on the system L(t) = {c(t), α(t)}. In an atlas construction, the deformation of the template complex is the last one (φ 1 ) of the flow of diffeomorphisms {φ t } t∈[0,1] . Among all the possible paths connecting φ 0 to φ 1 we use the geodesic one, which means the one that minimizes the total kinetic energy along the path:

1 0 ||v t || 2 V dt = 1 0 C p k=1 C p p=1 α k (t) T K V (c k (t), c p (t))α p (t)dt (16)
It has been demonstrated in [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF] that the extremal paths are such that the system L(t) satisfies:

       ċk (t) = C p p=1 K(c k (t), c p (t))α p (t) αk (t) = - C p p=1 α k (t) T α p (t)∇ 1 K(c k (t), c p (t)) (17)
Denoting L(0)=L 0 ={c(0), α(0)} the initial condition of the system, the previous set of ODEs can be rewritten as: L(t) = F[L(t)] with L(0) = L 0 . This equation shows that the entire flow of diffeomorphisms, and more precisely the last one, is completely determined by the initial state of the system L 0 . Thus, given c(0) and α(0), one first integrates Eq.17 obtaining the values of c(t) and α(t) at every time t. Then, every point x in the ambient space is deformed by integrating Eq. 15.

In geometrical terms, the value L 0 defines the initial velocity of the geodesic path in the tangent space of our group of diffeomorphisms at identity. The parameters {c(0), α(0)} can be used to perform the so-called tangent-space statistics [START_REF] Vaillant | Statistics on diffeomorphisms via tangent space representations[END_REF]).

Optimization procedure

This diffeomorphic setting can be easily included into the proposed Bayesian framework. We choose to employ a unique set of control points c(t) for all subjects (and groups) and N subject-specific sets of momenta {α i (t)} i=1,...,N . It is important to notice that initial momenta of different groups are defined on the same control points and therefore they can be compared among each other. Initial control points and momenta {c(0), α i (0)} = {c 0 , α i0 } represent the deformation parameters that warp the template complex towards the shape complex of subject i. Controls points c 0 are considered as parameters of our model since they are a fixed effect of the entire population. Initial momenta are instead considered as hidden variables (i.e., Z = {{α i0 } i=1,...,N }. Assuming that all random variables are independent, the expression of the cost function in Eq.9 does not change except for the second framed term where we simply substitute α i with α i0 : 1 2 N i=1 (α i0 ) T (Γ α ) -1 α i0 . The variables T, {α i 0 } and c 0 are minimised using an accelerated version of the line search gradient descent method based on the Nesterov's scheme [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate $O(1/kˆ2$)[END_REF]. Instead, due to the use of conjugate priors, the optimal values for {σ 2 j } and Γ α can be computed in closed form:

Γα = N i=1 (α i0 )(α i0 ) T + w α P T α (w α + N) σ2 j = N i=1 ||S i j -φ i (T j )|| 2 W * Λ j + w j P j (w j + NΛ j ) (18)
The first parameter Γα is equal to a weighted sum between the sample covariance matrix of the initial momenta α i0 and the prior P α . A good choice for the prior seems to be:

P α = K -1 V
, where K V is a block matrix whose blocks are 3D Gaussian kernels between two different control points. This is the kernel of the RKHS to which the velocity field v t belongs to. If the number of subjects N is small with respect to w α then Γα ∝ K -1

V and this means that the deformation regularity in Eq.9 becomes N i=1 (α i0 ) T K V α i0 , which is the sum of the geodesic distances between the template complex and all its transformations. This is the deformation regularity term used in previous atlas constructions not based on a statistical setting [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF]. Moreover, the use of this prior makes the inversion of Γ α possible even when the number of subjects N is smaller than the number of deformation parameters.

The second parameter σ2 j is equal to a weighted sum between the data-term of the j-th structure and the prior. It balances the importance of structure j with respect to the other structures and at the same time it is a trade-off between the data-term and the other parts of the cost function (Eq.9). The automatic estimate of this parameter is rather important as shown in Sec.3.2. The use of an inverse Wishart prior prevents over-fitting. In fact, using a maximum likelihood estimator, the minimization process might focus only on structure k reducing its residuals almost to zero and consequently log(σ 2 k ) →-∞. This could also result in almost ignoring all the other structures. The role played by the prior is to impose a minimum value to σ 2 k in order to avoid such a situation. It is important to notice that the closed-form solutions for both Γα and σ2 j in Eq.18 would be different if using the standard Inverse Wishart distributions instead than the ones proposed in Eq.7 and Eq.8. They would not be equal to a weighted average between prior and empirical estimates, making therefore more difficult their interpretation and use. The gradients of the cost function E (Eq.9) with respect to T,{α i 0 } and c 0 are respectively equal to:

∇ T k E = N i=1 1 2σ 2 k ∇ T k [D ik ] ∇ α s0 E = M j=1 1 2σ 2 j ∇ α s0 [D s j ] + Γ -1 α α s0 (19)
∇ c 0 E = N i=1 M j=1 1 2σ 2 j ∇ c 0 [D i j ] + w α 2 ∇ c 0 tr((Γ α ) -1 P α )
where

D i j =||S i j -φ i (T j ))|| 2 W * Λ j
. The differentiation of the prior P α with respect to c 0 is not taken into account since its norm is negligible with respect to the one of the data term.

In order to compute the three derivatives of the data-terms D i j , we need first to compute the deformed template complex φ i (T) where the deformation φ i is the final diffeomorphism of the flow generated by integrating forward in time Eq.17 and then Eq.15 which can be rewritten in matrix notation as: Ṫi (t)=G[T i (t), L i (t)] with T i (0)=T(0)=T. After that, we can compute the gradient of D i j with respect to the deformed template complex T i (1) = φ i (T) and then bring back this information to t = 0 in order to update the initial template points T, control points c 0 and momenta α i0 . This is done by integrating backward in time from t=1 to t=0 a set of linearised ODEs called adjoint equations like in [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF]. It results:

∇ T E = N i=1 θ i (0) ∇ α i0 E = ξ α i (0) + Γ -1 α α i0 (20)
∇ c 0 E = N i=1 ξ c i (0)
where the auxiliary variables ξ i (t) = {ξ α i (t), ξ c i (t)} and θ i (t) satisfy the linearised ODEs:

θi (t) = -∂ T i G[T i (t), L i (t)] T θ i (t) θ i (1) = 1 2σ 2 ∇ T i (1) [ D i ] (21) ξi (t) = -∂ L i G[T i (t), L i (t)] T θ i (t) + d L i F[L i (t)] T ξ i (t) ξ i (1) = 0 (22)
ODEs are integrated with the Euler's method using 10 steps. More details can be found in the Appendix of Durrleman et al. (2014). The computations of the gradients for the multi-population extension can be found in Appendix A.

Template initialisation

Since we use a gradient descent scheme to minimize the cost function, we need an initial template. This should not be the mesh of a subject, as otherwise it would inherit all the imperfections due to segmentation errors and data noise, biasing consequently the atlas construction.

For surfaces, we use a regular sphere. It is first centred with respect to the ensemble of subject meshes. Then, it is scaled to an ellipsoid using the three main modes of variability of the meshes considered as points cloud.

For bundles of curves, we propose instead to initialise the template as the "most relevant" curves in the population. We first gather the curves F of every subject bundle into a raw initial template T = ∪ P h=1 F h , where P is the total number of curves. Then, we select the set of H curves which most resembles the raw template T . Calling {F} k one of the possible P H combinations of H curves out of P, one writes:

argmin k ||T -{F} k || 2 W * k ∈ P H (23)
In many applications the number of curves of the raw template can be huge, thus we have developed a greedy algorithm which does not test all the possible combinations of curves but it is based on an iterative scheme. At each step, we divide T in subsets of Z curves, we look for the best one and we remove it from the raw template T saving it in a new template. This process is repeated until the new template has H curves, where H is the average number of curves in the population and Z is a submultiple of H. At each iteration, the best subset of Z curves is selected minimising:

argmin k ||T - (k+1)Z h=(kZ)+1 F h || 2 W * k = 0, ..., (H/Z) -1 (24)
after having randomized the order of the curves in T . The new template is the initial template used in the atlas construction. In Sec.3.4, we will show that the proposed initialisation method for templates of fiber bundles is stable and that it produces similar results even when varying the value of Z.

Since the number of points per curve and the number of curves in the bundles can be different among the population, we have implemented this algorithm with the metric of varifolds, but any metric could be used.

An example of template initialisation applied to three brain fiber bundles is shown in Fig. 2. The last row represents the initial templates used in the following experiments.

Experiments

In this section, we first present the dataset used in the following experiments. Then, we evaluate the robustness of our Bayesian framework with respect to the values of the hyperparameters, and we compare it with the robustness of a previous atlas procedure without automatic estimate of σ 2 j . Furthermore, we also analyse the stability of the Bayesian atlas construction, as a function of the number of points of the varifold grid Λ j and of the initial template of fiber bundles T j . Eventually, we present how the proposed multi-population atlas construction approach can be used to highlight morphological differences between two groups.

Materials

The dataset used throughout this article contains 40 subjects: 20 controls and 20 patients subject to Gilles de la Tourette syndrome. This neurodevelopmental disorder is thought to be associated with dysfunctions of the corticostriato-pallido-thalamic circuits which are composed of sub-cortical structures linked to the cortical surface by fiber bundles. We will investigate the morphological characteristics of these structures which distinguish the group of patients from the one of controls and which might be due to the Gilles de la Tourette syndrome.

For each subject, we consider three sub-cortical structures (left caudate, left putamen and left thalamus) and the fiber bundles connecting them to the left hemisphere of the cortical surface. The three sub-cortical structures are segmented with FSL (Patenaude et al., 2011) from T1-weighted images (3T). Since we want to analyse as a single structure nucleus accumbens and caudate, we merge these two segmentations together. We use the marching cubes algorithm to create the 3D meshes of all structures (AimsMesh function of Brainvisa 4.3.0). Fiber bundles result from a probabilistic tractography method using 8 seeds [START_REF] Perrin | Fiber Tracking in q-Ball Fields Using Regularized Particle Trajectories[END_REF] applied to diffusion-weighted images (50 directions, B-factor=1000) using an analytical Q-Ball model to estimate the local underlying orientation distribution function (ODF) [START_REF] Descoteaux | Regularized, fast, and robust analytical Q-ball imaging[END_REF]. Every subject bundle is downsampled by randomly selecting 10% of its fibers and it is then approximated as proposed in [START_REF] Durrleman | Statistical models of sets of curves and surfaces based on currents[END_REF] using a Matching Pursuit Algorithm. For more details about the acquisition and tractography, the reader is referred to [START_REF] Worbe | Altered structural connectivity of cortico-striato-pallido-thalamic networks in Gilles de la Tourette syndrome[END_REF].

Robustness with respect to the hyperparameters

We evaluate here the robustness of the proposed algorithm with respect to the hyperparameters {w j , P j , w α }. We do not consider P α since, as explained in Sec.2.5, we fix it to P α = K -1

V . We compute 18 different atlases changing every time only one of the hyperparameters and keeping fixed the others at a certain value. Both w j and w α are previously normalised in order to use the same range of values, namely w = w NΛ and w α = w α N, giving as result:

Γα = 1 N N i=1 (α i0 )(α i0 ) T + w α P T α (1 + w α ) σ2 j = 1 NΛ j N i=1 ||S i j -φ i (T j )|| 2 W * + w j P j (1 + w j ) (25)
In order to understand a plausible range of testing values, we compute the residuals 1

NΛ j N i=1 ||S i j -φ i (T j )|| 2
W * using the aforementioned dataset and initial templates and we notice that the maximum value considering either only the sub-cortical structures or only the fiber bundles is never above 10 and 100 respectively. This means that the product w j P j should not surpass this limit, otherwise a reduction of the data-term would impact very little the estimate of σ2 j and consequently the minimization of the cost function. Noting that, we choose as fixed value 10 -2 and as testing range for all hyperparameters the interval [10 -2 -10 3 ] for sub-cortical structures and [10 -2 -10 4 ] for fiber bundles.

In order to test the robustness of the results, we compute the norm of the difference between the resulting template complex of all 18 atlases and a reference template complex. Moreover, we also analyse the norm of the residuals obtained at the end of the atlas constructions.

To show the importance of an automatic estimate of σ 2 j and Γ α , we compare our results with the ones obtained using a previous atlas procedure (Fixed) where both σ 2 j and Γ α are fixed. In order to be consistent with the proposed Bayesian algorithm, we set

Γ α = K -1
V , which is the value of the prior P α , and we let σ 2 j change along the same range of values as the hyperparameters. This means building 6 other atlases changing only the value of σ 2 j . In all 24 atlas procedures, the other parameters, like λ W and λ V , are kept equal. Results using only left caudates or left caudate bundles of the 20 controls are shown in Fig. 3 where we have used as reference template the one obtained using the Fixed method with σ 2 j =10 -1 . In Fig. 4 we compare the robustness of the results using simultaneously two structures in each atlas construction: left caudate and left caudate bundle. This means that we need to fix two σ 2 j values with the Fixed method and 5 hyperparameters in the Bayesian framework (w α , two w j and two P j ). Using as range of values [10 -2 -10 3 ], the number of possible combinations for the parameters of the Fixed technique is 36 whereas for the hyperparameters of the Bayesian framework is 7776. This makes infeasible an extensive and complete analysis of the robustness of the results. Nevertheless, we compute 30 atlases (15 Fixed + 15 Bayesian) picking randomly the values of the hyperparameters and σ 2 j in the range [10 -2 -10 2 ]. The employed values for both the hyperparameters and σ 2 j are listed in Table B.3. The reference template is the one obtained at the end of the 9-th case using the Fixed method. It is possible to conclude that the Fixed method leads to much more variable results than using Bayesian priors. The choice of the hyperparameters in our Bayesian framework is therefore easier than the one of σ 2 j . We decide to set them at: ω j =0.01Λ j N and P j = 0.05R 0 j ω j where R 0 j is the initial data-term of structure j. This brings to

σ 2 j = R j +0.05R 0 j NΛ j (1+0.01)
which means that the minimum value of σ 2 j is equal to about 5% of R 0 j NΛ j . This choice reduces the risk of over-fitting since a value of R j smaller than 0.05R 0 j NΛ j would almost not affect the estimate of σ 2 j . We also fix ω α to 10 -3 . This makes the contribution of P α negligible for the estimate of Γα , especially when N is big, but it still permits the inversion of Γα .

Robustness with respect to the number of points Λ of the varifold grid

Another parameter of our algorithm is the number of points Λ of the varifold grid where template and shape complexes are projected to. We decide to make it dependent on λ W , which is the bandwidth of the varifold kernel, as explained in Sec.2.3. In Fig. 5, we show that the results are not so sensitive to a change of the number of grid points and that our choice seems a good one. We compute 9 different atlases using the left-caudate surfaces of the 20 controls. All atlas constructions share the same set of parameters except for Λ which is tested in the range: [1-10 5 ]. As before, we analyse the norm of the difference between a reference template (Λ=1) and the estimated templates at the end of the atlas procedures. We also show the norm of the residuals and the logarithm of the determinant of the covariance matrix of the initial momenta (log(|Γ α |)). As it is possible to notice, there is a minimum value of grid points ( 1000) above which results are very similar and also more satisfactory.

In order to explain these results, it can be noticed that by replacing σ 2 j with σ2 j into Eq.9, one obtains a new data-term which is equal to (w + ΛN) log( σ2 ). When Λ is small, the reduction of the residuals, and consequently of log( σ2 ), is not enough amplified by the multiplication with Λ to compensate a great increase in the regularity term. Thus, the estimated initial momenta keep a norm close to zero throughout the atlas construction. The deformations of the template are therefore insignificant and the covariance matrix Γ α depends only on the prior, which explains why log(|Γ α |) does not vary. The number of grid points obtained with the proposed method is equal to: 3456.

Robustness with respect to the initial template of fiber bundles

The proposed initialisation procedure for templates of fiber bundles is iterative. At each iteration, a raw template is randomly subdivided in subsets of Z fibers. The best subset is then removed from the raw template and added to the employed initial template, as explained in Sec.2.6. Here, we show that the initial templates obtained using different Z values are very similar to each other. Moreover, we also verify that the results of the atlas construction are rather stable when using different initial templates.

First, we compute 40 initial templates for each fiber bundle considered in this paper: left-caudate bundle (LCD), left-putamen bundle (LPU) and left-thalamus bundle (LTH). We use four different values of Z: 25,50,100 and 150. For each value, we compute 10 different initial templates. Then, we compare all the resulting initial templates T with the one shown in Fig. 2, called T 0 . We use two dissimilarity measures. The first one is the varifold norm of the difference between every T and T 0 divided by the varifold norm of T 0 . The second one is the absolute value of the difference between the data-terms of every template D(T ) and D(T 0 ) divided by D(T 0 ), where D(T ) : 1

N N i=1 ||S i -φ i (T )|| 2
W * . Results are shown in Fig. 6. They are very stable and the difference between D(T ) and D(T 0 ) is always very small. Furthermore, we also compute 4 multi-population atlases with Γ c α =Γ p α using the left caudates and left-caudate bundles of all 40 subjects (20 controls and 20 patients). For each atlas, we use a different initial template for the left-caudate bundles, one for each group of Z. We then compare the results of the atlas constructions with respect to the ones of the first atlas, chosen as reference. We use 5 dissimilarity measures shown in Table 1. The subscript 0 refers to the reference atlas. Results are very stable and therefore we can conclude that they are not very sensitive to a change of initial template of fiber bundles.

Average [%] Standard deviation [%]

||CP-CP 0 || 2 ||CP 0 || 2 0.50 0.04

|D(T c )-D(T c 0 )| D(T c 0 )
0.51 0.41

|D(T b )-D(T b 0 )| D(T b 0 )
0.41 0.02

|σ 2 c -σ 2 c0 | σ 2 c0 0.51 0.21 |σ 2 b -σ 2 b0 | σ 2 b0 1.50 0.61
Table 1: Comparison between the results of four multi-population atlases computed using left caudate surfaces and left caudate bundles. Each atlas is characterised by a different initial template for the left caudate bundle. The subscript 0 refers to the first atlas, which is used as reference.

The variable CP refers to control points, the superscript c to the left caudate surface, the superscript b to the left caudate bundle and D(T ) to the data-term. 6: Robustness analysis with respect to the initial template of fiber bundles. We compute 40 initial templates for each fiber bundle employed in this paper: left-caudate bundle (LCD), left-putamen bundle (LPU) and left-thalamus bundle (LTH). Experiments are divided into four groups depending on the value of Z: 25 (10 templates), 50 (10 templates), 100 (10 templates) and 150 (10 templates). We use two dissimilarity measures, as explained in Sec.3.4, highlighted here with a continuous and dotted line respectively and with a different color for each bundle. Multi-population Atlas. (Γ c α Γ p α ) Here, we assume that the initial momenta of patients and controls do not share the same covariance matrix. To test this assumption, we use the test proposed in Srivastava and Yanagihara (2010). The p-value (0.3461) is too big to reject the null hypothesis: Γ c α = Γ p α . This might be due to the fact that the sample size (20 subjects per group) is too small with respect to the number of deformation parameters (3240). Another possible explanation could be that the two groups share most of their morphological variability which means that the first modes of variation of both covariance matrices are similar. To test this hypothesis, we first compute a Principal Component Analysis (PCA) for both covariance matrices. Then, we calculate the angles between the modes of variation of the two groups. The results for the first 5 modes are shown in Table 2. As it is possible to notice, the first mode of the group of controls is almost parallel to the first mode of the group of patients. This means that they produce similar morphological changes in the template complex. Moreover, they explain almost 35% and 45% of the total variability of their groups respectively (see Fig. 12). This might explain why the previous test fails in rejecting the hypothesis of equality between the two covariance matrices.

The first mode of both groups indicates the shared morphological variability between controls and patients. The second mode should instead produce shape variations which are more characteristic of each group (angle ∼ 60 • ) and which contribute also substantially (∼ 10%) to the total variability. In Fig. 13, we show the second mode of variation of both groups. Colors refer to the displacement from the final common template complex, shown in gray in the middle of the figure. It is possible to notice that the dorso-lateral part of the putamen varies almost in the same way between the two groups. Instead, the posterior-lateral part of the caudate and the ventral part of the thalamus present morphological variations which are different between the two groups.

Discussion and Conclusions

This work provides a Bayesian framework to embed a multi-object atlas construction into a statistical setting. It is general and it can be applied to any parametric deformation framework and to all shape models with which it is possible to define probability density functions. We can automatically estimate important balancing parameters which were fixed by the user in previous methods, namely the noise variances of every structure. We demonstrated that these parameters influence more the results than the hyper-parameters introduced with the proposed Bayesian priors. This statistical setting makes therefore multi-object atlas constructions more feasible and reproducible since the user is not obliged to tune or fix balancing weights which grow with the number of analysed structures.

With the proposed method, it is also possible to estimate a well-conditioned covariance matrix of the deformation parameters. In previous works [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF][START_REF] Avants | Geodesic estimation for large deformation anatomical shape averaging and interpolation[END_REF][START_REF] Ma | A Bayesian Generative Model for Surface Template Estimation[END_REF], the authors proposed to use the sample covariance matrix computed at the end of the atlas procedure. Since the number of subjects is usually smaller than the number of deformation parameters, the sample covariance matrix might be ill-conditioned. A standard solution is to regularize it a posteriori using an identity matrix multiplied by a small scalar. This strategy is not satisfactory since it means that the estimates of the deformation parameters are based on a degenerate sample covariance matrix or on a wrong covariance matrix. Moreover, the regularization seems also too simplistic. In our approach, we directly estimate a well-conditioned covariance matrix using a more natural and coherent regularisation term given by the kernel of the RKHS to which belongs the velocity field used to compute the diffeomorphisms. Furthermore, the estimate of the covariance matrix does not change if the regularity term, and its related prior, are multiplied by a scalar. This means that we automatically take into account also the so-called regularity level.

We also extend the proposed statistical setting to a multi-population strategy, where we take into consideration the affiliation of a subject to a group. We can thus employ more pertinent models at no additional cost, meaning without fixing or tuning new parameters. One of the main advantages is that we can estimate different distributions of initial momenta, one for each group. We can then quantitatively compare their averages and modes of variations.

It is interesting to notice that one could also compute the sample covariance matrices of different groups at the end of a single-population atlas and then compare their modes. However, this would not be consistent with the fact that in a single-population atlas the initial momenta of different groups are estimated as belonging to the same Gaussian distribution.

We propose to model both curve and surface meshes as Gaussian random varifolds. They can be seen therefore as instances of the same mathematical object. Since the space of varifolds is of infinite dimension, we create a grid for every structure j on which we project both shapes and template. We decide to make the number of grid points Λ j dependent on the varifold bandwidth λ j W . This choice does not influence the results in a single-object atlas, as demonstrated in Sec.3.3. In a multi-object analysis, the importance of structure j is highly determined by λ j W . In fact, the contribution of structure j mainly depends on the residuals and on Λ j , which are both related to λ j W . On the one hand, this choice simplifies the atlas construction, since the user has to fix only one parameter. On the other hand, it makes the choice of λ j W even more crucial. In future works, we will investigate how to automatically estimate λ j W . Another important parameter fixed by the user is the diffeomorphic kernel bandwidth λ V . It is chosen based on the desired registration accuracy and it defines the number of control points. Its estimate, together with the one of λ j W , would make the atlas construction completely automatic whatever the number or kind of meshes. A possible strategy could be to select a range of suitable values of λ V and then use sparse multi-scale methods [START_REF] Sommer | Sparse Multi-Scale Diffeomorphic Registration: The Kernel Bundle Framework[END_REF]. Otherwise, one could also estimate the best deformation modules from a dictionary using sparse techniques in order to disentangle the single complicated diffeomorphism into interpretable transformations [START_REF] Gris | A sub-riemannian modular approach for diffeomorphic deformations[END_REF]. This would augment the computational load and execution time but it would also make the analysis more objective.

Another improvement for our method would be the employment of sampling algorithms, like MCMC, to ensure the convergence of the EM scheme whatever the quality of the data. As already said, these methods require a great computational burden. A possible solution would be to use a GPU implementation based on Multi-Graphics Processing Units (Multi-GPU) as in [START_REF] Ha | Fast Parallel Unbiased Diffeomorphic Atlas Construction on Multi-Graphics Processing Units[END_REF]. The authors showed that a Multi-GPU implementation gives a computational gain up to sixty times faster than a single CPU implementation. We tested the proposed algorithm comparing two populations of shape complexes, one of controls and one of Gilles de la Tourette patients. Shapes consisted of three left sub-cortical nuclei and the fiber bundles connecting them to the left cortical surface. Results proved the effectiveness of our method in detecting morphological differences between the two populations. They could indicate atypical connections resulting from abnormal brain development due to Gilles de la Tourette syndrome. Future works will aim to confirm this study by adding more structures, like the cortical surface, and possibly more subjects. Moreover, we plan to apply our algorithm to other human organs such as the heart or the liver.

  Some examples are Gorczowski et al. (2010); Cates et al. (2008); Bossa et al. (2011); Qiu et al. (2010); Durrleman et al. (2014) where the authors have shown that one better discriminates controls from patients when also considering the relative position between structures.

Figure 1 :

 1 Figure 1: Graphical description of the proposed model. The variables in square boxes and circles are constants and random variables respectively.The shape S i j is generated by a deformation φ i of the template T j parametrised by α i . The noise i j is parametrised by a scalar value σ 2 j and a constant matrix Σ j which depends on the computational model employed for S i j .

Figure 2 :

 2 Figure 2: Template initialisation for three brain fiber bundles: left caudate, left putamen and left thalamus bundle. The top row contains the three raw initial templates where the fibers of every subject present in the population have been gathered together (see 3.1). The second row presents the initial templates obtained at the end of the proposed initialisation process cut at the intersection with their respective sub-cortical initial templates, transparent and coloured in red.

Figure 3 : 2 j

 32 Figure 3: Analysis of the robustness of the results when changing the hyper-parameters of the proposed algorithm (Bayesian) and the value of σ 2 j in an atlas without automatic estimates (Fixed). Every dot represents the result of an atlas construction. In all Bayesian estimations, it has been changed only one of the hyperparameters, fixing the others to 0.01. "Surfaces" are the left caudates and "Bundles" are the left-caudate bundles. Figures on the left represent the norms of the differences between the templates obtained at the end of the atlas constructions and a reference template (Fixed method with σ 2 j =0.01). Figures on the right refer to the norms of the residuals obtained at the end of the atlas procedures.

Figure 4 :

 4 Figure 4: Analysis of the robustness of the results when changing the hyper-parameters of the proposed algorithm (Bayesian) and the value of σ 2 j in an atlas without automatic estimates (Fixed). Every dot represents the result of an atlas construction. In all atlas estimations the template complex is composed of two structures: left caudate and left caudate bundle. The values of the hyperparameters for the Bayesian estimation and the values of the two σ 2 j for the Fixed method are listed in Table B.3. Figures on the left represent the norms of the differences between the templates obtained at the end of the atlas constructions and a reference template (9-th case of the Fixed method). Figures on the right refer to the norms of the residuals obtained at the end of the atlas procedures.

Figure 5 :

 5 Figure 5: Robustness of the results with respect to the number of points Λ of the varifold grid. We compute 9 atlases using 20 left-caudate surfaces and changing only Λ within the range [1-10 5 ]. Top-left figure represents the norm of the difference between a reference template (Λ=1) and the templates obtained at the end of the atlas constructions. Top-right figure shows the norm of the residuals during the atlas constructions. Bottom figure is about the evolution of log(|Γ α |) during the atlas procedures. The maximum number of iterations in every atlas construction is 150. Colors refer to the different grid sizes.

Figure

  Figure6: Robustness analysis with respect to the initial template of fiber bundles. We compute 40 initial templates for each fiber bundle employed in this paper: left-caudate bundle (LCD), left-putamen bundle (LPU) and left-thalamus bundle (LTH). Experiments are divided into four groups depending on the value of Z: 25 (10 templates), 50 (10 templates), 100 (10 templates) and 150 (10 templates). We use two dissimilarity measures, as explained in Sec.3.4, highlighted here with a continuous and dotted line respectively and with a different color for each bundle.

Figure 7 :

 7 Figure 7: Atlas construction procedure. From left to right, we present first the initial template complex and the initial set of control points. Then, we show the final template complex obtained at the end of the atlas construction. The top and bottom arrows point to the final template complex deformed accordingly to the average initial momenta of controls µ (top) and patients -µ (bottom). The averages of initial momenta (µ and -µ) are shown respectively above and below the two arrows.

Figure 8 :

 8 Figure 8: Shape dissimilarities of the sub-cortical structures between the two groups. Colors refer to the absolute value of the difference between the displacements from the final template complex to the average configurations of patients and controls. The two average configurations are obtained by deforming the final template complex with µ and -µ. The four frames represent the same three structures from different points of view. Letters {c, p, t} refer to caudate, putamen and thalamus respectively.

Figure 9 :

 9 Figure9: Shape dissimilarities of the fiber bundles between the two groups. Colors refer to the absolute value of the difference between the displacements from the final template complex to the average configurations of patients and controls. The two average configurations are obtained by deforming the final template complex with µ and -µ.

Figure 10 :

 10 Figure 10: Statistically significant differences between the average initial velocities of the two groups (controls and patients) computed at each vertex of the templates of the sub-cortical structures. P-values result from a Hotelling's two-sample T-squared test. In the second row, p-values are corrected using a false discovery rate (FDR) method. Letters A and P refer to the Anterior and Posterior part of the brain respectively.

Figure 11 :

 11 Figure 11: Statistically significant differences between the average initial velocities of the two groups (controls and patients) computed at each vertex of the templates of fiber bundles. P-values result from a Hotelling's two-sample T-squared test. In the second row, p-values are corrected using a false discovery rate (FDR) method. Letters A and P refer to the Anterior and Posterior part of the brain respectively.

Figure 12 :

 12 Figure 12: PCA of the two covariance matrices of initial momenta based on the assumption Γ c α Γ p α . The magenta line represents the percentage of variability explained by the first mode. The red line shows the number of components needed to explain 90% of the total variability. The green line indicates the magnitude of the eigenvalues and the blue line shows the explained variability using all the previous components.

Figure 13 :

 13 Figure 13: Second mode of variation at ± 1.5 standard deviation (σ) of two PCAs based on the two estimated covariance matrices of initial momenta Γ c α and Γ p α . Colors refer to the displacement from the final common template complex shown in the middle of the figure in gray. Arrows indicate the areas of the template complex which vary in a different way between the two groups.

Table 2 :

 2

				Patients
		1	2	3	4	5
		1 25.71	82.49	87.10	87.77	83.32
		2 87.70	60.39	72.57	83.62	80.30
	Controls	3 80.60	62.59	67.17	84.39	79.76
		4 86.74	75.15	82.60	88.98	77.35
		5 83.70	81.24	89.51	85.39	83.80

Angles in degrees between the first 5 modes of variation of the group of controls and the first five modes of the group of patients. The average angle between all modes is 85.39 • .

Table B .

 B 3: Values of the hyperparameters for the Bayesian estimation and values of the two σ 2 j for the Fixed method used in the robustness analysis shown in Fig.4.

	Case	Fixed			Bayesian
		σ 2 sur f ace	σ 2 bundle	w α	w sur f ace	P sur f ace	w bundle	P bundle
		0.01	0.01	1	10	0.1	10	0.1
		0.01	0.1	0.01	0.01	0.01	0.01	0.01
		0.01	1	0.1	0.1	0.1	0.1	0.1
		0.01	10	1	1	1	1	1
		0.01	100	10	10	10	10	10
		100	0.01	100	100	100	100	100
		100	0.1	10	1	0.1	0.01	100
		100	1	0.1	10	1	10	0.1
		100	10	0.01	100	10	1	0.01
	10	100	100	100	0.1	0.1	10	1
	11	1	0.01	0.1	0.1	100	1	10
	12	1	0.1	0.1	10	10	1	1
	13	1	1	0.01	1	10	0.1	100
	14	1	10	10	10	0.1	100	1
	15	1	100	10	1	0.01	0.1	1

Acknowledgements

The research leading to these results has received funding from the program Investissements d'avenir ANR-10-IAIHU-06.

Population differences

In all the following experiments, we use the same parameters. All structures, except for the caudate, have a varifold parameter λ W equal to 2mm. For the caudate we choose a value of 3mm. About the bandwidth of the diffeomorphic kernel K V , we use a value of 7mm with a consequent number of control points equal to 1080. The maximum number of iterations is 120 and all the computations are performed on a Intel Xeon, 32 cores, CPU E5-2650, 2.60GHz with a graphic card NVIDIA Quadro 5000. All shape complexes are previously rigidly registered to a reference shape complex.

Multi-population Atlas. (Γ c α =Γ p α =Γ α ) We compute a multi-population atlas using the groups of controls and patients. The template complex and control points are shared between the two groups. Initial momenta are assumed to follow two Gaussian distributions with opposite mean and equal covariance matrix. The overall processing time is about 100 hours based on a single-CUDA implementation.

In Fig. 7 we show a scheme of the process. On the left, there are the initial template complex and the initial set of control points, which is initialised as an equally-distanced grid covering the entire ambient space. The distance between couples of initial control points is equal to the bandwidth of the diffeomorphic kernel: λ V . Momenta α i0 are initially set to zero. At the end of the atlas construction we obtain: a final template complex, an updated set of control points, optimised subject-specific initial momenta, the average of the initial momenta µ and their common covariance matrix Γ α . The final template complex shows the morphological features common to both groups. At the top and at the bottom, we present the average of the initial momenta for each group: µ and -µ. Then, on the right, we show the final common template complex deformed according to µ and -µ. The two complexes represent the anatomical configurations typical of each group and they can be directly compared since they stem from the same final template complex. We take advantage of that by computing the absolute value of the difference between the displacements from the final template complex along the two average directions (i.e. µ and -µ). This is shown in Fig. 8 for subcortical structures and in Fig. 9 for fiber bundles. From the first figure, it is clear that the main differences are in the dorso-lateral part of the three sub-cortical structures, especially for the caudate. About the bundles, the differences are mainly in the central part of the caudate-cortico and putamen-cortico bundles.

A question that naturally arises is whether the shape dissimilarities shown in Fig. 8 and Fig. 9 are statistically significant. To check that, we compute two statistical tests: a global and a local one. The former is a permutation test which assesses whether the initial momenta of the two groups (controls and patients) come from the same distribution. We employ as statistics the Mahalanobis norm of the difference between the average initial momenta of the two groups. Calling α 1 and α 2 the initial momenta of the two groups at each permutation, ᾱ1 and ᾱ2 their respective averages and Γ 12 their pooled covariance matrix (computed as in Eq.A.3), the employed statistics is equal to: ( ᾱ1 -ᾱ2 ) T Γ -1

12 ( ᾱ1 -ᾱ2 ). We perform 10.000 permutations and the resulting p-value is lower than 0.05. Thus, we can reject the null hypothesis of identical distribution. This shows that the average anatomical complexes of the two groups have a statistically significantly different global shape.

In addition to a multivariate test, we also perform mass univariate tests to produce local maps of significance. First, we compute the initial velocities at each vertex of the template complex (Eq.15) using the estimated initial momenta of both groups. For each vertex and for each group, we obtain 20 3D vectors which are normally distributed. This comes from the fact that initial momenta follow a prior Gaussian distribution. At each vertex, we can then use a Hotelling's two-sample T-squared test for assessing the statistical significance of the difference between the average initial velocities of the two groups. The ensemble of resulting p-values, one for each vertex, is corrected for multiple comparisons using a false discovery rate (FDR) method [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF]. Results are shown in Fig. 10 and in Fig. 11 for sub-cortical structures and fiber bundles respectively. After the FDR correction, the only structures which still have statistically significant differences at the 5% level are the caudate (dorsal part), the thalamus (ventral and dorsal part) and the caudate bundle (central part). This is probably due to the relatively small number of subjects employed in the experiments.

Appendix A. Multi-population atlas construction

Here we compute the cost function and gradients for the multi-population extension of the atlas construction.

Appendix A.1. Γ c α =Γ p α =Γ α Let N c and N p be the number of subjects of the two groups and N=N c +N p . Moreover, we define α c i0 ∼ N(µ, Γ α ) and α p i0 ∼ N(-µ, Γ α ). Eq.9 becomes:

where the norm || • || can be both the L 2 -norm and the varifold one. The gradients with respect to T, c 0 and {σ 2 j } are exactly the same as in Eq.19 and Eq.18. The one with respect to α i0 depends on the group and it is equal to:

There is a closed form solution for both µ and Γ α . They are equal to:

As before, let N c and N p be the number of subjects of the two groups and N=N c +N p . Moreover, we define α c i0 ∼ N(µ, Γ c α ) and α p i0 ∼ N(-µ, Γ p α ). Priors on Γ p α and Γ p α are both equal to W -1 (P α , w α ). In this case Eq.9 becomes:

where the norm || • || can be both the L 2 -norm and the varifold one. The gradients with respect to T, c 0 and {σ 2 j } are exactly the same as in Eq.19 and Eq.18. The one with respect to α i0 depends on the covariance matrix of the group and its formulation is equal to the one in Eq.A.2. Furthermore, there is a closed form solution for µ, Γ c α and Γ p α . For the two covariance matrices, it is equal to the one in Eq.18 with a sum over N c and N p respectively. Instead, the optimal value for the average of the initial momenta is:

(A.5)

Appendix B. Robustness Analysis

We present in Table B.3 the values of the parameters and hyperparameters used in the experiments shown in Fig. 4.