Stanley Durrleman

Marcel Prastawa

Nicolas Charon

Julie R Korenberg

Sarang Joshi

Guido Gerig

Alain Trouvé

Morphometry of anatomical shape complexes with dense deformations and sparse parameters

Keywords: morphometry, deformation, varifold, anatomy, shape, statistics

We propose a generic method for the statistical analysis of collections of anatomical shape complexes, namely sets of surfaces that were previously segmented and labeled in a group of subjects. The method estimates an anatomical model, the template complex, that is representative of the population under study. Its shape reflects anatomical invariants within the dataset. In addition, the method automatically places control points near the most variable parts of the template complex. Vectors attached to these points are parameters of deformations of the ambient 3D space. These deformations warp the template to each subject's complex in a way that preserves the organization of the anatomical structures. Multivariate statistical analysis is applied to these deformation parameters to test for group differences. Results of the statistical analysis are then expressed in terms of deformation patterns of the template complex, and can be visualized and interpreted. The user needs only to specify the topology of the template complex and the number of control points. The method then automatically estimates the shape of the template complex, the optimal position of control points and deformation parameters. The proposed approach is completely generic with respect to any type of application and well adapted to efficient use in clinical studies, in that it does not require point correspondence across surfaces and is robust to mesh imperfections such as holes, spikes, inconsistent orientation or irregular meshing.

The approach is illustrated with a neuroimaging study of Down syndrome (DS). Results demonstrate that the complex of deep brain structures shows a statistically significant shape difference between control and DS subjects. The deformation-based modeling is able to classify subjects with very high specificity and sensitivity, thus showing important generalization capability even given a low sample size. We show that results remain significant even if the number of control points, and hence the dimension of variables in the statistical model, are drastically reduced. The analysis may even suggest that parsimonious models have an increased statistical performance.

The method has been implemented in the software Deformetrica, which is publicly available at www.deformetrica.org.

Introduction

Non-invasive imaging methods such as Magnetic Resonance Imaging (MRI) enable analysis of anatomical phenotypic variations over large clinical data collections. For example, MRI is used to reveal and quantify effects of pathologies on anatomy, such as hippocampal atrophy in neurodegenerative diseases or change in neuronal connectivity in neurodevelopmental disorders. Subject-specific digital anatomical models are built from the segmentation and labeling of structures of interest in images. In neuroanatomy, these structures of interest are often volumes whose boundaries take the form of 3D surfaces. For a given individual, the set of such labeled surfaces, which we call an anatomical complex, is indicative of the shape of different brain objects and their relative position. Our goal is to perform statistics on a series of such anatomical complexes from subjects within a given population. We assume that the complex contains the same anatomical structures in each subject, so that interindividual differences are not due to the presence or absence of a structure or a split of one structure into two. The quantification of phenotypic variations across individuals or populations is crucial to find the anatomical substrate of neurologic diseases, for example to find an early biomarker of disease onset or to correlate phenotypes with functional or genotypic variables. Not only the quantification, but also the description of the significant anatomical differences are important in order to interpret the findings and drive the search for biological pathways leading to pathologies.

The core problem is the construction of a computational model for such shape complexes that would allow us to measure differences between them and to analyze the distribution across a series of complexes. Geometric morphometric meth-ods make use of the relative position of carefully defined homologous points on surfaces, called landmarks [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF][START_REF] Dryden | Statistical Shape Analysis[END_REF]. Landmark-free methods often use geometric characteristics of the surfaces. They therefore need to make strong assumptions about the topology of the surface, for example limiting analysis to genus zero surfaces [START_REF] Chung | Deformation-based surface morphometry applied to gray matter deformation[END_REF][START_REF] Boyer | Algorithms to automatically quantify the geometric similarity of anatomical surfaces[END_REF] or using medial representations [START_REF] Styner | Morphometric analysis of lateral ventricles in schizophrenia and healthy controls regarding genetic and disease-specific factors[END_REF][START_REF] Bouix | Hippocampal shape analysis using medial surfaces[END_REF][START_REF] Gorczowski | Multi-object analysis of volume, pose, and shape using statistical discrimination[END_REF] or Laplace-Beltrami eigenfunctions [START_REF] Reuter | Laplace-Beltrami spectra as 'Shape-DNA' of surfaces and solids[END_REF]. Such methods can rarely be applied to raw surface meshes resulting from segmentation algorithms since such meshes may include small holes, show irregular sampling or split objects into different parts.

More important, such methods analyze the intrinsic shape of each structure independently, therefore neglecting the fact that brain anatomy consist of an intricate arrangement of various structures with strong interrelationships. By contrast, we aim at measuring differences between shape complexes in a way that can account for both the differences in shape of the individual components and the relative position of the components within the complex. This goal cannot be achieved by concatenating the shape parameters of each component or by finding correlations between such parameters [START_REF] Tsai | A shape-based approach to the segmentation of medical imagery using level sets[END_REF][START_REF] Gorczowski | Multi-object analysis of volume, pose, and shape using statistical discrimination[END_REF], as such approaches do not take into account the fact that the organization of the shape complex would not change, and in particular, that different structures must not intersect.

One way to address this problem is to consider surfaces as embedded in 3D space and to measure shape variations induced by deformations of the underlying 3D space. This idea stems from Grenanders group theory for modeling objects [START_REF] Grenander | General Pattern Theory: a Mathematical Theory of Regular Structures[END_REF], which revisits morphometry by the use of 3D space deformations. The similarity between shape complexes is then quantified by the "amount" of deformation needed to warp one shape complex to another. Only smooth and invertible 3D deformations (i.e., diffeomorphisms) are used, so that the internal organization of the shape complex is preserved during deformation since neither surface intersection nor shearing may occur. The approach determines point correspondences over the whole 3D volume by using the fact that surfaces should match as a soft constraint. The method is therefore robust to segmentation errors in that exact correspondences among points lying on surfaces are not enforced. In this context, a diffeomorphism could be seen as a low-pass filter to smooth shape differences. In this paper, it is our goal to show that the deformation parameters capture the most relevant parts of the shape variations, namely the ones that would distinguish between normal and disease.

Here, we propose a method that builds on the implementation of Grenanders theory in the LDDMM framework [START_REF] Miller | Geodesic shooting for computational anatomy[END_REF][START_REF] Vaillant | Diffeomorphic metric surface mapping in subregion of the superior temporal gyrus[END_REF][START_REF] Mclachlan | Discrete mechanics and optimal control for image registration[END_REF]. The method has 3 components: (i) estimation of an average model of the shape complex, called the template complex, which is representative of the population under study; (ii) estimation of the 3D deformations that map the template complex to the complex of each subject; and (iii) statistical analysis of the deformation parameters and their interpretation in terms of variations of the template complex. The first two steps are estimated simultaneously in a combined optimization frame-work. The resulting template complex and set of deformations are now referred to as an atlas.

Previous attempts to estimate template shapes in this framework offered little control over the topology of the template, whether it consists in the superimposition of a multitude of surface sheets [START_REF] Glaunès | Template estimation from unlabeled point set data and surfaces for computational anatomy[END_REF] or a set of unconnected triangles [START_REF] Durrleman | Statistical models of sets of curves and surfaces based on currents[END_REF]. The topology of the template may be chosen as one of a given subject's complex [START_REF] Ma | Bayesian template estimation in computational anatomy[END_REF], but this topology then inherits the mesh imperfections that result from an individual segmentation. In this paper, we follow the approach initially suggested by [START_REF] Durrleman | Topology preserving atlas construction from shape data without correspondence using sparse parameters[END_REF], which leaves the choice of the topology of the template with number of connected components to the user. This method estimates the optimal position of the vertices so that the shape of the template complex is an average of the subjects complexes.

Here, we extend this approach in order to guarantee that no selfintersection could occur during the optimization.

The set of deformations that result from warping the template complex to each subjects complex captures the variability across subjects. The deformation parameters quantify how the subjects anatomy is different from the template, and can be used in a statistical analysis in the same spirit as in [START_REF] Vaillant | Statistics on diffeomorphisms via tangent space representations[END_REF] and Pennec (2006). We follow the approach initiated in [START_REF] Durrleman | Optimal data-driven sparse parameterization of diffeomorphisms for population analysis[END_REF][START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF], which uses control points to parameterize deformations. The number of control points is fixed by the user, and the method automatically adjusts their position near the most variable parts of the shape complex. The method therefore offers control over the dimension of the shape descriptor that is used in statistics, and thus avoids an unconstrained increase with the number of surfaces and their samplings [START_REF] Vaillant | Surface matching via currents[END_REF]. We show that statistical performance is not reduced by this finite-dimensional approximation and that the parameters can robustly detect subtle anatomical differences in a typical low sample size study. We postulate that in some scenarios, the statistical performance can even be increased, as the ratio between the number of subjects and the number of parameters becomes more favorable.

An important key element of the method is a similarity metric between pairs of surfaces. Such a metric is needed to optimize the deformation parameters that enable the best matching between shape complexes. We use the varifold metric that has been recently introduced in [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF]. It extends the metric on currents [START_REF] Vaillant | Surface matching via currents[END_REF] in that it considers the non-oriented normals of a surface instead of the oriented normals. The method is therefore robust to possible inconsistent orientation of the meshes. It also prevents the "canceling effect" of currents, which occurs if two surface elements with opposite orientation face each other, and which may cause the template surface to fold during optimization. Otherwise, the metric inherits the same properties as currents: it does not require point-correspondence between surfaces and is robust to mesh imperfections such as holes, spikes or irregular meshing [START_REF] Vaillant | Surface matching via currents[END_REF][START_REF] Durrleman | Statistical models of sets of curves and surfaces based on currents[END_REF].

This paper is structured as follows to give a self-contained presentation of the methodology and results. We first focus on the main steps of the atlas construction, while discussing the technical details of the theoretical derivations in the appendices.

We then present an application to neuroimage data of a Down syndrome brain morphology study. This part focuses on the new statistical analysis of deformations that becomes possible with the proposed framework, and it also presents visual representations that may support interpretation and findings in the context of the driving clinical problem. The analysis also includes an assessment of the robustness of the method in various settings.

Mathematical Framework

Kernel formulation of splines

In the spline framework, 3D deformations φ are of the form φ(x) = x + v(x), where v(x) is the displacement of any point x in the ambient 3D space, which is assumed to be the sum of radial basis functions K located at control point positions {c k } k=1,...,N cp :

v(x) = N cp X k=1 K(x, c k )↵ k .
(1)

Parameters ↵ 1 ,...,↵ N cp are vector weights, N cp the number of control points and K(x, y) is a scalar function that takes any pair of points (x, y) as inputs. In the applications, we will use the Gaussian kernel K(x, y) = exp(-|x -y| 2 /σ 2 V ), although other choices are possible such as the Cauchy kernel K(x, y) = 1/(1 + |x -y| 2 /σ 2 V ) for instance. It is beneficial to assume that K is a positive definite symmetric kernel, namely that K is continuous and that for any finite set of distinct points {c i } i and vectors {↵ i } i :

X i X j K(c i , c j )↵ i T ↵ j ≥ 0, (2)
the equality holding only if all ↵ i vanish. Translation invariant kernels are of particular interest. According to Bochner's theorem, functions of the form K(xy) are positive definite kernels if and only if their Fourier transform is a positive definite operator, in which case (2) becomes a discrete convolution. This theorem enables an easy check if the previous Gaussian function is indeed a positive-definite kernel, among other possible choices.

Assuming K is a kernel allows us to define the pre-Hilbert space V as the set of any finite sums of terms K(., c)↵ for vector weights ↵. Given two vector fields v 1 = P i K(., c i )↵ i and v 2 = P j K(., c 0 j )β j , (2) ensures that the bilinear map

hv 1 , v 2 i V = X i X j K(c i , c 0 j )↵ i T β j (3)
defines an inner-product on V. This expression also shows that any vector field v 2 V satisfies the reproducing property:

hv, K(., c)↵i V = v(c) T ↵, (4)
defined for any point c and weight ↵. The space of vector fields V could be "completed" into a Hilbert space by considering possible infinite sums of terms K(., c)↵, for which (4) still holds. Such spaces are called Reproducing Kernel Hilbert Spaces (RKHS) [START_REF] Zeidler | Applied Functional Analysis: Application to Mathematical Physics[END_REF].

Using matrix notations, we denote c and α (resp. c 0 and β) in R 3N (resp. R 3M ) the concatenation of the 3D vectors c i and ↵ i (resp. c 0 j and β j ), so that the dot product (3) writes hv 1 , v 2 i V = α T K(c, c 0 )β, where K(c, c 0 ) is the 3N ⇥ 3M matrix with entries K(c i , c 0 j )I 3⇥3 .

Flows of diffeomorphisms

The main drawback of such deformations is their noninvertibility, as soon as the magnitude of v(x) or its Jacobian is "too" large. The idea to build diffeomorphisms is to use the vector field v as an instantaneous velocity field instead of a displacement field. To this end, we make the control points c k and weights ↵ k to depend on a "time" t that plays the role of a variable of integration. Therefore, the velocity field at any time t 2 [0, 1] and space location x is written as:

v t (x) = N cp X k=1 K(x, c k (t))↵ k (t) (5)
for all t 2 [0, 1] A particle that is located at x 0 at t = 0 follows the integral curve of the following differential equation:

dx(t) dt = v t (x(t)), x(0) = x 0 , (6)
This equation of motion also applies for control points. Using matrix notations, their trajectories follow the integral curves of

ċ(t) = K(c(t), c(t))α(t), c(0) = c 0 . (7)
At this stage, point trajectories are entirely determined by timevarying vector weights ↵ k (t) and initial positions of control points c 0 .

For each time t, one may consider the mapping x 0 ! φ t (x 0 ), where φ t (x 0 ) is the position at time t of the particle that was at x 0 at time t = 0, namely the solution of (6). At time t = 0, φ 0 = Id R 3 (i.e., φ 0 (x 0 ) = x 0 ). At any later time t, the mapping is a 3D diffeomorphism. Indeed, it is shown in [START_REF] Miller | Geodesic shooting for computational anatomy[END_REF] that (6) has a solution for all time t > 0, provided that timevarying vectors ↵ k (t) are square integrable. It is also shown that these mappings are smooth, invertible and with smooth inverse. In particular, particles cannot collide, thus preventing self-intersection of shapes. At any space location x, one can find a particle that passes by this point at time t via backward integration, thus preventing shearing or tearing of the shapes embedded in the ambient space.

For a fixed set of initial control points c 0 , the time-varying vectors α(t) define a path (φ t ) t in a certain group of diffeomorphisms, which starts at the identity φ 0 = Id R 3 , and ends at φ 1 , the latter representing the deformation of interest. We aim to estimate such a path, so that the mapping φ 1 brings the template shapes as close as possible to the shapes of a given subject. The problem is that the vectors, which enable us to reach a given φ 1 from the identity, are not unique. It is natural to choose the vectors that minimize the integral of the kinetic energy along the path, namely 1 2

Z 1 0 kv t k 2 V dt = 1 2 Z 1 0 α(t) T K(c(t), c(t))α(t)dt. (8)
We show in Appendix A that the minimizing vectors α(t), considering c(0) and c(1) fixed, satisfy a set of differential equations. Together with the equations driving motion of control points (7), they are written as:

8 > > > > > > > > > > < > > > > > > > > > > : ċk (t) = N cp X p=1 K(c k (t), c p (t))↵ p (t) ↵k (t) = - N cp X p=1 ↵ k (t) T ↵ p (t)r 1 K(c k (t), c p (t)) (9) Denoting S(t) = c(t) α(t)
! the state of the system of control points at time t, (9) could be written in short as

Ṡ(t) = F(S(t)), S(0) = c 0 α 0 ! . ( 10
)
The flow of deformations is now entirely parameterized by initial positions of control points c 0 and initial vectors α 0 (called momenta in this context). Integration of (10) computes the position of control points c(t) and momenta α(t) at any time t from initial conditions. Control points and momenta define, in turn, a time-varying velocity field v t via (5). Any configuration of points in the ambient space, concatenated into a single vector X 0 , follows the trajectory X(t) that results from the integration of (6). Using matrix notation, this ODE is written as Ẋ(t) = v t (X(t)) = K(X(t), c(t))α(t) with X(0) = X 0 , which can be further shortened to:

Ẋ(t) = G(X(t), S(t)), X(0) = X 0 (11)
A given set of initial control points c 0 defines a sub-group of finite dimension of our group of diffeomorphisms. Paths of minimal energy, also called geodesic paths, are parameterized by initial momenta α 0 , which play the role of the logarithm of the deformation φ 1 in a Riemannian framework. Integration of (10) computes the exponential map. It is easy to check that kv t k V is constant along such geodesic paths. Therefore, the length of the geodesic path that connects φ 0 = Id R 3 to φ 1 (i.e., R 1 0 kv t k V dt) simply equals the norm of the initial velocity (i.e., kv 0 k V ).

Varifold metric between surfaces

Deformation parameters c 0 , α 0 will be estimated so as to minimize a criterion measuring the similarity between shape complexes. To this end, we define a distance between surface meshes in this section, and show how to use it for shape complexes in the next section. If the vertices in two meshes correspond, then the sum of squared differences between vertex positions could be used. However, finding such correspondences is a tedious task and is usually done by deforming an atlas to the meshes. This procedure leads to a circular definition, since we need this distance to find deformations between meshes! Among distances that are not based on point correspondences, we will use the distance on varifolds [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF]. In the varifold framework, meshes are embedded into a Hilbert space in which algebraic operations and distances are defined. In particular, the union of meshes translates to addition of varifolds. The inner-product between two meshes S and S 0 is given as:

⌦ S, S 0 ↵ W ⇤ = X p X q K W (c p , c 0 q ) ⇣ n T p n 0 q ⌘ 2 n p n 0 q (12)
where c p and n p (resp. c 0 q and n 0 q ) denotes the centers and normals of the faces of S (resp. S 0 ). The norm of the normals n p equals the area of the mesh cell. K W is a kernel, typically a Gaussian function with a fixed width σ W .

The distance between S and S 0 then simply writes:

d W (S, S 0 ) 2 = kS -S 0 k 2 W ⇤ = hS, Si W ⇤ + hS 0 , S 0 i W ⇤ -2 hS, S 0 i W ⇤ .
One notices that the inner-product, and hence the distance, does not require vertex correspondences. The distance measures shape differences in the difference in normals directions, by considering every pair of normals in a neighborhood of size σ W . It considers meshes as a cloud of undirected normals and therefore does not make any assumptions about the topology of the meshes; one mesh may consist of several surface sheets, have small holes or have irregular meshing. Differences in shape at a scale smaller than the kernel width σ W are smoothed, thus making the distance robust to spikes or noise that may occur during image segmentation. The inner-product resembles the one in the currents framework [START_REF] Vaillant | Surface matching via currents[END_REF][START_REF] Durrleman | Statistical models of sets of curves and surfaces based on currents[END_REF], except that (n

T p n 0 q ) 2 |np||n 0 q | now replaces ⇣ n T p n 0 q ⌘
. With this new expression, the distance is invariant if some normals are flipped. It does not require the meshes to have a consistent orientation. Contrary to other correspondence-free distance such as the Hausdorff distance, the gradient of this distance with respect to the vertex positions is easy to compute, which is particularly useful for optimization.

We explain now how (12) is obtained. In the varifold framework, one considers a rectifiable surface embedded in the ambient space as an (infinite) set of points with undirected unit vectors attached to them. The set of undirected unit vectors is defined as the quotient of the unit sphere in R 3 by the two elements group {±Id R 3 }, and is denoted ! S . We denote ! u the class of u 2 R 3 in ! S , meaning that u, u/ |u| and -u/ |u| are all considered as the same element: ! u . In a similar construction as the currents, we introduce square-integrable test fields ! which is function of space position x 2 R 3 and undirected unit vectors ! u 2 ! S . Any rectifiable surface could integrate such fields ! thanks to:

S(!) = Z Ω S !(x, -! n(x)) |n(x)| dx, ( 13
)
where x denotes a parameterization of the surface S over a domain Ω S , and where n(x) denotes the normal of S at point x.

This expression is invariant under surface re-parameterization. It shows that the surface is a linear form on the space of test fields W. The space of such linear forms, denoted W ⇤ the dual space of W, is the space of varifolds.

For the same computational reasons as for currents, we assume W to be a separable RKHS on R 3 ⇥ ! S with kernel K chosen as:

K ⇣⇣ x, ! u ⌘ , ⇣ y, ! v ⌘⌘ = K W (x, y) u T v |u||v| ! 2 . ( 14
)
It is the same kernel as currents for the spatial part K W , and a linear kernel for the set of undirected unit vectors.

The reproducing property (4) shows that:

!(x, -! n(x)) = ⌧ !, K ✓ (x, -! n(x)), (., .) ◆ W .
Plugging this equation in ( 14) leads to

S(!) = ⌧ !, R Ω S K ✓ (x, -! n(x)), (., .) ◆ |n(x)| dx W .
The second part of the inner-product could be then identified with the Riesz representant of the varifold S in W, denoted L -1 W (S). Therefore, the inner-product between two rectifiable surfaces S and S 0 is hS, S 0 i

W ⇤ = S(L -1 W (S 0 )) = Z Ω S Z Ω S 0 K W (x, x 0 ) n(x) T n(x 0 ) |n(x)||n(x 0 )| ! 2 |n(x)| n 0 (x) dxdx 0 (15)
The expression in ( 12) is nothing but the discretization of this last equation.

For S a rectifiable surface and φ a diffeomorphism, the surface φ(S) can still be seen as a varifold. Indeed, a change of variables shows that for [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF]. Therefore, the varifold metric can be used to search for the diffeomorphism φ that best matches S to S 0 by minimizing d W (φ(S), S 0 ) 2 = kφ(S) -S 0 k 2 W ⇤ . In practice, the deformed varifold is computed by moving the vertices of the mesh and leaving unchanged the connectivity matrix defining the mesh cells. This scheme amounts to an approximation of the deformation by a linear transform over each mesh cell. Therefore, the distance kφ(S) -S 0 k 2 W ⇤ is only a function of X(1), i.e. A(X( 1)), where we denote X 0 the concatenation of the vertices of the mesh S and X(1) the position of the vertices after deformation. Indeed, from the coordinates in X(1), we can compute centers and normals of faces of the deformed mesh that can be then plugged into (12) to compute the distance kφ(S) -S 0 k 2 W ⇤ . Note that the varifold framework extends to 1D mesh representing curves in the ambient space, by replacing normals by tangents. In its most general form, varifold is defined for submanifolds with tangent-space attached to each point and uses the concept of Grassmannian [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF].

! 2 W, φ(S)(!) = S(φ?! ) where φ?!(x, ! n ) = (d x φ -1 ) T n !(φ(x), ------! (d x φ -1 ) T n)

Distances between anatomical shape complexes

The above varifold distance between surface meshes extends to a distance between anatomical shape complexes. An anatomical complex O is the union of labeled surface meshes, each label corresponding to the name of an anatomical structure. Meshes are pooled according to their labels into S 1 ,...,S N , where each S k contains all vertices and edges sharing the same label k. Let O 0 = {S 0 1 ,...,S 0 N } be another shape complex with the same number N of anatomical structures, but where the number of vertices and connected components in each S 0 k may be different than in S k . The similarity measure between both shape complexes is then defined as the weighted sum of the varifold distance between pairs of homologous structures:

d W (O, O 0 ) 2 = N X k=1 1 2σ 2 k S k -S 0 k 2 W ⇤ (16)
The values of σ k balance the importance of each structure within the distance. They are set by the user.

This distance cannot be used 'as' in a statistical analysis, since it is too flexible and, by construction, does not penalize changes in the organization of shape complexes. The idea is to use the distance on diffeomorphisms as a proxy to measure distances between shape complexes, the distance on varifolds being used to find such diffeomorphisms. Let O and O 0 be two shape complexes and {φ t } t2[0,1] be a geodesic path connecting φ 0 = Id R 3 to φ 1 such that, φ 1 (O) = O 0 . We can then define the distance between O and O 0 as the length of this geodesic path, which equals the norm of the initial velocity field v 0 . Formally, we define:

d φ (O, O 0 ) 2 = kv 0 k 2 V = α T 0 K(c 0 , c 0 )α 0 , (17)
for a given set of initial control points c 0 and with α 0 such that φ α 0 1 (O) = O 0 . However, it is rarely possible to find such a diffeomorphism that exactly matches O and O 0 . It is even not desirable since such a matching will be likely to capture shape differences that are specific to these two shape complexes and that poorly generalize to other instances. We prefer to replace the expression in (17) with the following relaxed formulation:

d φ (O, O 0 ) 2 = α T 0 K(c 0 , c 0 )α 0 with α 0 = arg min α d W (φ α 1 (O), O 0 ). ( 18
)
In this expression, the distance between varifolds d W is used to find the deformed shape complex φ 1 (O) that is the closest to the target complex O 0 and the distance in the diffeomorphism group between O and φ 1 (O) quantifies how far the two shape complexes are. The minimizing α 0 gives the relative position of φ 1 (O) (which is similar to O 0 ) with respect to O 0 .

In the following, O will represent the template shape complex that will be a smooth mesh with a simple topology and regular meshing. By construction, the deformed template φ 1 (O) is as smooth and regular as the template itself, whereas the subjects' shape complex O 0 may have irregular meshing, small holes, spikes, etc.. On the one hand, d W is flexible and loose in the sense that it measures a global discrepancy between the deformed template φ 1 (O) and the observation O 0 , but does not provide an accurate and computable description of the shape differences. On the other hand, d φ captures only shape differences that are consistent with a smooth and invertible deformation of the shape complex O, leaving in the residual norm d W (φ 1 (O), O 0 ) all other differences including noise and such very small scale mesh deformations. Deformations can be seen as a smoothing operator that captures only certain kind of shape variations and encode them into a descriptor α 0 , which will be used in the statistical analysis. The varifold metric d W allows us to compute this distance d φ without the need to smooth meshes, to build single connected components, to control for mesh quality, etc.

Atlas construction method

We are now in a position to introduce the estimation of an atlas from a series of anatomical shape complexes segmented in a group of subjects. An atlas refers here to a prototype shape complex, called a template, a set of initial control points located near the most variable parts of the template and momenta parameterizing the deformation of the template to each subject's complex.

For N su subjects, let {O 1 ,...,O N su } be a set of N su surface complexes, each complex O i being made of labeled meshes S i,1 ,...,S i,N . We define the template shape complex, denoted O 0 , as a Fréchet mean, which is defined as the minimizer of the sample variance:

O 0 = arg min O P i d φ (O, O i ) 2 .
The computation of d φ in (18) requires the estimation of a diffeomorphism φ by minimizing the varifold metric d W (φ(O), O i ). The combination of the two minimization problems leads to the optimization of the single joint criterion:

E(O 0 , c 0 {α i 0 }) = N su X i=1 ( N X k=1 1 2σ 2 k d W (φ α i 0 1 (S 0,k ), S i,k ) 2 + (α i 0 ) T K(c 0 , c 0 )α i 0 ) . (19)
The sum

P N i=1 (α i 0 ) T K(c 0 , c 0 )α i 0 = P N i=1 v i 0 2 is the sample vari-
ance. This term attracts the template complex to the "mean" of the observations. The other term with the varifold metric acts on the deformation parameters so as to have the best matching possible between the template complex and each subject's complex. The weights σ k can be now interpreted as Lagrange multipliers. The momentum vectors α i 0 parameterize each templateto-subject deformation. We assume here that they are all attached to the same set of control points c 0 , thus allowing the comparison of the momentum vectors of different subjects in the statistical analysis.

We further assume that the topology of the template complex is given by the user, so that the criterion depends only on the positions of the vertices of the template meshes. The number of control points is also set by the user, so that the criterion depends only on the positions of such points. In practice, the user gives as input of the algorithm a set of N meshes (typically ellipsoid surface meshes) whose number of vertices and edges connecting the vertices are not be changed during optimization. The user also gives a regular lattice of control points as input of the algorithm. Optimization of (19) finds the optimal position of the vertices of the template meshes, the optimal position of the control points and the optimal momentum vectors.

Let S i 0 = {c 0,p , ↵ i 0,p } denote the parameters of v i 0 , and X 0 the vertices of every template surface concatenated into a single vector. The flow of diffeomorphisms results from the integration of N su differential equations, as in (10): Ṡi (t) = F(S i (t)) with S i (0) = S i 0 . As in (11), X 0 follows the integral curve of N su differential equations:

Ẋi (t) = G(X i (t), S i (t)) with X i (0) = X 0 . The final value X i (1) = φ v i 0 1 (X 0 )
gives the position of the vertices of the deformed template meshes, from which we can compute centers and normals of each face of the deformed meshes, pool them according to mesh labels and compute each term of the kind d W (φ

α i 0 1 (S 0,k ), S i,k
) 2 using the expression in ( 12). Therefore, the varifold term essentially depends on the vector X i (1) and is denoted A(X i (1)). By contrast, the norm of the initial velocity, α i 0 T K(c 0 , c 0 )α i 0 depends only on the initial conditions S i 0 and is written as L(S i 0 ). The criterion ( 19) can be rewritten now as:

E(X 0 , {S i 0 }) = N su X i=1 ⇣ A(X i (1)) + L(S i 0 ) ⌘ , s.t. ( Ṡi (t) = F(S i (t)) S i (0) = S i 0 Ẋi (t) = G(X i (t), S i (t)) X i (0) = X i 0 . (20)
We notice that the parameters to optimize are the initial conditions of a set of coupled ODEs and that the criterion depends on the solution at time t = 1 of these equations.

8 > > > > > < > > > > > :
r α i 0 E = ξ ↵,i (0) + r α i 0 L(S i 0 ) r c 0 E = N su X i=1 ⇣ ξ c,i (0) + r c 0 L(S i 0 ) ⌘ , r X 0 E = P N su i=1 θ i (0),
where the auxiliary variables ξ i (t) = {ξ c,i (t), ξ ↵,i (t)} (of the same size as S i (t)) and θ i (t) (of the same size as X 0 ) satisfy the linear ODEs (integrated backward in time):

8 > > > < > > > : θi (t) = - ⇣ @ 1 G(X i (t), S i (t)) ⌘ T θ i (t), θ i (1) = r X i (1) A ξi (t) = -@ 2 G(X i (t), S i (t) T θ i (t) -d S i (t) F T ξ i (t), ξ i (1) = 0 .
Data come into play only in the gradient of the varifold metric with respect to the position of the deformed template r X i (1) A (derivation is straightforward and given in Appendix C). This gradient indicates in which direction the vertices of the deformed template have to move to decrease the criterion. This decrease could be achieved in two ways, by optimizing the shape of the template complex or the deformations matching the template to each complex. The vector θ i transports the gradient back to t = 0 where it is used to update the position of the vertices of the template complex. The vector ξ i interpolates at the control points the information in θ i , which is located at the template points, and is used at t = 0 to update deformation parameters. A striking advantage of this formulation is that one single gradient descent optimizes simultaneously the shape of the template complex and deformation parameters.

By construction, only the positions of the vertices of the template shape complex are updated during optimization. The edges in the template mesh remain unchanged, so that no shearing or tearing could occur along the iterations. However, the method does not guarantee that the template meshes do not selfintersect after an iteration of the gradient descent. To prevent such self-intersection, we propose to use a Sobolev gradient instead of the current gradient, which was computed for the L 2 metric on template points X 0 . The Sobolev gradient for the metric given by a Gaussian kernel K X with width σ X , is simply computed from the L 2 gradient as:

r X x 0,k E = N su X i=1 N x X p=1 K X (x 0,k , x 0,p )✓ i p (0). ( 21
)
We show in Appendix D that this new gradient r X E is the restriction to X 0 of a smooth vector field u s . Denoting X 0 (s) the positions of the vertices of the template meshes at iteration s of the gradient descent, we have that X 0 (s) = s (X 0 (0)) where s is the family of diffeomorphisms integrating the flow of u s . At convergence, the template meshes, therefore, have the same topology as the initial meshes.

Eventually, the criterion is minimized using a line search gradient descent method. The algorithm is initialized with template surfaces given as ellipsoidal meshes, control points located at the nodes of a regular lattice and momenta vectors set to zero (i.e., no deformation). At convergence, the method yields the final atlas: a template shape complex, optimized positions of control points and deformation momenta.

Computational aspects 2.6.1. Numerical schemes

The criterion for atlas estimation is minimized using a line search gradient descent method combined with Nesterov's scheme [START_REF] Nesterov | Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements[END_REF]. Differential equations are integrated using a Euler scheme with prediction correction, also known as Heun's method, which has the same accuracy as the Runge-Kutta scheme of order 2. Sums over the control points or over template points are computed using projections on regular lattices and FFTs using the method in Durrleman (2010, Chap. 2).

The method has been implemented in a software called "Deformetrica", which can be downloaded freely at www. deformetrica.org.

Parameter setting

The method depends on the kernel width for the deformation σ V , for the varifolds σ W and for the gradient σ X , as well as the weights σ k that balance each data term against the sum of squared geodesic distances between template and observations. The kernel widths σ V and σ W compare with the shape sizes. The varifold kernel width σ W needs to be large enough to smooth noise and to be sensitive to differences in the relative position between meshes (Durrleman, 2010, Ch. 1); otherwise values that are too small tend to make the shapes orthogonal. However, too large values tend to make all shape alike and therefore alter matching accuracy. The deformation kernel width σ V compares with the scale of shape variations that one expects to capture. Deformations are built essentially by integrating small translations acting on the neighborhoods of radius σ V . With smaller values, the model considers more independent local variations and the information in larger anatomical regions is not well integrated. With larger values, the model is based on almost rigid deformations.

The value of σ X is essentially a fraction of σ V : σ V or 0.5σ V work well in practice. The weights σ k are chosen so that data terms have the same order of magnitude as the sum of squared geodesic lengths. Values that are too small over-weight the importance of the data term and prevent the template from converging to the "mean" of the shape set. Values that are too large alter matching accuracy and thus shape features captured by the model.

A reasonable sampling of control points is reached for a distance between two control points being equal to the deformation kernel width σ V . Finer sampling often induces a redundant parameterization of the velocity fields as shown in [START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF]. Nonetheless, coarser sampling also may be sufficient for the description of the observed variability, as shown in the next experiments.

Kernel widths are chosen after few trials to register a pair of shape complexes. The weights σ k were then assessed while building an atlas with 3 subjects. The initial distribution of the control points was always chosen as the nodes of a regular lattice with step σ V or a down-sampled version of it. We always keep σ X = 0.5σ V . A qualitative discussion about the effects of parameter settings can also be found in [START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF].

We will show that the method works well without fine parameter tuning and that statistical results are robust with respect to changes in parameter settings.

Application to a Down syndrome neuroimaging study

We evaluate our method on a dataset of 3 anatomical structures segmented from MRIs of 8 Down syndrome (DS) subjects and 8 control cases. The hippocampus, amygdala and putamen of the right hemisphere (respectively in green, cyan and orange in figures) form a complex of grey matter nuclei in the medial temporal lobe of the brain. This study aims to detect complex non-linear morphological differences between both groups, thus going beyond size analysis, which already showed DS subjects to have smaller brain structures than controls [START_REF] Korenberg | Down syndrome phenotypes: the consequences of chromosomal imbalance[END_REF][START_REF] Mullins | Dementia in Down's syndrome: an MRI comparison with Alzheimer's disease in the general population[END_REF]. Whereas our sample size is small in view of standard neuroimaging studies, the previous findings in neuroimaging of DS suggest large morphometric differences. We therefore hypothesize that such differences would also be reflected in the shapes of anatomical structures, so that the proposed method could demonstrate its Control points that were initially on a regular lattice move to the most variable place of the shape complex during optimization.

Arrows parameterize space deformations and are used as a shape descriptor of each subject in the statistical analysis.

strength to differentiate intra-group variability from inter-group differences. To discard any linear differences, including size, we co-register all shape complexes using affine transforms.

We then construct an atlas using all data, setting σ V = 10 mm, σ W = 5 mm, σ X = σ V /2 and σ k = σ V for all nuclei, and control points initially located at the nodes of regular lattice of step σ V , yielding a set of 105 points. Robustness of results with respect to these values is discussed in Sec. 3.6.

The resulting template shape complex (Fig. 1-a) averages the shape characteristics of every individual in the dataset. The position of each subject's anatomical configuration (either DS or controls) with respect to the template configuration is given by initial momentum vectors located at control point positions (arrows in Fig. 1). These momentum vectors lie in a finitedimensional vector space, whose dimension is 3 times the number of control points. Standard methods for multivariate statistics can be applied in this space. The resulting statistics are expressed in terms of a set of momentum vectors. The template shape complex can be deformed in the direction pointed by the statistics via the integration of the geodesic shooting equations (10) followed by the flow equations ( 11). This procedure, also known as tangent-space statistics, is a way to translate the statistics into deformation patterns, and hence eases the interpretation of the results.

In the following sections, we show how such statistics can be computed and visualized, using the Down syndrome data as a case study.

Group differences

The first step is to show the differences between healthy controls (HC) and DS subjects that have been captured by the atlas. We compute the sample mean of the momenta for each group separately:

α HC = 1 N HC su P i2HC α i and α DS = 1 N DS su P i2DS α i ,
where HC (resp. DS) denotes the set of indices corresponding to healthy controls (resp. DS subjects). We then deform the template complex in the direction of both means, thus showing anatomical configurations that are typical of each group (Fig. 2). The figure shows that nuclei of DS subjects are turned toward the left part of the brain, with another torque that pushes the hippocampus tail (its posterior part) toward the superior part of the brain, and the head toward the inferior part. These two torques are more pronounced near the hippocampus/amygdala boundary than in the hippocampus tail or upper putamen region. The DS subjects' amygdala also has lesser lateral extension than that of the controls.

We perform Linear Discriminant Analysis (LDA) to exhibit the most discriminative axis between both groups in the momenta space. For this purpose, we compute the initial velocities of the control points v i = K(c 0 , c 0 )α i . The sample covariance matrix of these velocities, assuming equal variance in both groups, is given by:

Σ = 1 N su 0 B B B B B @ X i2HC (v i -v HC ) ⇣ v i -v HC ⌘ T + X i2DS (v i -v DS )(v i -v DS ) T 1 C C C C C A .
The direction of the most discriminative axis in the velocity space is defined as

v LDA ± = v ± Σ -1 (v HC -v DS ) where v = 1 2 (v HC + v DS ).
The associated momentum vectors are given as:

α LDA ± = K(c 0 , c 0 ) -1 v LDA ± .
The anatomical configurations are generated deforming the template shape complex in the two directions α LDA ± . We normalize the directions, so that their norm equals the norms between the means:

α LDA ± W ⇤ = α HC -α DS W ⇤ .
Therefore, the sum of the geodesic distance between the template complex and each of the deformed complexes is twice the norm between the means.

Results in Fig. 3 reveal similar thinning effects and torques as in Fig. 2. The figure also shows that putamen structures of DS subjects are more bent than those of controls.

Remark 3.1. Note that if the number of observations is smaller than 3 times the number of control points, then Σ is not invertible, and we use instead the regularized matrix Σ + "I 3 . In practice, we use " = 10 -2 , which leads to a condition number of the covariance matrix of order 1000. Statistics are not altered if this number is increased to 0.1 and 1, for which the condition number become 100 and 10 (results not shown). Remark 3.2. Note that we perform the statistical analysis using the velocity field sampled at the control points v = K(c 0 , c 0 )α and the usual L 2 inner-product. However, it would seem more natural to use the RKHS metric on the momenta α instead. Using the RKHS metric amounts to using ṽ = K 1/2 α so that the inner-product becomes (ṽ i ) T ṽ j = α i T K(c 0 , c 0 )α j , which is the inner-product between the velocity fields in the RKHS V. One can easily check that without regularization (" = 0), the most discriminant axis is the same in both cases, as will be the LDA and ML classification criteria introduced in the sequel. Using the identity matrix as a regularizer for the sample covariance matrix above amounts to using the matrix K(c 0 , c 0 ) -1 as a regularizer in the RKHS space. More precisely, the matrix Σ + "I 3 becomes Σ + "K(c 0 , c 0 ) -1 where Σ is the sample covariance matrix of the ṽi 's. It is natural to use this regularizer, since the criterion for atlas construction precisely assumes the momentum vectors to be distributed with a zero-mean Gaussian distribution with covariance matrix K(c 0 , c 0 ) -1 (which leads to 19)). For this reason, the same matrix is used in [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF] as a prior in a Bayesian estimation framework.

Anterior Right

v i 0 2 V = α i 0 T K↵ i 0 in (

Statistical significance

We estimate the statistical significance of the above group differences using permutation tests in a multivariate setting. In our experiments, the number of subjects is always smaller than the dimension of the concatenated momentum vectors, which is 3 times the number of control points. In this case the distribution Hotelling T 2 statistics cannot be computed and we use permutations to give an estimate of this distribution.

Let (u k , λ 2 k ) be the eigenvectors and eigenvalues sorted in decreasing order of the sample covariance matrix Σ (without regularization, i.e., " = 0). We truncate the matrix up to the N modes largest eigenvalues that explain 95% of the variance: Σ = P N modes k=1 λ 2 k u k u T k . Its inverse is given by: Σ-1 =

P N modes k=1 1 λ 2 k u k u T k .
We then compute the T 2 Hotelling statistics as:

T 2 = N su -2 4 (v HC -v DS ) T Σ-1 (v HC -v DS ).
To estimate the distribution of the statistics under the null hypothesis of equal means, we compute the statistics for 10 5 permutations of the subjects' indices i. Each permutation changes the empirical means and within-class covariance matrices, and thus the selected subspace and the statistics. The resulting pvalue equals p = 2.6 10 -4 , thus showing that our shape descriptors are significantly different between DS and HC subjects at the usual 5% level. The anatomical differences highlighted in Fig. 2 and 3 are not due to chance.

Sensitivity and specificity using cross-validation

Over-fitting is a common problem of statistical estimations in a high dimension low sample size setting. We perform leave-out experiments to evaluate the generalization errors of our model, namely its sensitivity and specificity.

We compute an atlas with the same parameter setting and initial conditions but with one control and one DS subject data out, yielding 8 2 = 64 atlases. Note that this is a design choice since one does not necessarily need to have balanced groups to apply the method. For each experiment, we register the template shape complex to each of the left-out complex by minimizing (19) for N su = 1 and considering template and control points of the atlas fixed. The resulting momentum vectors are compared with those of the atlas. We classify them based on Maximum Likelihood (ML) ratios and LDA. Let α test be the initial momenta parameterizing the deformation of the template shape complex to a given left-out shape complex (seen as a test data), and v test = K(c 0 , c 0 )α test . In this section, v, v HC and v DS denotes the sample mean using only the training data (7 HC and 7 DS). In LDA, we write the classification criterion as:

C(v test ) = (v test -v) T Σ -1 (v HC -v DS ), ( 22
)
where Σ denotes the regularized sample covariance matrix of training data (for " = 10 -2 , see Remark 3.1). For a threshold η, the test data is classified as healthy control if C(v test ) > η and DS subject otherwise. ROC curves are built when the threshold η is varied. For estimating classification scores, we estimate the threshold η on the training dataset so that the best separating hyperplane (orthogonal to the most discriminative axis Σ -1 (v HCv DS )) is positioned at equal distance to the two classes. This threshold value is used for classifying the test data.

For classifying in a Maximum Likelihood framework, we compute the regularized sample covariance matrices

Σ DS = 1 N Ds su P i2DS (v i -v DS )(v i -v DS ) T and Σ HC = 1 N HC su P i2HC (v i -v HC )(v i - v HC ) T .
The classification criterion, also called the Mahalanobis distance, is given by:

C(v test ) = (v test -v DS ) T Σ -1 DS (v test -v DS ) -(v test -v HC ) T Σ -1 HC (v test -v HC ) (23)
and the classification rule remains the same.

The very high sensitivity and specificity reported in Table 1 (first row) show that the anatomical differences between DS and controls that were captured by the model are not specific to this particular dataset, but are likely to generalize well to independent datasets.

Shape complexes versus individual shapes

In this section, we aim to emphasize the differences between using a single model for the shape complex and using different models for each individual component of a shape complex.

We perform the same experiments as described above, but for each of the three structures independently. The atlas of each structure has its own set of control points and momentum vectors. The hypothesis of equal means for DS and control subjects is rejected with a probability of false positive of p = 3.5 10 -3 for the hippocampus, p = 4.7 10 -3 for the putamen and p = 1.2 10 -4 for the amygdala. The statistical significance is lower for the hippocampus and the putamen than for the shape complex (p = 2.6 10 -4 ), and higher for the amygdala. The classification scores reported in Table 1 (rows 2 to 4) show that none of the structures alone may predict the subject's status with the same performance as the shape complex. Although the model for the amygdala has a higher statistical significance, it has a lower specificity in the Maximum Likelihood approach.

For visualization of results from individual analyses, we deform each structure along its most discriminative axis. Because the 3 deformations are not combined into a single space deformation, intersections between surfaces occur (Fig. 4). The deformation of the amygdala, though highly significant, is not compatible with the deformation of the hippocampus. From an anatomical point of view, both parts of the amygdala/hippocampus boundary should vary together, since almost nothing separates the two structures at the image resolution.

The shape complex analysis in Fig 2 and3 showed that the most discriminative effects involve deformations of specific subregions, and in particular the most lower-anterior part of the complex where the amygdala is located. Therefore, it is not surprising that this structure shows higher statistical performance than the hippocampus and putamen in an independent analysis of each structure. However, the most discriminative deformations of each structure are not consistent among themselves, thus misleading the interpretation of the findings. By contrast, the shape complex analysis shows that the discriminative effect is not specific to the amygdala but to the whole lower anterior part of the medial temporal lobe with strong correlations between parts of the structures within this region. The shape complex model may be slightly less significant, but it highlights shape effects that can be interpreted in the context of anatomical deformations related to underlying neurobiological processes.

One could argue that independently analyzing each structure does not take into account the correlations among structures. To mimic what previously reported shape analysis methods do, we build a composite shape descriptor v i by concatenating the velocities of each individual atlas v i 1 , v i 2 and v i 3 (for each structure s = 1, 2, 3 and subject i, v i s = K(c 0,s , c 0,s )α i s where α i s 's are the initial momenta in each atlas). We use this composite descriptor to compute means, sample covariance matrices, most discriminative axis and classification scores as above. This approach achieves a classification nearly as good as with the single atlas method (Table 1, last row) with a very high statistical significance p < 10 -5 . The direction of the most discriminative axis v LDA takes into account the correlations between each structures. However, this vector does not parameterize a single diffeomorphism-only each of its three components does. To display these correlations, we compute the initial momentum vectors associated with each component:

α LDA s = K(c 0,s , c 0,s ) -1 v LDA s
for s = 1, 2, 3, and then deform each structure using a different diffeomorphism. Even in this case, surfaces intersect, thus showing that this way of taking into account correlations does not prevent generating anatomical configurations that are not compatible with the data (Inline Supplementary Figure S1). By contrast, the single atlas method proposed in this work integrates topology constraints into the analysis by the use of a single deformation of the underlying space, and therefore correctly measures correlations that preserve the internal organization of the anatomical complex.

Effects of dimensionality reduction

Our approach offers the possibility to control the dimension of the shape descriptor by choosing the number of control points given as input to the algorithm. In 3D, the dimension of the shape descriptor is 3 times the number of control points. In this section, we evaluate the impact of this dimensionality for atlas construction and statistical estimations given our low sample size setting.

We start with 105 control points on a regular lattice with spacing equal to the deformation kernel width σ V and then Table 2: Decrease of the data term during optimization for different number of control points and σ V = 10 mm successively down-sample this lattice. With only 8 points, the number of deformation parameters is decreased by more than one order of magnitude and the initial ellipsoidal shapes still converge to a similar template shape complex (Fig. 1-b). The main reason for it is that control points are able to move to the most strategic places, noticeably at the tail of the hippocampus and the anterior part of the amygdala where the variability is the greatest. Qualitatively, the most discriminant axis is stable when the dimension is varied (Inline Supplementary Figure S2), as is the spectrum of the sample covariance matrices of the momentum vectors (Inline Supplementary Figure S3). The method is able to optimize the "amount" of variability captured for a given dimension of deformation parameters. Nevertheless, the residual data term at convergence increases. The initial data term (i.e., varifold norm) decreases by 97.8% for 105 points, and only by 93.3 for 8 points, thus showing that the sparsest model captured less variability in the dataset (Table 2).

If there could be an infinite number of control points, their optimal locations would be on surface meshes themselves. Therefore, one might place one control point at each vertex [START_REF] Vaillant | Surface matching via currents[END_REF][START_REF] Ma | Bayesian template estimation in computational anatomy[END_REF]. In our case, such a parameterization would involve 23058 control points. Nonetheless, this number can be arbitrarily increased or decreased by up/down sampling of the initial ellipsoids, regardless of the variability in the dataset! We increase the number of control points to 650 and notice that the estimated template shapes are the same as with 105 control points (results not shown), and that the atlas explains the same proportion of the initial data term (Table 2). Therefore, increasing the number of control points does not allow us to capture more information, which is essentially determined by the deformation kernel width σ V , but distributes this information over a larger number of parameters. This conclusion is in line with [START_REF] Durrleman | Statistical models of sets of curves and surfaces based on currents[END_REF], who show that such high dimensional parameterizations are very redundant.

The statistical significance, as measured by the p-value associated with the Hotelling T 2 statistics, is not increased with higher dimensions (Fig. 5-b). It is even smaller than in small dimensions, the maximum being reached for 16 control points (p < 10 -5 ). Leave-2-out experiments give 100% specificity and sensitivity using the ML approach, regardless of the number of control points used. To highlight differences, we performed classification using the hippocampus shape only. Again, the performance of the classifier does not necessarily decrease with the number of control points (Table 3). ROC curves in Fig. 6 show that the atlases with 48 and 18 control points have poorer performance than atlases with 12 and 8 control points.

These results suggest that using atlases of small dimension could have even greater statistical power, especially in a small sample size setting. Nevertheless, two different dimensionality reduction techniques compete with each other in these experiments. The first is the use of a small set of control points, which is a built-in dimensionality reduction technique, which has the advantage to optimize simultaneously the information captured in the data and the encoding of this information in a space of fixed dimension. The second is a post-hoc dimensionality reduction using PCA when computing classification scores that project shape descriptors into the subspace, explaining 95% of the variance captured. The variation of the p-values, when the number of modes selected in the PCA is varied, shows that a number of modes optimizes the statistical significance, between 6 and 8 modes (Inline Supplementary Figure S5). For each number of modes, an optimal number of control points also maximizes significance, and this number is never greater than 105 when one control point is placed at every σ V .

It is difficult to distinguish the effects of the two techniques in such a low sample size setting. With 8 control points and a few dozen or more subjects, we could estimate full-rank covariance matrices and would not need the post-hoc dimensionality reduction techniques. A fair comparison between post-hoc and built-in dimensionality reduction would be then possible. Our hypothesis is that, in this regime, the trend of increased statistical significance when the number of control points is decreased would be amplified. Indeed, the ratio between the number of variables to estimate and the number of subjects is more favorable in such a scenario, thus making the statistical estimations more stable.

Effects of parameter settings

We assess the robustness of the results with respect to parameter settings. We change the values of the deformation and varifold kernel widths by ±50%, namely by setting σ V = 5, 10 or 15 mm and σ W = 2.5, 5, or 7.5 mm. Other settings are kept fixed, namely the weights σ k = 10 mm, the gradient kernel width σ X = 0.5σ V and the initial distance between control points, which always equals σ V . Classification scores are reported in estimates, noticeably for the ML method. We note a decrease in the specificity in the LDA classifier for the large deformation kernel width σ V = 15 mm. With large deformation kernel widths, the atlas captures more global shape variations, which might not be as discriminative as more local changes. This effect is more pronounced with increased varifold width σ W , as surface matching accuracy decreases, thus further reducing the variability captured in the atlas. These results show that the performance of the atlas is stable for a large range of reasonable values, and therefore that they are not due to fine parameter Table 3: Classification ratios based solely on hippocampus shape. LDA and ML classification are performed with a varying number of control points in the atlas. Ratios are in percentages. Reducing the number of control points to 12 or 8 may increase statistical performance.

tuning.

The shape of the template also depends on the parameter setting, and notably the deformation kernel width σ V . With larger values, the template shape captures more rigid variations, which translates into a smoother shape. With smaller values, the template captures finer details in the data (Inline Supplementary Figure S4)

The dimension of the atlas is intrinsically linked with the deformation kernel width. Deformations with smaller σ V need more control points to potentially deform every small region of the shape complex. Deformations with larger σ V have fewer degrees of freedom and could be decomposed using fewer control points. Placing one control point at the nodes of lattice of step σ V yields 15 control points for σ V = 15 mm, 105 control points for σ V = 10 mm and 650 control points for σ V = 5 mm. We build an atlas for each of these values of σ V and with down/up sampling the set of associated control points. All these atlases show a good significance level, far below the usual 0.05 threshold. On average, the statistical significance is decreased with increasing σ V , as the atlas represents a coarser and coarser description of the variability within the dataset (Fig. 5). With σ V = 15 mm (Fig. 5-c), the maximum significance is reached for 8 control points, and the significance is decreased with increasing dimensionality. With σ V = 5 mm (Fig. 5-a), the same trend is observed, except an unexpected increase in statistical significance at very high dimensions. These results show that the discussion about dimensionality reduction in the previous section does not depend on a particular choice of deformation kernel width.

We also assess the influence of the amount of regularization in the covariance matrices ", which otherwise are singular. We increase the value used in the previous experiment from " = 10 -2 to " = 0.1, " = 1 and " = 10. With these values, the condition number of the covariance matrix decreased from 1000 to 100, 10 and 1 respectively. A decrease in the sensitivity of the classifier was detected only for " = 10, that is when the regularization became of the same order as the largest eigenvalues of the matrix. The choice of this setting has, therefore, very little influence on the classification results.

It is clear that the weights σ k 's also should have been adjusted. As noted in [START_REF] Akin | You laid out the lands:" georeferencing the Chinese Yujitu [Map of the Tracks of Yu] of 1136[END_REF], adjusting the weights could increase matching accuracy, and possibly increase statistical performance. As explained in Sec. 2.6.2, these values were chosen so that the data term has the same order of magnitude as the sum of squared geodesic distances. However, it is clear from a statistical point of view that these values mea- 4: Classification scores when deformation and varifold kernel widths are varied. Regularization of the covariance matrices " = 10 -2 . Results are overall very stable when settings are varied. Very large kernel widths penalize the matching accuracy between the template and the subject shape complexes, thus eventually altering classification performance. sure noise variance, and therefore should be estimated from the data and not fixed by the user. This estimation could be done in a Bayesian framework by adapting to varifolds the method proposed in [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF][START_REF] Allassonnière | Construction of bayesian deformable models via a stochastic approximation algorithm: A convergence study[END_REF] for images.

Overall, these experiments demonstrate the reproducibility of our results under various parameter settings. They show that the method could be applied in real cases without fine parameter tuning.

Discussion and Conclusion

This paper presents a comprehensive framework for the statistical analysis of shape complexes extracted from 3D anatomical images. The method can deal with raw surfaces resulting from nearly any segmentation methods thanks to its robustness to noise, mesh imperfections and inconsistencies in mesh orientation. The scheme estimates a template shape complex with a fixed topology that is representative of the anatomy, and computes modes of deformation that preserve template structure and capture variability in data. Such topology constraints lead to modes of variations that are anatomically realistic and interpretable. The proposed approach therefore contrasts with the study of correlations between shape models that are estimated independently for each component within a shape complex. Given a typical neuroimaging study of a complex of deep brain structures in Down syndrome subjects, the method can find discriminative anatomical features with high statistical significance and small generalization errors, even with a limited number of observations. We show the robustness of these results in various experimental settings, demonstrating the effectiveness of the method without fine parameter tuning. The scientific community can evaluate the method by downloading the software Deformetrica, which is freely available at www.deformetrica.org.

The statistical analysis on deformations that we proposed is essentially multivariate. Statistics show the correlations between the deformation patterns in every region of the brain. The visualization of the deformations gives a comprehensive view of how these local deformations are combined into a consistent deformation of the underlying tissue. This analysis is therefore in strong contrast with voxel-based methods, which test at every voxel the difference in image intensities [START_REF] Ashburner | Voxel-based morphometry-the methods[END_REF] or the difference in the Jacobian determinant of the template-to-subject deformations [START_REF] Thompson | Growth patterns in the developing human brain detected by using continuum-mechanical tensor maps[END_REF]. In particular, the analysis of the Jacobian determinant only indicates local contraction or expansion of the tissue, while ignoring more complex deformations patterns such as torques or a shift between two structures. Such cofounding effects may be misleading when interpreting the results.

In contrast to such mass-univariate methods, our multivariate approach also avoids the problem of correction for multiple comparisons. The dimension of the variables used in the statistical analysis is essentially determined by the deformation kernel width σ V and therefore by the scale of anatomical variants that are captured by the model. In the current scheme, the choice of the number of control points is left to the user, using a practical heuristics that consists in placing one point for every deformation kernel width σ V . We show that this number could be even drastically reduced without altering statistical significance and generalization ability of the model. This built-in dimensionality reduction may lead to increased statistical performance as suggested by our results, although our initial results need to be confirmed and supported using more subjects and different datasets. The fact that the dimension is determined by the user before any experiments allows one to adjust the scale σ V according to the number of available subjects, and also eases the power calculations and sample sizes estimates required in clinical trials. This finite-dimensional setting also paves the way for estimating mean and covariance matrices during the optimization in a Bayesian framework, following research by [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF] and [START_REF] Allassonnière | Stochastic algorithm for bayesian mixture effect template estimation[END_REF]. Constraining statistical inference to take place in a small dimensional space is likely to increase the convergence speed of the statistical estimates, as compared to performing the inference in very high dimensions and then performing post-hoc dimensionality reduction, using PCA for instance.

Cross-validation showed the very good prediction capability of our model. The prediction of Down syndrome based on neuroimaging data has little clinical interest, since subjects are characterized by their genotype and especially the copy number of chromosome 21, which is known with very high confidence. However, the shape deformation studies as shown here may give new insights into anatomical changes linked to genetics, and associations between morphologic differences and cognitive and behavioral scores. Nonetheless, our model is completely generic and can be applied to different pathologies for which the clinical status may be more difficult to assess. This prediction capability of the method demonstrates its potential in computer-aided diagnosis or prognosis in studies where a subject's status is based only on clinical diagnosis with limited reproducibility, such as in neurodegenerative diseases, or for pre-diagnostic prediction of disease onset based on image data. Shape descriptors, which encode the joint shape variability of sets of anatomical structures with a small number of parameters, would be preferable to study correlations between anatomical phenotypes and genotype, in the spirit of [START_REF] Korbel | The genetic architecture of down syndrome phenotypes revealed by high-resolution analysis of human segmental trisomies[END_REF], where these image-derived parameters can take the place of clinical variables. K(c, c)α is a 3N cp vector, whose kth coordinate is the 3D vector:

P N x p=1 K(c k , c p )↵ p . Therefore, d c i (K(c, c)α) k = N cp X p=1 ↵ p r 1 K(c k , c p ) T δ(i -k) + ↵ i r 2 K(c k , c i ) T (30) Using the fact that K is symmetric (hence r 1 K(x, y) = r 2 K(y, x)) we have: ↵i = - 1 2 N cp X k=1 (d c i (K(c, c)α) k ) T ↵ k = - 0 B B B B B B @ N cp X k=1 r 1 K(c i , c k )↵ T k 1 C C C C C C A ↵ i (31)
Appendix B Gradient of the atlas criterion

We provide here the differentiation of the criterion for atlas construction:

E(X 0 , c 0 , {α i 0 }) = N su X i=1 ⇣ A(X i (1)) + L(S i 0 ) ⌘ (32)
subject to:

8 > > < > > : Ṡi (t) = F(S i (t)) S i (0) = {c 0 , α i 0 } Ẋi (t) = G(X i (t), S i (t)) X i (0) = X 0 (33)
where

L(S i 0 ) = α i 0 T K(c 0 , c 0 )α i 0 (34)
X is a vector of length 3N x , where N x is the number of points in the template shape, c and α are two vectors of length 3N cp each, where N cp is the number of control points, so that S is a vector of length 6N cp .

F(S) = F c (c, α) F ↵ (c, α)
! is a vector of length 6N cp , which is decomposed into two vectors of size 3N cp . The kth coordinate (among N cp ) of F c and F ↵ is the 3D vector:

F c (S) k = N cp X p=1 K(c k (t), c p (t))↵ p (t) F ↵ (S) k = - N cp X p=1 ↵ k (t) T ↵ p (t)r 1 K(c k (t), c p (t)) (35)
G(X, S) is a vector of size 3N x . Its kth coordinate (among N x ) is the 3D vector:

G(X, S) k = N cp X p=1 K(x k (t), c p (t))↵ p (t) (36) Similarly, L(S i 0 ) = N cp X p=1 N cp X q=1 ↵ i 0,p T K(c 0,p , c 0,q )↵ i 0,q (37)

B.1 Gradient in matrix form

The differentiation of the criterion can be done for each subject i independently. Therefore, we differentiate only one term of the sum in (32) for a generic subject's index i that we omit in the rest of this section for clarity purposes.

A small perturbation δS 0 induces a perturbation of the motion of the control points and momenta δS(t), which, in turn, induces a perturbation of the template points' trajectory δX(t) and then of the criterion δE, which we write, thanks to the chain rule δE = r X(1) A T δX(1) + r S 0 L T δS 0 .

(38)

According to (33), the perturbations δS(t) and δX(t) satisfy the linearized ODEs:

δS(t) = d S (t) FδS(t) δS(0) = δS 0 δ X(t) = @ 1 GδX(t) + @ 2 GδS(t) δX(0) = δX 0
The first ODE is linear. Its solution is given by:

δS(t) = exp Z t 0 d S(u) Fdu ! δS 0 . (39)
The second ODE is linear with source term. Its solution is given by:

δX(t) = Z t 0 exp Z t u @ 1 Gds ! @ 2 G(u)δS(u)du + exp Z t 0 @ 1 G(s)ds ! δX 0 (40)
Plugging ( 39) into ( 40) and then into (38) leads to:

8 > > > > < > > > > : r S 0 E = Z 1 0 ⇣ R 0t T @ 2 G(X(t), S(t)) T V t1 T r X(1) A ⌘ dt + r S 0 L r X 0 E = V 01 T r X(1) A , ( 41
) where we denoted R st = exp ⇣R t s d S(u) Fdu ⌘ and V st = exp ⇣R t s @ 1 G(X(u), S (u))du ⌘ . Let us denote θ(t) = V t1 T r X(1) A, g(t) = @ 2 G(t) T θ(t) and ξ(t) = R 1 t R ts T g(s)
ds, so that the gradient (41) can be rewritten as:

8 > > > > < > > > > : r S 0 E = Z 1 0 R 0s T g(s)ds + r S 0 L = ξ(0) + r S 0 L r X 0 E = θ(0)
. Now, we need to make explicit the computation of the auxiliary variables θ(t) and ξ(t). By definition of V t1 , we have V 11 = Id and dV t1 /dt = V t1 @ 1 G(t), which implies that θ(1) = r X(1) A and θ(t) = -@ 1 G(t) T θ(t).

For ξ(t), we notice that R ts = Id -R s t dR us du du = Id + R s t R us d S(u) F(u)du. Therefore, using Fubini's theorem, we get:

ξ(t) = Z 1 t R ts T g(s)ds = Z 1 t g(s) + d S(s) F T Z 1 s R su T g(u)du ! ds = Z 1 t ⇣ g(s) + d S(s) F T ξ(s) ⌘ ds.
This last equation is nothing but the integral form of the ODE given in the main text.

Given the actual values of S 0 and X 0 , one needs to integrate the geodesic shooting equations and the flow equation in (33) forward in time to give the full path of parameters S(t) and template shape points X(t). Then, one needs to compute the gradient of the data term r X(1) A, which is given in Appendix C. This term indicates in which direction one has to move the vertices of the deformed template shape in order to better match the observations. This term is transported back to time t = 0 by the coupled linear equations satisfied by ξ and θ. The values of time t = 0 of these auxiliary variables are used to update the deformation parameters (position of control points and momenta) and the position of the vertices of the template surfaces.

B.2 Gradient in coordinates

Expanding the variables S

i (t) = {c 0,k (t), ↵ i 0,k (t)}, X i (t) = {X i k (t)}, θ i (t) = {✓ i k (t)} and ξ i (t) = {ξ c,i k (t), ξ α,i k (t)}, we have r c 0,k E = N su X i=1 ξ c,i k (0) + r c 0,k L(S i 0 ) r α i 0,k E = N su X i=1 ξ α,i k (0) + r α i k L(S i 0 ) r x 0,p E = N su X i=1 θ i p (0)
where the gradient of L is given as (from now on, we omit the subject's index i for clarity purposes): The term ∂ 1 G(X(t), S(t)) is a block-matrix of size 3N cp ⇥ 3N x whose (k, p)th 3 ⇥ 3 block is given as:

r α 0,k L = 2
d X k G(X(t), S(t)) p = N cp X j=1
α j (t)r 1 K(X p (t), c j (t)) T δ(pk) so that the vector θ(t) is updated according to:

-θk (t) = N cp X p=1 α p (t) T θ k (t)r 1 K(X k (t), c p (t)) (42)
The terms ∂ c g G(X(t), S(t)) and ∂ α G(X(t), S(t)) are both matrices of size 3N x ⇥ 3N cp , whose (k, p) block is given, respectively, by:

d c k G p = α k ⇣ r 1 K(c k , X p ) ⌘ T d α k G p = K(c k , X p )I 3
The differential of the function F(S) = F c (c, α) F α (c, α) ! can be decomposed into 4 blocks as follows:

d S(t) F = ∂ c F c ∂ α F c ∂ c F α ∂ α F α ! (43)
Therefore, the update rules for the auxiliary variables ξ c (t) and ξ α (t) are given as:

8 > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > : -ξc k (t) = N x X p=1
α k (t) T θ p (t)r 1 K(c k (t), X p (t))

+ (∂ c F c ) T ξ c (t) k + (∂ c F α ) T ξ α (t) k -ξα k (t) = N x X p=1
K(c k (t), X p (t))θ p (t)

+ (∂ α F c ) T ξ c (t) k + (∂ α F α ) T ξ α (t) k with (∂ c F c ) T ξ c (t) k = N cp X p=1 ⇣ α p (t) T ξ c k (t) + α k (t) T ξ c p (t) ⌘ r 1 K(c k (t), c p (t)) (∂ c F α ) T ξ α (t) k = N cp X p=1 α k (t) T α p (t)r 1,1 K(c k (t), c p (t)) T ⇣ ξ α p (t) -ξ α k (t) ⌘ (∂ α F c ) T ξ c (t) k = N cp X p=1 K(c k (t), c p (t))ξ c j (t) (∂ α F α ) T ξ α (t) k = N cp X p=1 r 1 K(c k (t), c p (t)) T ⇣ ξ α p (t) -ξ α k (t) ⌘ α p (t)
In these equations, we supposed the kernel symmetric: K(x, y) = K(y, x). If the kernel is a scalar isotropic kernel of the form K = f (kxyk 2 )I 3 , then we have: r 1 K(x, y) = 2 f 0 (kxyk 2 )(xy) r 1,1 K(x, y) = 4 f 00 (kxyk 2 )(xy)(xy) T + 2 f 0 (kxyk 2 )I 3

Appendix C Gradient of the varifold metric for meshes

We derive here the gradient of the varifold metric with respect to the position of the vertex of the mesh. Let S be a triangular mesh. For each face f k , we denote n k its normal, p k its center and u k = n k / |n k | 1/2 . Let T be another triangular mesh, m k its normal, q k its center and v k = m k / |m k | 1/2 . Our goal is to compute the gradient of d(S, T ) 2 with respect to x i , a given vertex of S. The chain rule gives: where we sum over all the faces that have x i among their vertices.

r x i d(S, T ) 2 = X f k 3x i (d x i n k ) T (d n k u k ) T
Given the inner-product between varifolds (see main text), we have:

r u k d(S, T ) 2 = 4 0 B B B B B B @ N S X i=1 K W (p k , p i )u i u T i - N T X j=1 K W (p k , q j )v j v T j 1 C C C C C C A u k , (
N S X i=1 @K W (p k , p i ) @ p k,d u i u T i - N T X j=1 @K W (p k , q j ) @ p k,d v j v T j 1 C C C C C C A u k (46)
Finally, for a face f k , we have n k = 1 2 (X 1 -X 0 ) ⇥ (X 2 -X 0 ) and p k = 1 3 (X 0 + X 1 + X 2 ), where we denote X 0 , X 1 , and X 2 the vertices of the face. If we denote e the edge opposite to the vertex x i (i.e., e = X 2 -X 1 if x i = X 0 ), we have for a generic 3D-vector V:

(d x i n k ) T V = 1 2 e ⇥ V and (d x i p k ) T V = 1 3 V. ( 47
)
and since

u k = n k / |n k | 1/2 , d n k u k = 1 |n k | 1/2 0 B B B B @ I 3 - 1 2 n k n T k |n k | 2 1 C C C C A = 1 |u k | 0 B B B B @ I 3 - 1 2 u k u T k |u k | 2 1 C C C C A (48)
The gradient is computed by plugging ( 45), ( 46), ( 47) and ( 48) into (44). The gradient is computed by scanning each face of the mesh S and adding the contribution of this face to each of its vertices.

One can easily verify that (44) is independent of the ordering of the vertices, thus showing its invariance with respect to the local orientation of the mesh. dashed) lines correspond to the 10% (resp 5%) significance levels, respectively. For a given number of modes, the best p-value is never achieved for the largest number of control points, showing the interest of small-dimensional models. It seems also that there is an optimal number of modes to be selected, for which the statistical power is overall increased (between 6 and 8 modes). With a few subjects more, we could estimate a full-rank covariance matrix and make the method less and less sensitive to the number of modes selected. We hypothesize that the effect of the number of control points will be more pronounced in this regime (Note that the Fig. 5-b is built from these plots: for each number of control points, we picked the p-values that correspond to the number of modes explaining 95% of the variance, which was always either 8 or 9.)

Figure 1 :

 1 Figure 1: Atlas estimated from different initial conditions. Left: 105 control points with initial spacing equal to the deformation kernel width σ V = 10 mm, Right: 8 control points. Arrows are the momentum vectors of DS subjects (red) and controls (blue).Control points that were initially on a regular lattice move to the most variable place of the shape complex during optimization. Arrows parameterize space deformations and are used as a shape descriptor of each subject in the statistical analysis.

Figure 2 :

 2 Figure2: Template complex deformed using the mean deformation of controls (transparent shapes) and DS subjects (opaque shapes), which illustrates the anatomical differences that were found between both groups.

Figure 3 :

 3 Figure 3: Most discriminative deformation axis showing the anatomical features that are the most specific to the DS subjects as compared to the controls. Differences are amplified, since the distance between the two configurations is twice the distance between the means (black grids are mapped to the surface for visualization only)

Figure 4 :

 4 Figure 4: Most discriminative deformation axis computed for each structure independently. Surface intersection occurs in the absence of a global diffeomorphic constraint. (black grids are mapped to the surface for visualization only)

Figure 5 :

 5 Figure5: Statistical significance of the group means difference for a varying number of control points. The solid (resp. dashed) lines correspond to the 0.1 (resp. 0.05) significance thresholds, respectively. The ability of the classifier to separate DS subjects to controls is little altered by the deformation kernel width σ V . Increasing the number of control points, and hence the dimensionality of the atlas, does not necessarily increase statistical performance.

Figure 6 :

 6 Figure 6: ROC curves for hippocampus classification using a different number of control points in the atlas and ML classifier. Atlases with 48 and 18 control points exhibit poorer performance than those with 12 and 8 control points.

  0,k , c 0,p )α 0,p r c 0,k L = 2 N cp X p=1 α 0,p T α 0,k r 1 K(c 0,k , c 0,p )

  r u k d(S, T ) 2 + (d x i p k ) T r p k d(S, T ) 2 ,(44)

  45) and denoting p k,d the dth coordinate of the 3D vector p k ,

  FigureS5: P-values computed for a different number of control points and a different number of selected modes. Solid (resp. dashed) lines correspond to the 10% (resp 5%) significance levels, respectively. For a given number of modes, the best p-value is never achieved for the largest number of control points, showing the interest of small-dimensional models. It seems also that there is an optimal number of modes to be selected, for which the statistical power is overall increased (between 6 and 8 modes). With a few subjects more, we could estimate a full-rank covariance matrix and make the method less and less sensitive to the number of modes selected. We hypothesize that the effect of the number of control points will be more pronounced in this regime (Note that the Fig.5-b is built from these plots: for each number of control points, we picked the p-values that correspond to the number of modes explaining 95% of the variance, which was always either 8 or 9.)

False Positive Rate (in %) True Positive Rate (in %)

  Table 4 and show a great robustness of the statistical

	100
	97.5
	95
	92.5
	90
	87.5
	85
							48 Control Points
	82.5						18 Control Points
							12 Control Points
							8 Control Points
	80						4 Control Points
	0	2.5	5	7.5	10	12.5	15	17.5	20

Table

Acknowledgments. We thank Christine Pickett for her careful proofreading of the manuscript. This work has been partly funded by the program "Investissements d'avenir" ANR-10-IAIHU-06 and by the NIH grants U54 EB005149 (NA-MIC), 1R01 HD067731, 5R01 EB007688 and 2P41 RR0112553-12.

Appendix A Geodesic equations

We derive here the minimum action principle of Lagrangian mechanics. A variation δα(t) of the time-varying momentum vectors α(t) induces a variation of the control point positions δc(t), which in turn induces a variation δE of the quantity E = R 1 0 α(t) T K(c(t), c(t))α(t)dt.

Since ċ = K(c, c)α, we have

and

Therefore, we have

and

Assuming δc(0) = δc(1) = 0, integration by parts yields:

The linear ODE with source term (24) shows that there is a one-to-one relationship between δc and δα. Since δα is arbitrary, so is δc and

along extremal paths.

Appendix D Diffeomorphic template evolution

The purpose of this section is to prove that no selfintersection may occur during the optimization of the template shape, by showing that the updates of the template follow a geodesic flow of diffeomorphisms. Using notations of the main text, rE x 0,p is the gradient of the criterion with respect to the position of the vertex x 0,p of the current template using the L 2 metric, and r X E x 0,p its smoothed version using a metric given by a Gaussian kernel with width σ X > 0, K X , so that:

where u s is a vector field in V X , the RKHS associated with the Gaussian kernel K X . In particular, if s is the flow associated with integration of u s , we get X 0 (s) = s (X 0 (0)). An important point to be verified here is that this flow exists and generates a continuous curve s ! s of C 1 diffeomorphisms so that the template components cannot degenerate or self-intersect.

Let Ω X be the open set of the configurations X 0 such that all the mesh faces associated with X 0 are non-degenerated (positive area) and that any pairs of distinct vertices do not coincide in space. The total energy E(X 0 , {S i 0 }) is C 1 on an open set Ω X ⇥ R N S so that the local existence of the gradient descent follows from the Cauchy-Lipschitz theorem. Now, if we consider a maximal solution on [0, s f [, we will prove below (and this is the key estimate) that

so that the flow s is a flow of C 1 diffeomorphisms staying at a bounded distance d X (Id, s )  p E 0 from the identity and X 0 (s) = s (X 0 (0)) stays in a compact set of Ω X . In particular, since the differential d and d -1 can be controlled uniformly by d X (Id, ), we get that no face can degenerate during the gradient descent, that the distance between two distinct vertices or two surface patches (up to the continuous limit) cannot vanish. Now, we prove (49). From the RKHS property of the kernel we get

Supplementary Data

Right … towards … controls

DS subjects

Figure S1: Most Discriminative Axis computed using a composite descriptor. The direction takes into account the correlations among the three structures. However, it does not parameterize a single space deformation, but three of them, and intersections between surfaces occur. Moreover the patterns of shape variations are rather different from the results using a single atlas of shape complex, in particular the relative position of the amygdala (in blue) with respect to the hippocampus (in green)

Anterior Right … towards … controls DS subjects