B Charlier

N Charon

A Trouvé

The fshape framework for the variability analysis of functional shapes

This article introduces a full mathematical and numerical framework for treating functional shapes (or fshapes) following the landmarks of shape spaces and shape analysis. Functional shapes can be described as signal functions supported on varying geometrical supports. Analysing variability of fshapes' ensembles require the modelling and quantification of joint variations in geometry and signal, which have been treated separately in previous approaches. Instead, building on the ideas of shape spaces for purely geometrical objects, we propose the extended concept of fshape bundles and define Riemannian metrics for fshape metamorphoses to model geometrico-functional transformations within these bundles. We also generalize previous works on data attachment terms based on the notion of varifolds and demonstrate the utility of these distances. Based on these, we propose variational formulations of the atlas estimation problem on populations of fshapes and prove existence of solutions for the different models. The second part of the article examines thoroughly the numerical implementation of the tangential simplified metamorphosis model by detailing discrete expressions for the metrics and gradients and proposing an optimization scheme for the atlas estimation problem. We present a few results of the methodology on a synthetic dataset as well as on a population of retinal membranes with thickness maps.

group actions of diffeomorphisms [START_REF] Trouvé | Diffeomorphisms groups and pattern matching in image analysis[END_REF][START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF][START_REF] Miller | On the metrics and euler-lagrange equations of computational anatomy[END_REF][START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Miller | Geodesic Shooting for Computational Anatomy[END_REF][START_REF] Trouvé | Handbook of Mathematical Imaging, chapter Shape spaces[END_REF][START_REF] Younes | Shapes and diffeomorphisms[END_REF][START_REF] Michor | A zoo of diffeomorphism groups on R n[END_REF] that should be rooted to the seminal work of V. Arnold [START_REF] Arnold | Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits[END_REF] are powerful vehicles to build a full theoretical as well as computational framework for a comprehensive quantitative analysis of shape variability recently coined as diffeomorphometry and successfully applied in computational anatomy [START_REF] Miller | Diffeomorphometry and geodesic positioning systems for human anatomy[END_REF].

The main purpose of the paper is to develop a framework embedding the situation of geometrical shapes carrying functional information that we call here fshapes for functional shapes. The idea of functional shape is quite natural since in many scientific settings, a geometrical shape X (basically a submanifold) has associated with it a scalar field f : X → R attached to every geometrical point of the support X. What is much more unusual however is to consider a fshape i.e. the pair (X, f ) as a single object that should live in some well defined and well structured ensemble F of fshapes from which further processing can be derived. We believe this approach of putting the effort on a well grounded definition of shape spaces that play a similar role to functional spaces in modern analysis should be pushed forward to the setting of fshapes. This global point of view leading to the notion of geometry of shape spaces is now well established with fruitful development in many regards [START_REF] Micheli | Sobolev metrics on diffeomorphism groups and the derived geometry of spaces of submanifolds[END_REF] for purely geometrical shapes but as far as we know the extension of such approach to functional shapes has only started very recently in [START_REF] Charon | Analysis of geometric and functional shapes with extensions of currents. Application to registration and atlas estimation[END_REF].

A core issue in working with fshapes is their mixed geometrical and functional nature so that smooth infinitesimal transformations of a given fshape (X, f ) should combine a geometrical and functional component simultaneously. The geometrical part should transport the supporting manifold X, and the functional one should modify the functional signal f . This notion of combined geometrical and functional infinitesimal transformation has been introduced previously in [START_REF] Trouvé | An approach of pattern recognition through infinite dimensional group action[END_REF] in the setting of images on a fixed support X (where X in the unit square or a flat torus). In that case the infinitesimal signal evolution δf is the combination of two factors: on the one hand, the variation of the signal due to the geometric transport of pixel values, on the other hand a purely additive perturbation of the signal. This has been further studied in [START_REF] Trouvé | Local geometry of deformable templates[END_REF] and conceptualize as the metamorphosis framework in [START_REF] Trouvé | Metamorphoses through lie group action[END_REF]. The situation we are looking at here is more general since now the support is a submanifold that can freely evolve during the metamorphosis (as in the example of Figure 1). A first global outcome is that the associated space of fshapes is not only an orbit G.X 0 under the action of a subgroup of smooth diffeomorphisms G of the ambient space on a geometrical template X 0 : the orbit G.X 0 is now the base space of a vector bundle and above each manifold X in the orbit we consider the full vector space L 2 (X) of square integrable signals so that we end up with a vector bundle F over the orbit (see Figure 2). The precise construction of the fshape bundle F is described in Section 2 as well as the associated metamorphosis based Riemannian metric for which we can prove two key results: the existence of geodesic between any two fshapes (Theorem 2) and the existence of Karcher means for a population of such fshapes in F (Theorem 3).

Obviously there is no hope to build a fshape bundle F that can contain any possible pair (X, f ) in particular since any two fshapes within F will have diffeomorphic supports. Every fshape in F should be understood as an ideal model for truly noisy observed fshapes that generically do not belong to F. As a consequence, a second core issue is to build proper smooth data attachment terms or dissimilarity measures that can be defined between any arbitrary pair of fshapes with possibly non diffeomorphic supports. The powerful setting of mathematical currents equipped with dual norms for purely geometrical shapes [START_REF] Glaunès | Surface matching via currents[END_REF] has been successfully extended to the situation of fshapes in [START_REF] Charon | Functional currents : a new mathematical tool to model and analyse functional shapes[END_REF]. However, we develop here (Section 4) a different approach encoding non oriented tangential space information that is based on the concept of varifold along the lines of [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF] and extended here to the new situation of functional varifolds. We believe that this new type of dissimilarity measure, which is able to compare fshapes on smooth manifolds as well as on polyhedral meshes, is well suited in combination with metamorphoses distances to define a theoretical and computational framework for fshapes. In particular, we establish that useful functional varifolds metrics can be induced by smooth embeddings of functional varifolds into Reproducible Kernel Hilbert Spaces (RKHS) associated with computable kernels. We then establish that the resulting metric has several key regularity properties with respect to smooth geometrical and functional variations of fshapes (Theorem 5). Such regularity results open the way to various smooth and computationally tractable variational problems on fshapes involving observed fshapes. We concentrate on the central problem of template (or atlas in the context of computational anatomy) estimation which is the extension of the Karcher mean problem to the much broader situation where the observed fshapes do not belong to any a priori defined fshape bundle F. A much more studied situation is the case of purely image data that has been pushed forward on practical ground in the framework of actions of diffeomorphisms in [START_REF] Joshi | Unbiased diffeomorphic atlas construction for computational anatomy[END_REF]. Here the goal is to obtain, via a coherent fshape framework, a mean template (X, f ) encoding simultaneously the mean geometrical and functional information from a dataset (X i , f i ) of fshapes as the minimizer of a well posed variational problem. This variational problem is basically the sum of the square of the metamorphosis distances between the template (X, f ) and the approximating models ( Xi , f i ) of the observations (X i , f i ) within the fshape bundle F (associated to (X, f )) with the addition of the square of the fvarifold discrepancy measure for the difference between the models ( Xi , f i ) and the true observation (X i , f i ). We prove that introducing a hypertemplate in the spirit of [START_REF] Ma | Bayesian template estimation in computational anatomy[END_REF][START_REF] Ma | A bayesian generative model for surface template estimation[END_REF] and restricting the template to belong to the fshape bundle F 0 of a hypertemplate (X 0 , f 0 ) is sufficient to establish a well posed variational problem (Section 4). The proof of existence of an optimal template (X, f ) as the minimizer within the fshape bundle F 0 is established on fairly general assumptions that cover most of the practical situations (Theorem 6). Note that even restricted to the situation of purely geometrical shapes, as far as we know, there is no available results in the literature providing a rigorous proof of existence of an average shape from a population of noisy geometrical shapes. Available rigorous results have usually focused on the situation where the dataset is sampled from a Riemannian manifold (usually finite dimensional). When the functional part is also involved, the situation is more complex due to the interactions between function and geometry inside the Riemannian metamorphosis distance but also the fvarifold metric defining the dissimilarity term. We establish similar existence results in a slightly simpler situation, called the tangential model for the metamorphosis metric, where the infinitesimal signal variations along a path are computed with respect of a frozen L 2 (X) metric on the geometrical support X of the initial fshape.

Since the paper aims at presenting both a theoretical and a computational framework, we extensively explore the material needed to bridge the gap between the theoretical framework and its derivation into a computational and algorithmic one in the latter part of this paper. We believe that this part will give also more practical insights into what goals can be achieved with the fshape framework. Moreover, we provide the details necessary to conform to the minimal standard of reproducible research. Just as the absence of a well posed theoretical framework may often produce inconsistent algorithms, we believe that the lack of precise descriptions of the underlying numerical schemes and hidden "tricks" may also hinder the development of incrementally better solutions and fair comparison by other researchers. We start in Section 6 with the derivation of a discrete framework based on polyhedral representation of fshapes leading to particle based approximations of the fvarifold dissimilarity metric and of the Hamiltonian underlying the optimal control formulation of the geodesic trajectories. The actual implemented framework is based on the tangential model for the metamorphosis metric that leads to a slightly more straightforward variational problem. The optimization scheme is carefully described in Section 7 and numerical illustrations are provided in Section 8 on real and synthetic datasets. A special attention is given in Section 9 on various numerical issues and pitfalls that may affect the computation process. In particular we discuss the very important problem of open surfaces that introduce free boundaries evolutions during the optimisation process. The tuning of free parameters are also discussed.

As discussed in this introduction, this paper describes on the one hand a new theoretical shape space framework to work with functional shapes with an emphasis on precise mathematical statements and proofs. The more technical aspects of the novel framework have been collected and presented in Appendix to ease the reading of the paper. The reader more interested in applying the framework could read Section 2 together with the first two Subsections of Section 3, skip Section 4 and 5 and jump over directly to Section 6 to 9 dealing with the numerical part of the paper.

Riemannian metamorphosis framework for fshapes

In the classical Grenander's setting, shape spaces are modelled as sets of shapes homogeneous under the action of a group of space transformations. Metrics between shapes are then induced from right-invariant Riemannian metrics on the group itself. The goal of this section is to propose a similar but extended framework for the case of functional shapes. We show, in the first place, that sets of functional shapes can be structured naturally as vector bundles and we then define a Riemannian setting to model and quantify transformations within those spaces.

Fshape bundles

Let's introduce a finite dimensional vector space E that shall be the embedding space of all shapes. In the large deformations' model (LDDMM, cf [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF]), diffeomorphisms are constructed as flows of time-varying vector fields. The basic ingredient is a reproducing kernel Hilbert space (RKHS) of vector fields on E that is denoted V and which is continuously embedded in C 1 0 (E, E) the space of continuously differentiable functions from E to E vanishing at the infinity. Let G = G V be the associated group of diffeomorphisms obtained by flowing vector fields in L 2 ([0, 1], V ) (see [START_REF] Trouvé | Handbook of Mathematical Imaging, chapter Shape spaces[END_REF][START_REF] Younes | Shapes and diffeomorphisms[END_REF][START_REF] Arguillere | Shape deformation analysis from the optimal control viewpoint[END_REF]). Let X 0 be a homogeneous space generated by a finite volume d-dimensional rectifiable compact subset X 0 (see [START_REF] Morgan | Geometric measure theory, a beginner's guide[END_REF]) i.e. X 0 = G.X 0 = {φ(X 0 ) | φ ∈ G}. Now we consider the space

F . = { (X, f ) | X ∈ X 0 and f ∈ L 2 (X) } , (1)
where L 2 (X) is the space of square integrable functions on X, i.e the set of functions f :

X → R such that X f 2 (x)dH d (x) < +∞
for H d the d-dimensional Hausdorff (or volume) measure. F can be considered as a vector bundle with fiber L 2 (X 0 ) (here we will not try to define any explicit differentiable structure on it so that the last statement should be considered as formal). Each element of F is thus given as a couple of a geometrical shape X and a signal function living on X, that we call a functional shape (see [START_REF] Charon | Functional currents : a new mathematical tool to model and analyse functional shapes[END_REF]) or fshape.

Metamorphosis distance on fshapes bundles

For a purely geometrical deformation φ ∈ G V , a natural transport of a fshape (X, f ) would be given by φ.(X, f ) = (φ(X), f • φ -1 ), which corresponds to deforming the support X and transporting the values of the signal f onto the deformed shape, which is the exact generalization of image deformation. Such a model, originally presented in [START_REF] Charon | Functional currents : a new mathematical tool to model and analyse functional shapes[END_REF], remains insufficient to account for variations of signals within fibers themselves. This justifies the following efforts to propose a model of joint geometric and functional transformations. Let's consider (X, f ) ∈ F and instantaneous velocities

X v t .X t X X 1 = φ v 1 (X) L 2 (X) f h t f + ζ h 1 L 2 (X 1 )
f 1 = (f + ζ h 1 ) • (φ v 1 ) -1
(v, h) ∈ L 2 ([0, 1], V × L 2 (X)). We can define from (v, h) a path (φ v , ζ) t∈[0,1]
where φ v t is the usual flow of v starting from the identity and t → ζ t = t 0 h s ds is the path in L 2 (X) with instantaneous speed given by h. From that path, we get a path t → (X t , f t ) in F defined by

(X t , f t ) . = (φ v t .X, (f + ζ h t ) • (φ v t ) -1 ) . (2)
where

φ v t .X . = φ v t (X)
is the natural diffeomorphisms transport action on rectifiable subsets of E. We denote (φ v 1 , ζ h 1 ) the end point value of the path. Now we denote for γ V , γ f > 0

E X (v, h) . = γ V 2 1 0 |v t | 2 V + γ f 2 1 0 X |h t | 2 (x)|d x φ v t ↾ TxX |dH d (x) (3)
where T x X is the tangent space to X at point x (defined H d -almost everywhere on X if X is rectifiable, cf [13] 3.2.19) and, for any vector space U , |d x φ v t ↾ U | denotes the Jacobian of φ v t at x restricted to U , which equals

|d x φ v t (u 1 ) ∧ ... ∧ d x φ v t (u d )| = (det( d x φ v t (u i ), d x φ v t (u j ) ) i,j ) 1/2 if (u 1 , ..., u d ) is an orthonormal basis of U . It is clear that E X (v, h) < ∞ since by definition 1 0 |v t | 2
V dt < ∞ and since the usual controls on φ t along finite energy paths give a uniform control in (x, t) of d x φ t and the existence of an increasing function C : R + → R + (independent of X) such that

E X (v, h) ≤ γ V 2 1 0 |v t | 2 V dt + C( 1 0 |v t | 2 V dt) γ f 2 1 0 |h t | 2 L 2 (X) dt < ∞ . (4)
Now for any (X, f ) and (X ′ , f ′ ) ∈ F we define

d F ((X, f ), (X ′ , f ′ )) . = (inf{E X (v, h) | φ v 1 .X = X ′ , f ′ = (f + ζ h 1 ) • (φ v 1 ) -1 }) 1/2 (5)
which is a true distance on F as stated in the following Theorem:

Theorem 1. The function d F given by ( 5) defines a distance on F ie. is symmetric, satisfies the triangle inequality and

d F ((X, f ), (X ′ , f ′ )) = 0 iff X = X ′ and f = f ′ ∈ L 2 (X).
Moreover, the distance is finite everywhere on F × F.

Proof. The symmetry comes from a usual process of time reversal.

Let (X, f ), (X ′ , f ′ ) ∈ F and (v, h) ∈ L 2 ([0, 1], V × L 2 (X)) such that X ′ = φ v 1 .X, f ′ = (f + ζ h 1 ) • (φ v 1 ) -1 . (6)
If we define

ṽt . = -v 1-t and ht . = -h 1-t • (φ v 1 ) -1 for any t ∈ [0, 1], then for X t . = φ v t .X (so that X 0 = X and X 1 = X ′ ) we have φ v t .X 0 = X t = φ ṽ 1-t .X 1 and (ṽ, h) ∈ L 2 ([0, 1], V × L 2 (X ′ )) with E X (v, h) = E X ′ (ṽ, h) . (7)
Since one easily checks that

ζ h 1 = -ζ h 1 • (φ v 1 ) -1 so that if f ′ = (f + ζ h 1 ) • (φ v 1 ) -1 we have f = (f ′ + ζ h 1 )
• (φ ṽ 1 ) -1 with X = φ ṽ 1 .X ′ This gives immediately the symmetry. Concerning the triangle inequality, is comes from a usual process of path concatenation. Let (X, f ), (X, f ′ ) and (X ′′ , f ′′ ) be three fshapes in

F such that d F ((X, f ), (X ′ , f ′ )) > 0 and d F ((X ′ , f ′ ), (X ′′ , f ′′ )) > 0. One easily checks that for (v, h) ∈ L 2 ([0, 1], V × L 2 (X)) and (v ′ , h ′ ) ∈ L 2 ([0, 1], V × L 2 (X ′ )) with (X ′ , f ′ ) = (φ v 1 .X, (f + ζ h 1 ) • (φ v 1 ) -1 ) and (X ′′ , f ′′ ) = (φ v ′ 1 .X ′ , (f ′ + ζ h ′ 1 ) • (φ v ′ 1 ) -1 ) then denoting for α, β > 1 such that 1/α + 1/β = 1 and s ∈ [0, 1] Cat α ((v ′ , h ′ ), (v, h)) s . = β(v ′ β(s-1/α) , h ′ β(s-1/α) • φ)1 s≥1/α + α(v αs , h αs )1 0≤s<1/α (8)
with φ = φ v 1 , we have for (ṽ, h)

. = Cat α ((v ′ , h ′ ), (v, h)) that (ṽ, h) ∈ L 2 ([0, 1], V × L 2 (X)) and (X ′′ , f ′′ ) = (φ ṽ 1 .X, (f + ζ h 1 ) • (φ ṽ 1 ) -1 ) (9) so that d F ((X, f ), (X ′′ , f ′′ )) ≤ E X (ṽ, h) 1/2 . ( 10
)
However, for α

* = (E X (v, h) 1/2 + E X ′ (v ′ , h ′ ) 1/2 )/E X (v, h) 1/2 , we check easily that E X (ṽ, h) 1/2 = E X (v, h) 1/2 + E(v ′ , h ′ ) 1/2 (11)
and the triangle inequality follows immediately. Now, if d F ((X, f ), (X ′ , f ′ )) = 0 then there exists a sequence φ n such that φ n .X = φ n (X) = X ′ and φ n → Id on X. In particular, X ′ is dense in X. As X ′ is compact (since X is), X ⊂ X ′ . By symmetry we get X = X ′ . Moreover, there exists also a sequence

ζ n ∈ L 2 (X) such that ζ → 0 in L 2 (X) and f ′ = (f + ζ n ) • φ -1 n with φ n (X) = X. We get X (f ′ -f ) 2 ≤ 2 X (f • φ -1 n -f ) 2 + 2 X (ζ n • φ -1 n ) 2 → 0
and the result is proved.

A last point to be verified is that the distance d F is finite between any pair of fshapes. Indeed, by definition there exists between any X and

X ′ ∈ X 0 a path t → φ v t .X with v ∈ L 2 ([0, 1], V ). Now, if h t = (f ′ • φ v 1 -f ), we get (v, h) ∈ L 2 ([0, 1], V × L 2 (X)) and f ′ = (f + ζ v 1 ) • (φ v 1 ) -1 so that d F ((X, f ), (X ′ , f ′ )) ≤ E X (v, h) 1/2 < ∞.
As stated in the introduction, the framework extends metamorphosis on images [START_REF] Trouvé | Metamorphoses through lie group action[END_REF] which corresponds to the case where X 0 is E itself or is the unit cube of E (with the extra condition that φ.X 0 = X 0 for any φ ∈ G). Here the support can be d-dimensional and is not fixed.

Existence of geodesics

A natural question is the existence of a minimizing geodesic between fshapes in F. We have the following Theorem:

Theorem 2. For any (X, f ), and (X ′ , f ′ ) in F, there exists (v, h) ∈ L 2 ([0, 1], V × L 2 (X)) such that E X (v, h) 1/2 = d F ((X, f ), (X ′ , f ′ )).
In particular, if

(X t , f t ) . = (φ v t .X, (f + t 0 h s ds) • (φ v t ) -1 ) (12)
the path t → (X t , f t ) can be considered as a minimizing geodesic between (X, f ) and (X ′ , f ′ ).

Proof. The proof extends the results for metamorphosis. It is sufficient to show that (v, h) → E X (v, h) is lower semi-continuous for the weak convergence of the space L 2 ([0, 1], V × L 2 (X)). Indeed, if this is the case, then from any minimizing sequence 2 we deduce that v n is bounded on L 2 ([0, 1], V ) and using the inequality

(v n , h n ) such that E X (v n , h n ) → d F ((X, f ), (X ′ , f ′ ))
1 0 X |h n,t (x)| 2 dH d (x)dt ≤ 1 0 (sup x∈X |d x φ vn t ↾ TxX | -1 ) X |h n,t (x)| 2 |d x φ vn t ↾ TxX |dH d (x) ≤ C( 1 0 |v n,t | 2 V dt)E X (v n , h n ) ( 13
)
where C is a increasing function depending only on V , we get that the sequence

(h n ) is bounded in L 2 ([0, 1], L 2 (X)). Hence, by weak compactness of strong balls in L 2 ([0, 1], V × L 2 (X)), we can assume that, up to the extraction of a sub-sequence, that (v n , h n ) weakly converges to- wards (v ∞ , h ∞ ) ∈ L 2 ([0, 1], V × L 2 (X))
and by lower semi-continuity of E X , we deduce that

E X (v ∞ , h ∞ ) ≤ d F ((X, f ), (X ′ , f ′ )) 2 . We only need to check that (X ′ , f ′ ) = (φ v∞ 1 .X, (f + ζ h∞ 1 ) • (φ v∞ 1 ) -1
). This last result follows from the fact that, if

(v n , h n ) weakly converges to (v ∞ , h ∞ ), then φ vn 1 → φ v∞ 1
uniformly on any compact sets (which is a well known result) and

(f + ζ hn 1 ) • (φ vn 1 ) -1 weakly converges to (f + ζ h∞ 1 ) • (φ v∞ 1 ) -1 .
The last weak convergence is straightforward.

The proof of the lower semi-continuity itself is done now. We know from a classical result on weak convergence that

1 0 |v ∞,t | 2 V dt ≤ lim inf n→∞ 1 0 |v n,t | 2 V dt. Moreover, we have 1 0 X |h ∞,t (x)| 2 |d x φ v∞ t ↾ TxX |dH d (x) = lim n→∞ 1 0 X h n,t (x)h ∞,t (x)|d x φ v∞ t ↾ TxX |dH d (x) (14)
by definition of the weak convergence of (h n ). Now, since

d x φ vn t → d x φ v∞ t uniformly on t ∈ [0, 1] and x ∈ X, we get 1 0 X |h ∞,t (x)| 2 |d x φ v∞ t ↾ TxX |dH d (x) = lim n→∞ 1 0 X h n,t (x)h ∞,t (x)(|d x φ vn t ↾ TxX |) 1/2 |d x φ v∞ t ↾ TxX | 1/2 dH d (x) ≤ lim inf n→∞ 1 0 X |h n,t (x)| 2 |d x φ vn t ↾ TxX |dH d (x) 1/2 1 0 X |h ∞,t (x)| 2 |d x φ v∞ t ↾ TxX |dH d (x) 1/2 (15) so that 1 0 X |h ∞,t (x)| 2 |d x φ v∞ t ↾ TxX |dH d (x) ≤ lim inf n→∞ 1 0 X |h n,t (x)| 2 |d x φ vn t ↾ TxX |dH d (x) . (16)

Karcher means on fshape bundles

We consider here the existence problem of Karcher mean on the fshape bundle F. The problem can be stated as follows: let (X i , f i ) 1≤i≤N be a family of fshapes in F. Does there exist (X * , f * ) ∈ F minimizing the sum of the square distances to each (X i , f i ), i.e. solution of the minimizing problem:

min (X,f )∈F N i=1 d F ((X, f ), (X i , f i )) 2 . ( 17
)
Note that this problem does not make sense on practical ground since usually there is no reason for a family of observed fshapes to belong to a pre-defined fshape bundle F. However, on a theoretical perspective the existence of such Karcher mean is a quite important point in addition to the existence of geodesic between any two fshapes in F. Such a result shows that F has basic important properties for further statistical analysis.

Our main result is the following: Theorem 3. For any family (X i , f i ) 1≤i≤N of fshapes in F, there exists at least a Karcher mean (X * , f * ) in F i.e. a solution of problem [START_REF] Joshi | Unbiased diffeomorphic atlas construction for computational anatomy[END_REF].

Proof. The first step of the proof is to reparametrize problem [START_REF] Joshi | Unbiased diffeomorphic atlas construction for computational anatomy[END_REF] in the Hilbert space

W . = N i=1 L 2 ([0, 1], V × L 2 (X i )) .
Indeed, using the symmetry of the problem, we can exchange the role of (X, f ) and (X i , f i ) and look for a family of minimizing paths t → (X i t , f i t ), starting from every fshape (X i , f i ) and ending on a common fshape (X, f ) ∈ F, each path being parametrized by

(v i , h i ) ∈ L 2 ([0, 1], V × L 2 (X i )). Hence, introducing W 0 . = {α = (v i , h i ) 1≤i≤N ∈ W | (X i 1 , f i 1 ) = (X 1 1 , f 1 1 ), ∀1 ≤ i ≤ N } (where (X i 1 , f i 1 ) = (φ v i 1 .X i , (f i + ζ h i 1 ) • (φ v i 1 ) -1
) we need to prove that if we define

J(α) . = N i=1 E X i (v i , h i ) (18)
then W 0 is closed and J is lower semi continuous (l.s.c.) for the weak topology on W . Since the fact that J is l.s.c. was established inside the proof of Theorem 2, we just need to check that W 0 is closed for the weak topology. So, let (α n ) n≥0 be a sequence in W 0 weakly converging to α

∞ . = (v i ∞ , h i ∞ ) 1≤i≤N ∈ W . We get that φ v i n 1 → φ v i ∞ 1
and (φ

v i n 1 ) -1 → (φ v i ∞ 1
) -1 for the uniform C 1 topology on compact sets so that in particular, if

X * = φ v 1 ∞ 1 .X 1 then since by construction φ v i n 1 .X i = φ v 1 n 1 .X 1 we get that d H (φ v i n 1 .X i , X * ) → 0 where d H denotes the Hausdorff distance 1 and φ v i ∞ 1 .X i = X * . (19)
Now, for any

1 ≤ i ≤ N , if we denote f i 1,n . = (f i + ζ h i n 1 ) and ψ i n . = (φ v i n 1 ) -1 • φ v 1 n 1 : X 1 → X i , we have (since α n ∈ W 0 and α n ⇀ α ∞ ) f i 1,n • ψ i n = f 1 1,n and f 1 1,n ⇀ f 1 1,∞ (20)
where the last statement is straightforward. Let us check that

f i 1,n • ψ i n ⇀ f i 1,∞ • ψ i ∞ so that we will get f i 1,∞ • ψ i ∞ = f 1 1,
∞ and since we have [START_REF] Lee | Optic nerve head and peripapillary morphometrics in myopic glaucoma[END_REF] we will get α ∞ ∈ W 0 and W 0 weakly closed. For that, let g : X 1 → R be a Lipschitz mapping on X 1 and denote now f, f ′

X 1 . = X 1 f (x)f ′ (x)dH d (x) the usual dot product on L 2 (X 1 ). Since X 1 is compact and H d (X 1 ) < ∞, we have that H d ↾ X 1
is a Radon measure on the compact metric space X 1 so that such Lipschitz mapping g are dense in L2 (X 1 ) and we just need to check that

f i 1,n • ψ i n -f i 1,∞ • ψ i ∞ , g X 1 → 0 to prove the weak convergence. We have f i 1,n • ψ i n -f i 1,∞ • ψ i ∞ , g X 1 = f i 1,n • ψ i n -f i 1,n • ψ i ∞ , g X 1 An + f i 1,n • ψ i ∞ -f i 1,∞ • ψ i ∞ , g X 1 Bn . (21)
Concerning the B term, by change of variable and using the fact that the Jacobian d(ψ i ∞ ) -1 is bounded on X 1 we get

B n = f i 1,n -f i 1,∞ , g • (ψ i ∞ ) -1 |d(ψ i ∞ ) -1 .ξ| X i → 0 , ( 22
)
where ξ is the unit d-vector representing the tangent space defined H d ↾ X i -a.e. on X i . Concerning the A term, we have again by change of variable

A n = f i 1,n , g • (ψ i n ) -1 (|d(ψ i n ) -1 .ξ| -|d(ψ i ∞ ) -1 .ξ|) X i + f i 1,n , (g • (ψ i n ) -1 -g • (ψ i ∞ ) -1 )|d(ψ i ∞ ) -1 .ξ| X i . (23)
The first term of ( 23) is bounded by

| f i 1,n | L 2 (X i ) |g| L ∞ (X 1 ) |d(ψ i n ) -1 .ξ| -|d(ψ i ∞ ) -1 .ξ| L ∞ (X i ) → 0 (24)
since the weak convergence of

f i 1,n implies that | f i 1,n | L 2 (X i ) stays bounded.
The second term of ( 23) is bounded by

k g | f i 1,n | L 2 (X i ) (ψ i n ) -1 -(ψ i ∞ ) -1 L ∞ (X i ) |d(ψ i ∞ ) -1 | L ∞ (X i ) → 0 ( 25
)
where k g is the Lipschitz constant of g on X i .

The tangential model

The previous metamorphosis metrics can be also approximated in a simpler setting, which becomes closer to a natural extension of the classical image framework in the context of fshapes. We shall call it the tangential model. Instead of computing the cost of the functional evolution t → ζ t along a time dependent L 2 (X t ) metric as

1 0 X | ζt (x)| 2 |d x φ v t ↾ TxX |dH d (x)dt = 1 0 Xt | ζt • (φ v t ) -1 | 2 (y)dH d (y)dt (26)
one can freeze the metric to its initial value X 0 = X neglecting the change of measure weight:

1 0 X | ζt (x)| 2 |d x φ v t ↾ TxX |dH d (x)dt ≈ 1 0 X | ζt (x)| 2 dH d (x)dt (27)
which gives after optimization with fixed end points condition ζ 1 the usual L 2 (X) cost

|ζ 1 | 2 X . = X |ζ 1 (x)| 2 dH d (x) (28)
which can be interpreted naturally as the log-likelihood of a Gaussian noise in the statistical framework.

Note as we say in introduction of this section, that this is the usual noise term that appears in the classical situation of image matching and this can be considered as some kind of tangential metric for the functional part. Concerning the geometrical part, we can keep the usual

L 2 ([0, 1], V ) penalization to define ẼX (v, ζ) . = γ V 2 1 0 |v t | 2 V dt + γ f 2 |ζ| 2 X ( 29
)
and consider

dF ((X, f ), (X ′ , f ′ )) . = (inf{ ẼX (v, ζ) | X ′ = φ v 1 .X, f ′ = (f + ζ) • (φ v 1 ) -1 }) 1/2 . ( 30
)
With that definition, dF is no more a distance satisfying the symmetry and triangle inequality as previously, but the minimization of ẼX can be seen as the minimization of the log-likelihood of (v, ζ) given (X, f ) and (X ′ , f ′ ) for a quite natural statistical model.

In this framework, one can still obviously consider the problem of existence for two given fshapes (X, f ) and

(X ′ , f ′ ) in F of an optimal (v * , ζ * ) such that            (X ′ , f ′ ) = (φ v * 1 .X, (f + ζ * ) • (φ v * 1 ) -1 ) and ẼX (v * ; ζ * ) = min{ ẼX (v, ζ) | (v, ζ) ∈ F (X ′ ,f ′ ) } where F (X ′ ,f ′ ) . = {(v, ζ) ∈ L 2 ([0, 1], V ) × L 2 (X) | (X ′ , f ′ ) = (φ v 1 .X, (f + ζ) • (φ v 1 ) -1 )} (31)
The proof existence of (v * , ζ * ) can be done along the very same lines that the proof of existence of geodesics and is omitted. The proof is even more simpler since the functional metric is frozen and does not change along the path. We get eventually the following result: Theorem 4. For (X, f ) and (X ′ , f ′ ) in F, there exists an optimal solution of the exact matching problem between (X, f ) and (X ′ , f ′ ) i.e. there exists

(v * , ζ * ) in L 2 ([0, 1], V ) × L 2 (X) such that (31) holds.
3 Dissimilarity measure between fshapes

Dissimilarity measures : state of the art

The previous framework has been focusing so far on the comparison of fshapes within one given common bundle F. It is clear that it becomes inoperant in most practical situations where the elements of a dataset cannot be assumed to belong to the same bundle. The reason is that the geometrical supports of two subjects need not be obtained from one the other by a deformation belonging to the group G. In more usual cases of shape spaces such as sets of images, landmarks, curves or surfaces, the general approach has been to introduce dissimilarity measures between the subjects, which provides data attachment terms to perform inexact registration and atlas estimation in practice.

For images and landmarks, L 2 distances are natural and have been commonly used to compare such objects. The case of curves, surfaces and submanifolds is theoretically more involved and has drawn consistent attention in several fields of mathematics. In computational anatomy, one possible and very powerful setting was proposed through the adaptation of the concept of currents [START_REF] Glaunès | Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF][START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF]. The dissimilarity measures between oriented submanifolds of given dimension defined in these frameworks derive from Hilbert metrics on spaces of differential forms and combine the advantages of not relying on parametrizations or point to point matchings, of being easily computable for discrete shapes and robust to shape sampling. One potential drawback is the issue of consistent orientation that is required in the currents' representations. More recently, an alternative methodology based on varifolds was introduced and implemented in [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF], which has the interest of being close to the idea of currents while overcoming the problem of orientation. Now, fshapes do pose additional difficulties with respect to the definition of dissimilarity measures. This problem has been only addressed very recently in [START_REF] Charon | Functional currents : a new mathematical tool to model and analyse functional shapes[END_REF] where the authors define the extended notion of functional currents that, similarly to usual currents, embeds fshapes in some dual of spaces of differential forms. The Hilbert metrics that are defined on functional currents can be then used again as dissimilarity measures to perform inexact registration between fshapes, as exposed in the article. Functional currents could be thus very well used in the generalized problem of atlas estimation that we focus on in this paper. This has been touched upon, essentially from the numerical point of view, in [START_REF] Charon | Analysis of geometric and functional shapes with extensions of currents. Application to registration and atlas estimation[END_REF] (chapter 4). Although the results and algorithms presented in this paper could be transposed quite easily to fshape dissimilarities provided by functional currents, we shall work instead with the slightly modified approach of varifolds. As aforementioned, this is a way to get rid of shape orientation which is usually quite desirable in applied situations. In the following section, we briefly present the functional varifolds' mathematical setting, that stands for a natural generalization of varifolds to fshapes. The more technical and discrete computations shall be detailed in Section 6 dedicated to the algorithmic part itself.

Functional varifolds

The mathematical concept of varifold goes back to geometric measure theory and the original exposition of F. Almgren [START_REF] Almgren | Plateau's Problem: An Invitation to Varifold Geometry[END_REF] which was further developed considerably by W. Allard in [START_REF] Allard | On the first variation of a varifold[END_REF] subsequently. The link to computational anatomy's problems is fairly more recent and is in most part presented in [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF] or [START_REF] Charon | Analysis of geometric and functional shapes with extensions of currents. Application to registration and atlas estimation[END_REF]. We shall frequently refer the reader to these references for additional mathematical details. What we present in this section is an extension of varifolds to represent functional shapes, which is much related to the functional current idea of [START_REF] Charon | Functional currents : a new mathematical tool to model and analyse functional shapes[END_REF].

Representation of fshapes as functional varifolds

As previously, we call E the n-dimensional vector space embedding all geometrical supports of fshapes. For any integer 1 ≤ d ≤ n, we will write G d (E) for the Grassmann manifold of all ddimensional (non-oriented) subspaces of E. G d (E) is a compact manifold that can be embedded trivially in the space L(E) by identifying any V ∈ G d (E) with the orthogonal projector p V on V . The usual space of d-dimensional varifolds is defined intuitively as the space of spatially spread distributions of Grassmannians or, more rigorously, as finite Borel measures on the product E × G d (E) (cf [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF]). Now, to account for the existence of signals at each point, one can simply augment the varifold representation with a signal component, which we formalize by the following definition: Definition 1. We say that µ is a d-dimensional functional varifold (fvarifold in short) if µ is a Borel finite measure on the space

E × G d (E) × R or equivalently if µ ∈ C 0 (E × G d (E) × R) ′ .
Note that we consider only the case of real-valued signals here, but this could be extended to different situations, as vector or tensor-valued signals (cf [START_REF] Charon | Analysis of geometric and functional shapes with extensions of currents. Application to registration and atlas estimation[END_REF]). Now, any functional shape (X, f ), with X a d-dimensional rectifiable subset, can be represented by a functional varifold µ (X,f ) defined by:

µ (X,f ) (ω) = X ω(x, T x X, f (x))dH d (x) (32)
where H d is the d-dimensional Hausdorff measure on E. Some particularly simple functional varifolds are the Diracs that, in this context, can be written in the form δ (x,V,f ) with x ∈ E, V ∈ G d (E) and f ∈ R, and act on any function

ω ∈ C 0 (E × G d (E) × R) by the relation: δ (x,V,f ) (ω) = ω(x, V, f ) . (33)
We shall detail, in Section 6, how to approximate polyhedral fshapes by a finite sums of Dirac in order to perform practical computations. Now, an important point is to express the way that geometrico-functional transformations act on functional varifolds, in such a way that this action is consistent with the transport of fshapes that we have been considering in section 2.2. In a very similar fashion as with regular varifolds (cf [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF]), one can express it by usual pull-back and push-forward operations. Let φ ∈ Diff(E) and ζ a measurable function on E, this is given by the following set of equations:

   ∀ω ∈ C 0 (E × G d (E) × R), ((φ, ζ) * µ) (ω) = µ ((φ, ζ) * ω) where ((φ, ζ) * ω) (x, V, f ) = |d x φ ↾ V | ω(φ(x), d x φ(V ), f + ζ) (34)
where for V ∈ G d (E), |d x φ ↾ V | denotes as previously the Jacobian of φ along subspace V (i.e. the volume change along V at point x) and d x φ(V ) is the image of V by the invertible linear application d x φ. It is then a simple verification that one has the following consistency property:

Proposition 1. Let (X, f ) be a fshape where X is a d-dimensional rectifiable subset and f a L 2 function on X. If φ ∈ Diff(E) and ζ ∈ L 2 (X) then: (φ, ζ) * µ (X,f ) = µ (φ(X),(f +ζ)•φ -1 ) .
The proof does not involve any additional difficulty than in the case of usual varifolds, for which we refer the reader to [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF].

Thus, functional varifolds enable all possible fshapes of given dimension to be embedded into a common space of distributions. The next step then is to equip such a space with a metric that shall induce a dissimilarity measure on the set of fshapes. This can be efficiently addressed by introducing reproducing kernels on the product E × G d (E) × R, as we explain in the following text.

RKHS of functional varifolds

The use of reproducing kernel in the context of currents and varifold spaces has been argued and implemented many times, for instance in [START_REF] Glaunès | Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF][START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF][START_REF] Charon | Analysis of geometric and functional shapes with extensions of currents. Application to registration and atlas estimation[END_REF]. The general scheme is to consider Hilbert spaces embedded in C 0 (E × G d (E) × R) which are built as the reproducing kernel Hilbert space (RKHS) associated to a certain positive kernel on E × G d (E) × R. A natural and convenient (but not exhaustive) way to define kernels on product spaces is to consider tensor products of kernels, the classical result from kernel theory being that: Lemma 1. Let A and B be two sets and k A , k B positive kernels respectively on A and B. Then

k A ⊗ k B defined for all a 1 , a 2 ∈ A and b 1 , b 2 ∈ B by: k A ⊗ k B ((a 1 , b 1 ), (a 2 , b 2 )) . = k A (a 1 , a 2 )k B (b 1 , b 2 ) is a positive kernel on A × B.
Now going back to the case of functional varifolds itself, one has the following property:

Proposition 2.
Let k e be a positive kernel on the space E such that k e is continuous, bounded and for all x ∈ E, the function k e (x, .) vanishes at infinity, k t a kernel on the manifold G d (E) that is also continuous, and k f a kernel on R continuous, bounded and such that k f (f, .) vanishes at infinity for all f ∈ R. Then W , the RKHS associated to the positive kernel k

. = k e ⊗ k t ⊗ k f , is continuously embedded into the space C 0 (E × G d (E) × R).
Proof. The proof is very similar to the one of [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF] in the case of varifolds. By lemma 1, k is indeed a positive kernel on E × G d (E) × R and by definition:

k (x, V, f ), (x, Ṽ , f ) = k e (x, x) k t (V, Ṽ ) k f (f, f ) (35)
and so, thanks to the assumptions on the kernels, k((x, V, f ), .)

is continuous on E × G d (E) × R and belongs to C 0 (E × G d (E) × R).
The vector space W 0 generated by these functions is thus included in C 0 (E × G d (E) × R). Moreover, if ω ∈ W 0 , by the reproducing kernel property, we have that:

ω(x, V, f ) = δ (x,V,f ) (ω) = k((x, V, f ), .), ω W . With Cauchy-Schwarz inequality: |ω(x, V, f )| ≤ k((x, V, f ), .) W . ω W . In addition, k((x, V, f ), .) W = k((x, V, f ), (x, V, f ))
and all three kernels k e , k t and k f are bounded so that k is also bounded. We conclude that |ω| ∞ ≤ |k| ∞ . ω W . Thus Cauchy sequences in W 0 for the W -norm are also Cauchy sequences for the infinite norm. It results that all their limits belong to

C 0 (E × G d (E) × R) and therefore W is included in C 0 (E × G d (E) × R).
The previous inequality then holds for all ω ∈ W , which shows that the inclusion embedding ı :

W ֒→ C 0 (E × G d (E)) is indeed continuous.
Consequently, there exists a continuous mapping i * of the space of fvarifolds C 0 (E × G d (E) × R) ′ into the dual of W . This induces a pseudo-distance on fvarifolds resulting from the Hilbert structure of the RKHS. If we introduce the isometry K W : W ′ → W defined by K W µ, ω W = µ(ω) for all µ ∈ W ′ and ω ∈ W , by the reproducing kernel property, we know that

K W δ (x,V,f ) = k((x, V, f ), .). Then for all x 1 , x 2 ∈ E, V 1 , V 2 ∈ G d (E) and f 1 , f 2 ∈ R, δ (x 1 ,V 1 ,f 1 ) , δ (x 2 ,V 2 ,f 2 ) W ′ = K W δ (x 1 ,V 1 ,f 1 ) , K W δ (x 2 ,V 2 ,f 2 ) W = K W δ (x 2 ,V 2 ,f 2 ) (x 1 , V 1 , f 1 )
and thus we have the following expression for the inner product between two Diracs:

δ (x 1 ,V 1 ,f 1 ) , δ (x 2 ,V 2 ,f 2 ) W ′ = k e (x 1 , x 2 ) k t (V 1 , V 2 ) k f (f 1 , f 2 ) . (36)
Now, if (X, f ) and (Y, g) are two fshapes (of dimension d) and µ (X,f ) and µ (Y,g) the associated fvarifolds (as defined previously), one can prove easily the following formula:

µ (X,f ) , µ (Y,g) W ′ = X Y k e (x, y) k t (T x X, T y Y ) k f (f (x), g(y))dH d (x)dH d (y) . (37)
One can thus propose a dissimilarity measure between any two fshapes that is simply the norm of the difference in the space of fvarifolds i.e. µ (X,f ) -µ (Y,g) 2 W ′ , which can be computed easily using (37). We shall detail more specifically the discrete expressions in Section 6.

Yet, as we mentioned earlier, this does not necessarily result from a real distance on the space of varifolds because the dual application ı * does not need to be an embedding. This is actually the case if and only if the RKHS W is dense in C 0 (E × G d (E) × R), in which case the kernel k is said to be C 0 -universal. In our previous construction, this holds in particular if all kernels k e , k t and k f are themselves C 0 -universal. In [START_REF] Carmeli | Vector valued reproducing kernel hilbert spaces and universality[END_REF], authors study thoroughly construction of C 0 -universal kernels on vector spaces. Notably, it is proven that all Gaussian kernels satisfy this property. Thus, such kernels are easy to provide in the case of kernels k e and k f defined on the spaces E and R. As for kernel k t on the Grassmann manifold, kernels can be defined by using the embedding of G d (E) in L(E) (we refer to [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF] for more details), and therefore one can obtain a similar notion of Gaussian kernels that can be also shown to verify the C 0 -universality property.

To conclude this section, let us mention a generalization of proposition 2 to higher-order regularities that shall be useful in the rest of the paper. Proposition 3. Let k be a positive kernel on the product space E × G d (E) × R such that k is continuously differentiable of order 2p and such that all the derivatives of k up to the order p are bounded. Assume in addition that for any

(x, V, f ) ∈ E × G d (E) × R, the function k((x, V, f ), .
) and all its derivatives up to order p vanish at infinity. Then, the RKHS associated to k is continuously embedded into

C p 0 (E × G d (E) × R).
The proof follows the same pattern as the one of proposition 2 and can be also adapted from the one given in [START_REF] Glaunès | Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF] (chapter 2).

Properties of the metrics

Control results

We now show a few control results on fvarifolds norm that shall be useful for the following. As a measure on E × G d (E) × R, we will call the total variation norm of a given fvarifold µ the quantity µ(E × G d (E) × R). In particular, one can easily check that the total variation norm of a rectifiable fvarifold µ (X,f ) equals H d (X). We have the following control of W ′ -norms:

Proposition 4. For a RKHS W continuously embedded into C 0 (E × G d (E) × R), there exists a constant c W > 0 such that for all µ ∈ C 0 (E × G d (E) × R) ′ : µ W ′ ≤ c W .µ(E × G d (E) × R)
Proof. From the continuous embedding property, we have the existence of c W > 0 such that for all ω ∈ W ′ , |ω| ∞ ≤ c W ω W . Now, by definition of the kernel metric µ 2 W ′ = µ(K W µ) for K W the Riesz isometry between W ′ and W and thus:

µ 2 W ′ = µ(K W µ) = K W µ(x, V, f )dµ(x, V, F ) ≤ |K W µ| ∞ µ(E × G d (E) × R) ≤ c W K W µ W µ(E × G d (E) × R) ≤ c W µ W ′ µ(E × G d (E) × R)
so that dividing on both sides by µ W ′ , we obtain the result.

A direct corollary is that the RKHS norm of a rectifiable fvarifold is controlled by the d-volume of its geometrical support.

A second important property to establish is the continuity of RKHS norms with respect to deformations of geometrical supports. We first consider general C 1 -diffeomorphisms of E such that φ and dφ tends to Id at infinity and we will denote

φ -Id 1,∞ = |φ -Id| ∞ + |dφ -Id| ∞ .
The result we show is the following: Proposition 5. Let X be a d-dimensional rectifiable subset of E of finite volume and f ∈ L 2 (X). Assume that W is constructed as in proposition 2 and is continuously embedded into

C 1 0 (E × G d (E) × R). Then: µ (X,f ) -µ (φ(X),f •φ -1 ) W ′ → 0 as φ -Id 1,∞ → 0.
Proof. We start by writing:

µ (X,f ) = X δ (x,TxX,f (x)) dH d (x)
and, in the same way,

µ (φ(X),f •φ -1 ) = φ(X) δ (y,Tyφ(X),f •φ -1 (y)) dH d (x) .
Now, applying the area formula (corollary 3.2.20 in [START_REF] Federer | Geometric measure theory[END_REF]), we have:

µ (φ(X),f •φ -1 ) = X δ (φ(x),dxφ(TxX),f (x)) |d x φ ↾ TxX |dH d (x) .
where |d x φ ↾ TxX | is the local volume change at x along T x X. We obtain therefore:

µ (X,f ) -µ (φ(X),f •φ -1 ) W ′ ≤ X |d x φ ↾ TxX |.δ (φ(x),dxφ(TxX),f (x)) -δ (x,TxX,f (x)) W ′ dH d (x) . (38)
Now, we assumed that the kernel k of W is built as a tensor product k = k e ⊗ k t ⊗ k f . We can then write W = W g ⊗ W f where W g is the RKHS of kernel k e ⊗ k t and W f the one of kernel k f . It is then straightforward to check that:

|d x φ ↾ TxX |.δ (φ(x),dxφ(TxX),f (x)) -δ (x,TxX,f (x)) W ′ = δ f (x) W ′ f . |d x φ ↾ TxX |.δ (φ(x),dxφ(TxX)) -δ (x,TxX) W ′ g Since δ f (x) W ′ f = k f (f (x), f (x)) and k f is bounded, we have δ f (x) W ′ f ≤ Cte. Moreover: |d x φ ↾ TxX |.δ (φ(x),dxφ(TxX)) -δ (x,TxX) W ′ g ≤ |d x φ ↾ TxX | -1 . δ (φ(x),dxφ(TxX)) W ′ g + δ (φ(x),dxφ(TxX)) -δ (x,TxX) W ′ g . (39)
Focusing on the first term, we have again, since k e and k t are bounded, that δ (φ(x),dxφ(TxX)) W ′ g is uniformly bounded on X. In addition, the volume variation |d x φ ↾ TxX | -1 converges to 0 uniformly on X whenever φ -Id 1,∞ → 0. As for the second term in the sum of (39), we know that δ (φ(x),dxφ(TxX)) -δ (x,TxX) W ′ g = sup ω Wg =1 (δ (φ(x),dxφ(TxX)) -δ (x,TxX) )(ω) and for any ω ∈ W g , we have:

(δ (φ(x),dxφ(TxX)) -δ (x,TxX) )(ω) = |ω(φ(x), d x φ(T x X)) -ω(x, T x X)| ≤ |dω| ∞ .d E×G d (E) ((φ(x), d x φ(T x X)), (x, T x X)) and d E×G d (E) ((φ(x), d x φ(T x X)), (x, T x X)
) is obviously uniformly upper bounded for x ∈ X by Cte. φ -Id 1,∞ . In addition, with the assumptions on the kernels, we also have that |dω| ∞ ≤ Cte ω Wg . It results eventually that δ (φ(x),dxφ(TxX)) -δ (x,TxX) W ′ g is bounded above by Cte. φ-Id 1,∞ . Now, since both terms in (39) converge to 0 uniformly for x ∈ X, it results the same behavior for the integral of (38), which completes the proof. Now, going back to the previous model of deformations, we consider diffeomorphisms φ = φ v 1 obtained as the flow at time 1 of a time varying vector fields v ∈ L 2 ([0, 1], V ), where V is a given Hilbert space of vector fields on E. The result of proposition 5 implies the following important corollary:

Corollary 1. If V is continuously embedded into C 2
0 (E, E) then, for any fixed and bounded rectifiable set X of finite volume and

L 2 signal f on X, the application v → µ (φ v 1 (X),f •(φ v 1 ) -1 ) is weakly continuous from L 2 ([0, 1], V ) to W ′ .
Proof. This essentially relies on classical results on differential equations and flows that can be found in [START_REF] Glaunès | Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF][START_REF] Younes | Shapes and diffeomorphisms[END_REF]. It is shown in particular that, since V ֒→ C 2 0 (E, E), diffeomorphisms φ v 1 and their differentials tend to Id at infinity. In addition, if v n is a sequence that weakly converges to v in L 2 ([0, 1], V ), then φ vn 1 and dφ vn 1 converge respectively to φ v 1 and dφ v 1 uniformly on every compact subset of E. Then, thanks to proposition 5, we are allowed to conclude that

µ (φ vn 1 (X),f •(φ vn 1 ) -1 ) W ′ --→ µ (φ v 1 (X),f •(φ v 1 ) -1
) , which proves the weak continuity.

Variation formula for fvarifold metrics

We now extend the variation formula shown in [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF] to the present setting of functional varifolds.

The aim is to have a theoretical description of the variations of fvarifold metrics with respect to variations of a functional shape in both its geometrical support and signal. As we shall see, the behaviour is qualitatively similar, except for the orientation, to the setting of fcurrents that was examined in [START_REF] Charon | Analysis of geometric and functional shapes with extensions of currents. Application to registration and atlas estimation[END_REF]. Fixing a functional shape (X, f ) with X a compact submanifold and f : X → R a C 1 signal on X, we wish to compute variations of terms like µ (X,f ) , µ W ′ for any µ ∈ W ′ . Using the previous isometry K W , we know that µ (X,f ) , µ

W ′ = K W µ (X,f ) , K W µ W = µ (X,f ) (ω)
with ω = K W µ thanks to the reproducing kernel property and thus one is led to consider variation of terms µ (X,f ) (ω) for ω ∈ W . Rigorously, such variations can be expressed with respect to the action of infinitesimal deformations in the geometrico-functional domain. One can define such infinitesimal geometrico-functional deformations by considering, exactly as in the previous tangential model, a C 1 compactly supported vector field v on E and a L 2 function h on X. Now we shall denote by φ t the flow of v at time t and by ψ t the application defined on X × R by ψ t (x, m) = (φ t (x), m + th(x)). In addition, ψ can be extended straightforwardly to the whole space E × R by setting

ψ t (x, m) = (φ t (x), m + t h(x))
where h↾ X = h and h = 0 elsewhere. The function h is then measurable on E and the push-forward action ψ * t ω is well-defined by equation [START_REF] Trouvé | Handbook of Mathematical Imaging, chapter Shape spaces[END_REF]. The variation we wish to compute can be expressed by:

d dt ↾ t=0 µ (φt(X),(f +th)•φ -1 ) (ω) = d dt ↾ t=0 X ψ * t ω(x, T x X, f (x))dH d (x) = X d dt ↾ t=0 ψ * t ω(x, T x X, f (x))dH d (x) = X (£ (v,h) ω)(x, T x X, f (x))dH d (x)
where

£ (v,h) ω . = d dt ↾ t=0 ψ *
t ω is a notation for the derivative of function ω in the direction of (v, h). In the following, we shall adopt the shortcut notation X g to denote the integral X g(x)dH d (x). Now, the result is the following: Theorem 5. Let X be a compact submanifold and f a C 1 function on X. Let v be any C 1 vector field with compact support defined on E and h a L 2 function on X. We denote by ψ t (x, m) . = (φ t (x), m+ th(x)). Then, we have for any

C 1 function (x, V, m) → ω(x, V, m) on E × G d (E) × R: X (£ (v,h) ω)(x, T x X, f (x))dH d (x) = X ∂ω ∂x -div X ∂ω ∂V -ωH X |v ⊥ + ∂ω ∂m .(h -∇f, v ⊤ ) + ∂X ν, ωv ⊤ + ∂ω ∂V |v ⊥ (40)
where v ⊤ and v ⊥ denote the tangential and normal part of v along X, ν is the unit outward normal along ∂X, and H X the mean curvature vector to X.

The proof is given in appendix A. This formula is interesting at several levels because it gives qualitative information on the gradients of the kernel metrics with respect to points and signal values of the fshape (X, f ). In particular, as already pointed out in [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF], for a constant signal f , we see that, for points in the interior of X, the variation with respect to the geometrical support only involves the orthogonal component of v which means that the gradient of the metric is orthogonal to the shape. This is somehow quite natural since tangential components of v do not change the shape itself in that situation. However, we see that for non-constant signals, a term involving the tangential component v ⊤ appears in the variation and is concentrated on regions of important gradient of the signal f . The second important consequence to mention is the presence of special terms on the boundary of X. On points located on ∂X, the gradient is orthogonal to the boundary (but not necessarily to X). The presence of these boundary singularities have important consequences on the numerical behaviour of the gradients that we shall address with more details in section 9.

Mathematical formulation of atlas estimation

To go beyond Karcher means within a single bundle (as exposed in Section 2.4) and address the problem of atlas computation, we have to consider now a family of (possibly noisy) observed fshapes (X i , f i ) 1≤i≤N . The main point here is the optimization of a template (X, f ) given the observations. We propose, in this section, a way to formalize the atlas estimation as a variational problem. The existence of solutions to such problems will be examined in the next section.

Template space

Ideally, the template should be optimized among all the possible templates to avoid any bias effect towards a particular configuration but this looks as a quite badly posed optimization problem. We follow here a more simple and secure route by introducing a priori a large but restricted space of possible templates as another vector bundle F 0 generated by a hypertemplate (X 0 , f 0 ). As presented previously, F 0 is defined from the (X 0 , f 0 ) as the set of every fshape defined from the mapping

(φ, ζ) → (φ.X 0 , (f 0 + ζ) • φ -1 ) for (φ, ζ) ∈ G 0 × L 2 (X 0 ) with G 0 .
= G V 0 is the group of diffeomorphisms generated by an additional RKHS space V 0 of vector fields on E:

F 0 . = { (X, f ) | (X, f ) = (φ.X 0 , (f 0 + ζ) • φ -1 ), (φ, ζ) ∈ G 0 × L 2 (X 0 )} . ( 41
)
One could take V 0 = V but we may want less regularity on V 0 than on V to generate a larger space and get closer to the ideal case. Now, given any (X, f ) ∈ F 0 , we should consider the associated fshape bundle F X

F X . = { (X ′ , f ′ ) | (X ′ , f ′ ) = (φ.X, (f + ζ) • φ -1 ), (φ, ζ) ∈ G × L 2 (X)} (42)
corresponding to the previously introduced space F.

Variational formulation (metamorphosis riemannian setting)

In the Riemannian setting, the problem of atlas estimation can be addressed as the minimization of

J riem 1 ((X, f ), (X i , f i ) i ) . = d F 0 ((X 0 , f 0 ); (X, f )) 2 + N i=1 d F X ((X, f ); ( Xi , f i )) 2 + g( Xi , f i ; X i , f i ) (43) on (X, f ) ∈ F 0 and ( Xi , f i ) 1≤i≤N ∈ F N X ( 44
)
where g is the "out-of-orbit" or previous data attachment term, which, in that case, can be defined, as we explained in section 3.2 by fvarifold Hilbertian metrics. The problem can be then reformulated equivalently as a minimization problem on

       (v 0 , h 0 ) ∈ L 2 ([0, 1], V 0 × L 2 (X 0 )) (v i , h i ) ∈ L 2 ([0, 1], V × L 2 (X)) with X = φ v 0 1 .X 0 (45) of J riem 2 ((v 0 , h 0 ), (v i , h i ) i ) . = γ V 0 2 1 0 |v 0 t | 2 V 0 dt + γ 0 f 2 1 0 X 0 |h 0 t | 2 |d x φ v 0 t ↾ TxX |dH d (x) + N i=1 γ V 2 1 0 |v i t | 2 V dt + γ f 2 1 0 X |h i t | 2 |d x φ v i t ↾ TxX |dH d (x) + γ W 2 g(( Xi , f i ), (X i , f i )) (46)
where   

f i = (f + ζ h i 1 ) • (φ v i 1 ) -1 and f = (f 0 + ζ h 0 1 ) • (φ v 0 1 ) -1 .
(47)

Variational formulation (tangential setting)

The atlas estimation problem can also be formulated using the tangential metric and with a statistical flavor. It seems quite natural to choose the template as the pivotal point for the linearization towards the hypertemplate (X 0 , f 0 ) and the observations (X i , f i ) that is to freeze the functional metric on the L 2 (X) space.

The derivation is then quite straightforward and we jump directly to the associated minimization problem which is the minimization of the functional

J tan ((v 0 , ζ 0 ), (v i , ζ i ) i ) . = γ V 0 2 1 0 |v 0 t | 2 V 0 dt + γ 0 f 2 |ζ 0 | 2 X + N i=1 γ V 2 1 0 |v i t | 2 V dt + γ f 2 |ζ i | 2 X + γ W 2 g(( Xi , f i ), (X i , f i )) (48) in the variables        (v 0 , ζ 0 ) ∈ L 2 ([0, 1], V 0 ) × L 2 (X) (v i , ζ i ) ∈ L 2 ([0, 1], V ) × L 2 (X) with X = φ v 0 1 .X 0 (49)
where we have

   f i = (f + ζ i ) • (φ v i 1 ) -1 and f = f 0 • (φ v 0 1 ) -1 + ζ 0 . (50)
The signals (ζ i ) i=1,..,N will be named the residual functions of the atlas. From this general setting, we can consider an even simpler situation when f 0 ≡ 0. This new setting is equivalent to the minimization of

J tan 0 ((v 0 , f ), (v i , ζ i ) i ) . = γ V 0 2 1 0 |v 0 t | 2 V 0 dt + γ 0 f 2 |f | 2 X + N i=1 γ V 2 1 0 |v i t | 2 V dt + γ f 2 |ζ i | 2 X + γ W 2 g(( Xi , f i ), (X i , f i )) (51) in the variables        (v 0 , f ) ∈ L 2 ([0, 1], V 0 ) × L 2 (X) (v i , ζ i ) ∈ L 2 ([0, 1], V ) × L 2 (X) with X = φ v 0 1 .X 0 (52)
where we have

f i = (f + ζ i ) • (φ v i 1 ) -1 . ( 53
)
5 Existence results for fshape atlases

Introduction

The existence of an atlas for a population of fshapes with fcurrents or fvarifold based data term does not follow from the same arguments than in the more well-known pure geometrical case. The optimization of a signal f on the geometrical template and for the more sophisticated model of additional residuals ζ i to match each observations (X i , f i ) introduces new difficulties.

The main one is that the fvarifold data term µ (X,f i ) -µ (X,f )

2

W is not lower semi-continuous for the weak convergence in L 2 (X) as a function of f for X fixed. This comes from the nonlinearities appearing in the various kernels defining the RKHS norm on W . Yet, interestingly, it remains continuous for the weak convergence on the space of measures. We prove that minimizing measure sequences are tight and that existence results can be established under quite general assumption for the extended functional.

However, the existence of strong solution (i.e. a true fshape template (X, f )) does not seem to be guaranteed in full generality but is true as soon as the regularization on the L 2 penalization on f (and on the residuals ζ i when introduced) is strong enough. This condition, although imposed theoretically for the existence of a fshape solution, is not clearly required for practical fshape atlases that we have computed so far. We show also that the computed template function and residuals are smooth proving a regularization effect of fvarifold attachment terms. The result holds in both the tangential and metamorphosis setting for fshapes, which is summed up by Theorem 6 (tangential model) and Theorem 7 (metamorphosis model), whose proofs are the core issues of this section.

The most technical part of the proof is the existence of the fshape template (X, f ) on which we will focus in the first place. The proof we present follows basically the so-called direct method of geometric measure theory where we first show the existence of a minimizer in a larger space of varifolds and then show that this solution does indeed result from a true fshape.

Existence of the template fshape

We shall focus first on the tangential setting of 4.3. In this part, we will also consider the simpler situation where all variables v i and ζ i in equation ( 48) are frozen and show that a minimum of the functional exists in the variables X and f . We will start by assuming that, in addition, X is also a fixed bounded and finite volume d-dimensional rectifiable subset of E.

Existence with X fixed

This subsection is dedicated to the proof of the following proposition:

Proposition 6. Assume that W is continuously embedded in C 2 0 (E × G d (E) × R)
, that X and (X i ) 1≤i≤N are finite volume bounded d-dimensional rectifiable subsets and that

f i ∈ L 2 (X i ) for 1 ≤ i ≤ N . Assume that γ f /γ W is large enough with γ f , γ W ≥ 0. Then J X (f ) . = γ f 2 X |f (x)| 2 dH d (x) + 1 2 N i=1 γ W µ (X i ,f i ) -µ (X,f ) 2 W ′
achieves its minimum on L 2 (X) and any minimizer

f * is such that f * ∈ L ∞ (X). Moreover, if X is also a C p submanifold and W ֒→ C m 0 (E ×G d (E)×R) with m ≥ max{p, 2}, then f * ∈ C p-1 (X).
The proof proceeds in two steps. The first one is to show the existence of a minimizer in the larger space of fvarifolds: namely, for X fixed as previously, we consider the set M X of all Borel finite measures ν on E ×G d (E)×R such that:

h(x, V )dν(x, V, f ) = X h(x, T x X)dH d (x) . (54)
for all continuous and bounded applications h on E × G d (E). These are fvarifolds that have a particular marginal on E × G d (E). Note that any varifold µ (X,f ) for f a L 2 function on X belongs to M X . In addition, M X is a closed subset for the usual weak convergence of measures defined by

ν n ⇀ ν ∞ if for any ω ∈ C b (E × G d (E) × R) ν n (ω) → ν ∞ (ω) . (55)
Then, we consider the extended functional J defined on measures ν ∈ M X by

J(ν) . = γ f 2 ν(|f | 2 ) + γ W 2 N i=1 ν -µ (X i ,f i ) 2 W ′ (56)
where ν(|f | 2 ) is the notation we shall use for |f | 2 dν. Then, Lemma 2. There exists ν * ∈ M X that minimizes the functional J.

The second step consists in proving that ν * can be actually expressed as a fvarifold associated to a true fshape (X, f * ) with f * ∈ L ∞ (X). We have detailed the full proofs of Lemma 2 and Proposition 6 in Appendix B.

Existence with non-fixed X

We now consider the existence of a template when X is no more fixed and is estimated as well. In this case, as we mentioned earlier, it is enough to introduce a RKHS Hilbert space V 0 continuously embedded in C 2 0 (E, E), an initial hypertemplate X 0 and consider also an optimization of the template X in the orbit of X 0 under the action of φ 0 ∈ G 0 , the group of diffeomorphisms associated with V 0 i.e. X = φ 0 .X 0 . To prove existence result, we will need to introduce a penalization depending on the distance between X 0 and X i.e. on d G 0 (Id, φ 0 ). A typical functional would be, if X 0 is the orbit of X 0 under G 0 and d X 0 (X, X ′ ) . = inf φ∈G 0 ,φ.X=X ′ d G 0 (Id, φ) is the induced distance between two templates in the orbit,

J 1 (X, f ) = γ V 0 2 d X 0 (X 0 , X) 2 + γ f 2 X |f (x)| 2 dH d (x) + γ W 2 N i=1 µ (X i ,f i ) -µ (X,f ) 2 W ′ (57)
Since a diffeomorphism φ 0 ∈ G 0 is obtained as the flow of a time-varying vector field v 0 ∈ L 2 ([0, 1], V 0 ), we consider the minimization of the following functional:

J 2 (v 0 , f ) = γ V 0 2 v 0 2 L 2 ([0,1],V 0 ) + γ f 2 X |f (x)| 2 dH d (x) + γ W 2 N i=1 µ (X i ,f i ) -µ (X,f ) 2 W ′ for v 0 ∈ L 2 ([0, 1], V 0 ), f ∈ L 2 (X) with X = φ 0 .X 0 . However, it is more convenient for X = φ 0 .X 0 to consider the change of variable f → f 0 = f •φ 0 from L 2 (X) → L 2 (X 0
) so that we keep working in a fixed space L 2 (X 0 ). Hence we end up with functional

J 3 (v 0 , f 0 ) = γ V 0 2 v 0 2 L 2 ([0,1],V 0 ) penalization on X + γ f 2 X 0 |f 0 (x)| 2 |d x φ 0 ↾ TxX |dH d (x) = X |f (x)| 2 dH d (x) for f =f 0 •(φ 0 ) -1 ∈L 2 (X) + γ W 2 N i=1 µ (X i ,f i ) -φ 0 .µ (X 0 ,f 0 ) 2 W ′ (58)
with f 0 ∈ L 2 (X 0 ) and φ.µ ∈ W ′ denoting for a fvarifold µ its diffeormorphic transport by a C 1 diffeormorphism φ defined by (φ.µ)(ω)

. = |d x φ ↾ V |ω(φ(x), d x φ(V ), f )dµ(x, V, f ) for ω ∈ W .
The existence of a minimizer (v 0 * , f 0, * ) for J 3 gives immediately the existence of a minimizer (X * , f * )

for J 1 with X * . = φ v 0 * 1 .X 0 , f * . = f 0, * • (φ v 0 * 1 ) -1 ∈ L 2 (X *
). The existence result then becomes:

Proposition 7. Assume that W is continuously embedded in C 2 0 (E × G d (E) × R)
, that X 0 and (X i ) 1≤i≤N are finite volume and bounded d-dimensional rectifiable subsets and that

f i ∈ L 2 (X i ) for 1 ≤ i ≤ N . Assume γ V 0 > 0 and γ f /γ W is large enough with γ f , γ W ≥ 0. Then • J 1 given by equation (57) achieves its minimum on {(X, f ) | X ∈ X 0 , f ∈ L 2 (X)}; • any minimizer (X * , f * ) is such that f * ∈ L ∞ (X * ); • if X 0 is also a C p submanifold and W ֒→ C m 0 (E × G d (E) × R) with m ≥ max{p, 2}, f * ∈ C p-1 (X * ).
The proof of Proposition 6 can in fact be easily adapted to this new situation. We can consider the formulation of equation (58) with J 3 a functional on the vector field v 0 and the function f ∈ L 2 (X 0 ). With respect to v 0 , thanks to the penalization in (58), we can restrict the search of a minimum on a closed ball B 0 of given radius b in L 2 ([0, 1], V 0 ), which guarantees at the same time that the Jacobians |d x φ 0 ↾ TxX | are uniformly lower bounded. This closed ball is also compact for the weak topology in L 2 ([0, 1], V 0 ). In addition, it follows from Corollary 1 that

v 0 → N i=1 µ (X i ,f i ) -φ 0 * µ (X 0 ,f 0 ) 2 W ′ is weakly continuous on L 2 ([0, 1], V 0 ) and it is also classical that v 0 → v 0 L 2 ([0,1],V 0 )
is lower semicontinuous for the weak convergence topology. Therefore, for all fixed f 0 ∈ L 2 (X 0 ), v 0 → J 3 (v 0 , f 0 ) is weakly lower semicontinuous on L 2 ([0, 1], V 0 ). It results that we obtain existence of a minimizing vector field v 0 and, reasoning as in the previous subsection and Appendix B, one deduces easily the claim of Proposition 7.

Existence of full fshape atlases (tangential setting)

We now generalize the previous results to the existence of complete atlases of fshapes' datasets. In addition to the template (X, f ), one wants to simultaneously estimate transformations from the template to each subject. In the tangential model of Section 4.3, these consist in deformations (φ i ) i=1,..,N obtained as flows of time varying vector fields (v i ) and residual signals (ζ i ) that are L 2 functions on X. In that situation, following Section 4.3, the optimization functional for atlas estimation writes:

J(X, f, (ζ i ), (v i )) . = γ V 0 2 d X 0 (X 0 , X) 2 + γ f 2 X |f (x)| 2 dH d (x) + 1 2 N i=1 γ V v i 2 L 2 ([0,1],V ) + γ ζ X |ζ i (x)| 2 dH d (x) + γ W µ (X i ,f i ) -µ (φ v i 1 (X),(f +ζ i )•(φ v i 1 ) -1 ) 2 W ′ (59)
The main existence result is the following:

Theorem 6. Assume that W is continuously embedded in C 2 0 (E × G d (E) × R)
, that X 0 and (X i ) 1≤i≤N are finite volume bounded d-dimensional rectifiable subsets and that

f i ∈ L 2 (X i ) for 1 ≤ i ≤ N . Assume γ V 0 , γ V > 0 and γ f /γ W and γ ζ /γ W are large enough with γ f , γ W , γ ζ ≥ 0. Then • J given by equation (59) achieves its minimum on {(X, f, ζ i , (v i )) | X ∈ X 0 , f ∈ L 2 (X), ζ = (ζ i ) ∈ L 2 (X) N , (v i ) ∈ L 2 ([0, 1], V ) N }; • any minimizer (X * , f * , (ζ i * ), (v i * )) is such that f * and ζ i * for 1 ≤ i ≤ N are in L ∞ (X * ); • if X 0 is also a C p submanifold and W ֒→ C m 0 (E × G d (E) × R) with m ≥ max{p, 2}, f * and the ζ i * 's are in C p-1 (X * ).
The proof leans essentially on the same arguments as detailed in the previous subsections: we refer the reader to Appendix C. This result is of fundamental importance for the rest of the paper since it ensures that the atlas estimation problems that we shall study numerically in the next sections do have at least a solution, and we see that this holds with only L 2 regularity assumptions on the signals.

Existence in the metamorphosis framework

Interestingly, the previous existence of solutions in the simplified tangential setting can be also used to show existence of solutions to the corresponding variational problem in the metamorphosis framework given by (43). In that model, the subjects are obtained approximately as metamorphoses of the template (X, f ), which is itself a metamorphosis of an hypertemplate (X 0 , f 0 ) where X 0 is a bounded finite-volume rectifiable subset of E and f 0 ∈ L 2 (X 0 ). The important lemma that bridges both approaches is the following:

Lemma 3. For v fixed and ζ 1 ∈ L 2 (X) fixed, the infimum over h ∈ L 2 ([0, 1], L 2 (X)) of γ f 2 1 0 X |h t (x)| 2 |d x φ v t ↾ TxX |dH d (x)
under the constraint that

1 0 h t dt = ζ 1 is reached on a unique point h * ∈ L 2 ([0, 1], L 2 (X))
given by

h * t (x) . = C(x) ζ 1 (x) |d x φ v t ↾ TxX | (60)
where C(x) . = (

1 0 1 |dxφ v s ↾ TxX | ds) -1 .
Proof. The proof could be deduced from general optimization results but for sake of completeness we give here a short proof. Let us first notice that over (x, t) ∈ X × [0, 1] we have 0

< inf |d x φ v t ↾ TxX | ≤ sup |d x φ v t ↾ TxX | < +∞. Now, if we denote α t (x) . = |d x φ v t ↾ TxX | -1/2
and consider the change of variable g t (x) = h t (x)/α t (x), the problem is equivalent (after the use of Fubini-Tonelli Theorem) to the minimization on

g ∈ L 2 ([0, 1], L 2 (X)) of γ f 2 X 1 0 |g t (x)| 2 dt dH d (x)
under the constraint 1 0 α s (x)g s (x)ds = ζ 1 (x) for a.e. x ∈ X. However, this is a separable problem and for any fixed x ∈ X we have to consider on

H . = L 2 ([0, 1]) the very simple quadratic minimization problem: inf |u| 2 H under the constraint α . (x), u H = ζ 1 (x) that has the unique solution u = λ(x)α . (x) with λ(x) = ζ 1 (x)/|α . (x)| 2 H . Hence, denoting h * t (x) = ζ 1 (x) |αt(x)| 2 |α.(x)| 2 H
, we get the result. Now, introducing, as in 2.2, the Riemannian energies

E X 0 (v 0 , h 0 ) . = γ V 0 2 1 0 |v t | 2 V 0 dt + γ f 0 2 1 0 X |h 0 t | 2 (x)|d x φ v 0 t ↾ TxX |dH d (x)
and

E X (v, h) . = γ V 2 1 0 |v t | 2 V dt + γ f 2 1 0 X |h t | 2 (x)|d x φ v t ↾ TxX |dH d (x)
the atlas estimation functional becomes:

J((v 0 , h 0 ), (v i , h i )) . = E X 0 (v 0 , h 0 ) + N i=1 E X (v i , h i ) + γ W 2 µ (X i ,f i ) -µ (φ v i 1 (X),(f +ζ h i 1 )•(φ v i 1 ) -1 ) 2 W ′ (61)
where

(X, f ) = (φ v 0 1 (X 0 ), (f 0 + ζ h 0 1 ) • (φ v 0 1 ) -1
). Theorem 7. Let W be continuously embedded in C 2 0 (E × G d (E) × R), X 0 and (X i ) 1≤i≤N finite volume bounded d-dimensional rectifiable subsets and f i ∈ L 2 (X i ) for 1 ≤ i ≤ N . Assume that γ V 0 , γ V > 0 and that γ f /γ W and γ ζ /γ W are large enough with γ f , γ W , γ ζ ≥ 0. Under these assumptions,

• J defined by equation (61) achieves its minimum on

v 0 ∈ L 2 ([0, 1], V 0 ), h 0 ∈ L 2 ([0, 1], L 2 (X 0 )), (v i ) ∈ L 2 ([0, 1], V ) N and (h i ) ∈ L 2 ([0, 1], L 2 (X)) N ; • any minimizer (X * , f * , (ζ i * ), (v i * )) with f * = (f 0 + ζ h 0 * 1 ) • (φ v 0 * 1 ) -1 and ζ i * = ζ h i * 1 • (φ v i * 1 ) -1 is such that f * and ζ i * for 1 ≤ i ≤ N are in L 2 (X * ); • if X 0 is also a C p submanifold, f 0 a C p-1 function on X and W ֒→ C m 0 (E × G d (E) × R) with m ≥ max{p, 2}, then f * and the ζ i * 's are in C p-1 (X * ).
The proof can be found in appendix D.

A discrete framework for fshape

The previous sections mainly dealt with the setting of a well posed theoretical framework to work with fshapes. Here we would like to tackle the problem of the actual processing of fshapes. We believe that the description of an explicit numerical scheme is of the utmost importance for the fshape setting to be of practical utility. Our main motivations come from medical imaging where acquisition is done in digital forms. It means that data are discrete and in very high dimensions (tens of thousands points). Our goal is to demonstrate that the fshape framework may be used to handle real data and not just low dimensional examples.

Fvarifold norm of polyhedral meshes

We assume hereafter that data at hand are finite functional polyhedral meshes composed by T cells of dimension d immersed in E = R n . Most common examples are piecewise linear planar fcurves corresponding to the case n = 2 and d = 1, piecewise linear curves in space corresponding to n = 3 and d = 1, and piecewise triangular surfaces in the space corresponding to n = 3 and d = 2. A finite polyhedral fshape (X, f ) is fully described by three matrices: a first P × n matrix x contains the coordinates of the P vertices x k in E (i.e. x = (x k ) 1≤k≤P ), a second P × 1 column vector f contains the P values f k . = f (x k ) ∈ R of the signal (i.e. f = (f k ) 1≤k≤P ) and a third T × (d + 1) matrix contains the list of edges for each cell so that each row contains (d + 1) integers corresponding to the indices of the vertices of the cell.

In the sequel, it will be convenient to introduce the corresponding geometrical diffeomorphic action in the discrete setting: for a diffeomorphism φ ∈ G and a discretized fshape

(x, f ) = (x k , f k ) 1≤k≤P we denote φ.(x, f ) . = (φ.x, f ) where φ.x . = (φ(x k )) 1≤k≤P (62)
Note that the discrete action is more straightforward than its continuous version φ.(X, f ) = (φ.X, f • φ -1 ) since along a flowing particle

x k we have f • φ -1 (φ(x k )) = f (x k ) = f k (Lagrangian point of view).
Note also that the connectivity matrix (i.e. list of edges) remains unchanged through the deformation process.

There is a rather natural way to approximately represent polyhedral meshes in the fvarifold space. The method is very similar to the one described in [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF] in the case of fcurrents space. Using formula [START_REF] Trouvé | Metamorphoses through lie group action[END_REF], each cell is coded by a Dirac δ (x,V, f ) so that a discrete polyhedral fshape is viewed as a distribution of Dirac spread in the space R n . Therefore, the measure µ (X,f ) of equation [START_REF] Trouvé | Local geometry of deformable templates[END_REF] which is associated to (X, f ) is approximated by the following finite sum of Dirac

µ (X,f ) ≈ µ (x,f ) . = T ℓ=1 r ℓ δ (x ℓ ,V ℓ , fℓ ) (63)
where xℓ ∈ R n is the center of the ℓ-th polyhedral cell, V ℓ ∈ G d (R n ) gives the direction of the cell, fℓ ∈ R is the average value of the signal on this cell and finally r ℓ ∈ R is equal to the d-volume of the cell. Note that for the case of surfaces (resp. curves), V ℓ will be simply a unit normal vector (resp. unit tangent) and the orientation will be removed directly "by the kernel": the formula of the kernel k t used to compute the scalar product in the fvarifold space will be invariant with respect to changes in the orientation as in equation (65). As an illustration, let us precisely describe how to get the representation of the ℓ-th triangle T ℓ belonging to a polyhedral functional surface (X, f ). If we assume that the vertices of T ℓ (read in the connectivity matrix) are

(x k 1 , x k 2 , x k 3 ) ∈ (R 3 ) 3 for some k 1 , k 2 , k 3 = 1, • • • , P , we have            xℓ = (x k 1 +x k 2 +x k 3 ) /3, V ℓ = (x k 2 -x k 1 )∧(x k 3 -x k 1 ) / (x k 2 -x k 1 )∧(x k 3 -x k 1 ) , fℓ = (f k 1 +f k 2 +f k 3 ) /3, r ℓ = (x k 2 -x k 1 ) ∧ (x k 3 -x k 1 ) . ( 64
)
where . is the standard Euclidean norm in R 3 and ∧ correspond to the cross product. To compute the fvarifold norm of a functional surface we use in our numerical experiments Gaussian kernels satisfying assumptions of Proposition 2. Let us define:

     k e (x 1 , x2 ) = exp -x1 -x 2 2 /σ 2 e , k t (V 1 , V 2 ) = exp -2 σ 2 t (1 -V 1 , V 2 2 ) , k f ( f1 , f2 ) = exp -| f1 -f2 | 2 /σ 2 f . (65)
where σ e , σ t , σ f > 0 are scale parameters fixed by the practitioner. Thence the scalar product of two fvarifolds µ (X,f ) and µ (Y,g) associated to two polyhedral triangular functional surfaces (X, f ) and (Y, g), and approximated respectively by Tx ℓ=1 r ℓ δ xℓ ,V ℓ , fℓ and Ty ℓ ′ =1 s ℓ ′ δ ŷℓ ′ ,W ℓ ′ ,ĝ ℓ ′ , may be written:

µ (X,f ) , µ (Y,g) W ′ ≈ µ (x,f ) , µ (y,g) W ′ = Tx ℓ=1 Ty ℓ ′ =1 r ℓ s ℓ ′ k e (x ℓ , ŷℓ ′ )k t (V ℓ , W ℓ ′ )k f ( fℓ , ĝℓ ′ ). (66)
We discuss in Section 9.2 some numerical aspects of the fvarifolds.

Discrete approximations of continuous fshapes.

A continuous fshape (X c , f c ) can be approximated by a finite polyhedral fshape in the fvarifold space as follows. First, extract P points and their corresponding signal from (X c , f c ). Then, compute a mesh with the extracted points to define a polyhedral fshape (X, f ). Finally, use (63) to approximate µ (X,f ) by a finite sum µ (x,f ) of Diracs.

Although this method of approximation seems reasonable in many practical cases, we do not provide explicit conditions to ensure the convergence of µ (x,f ) toward µ (Xc,fc) when the number P of extracted vertices tends to infinity. To the best of our knowledge, this is still an open problem and a famous illustration of pathological cases is the Schwarz polyhedron. Nevertheless, there exists various theoretical results concerning the polyhedral approximation of continuous surfaces for the currents norms and the varifold norms, see e.g. [START_REF] Simon | Lecture notes on geometric measure theory[END_REF]. Moreover, in [START_REF] Thibert | Sur l'approximation géométrique d'une surface lisse[END_REF], the author studies sufficient conditions to ensure the convergence of the area of triangular meshes toward the area of the continuous surface. A general result extending the aforementioned works to our geometrico-functional framework involves some materials beyond the scope of this paper.

Shooting equations for the metamorphosis Riemannian framework

In the large deformation setting for usual shapes, it was shown (see for instance [START_REF] Miller | Geodesic Shooting for Computational Anatomy[END_REF]) that the dynamic of optimal vector fields giving geodesics in groups of diffeomorphisms can be described through Hamiltonian systems of equations, called the forward equations. In the discrete situation of a finite set of points, these actually reduce to a coupled evolution of the position of the particles and extra variables called the momenta attached to every particle. In addition, the Hamiltonian structure implies that all geodesic trajectories are eventually parametrized only by the initial positions and momenta, which is the principle underlying geodesic shooting algorithms for shape matching (see [START_REF] Arguillere | Shape deformation analysis from the optimal control viewpoint[END_REF] for an extensive presentation of the Hamiltonian setting for shape deformation analysis).

We now describe how to obtain similar shooting frameworks in our more general situation of fshapes, first for the fshape metamorphosis model. As discussed above, any element (X, f ) in F will be discretized as a family (x, f ) = (x k , f k ) 1≤k≤P of points x with signal value f . In that setting, for a proper weighting matrix (the precise definition of this matrix depends on actual choices for the approximation of the L 2 norm), we can discretize the L 2 (X) dot product as (D(x)h|h) for any h = (h k ) ∈ R P . The continuous problem (3) is then approximated by a simple discrete control problem

min γ V 2 1 0 |v t | 2 V dt + γ f 2 1 0 (D(x t )h t |h t )dt (67)
for fixed end point conditions (x t , f t ) ↾ t=0 and (x t , f t ) ↾ t=1 and controlled dynamic in V × R P given by ẋt

= v t .x t ḟ t = h t , (68)
where v.x = (v(x k )) 1≤k≤P .

Remark 1. An important point is to note that the discrete evolution is based on a Lagrangian particle based representation and to understand where the equations (68) are coming from. In the discrete setting, t → x t = (x t,k ) 1≤k≤P is the evolution the vertices of a polyhedral mesh i.e. x t,k = φ v t (x k ) where φ v t is the flow of the time dependent vector field t → v t and x = (x k ) are the vertices on the initial mesh X so that ẋt,k = v t (φ v t (x k )) = v t (x t,k ) which gives the first equation of (68). Moreover, f t = (f t,k ) 1≤k≤P are the signal values attached to X t at the vertices positions x t = (x t,k ) 1≤k≤P so that f t,k is the discretization of the continuous signal

f t : X t → R in Lagrangian coordinates i.e. f t,k = f t (φ v t (x k )) = f t (x t,k ). In particular, one have ḟt,k = dft(φ v t (x k )) dt = d t 0 hs(x k )ds dt = h t (x k ) = h t,k where h t = (h t,k ) = (h t (x k ))
is the discretization of h t : X → R on the vertices (x k ) which gives the second equation of (68).

In that situation, introducing co-states or momenta (p, p f ) = (p k , p f k ) 1≤k≤P we get the associated Hamiltonian

H((x, f ), (p, p f ), (v, h)) = (p|v.x) + (p f |h) - γ V 2 |v| 2 V - γ f 2 (D(x)h|h)
so that we deduce from the Pontryagin's Maximum Principle that the optimal controls satisfy

v(•) = 1 γ V P k=1 K V (•, x k )p k and h = 1 γ f D -1 (x)p f (69)
where K V is the kernel associated with the RKHS space V . Then, plugging the optimal control into the Hamiltonian H we get the reduced Hamiltonian given as

H r ((x, f ), (p, p f )) = 1 2γ V (K x,x p|p) + 1 2γ f (D -1 (x)p f |p f ) (70) with K x,x . = (K V (x k , x k ′ )) 1≤k,k ′ ≤P .
Note that the weighting matrix D(x) may be chosen to be diagonal ("mass lumping") so that the computation of D -1 (x) is straightforward. From H r , we can derive the forward equation given by the Hamiltonian dynamic:

ẋ = ∂ p H r (x, p) ṗ = -∂ x H r (x, p) . ( 71
)
As usual, since the Hamiltonian is not depending on time, it is a conserved quantity during the geodesic evolution so that we get from (69) and (70) that

γ V 2 1 0 |v t | 2 V dt + γ f 2 1 0 (D(x t )h t |h t )dt = H r ((x, f ), (p, p f )) ↾ t=0 . (72)
We see also that the reduced Hamiltonian is a perturbation of the more familiar Hamiltonian for the pure geometrical case mentioned earlier,

H geo r (x, p) = 1 2γ V (K x,x p|p) = γ V 2 |v| 2 V with a geometrico-functional term H geo-f un r (x, f , p f ) . = 1 2γ f (D -1 (x)p f |p f ) . ( 73
)
Hence the shooting equations contains a new source term which is ∂ x H geo-f un r in the dynamic of the momenta p:

ṗ = -∂ x H geo r (x, p) -∂ x H geo-f un r (x, f , p f ) new source term (74)
and since the H geo-f un r does not depend on f , we have

ṗf = -∂ f H geo-f un r (x, f , p f ) = 0 . (75)
In particular, the functional speed ḟk is only modulated by the evolution of the local weight D(x) k,k which depends in the continuous limit both on the divergence of the "tangential" part of v to the manifold spanned by the x k 's and on its normal part when the mean curvature is non vanishing.

Shooting equations for the tangential setting

In the tangential model of section 2.5, equations are simplified once again. Indeed, the L 2 metric being frozen to the initial position of X, the second term in (67) becomes γ f 2 1 0 (D(x 0 )h t |h t )dt where D(x 0 ) is now a constant weighting matrix and the optimal functions h t for fixed time 1 end point condition f 1 are simply given by h

t ≡ ζ = f 1 -f , i.e. f t = f + tζ = tf 1 + (1 -t)f where ζ = (ζ(x k ))
1≤k≤P and ζ is the residual introduced in Section 4.3. In that situation, the new source term disappears in (74) and the dynamic of (x, p) is described by the usual Hamiltonian equations in the purely geometrical case, which gives:

   ẋ = ∂ p H geo r (x, p) ṗ = -∂ x H geo r (x, p) f t = f + tζ (76)
These are the forward equations in the tangential model. Compared to metamorphoses, the dynamic on the signal part can be expressed in closed form while the evolution in spatial position is described by the usual Hamiltonian system without extra term. This makes it particularly simple for implementation, which shall be exploited in Section 7. Optimizing over trajectories can be reduced to optimization with respect to initial momentum p 0 and signal f 1 (or the residual

ζ) since again v(•) = 1 γ V P k=1 K V (•, x k )p k , h = ζ and γ V 2 1 0 |v t | 2 V dt + γ f 2 1 0 (D(x 0 )h t |h t )dt = H geo r ((x, p)) ↾ t=0 + γ f 2 (D(x 0 )ζ|ζ) . (77)
Variation of functionals with respect to p 0 can be obtained by the usual backward integration of the adjoint Hamiltonian system (c.f. [START_REF] Charon | Analysis of geometric and functional shapes with extensions of currents. Application to registration and atlas estimation[END_REF] Chapter 1 or [3] Section 4).

Algorithms to compute mean template of fshapes

We present two different algorithms to compute a mean template from a sample of N discrete fshapes. Both methods consist in solving a variational problem via an adaptive gradient descent algorithm presented in Section 7.3. As the optimization method is the same in both cases, we just need to describe how to compute the functional and its gradient. These methods are implemented in Matlab, Cuda and C and some numerical experiments are shown Section 8 and 9.

Hypertemplate and tangential setting

Algorithms 1, 2 and 5 give a way to implement the atlas estimation in the tangential setting described in Section 4.3 when f 0 = 0. The principle of the method is illustrated by Figure 3. The method consists in computing a reduced version based on geodesic shooting and initial momenta of the functional J tan 0 of equation ( 51) with Algorithm 1 and its gradient with Algorithm 2 and then plugging this in the optimization box given by Algorithm 5.

(x 0 , 0) (x, f ) (x 1 , f 1 ) (x 2 , f 2 ) (x 2 , f 2 ) (x 3 , f 3 ) (x 4 , f 4 ) p 0 p 2
Figure 3: The hypertemplate method: the estimated mean template (x, f ) = (φ v p 0 1 (x 0 ), f ) is a deformed version of the hypertemplate (x 0 , 0).

More precisely, the inputs are N observations (x i , f i ), a hypertemplate (x 0 , f 0 ) with f 0 = 0, the momenta p 0 and p 1 , • • • , p N and the functionals f and ζ 1 , • • • , ζ N . All the fshapes should be provided with their corresponding connectivity matrix. The optimization is conducted on the reduced version of J tan 0 defined as

J tan 0 (p 0 , f , (p i , f i ) 1≤i≤N ) . = γ V 0 2 |v p 0 | 2 V 0 + γ 0 f 2 |f | 2 x + N i=1 γ V 2 |v p i | 2 V + γ f 2 |ζ i | 2 x + γ W 2 g((x i , f i ), (x i , f i )) (78) in the variables            (p 0 , f ) ∈ E P × R P (p i , ζ i ) ∈ E P × R P with x = φ v p 0 1 .x 0 , xi = φ v p i 1 .x, f i k = f + ζ i (79)
where |l| 2 x = (D(x)l|l) is the discretization of the L 2 (X) norm for l = (l k ) ∈ R P defined on the vertices x ∈ E P as introduced in Section 6.3.

The hypertemplate may be difficult to define and the choice of this initial guess may obviously affect the quality of the atlas estimation. Nevertheless, the outputs of the algorithms are found to be stable even for the non-trivial experiments shown in Figure 8. We recommend generating a simple and smooth fshape with a meshing program. This hypertemplate should be topologically equivalent to the observations as much as possible. Using one of the observations as the hypertemplate could be tempting but it may induce some bias in the estimated atlas as the final mean template may contain the same specific features as the chosen observation (recall that the template is a diffeomorphic deformation of the hypertemplate).

At the end of the minimization procedure, the outputs are: an estimation of the mean template (x, f ) = (φ v p 0 1 .x 0 , f ), the momenta (p i ) 1≤i≤N and the functional residuals

(ζ i ) 1≤i≤N so that (x i , f i ) = (φ v p i 1 .x, f + ζ i ) is close to (x i , f i ) for any i = 1, • • • , N .
Algorithm 1 Computation of the energy J tan 0 with formula (78)

1: Require: A hypertemplate (x 0 , f 0 = 0) and N fshapes (x i , f i ). Deformation: compute x by forward integration of (x 0 , p 0 ).

5:

for i = 1 to N do 6:

Deformation: compute xi by forward integration of (x, p i ); compute the signal

f i ← f + ζ i . 7:
Fvarifold norm: compute the fvarifold representation µ (x i , f i ) and µ (x i ,f i ) ; compute

g i ← µ (x i , f i ) -µ (x i ,f i ) 2 W ′ . 8:
Penalty terms: compute

|v p i | 2
V and |ζ i | 2

x .

9:

end for 10:

Penalty terms: compute

|v p 0 | 2 V 0 and |f | 2
x . 11: End 12: Outputs:

γ V 0 2 |v p 0 | 2 V 0 + γ 0 f 2 |f | 2 x + i γ f 2 |ζ i | 2 x + γ V 2 |v p i | 2 V + γ W 2 g i .
Algorithm 2 Computation of the gradient ∇J tan 0 1: Require: A hypertemplate (x 0 f 0 ) and N fshapes (x i , f i ). Deformation: compute x 0 by forward integration of (x 0 , p 0 ).

5:

for i = 1 to N do 6:

Deformation: compute x by forward integration of (x, p i ); compute the signal

f i ← f + ζ i . 7:
Gradient of g i wrt f , ζ i , p i and x: compute directly ∇ f g i and ∇ ζ i g i ; compute (∇ x g i , ∇ p i g i ) by backward integration of (∇ xi g i , 0).

8:

Gradient of penalty terms: compute directly

∇ x |v p i | 2 V , ∇ p i |v p i | 2 V , ∇ x |ζ i | 2 x and ∇ ζ i |ζ i | 2
x .

9:

end for 10:

Gradient of J tan 0 wrt p 0 : compute directly ∇ p 0 |v p 0 | 2 V 0 ; compute ∇ p 0 i ( γ f 2 |ζ i | 2 x + γ V 2 |v p i | 2 V + γ W 2 g i ) by backward integration of i ∇ x ( γ f 2 |ζ i | 2 x + γ V 2 |v p i | 2 V ) + γ W 2 g i , 0 . 11:
Gradient of penalty term: compute

∇ f |f | 2 x . 12: End 13: Outputs: ∇ p 0 J tan 0 = γ V 0 2 ∇ p 0 |v p 0 | 2 V 0 + ∇ p 0 i ( γ f 2 |ζ i | 2 x + γ V 2 |v p i | 2 V + γ W 2 g i ); ∇ f J tan 0 = γ 0 f 2 ∇ f |f | 2 x + γ W 2 ∇ f i g i ; ∇ p i J tan 0 = γ V 2 ∇ p i |v p i | 2 V + γ W 2 ∇ p i g i for i = 1, ..N ; ∇ ζ i J tan 0 = γ f 2 ∇ ζ i |ζ i | 2 x + γ W 2 ∇ ζ i g i for i = 1, ..N .

"Free" mean fshape and tangential model

Algorithms 3, 4 and 5 use a more direct approach to compute a mean fshape compared to the hypertemplate method of Section 7.1. The general framework is similar but the mean template is no longer defined as a deformed version of an hypertemplate. The coordinates of the points composing the mean template will be directly updated along the gradient descent. Notice that there is no extra cost as, in Algorithm 2, we compute the gradient ∇ x g i in the backward integration step (see also Section 9.2). The aim of this method is simply to use this gradient directly to update the mean fshape thereby making the evolution "free" along the optimization process. Figure 4 gives an illustration of the method.

(x init , f init ) (x, f ) (x 1 , f 1 ) (x 2 , f 2 ) (x 2 , f 2 ) (x 3 , f 3 ) (x 4 , f 4 ) (x t , f t ) p 2
Figure 4: The "free" mean template method: the estimated mean template (x, f ) is a gradientdescent-based update of an initial fshape (x init , f init ). The notation (x t , f t ) symbolizes the state of the mean template at iteration t.

The free mean fshape method is a gradient descent on the following functional:

J tan f ree ((x, f ), (p i , f i ) 1≤i≤N ) . = γ 0 f 2 |f | 2 x + N i=1 γ V 2 |v p i | 2 V + γ f 2 |ζ i | 2 x + γ W 2 g((x i , f i ), (x i , f i )) (80) in the variables            (x, f ) ∈ E P × R P , (p i , ζ i ) ∈ E P × R P with xi = φ v p i 1 .x, f i = f + ζ i (81)
This algorithm does not have yet the same theoretical warranties as those given by Theorem 6 for the hypertemplate version. We are not able to show, with the same arguments, that there exists proper minimizers of the continuous counterpart of (80) as the fshape (X, f ) is not constrained to belong to a single fshape bundle. However, we use in practice a regularized gradient descent to keep the evolution of the template smooth along the minimization. Therefore, if the size of the updates are small enough, the output of the free mean fshape algorithm may be considered as a discretization of a smooth deformation of (X init , f init ) as Property 4.2.1 of [START_REF] Charon | Analysis of geometric and functional shapes with extensions of currents. Application to registration and atlas estimation[END_REF] shows. This trick also prevents odd behaviours of the algorithm as discussed in Section 9.2. Compared to the hypertemplate algorithm of previous section, this numerical scheme still reduces the constraint imposed on the template evolution which is likely to provide less dependency in the choice of initialization.

The inputs are N observations (x i , f i ), an initial fshape (x init , f init ), the momenta p At the end of the minimization procedure, the outputs are: an estimation of the mean template (x, f ), the momenta (p i ) 1≤i≤N and the functional residuals

(ζ i ) 1≤i≤N so that (x i , f i ) = (φ v p i 1 .x, f + ζ i ) is close to (x i , f i ) for any i = 1, • • • , N .
Algorithm 3 Computation of the energy J tan f ree of formula (80) for i = 1 to N do 5:

Deformation: compute xi by forward integration of (x, p i ); compute the signal

f i ← f + ζ i . 6:
Fvarifold norm: compute the fvarifold representation µ (x i , f i ) and µ (x i ,f i ) ; compute

g i ← µ (x i , f i ) -µ (x i ,f i ) 2 W ′ . 7:
Penalty terms: compute

|v p i | 2
V and |ζ i | 2

x .

8:

end for 9:

Penalty term: compute |f | 2 x . 10: End 11: Outputs:

γ 0 f 2 |f | 2 x + i γ V 2 |v p i | 2 V + γ f 2 |ζ i | 2 x + γ W 2 g i .
Algorithm 4 Computation of the gradient ∇J tan f ree 1: Require: N fshapes (x i , f i ). for i = 1 to N do 5:

Deformation: compute x i by forward integration of (x, p i ); compute the signal

f i ← f + ζ i 6:
Gradient of g i wrt f , ζ i , p i and x: Compute directly ∇ f g i and ∇ ζ i g i ; compute (∇ x g i , ∇ p i g i ) by backward integration of (∇ xi g i , 0).

7:

Gradient of penalty terms: compute directly

∇ x |v p i | 2 V , ∇ p i |v p i | 2 V , ∇ x |ζ i | 2 x and ∇ ζ i |ζ i | 2
x .

8:

end for

9:

Gradient of penalty term: compute directly

∇ f |f | 2 x and ∇ x |f | 2 x . 10: End 11: Outputs: ∇ x J tan f ree = γ 0 f 2 ∇ x |f | 2 x + i ∇ x ( γ f 2 |ζ i | 2 x + γ V 2 |v p i | 2 V + γ W 2 g i ); ∇ f J tan f ree = γ 0 f 2 ∇ f |f | 2 x + γ W 2 ∇ f i g i ; ∇ p i J tan f ree = γ V 2 ∇ p i |v p i | 2 V + γ W 2 ∇ p i g i for i = 1, • • • , N ; ∇ ζ i J tan f ree = γ f 2 ∇ ζ i |ζ i | 2 x + γ W 2 ∇ ζ i g i for i = 1, • • • , N .

Optimization

An atlas estimation is a smooth, non-convex and high dimensional global optimization problem.

We use an adaptive gradient descent algorithm to solve it numerically. This method is popular in geometric atlas estimation as it is quite robust and performs relatively well on various real examples, see [START_REF] Glaunès | Diffeomorphic matching of distributions: A new approach for unlabelled point-sets and sub-manifolds matching[END_REF][START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF]. In our geometrico-functional framework, one main issue comes from the fact that the functionals to minimize depend on various types of variables living in different spaces: momenta, functional values or points coordinates. Adding the functional part in the optimization process makes it harder, compared to a pure geometrical approach, as some non desirable phenomenon may hold. For a thorough discussion, please see Section 9.1.

Let us assume hereafter that we want to minimize a functional J (u 1 , • • • , u K ) depending on K ≥ 1 types of variables. For instance, in an atlas estimation of fshapes with a hypertemplate (see Section 7.1), we have J = J tan 0 and K = 4 (namely the variables f , p, f i and p i )).

If the K types of variables are non homogeneous, we note that in practice, the different types of variables should not be updated at the same speed along the optimization process. It means that the gradient is not the proper descent direction to follow in order to reach a reasonable solution. As Algorithm 5 shows, we use K different steps sizes denoted (δ i ) 1≤i≤K and we adapt them separately. The quality of the final result will depend on a fine balance between the δ i 's.

Algorithm 5 Adaptive gradient descent 1: Require: Coefficients used to adapt the step sizes 0 < s -< 1 < s + . 2: Inputs: A functional J and its gradient ∇J = (∇

u i J ) K i=1 . Points (u init 1 , • • • , u init K ) and step sizes (δ init 1 , • • • , δ init K ). 3: Begin 4: Initialize : (u 1 , • • • , u K ) ← (u init 1 , • • • , u init K ) and (δ 1 , • • • , δ K ) ← (δ init 1 , • • • , δ init K ).
5:

repeat 6: Compute J cur ← J (u 1 , • • • , u K ) and ∇ u i J cur ← ∇ u i J (u 1 , • • • , u K ). 7:
Update all the variables simultaneously : compute

J new ← J (u new 1 , • • • , u new K ) where u new i ← u i + δ i ∇ u i J cur for i = 1, • • • , K. 8: Adapt steps : if J new < J cur then s ← s + else s ← s -end if 9:
evolution algorithm of Section 7.2 after the dataset is preregistered with respect to translations. The template is initialized with a prototype fshape with zero-valued signal everywhere in order to avoid as much as possible bias toward one particular individual. Figure 5 shows the current template at several intermediate steps of the estimation. We use all Gaussians for the kernels k e , k t and k f defining the data attachment terms as in equation ( 65). The template is obtained by refining twice the scales σ e and σ f throughout the process in order to have a coarse estimation of the atlas in the first place before being able to retrieve finer details from the dataset. The resulting template captures the average shape and signal behavior of the set of subjects: in particular, the algorithm tends to recover the signal patterns that are the most shared among the population, as the uniform coloring of the head for instance.

Figure 5: Template evolution along the gradient descent steps of the algorithm. On the left is the prototype fshape which serves as the initialization for the template. We use a multiscale approach that decreases, after a certain number of iterations, the characteristic sizes of the kernels on geometry and signal.

In addition to the template, we obtain, at convergence, a set of residual signals as well as deformation momenta that map the template on the different subjects both in terms of shape and texture. The kernel for the space V of deformation fields is taken as a sum of two Gaussian kernels, following the approach of [START_REF] Bruveris | Mixture of Kernels and Iterated Semidirect Product of Diffeomorphisms Groups[END_REF] which allows to introduce multiscale deformations. We show two of such mappings in Figure 6. The whole set of shape and texture matchings compared to the original subjects is finally summed up in Figure 7. In terms of numerics, the subjects of the dataset are unequally sampled between 2800 and 6500 points while the template of Figure 5 has 7000 points. Using our GPU implementation for kernel computations, each iteration of the algorithm takes approximately 38 secs for a total atlas estimation time of 3 hours (300 iterations) on a server equipped with a Nvidia GTX 555.

The dependency in the initialization of the template is also an important issue. Section 9.3 explains how to make the previous procedures independent of rescaling or sampling changes of the fshapes. But obviously, in both the hypertemplate and 'free' template evolution algorithm, the shape of the template is still constrained to live in the diffeomorphic orbit of the initial one. It results that one cannot expect to remove completely the bias resulting from the initialization's choice. In Figure 8, we show the template obtained after convergence for different initializations. Although these results do demonstrate some variations for the estimated template, they still show a quite remarkable stability in most of the important geometric and functional features of the dataset.

Initialization Estimated template

Figure 8: Variations of the final estimated template with respect to the initialization, respectively (from left to right) with the same prototype shape as in figure 5 but with constant 0.5 signal, with the shape of subject 4 and with subject 2.

OCT dataset

We now present a template estimation performed on a real dataset. The dataset contains surfaces segmented from volumetric images of the retina acquired by optical coherence tomography (OCT) as described in [START_REF] Lee | Comparative analysis of repeatability of manual and automated choroidal thickness measurements in nonneovascular age-related macular degeneration[END_REF][START_REF] Lee | Optic nerve head and peripapillary morphometrics in myopic glaucoma[END_REF]. The aim of these measurements is to detect early glaucoma by analysing changes of conformations of the inner layers of the retina around the optic disc (where the optic nerve slots into the retina). Data are fshapes: the surfaces represent the lower boundary of the nerves fiber layer (NFL) and the signals represent the thickness of the NFL. A loss of thickness may be an indicator of glaucoma. We depict below two typical observations in two different views: Data 1 is a normal subject in Figures 10b-10c and Data 2 is a glaucoma subject in Figure 11b-11c. The "view 1" has the same scale for the 3 axis and "view 2" is a flipped version of "view 1" where the depth axis has been scaled (×4) to better represent fine reliefs. The typical size of square boundary ranges from 5 to 7mm and the signal ranges from 0 to 0.3mm. The color scale for the signal is the same for all the pictures of this Section. The overall geometry of the NFL boundary does not seem to be challenging but these surfaces contains many boundaries making the problem hard to handle in practice. We discuss in Section 9.2 how we manage problems arising from boundary effects. Moreover, the difficulty is increased by the fact that some observations are not centered and the opening may be closed to the boundary of the acquisition area as in Figure 11. We ran our code, using the hypertemplate method described Section 7.1 on a dataset containing 51 observations (19 normal, 25 glaucomatous and 7 suspects). The observations were preregistered by hand with respect to translations so that the center of the opening is at the origin. We use downsampled data (raw data contains more than 130000 points and 270000 triangles each) and each observation contains about 5000 points and 9850 triangles. The hypertemplate is the flat rectangle with a hole and with a null signal depicted Figure 9a. It contains 5700 points and 11100 triangles. Computations of the mean template (x, f ) (Figure 9b and9c) and the deformations took 7 hours (120 iterations) using a server equipped with a Nvidia GTX 555 graphical processor unit. The estimated template is depicted Figures 9b and9c. The average shape is a bended version of the flat hypertemplate which is a reasonable guess. The baseline signal f captures a light but typical "C"-shaped signal. This blurred signal is due to the rather high variability of signals across the dataset. The results of the deformations and functional residuals estimations for Data 1 and 2 are given by Figures 10 and11 respectively. The functional part of all the 51 observations is well reconstructed and we are now able to compare these signals as they are all defined on the mean template, see Figure 10d and11d. The deformations are also satisfying although most of the energy of the deformations is spent to match the (non informative) outer boundaries as data are misaligned, compare Figures 10a and11a

Numerical pitfalls

In this section, we discuss some issues affecting the quality of the mean template estimation. We illustrate these problems by numerical experiments and propose some solutions to fix them.

Mass cancellation

The varifold norm was first introduced in the pure geometrical setting to avoid mass cancellation appearing with the classical current norms. As orientation matters in the currents' setting, a surface with positive area sufficiently crumpled may have an arbitrarily small current norm. This phenomenon is common in practical applications and particularly during an atlas estimation. This is described in detail in the introduction of [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF] or in the Chapter 3 of [START_REF] Charon | Analysis of geometric and functional shapes with extensions of currents. Application to registration and atlas estimation[END_REF]. In the (pure geometrical) varifold setting, the norm of a surface cannot decrease too much when two pieces of surface are folded: Theorem 3.4.1 in [START_REF] Charon | Analysis of geometric and functional shapes with extensions of currents. Application to registration and atlas estimation[END_REF] is in some sense a reciprocal inequality to Proposition 4. Unfortunately, there is no such a control for the fvarifold norm as one can exhibit a functional surface of positive area and non-zero signal with small fvarifold norm.

We follow here the notations introduced in Section 6. Recall that, in formula (66), the inner product of two fvarifolds was approximated by a double sum of inner products of Diracs. We assume hereafter that the kernel fvarifold inner product satisfies the assumptions of Proposition 2. Therefore, by formula (36) we have, for any Diracs δ (x 1 ,V 1 ,f 1 ) and δ (x 2 ,V 2 ,f 2 ) :

δ (x 1 ,V 1 ,f 1 ) , δ (x 2 ,V 2 ,f 2 ) W ′ ≤ |k e (•, •)| ∞ |k t (•, •)| ∞ k f (f 1 , f 2 ) . ( 82
)
The right hand side of (82) may be arbitrarily small when |f 1 -f 2 | is large as k f (f 1 , •) is continuous and vanishes at infinity. It means, in particular, that sufficiently high variations in the signal make any two Diracs orthogonal with respect to the fvarifold inner product. Now, let (x, f ) be a finite polygonal mesh approximated by a finite sum of Diracs T ℓ=1 r ℓ δ (x ℓ ,V ℓ , fℓ ) . Equations ( 37) and (82) yield

µ (x,f ) 2 W ′ ≤ |k e (•, •)| ∞ |k t (•, •)| ∞ T ℓ=1 r 2 ℓ k f ( fℓ , fℓ ) + T ℓ,ℓ ′ =1,ℓ =ℓ ′ r ℓ r ℓ ′ k f ( fℓ , fℓ ′ ) . (83)
From here on, we consider the case where k e , k t and k f are the Gaussian kernels of equation ( 65). Let 0 < ε < 1 and assume that fℓ = Kσ f ℓ for any ℓ = 1, . . . , T and K = √ ln ε. Inequality (83) becomes, in that case :

µ (x,f ) 2 W ′ ≤ T ℓ=1 r 2 ℓ + T ℓ,ℓ ′ =1,ℓ =ℓ ′ r ℓ r ℓ ′ exp(-K 2 |ℓ -ℓ ′ | 2 ) ≤ C 1 (T ) + εC 2 (T ),
where C 1 (T ) = T ℓ=1 r 2 ℓ and C 2 (T ) = T ℓ,ℓ ′ =1,ℓ =ℓ ′ r ℓ r ℓ ′ . Thence, if T remains fixed and K → ∞ the fvarifold norm of (x, f ) is less of equal to C 1 (T ). Assuming that the fshape (x, f ) is regularly discretized with P points, we may consider that r ℓ ≈ P -1 and C 1 (T ) ≈ P -1 . When the discretization of a fshape becomes finer, the number T of triangles increases and the fvarifold norm of (x, f ) may be arbitrarily small.

To illustrate this, we consider the discrete version of the variational problem of Proposition 6 in the simple case where N = 1 and γ W = 2. It gives:

       min f ∈R Px J f un (f )
where Figure 12: The source and target fshapes are respectively a big and a small flat square both lying in the yz-plane. The initial source signal is 0 and the target signal is 1. Figure 12c simply depicts the two fshapes together.

J f un (f ) = γ f 2 |f | 2 x + µ (x,f ) -µ (y,g) 2 W ′ (
In Figure 13, we show three outputs of an adaptive gradient descent on J f un with various choices of penalty parameter γ f . We denote f * ∈ R Px the solution found after 800 iterations. Figure 13a shows the result with no penalty term, i.e γ f = 0. The signal f * is almost equals to 1 on the overlapping square whereas f * oscillates dramatically on the complementary part to decrease the fvarifold norm of the remaining part of x. Figure 13b shows that with a small penalty, we are able to recover a signal close to 1 on the central square but oscillations are still present though weaker compared to Figure 13a. Finally, Figure 13c shows the solution f * found with a larger γ f . The oscillations have almost disappeared and the signal remains small in the non-overlapping part (in blue). The price to pay is a lower intensity in the central square (around 0.25). Note that this example is an extreme case but such a problem arises in real datasets. In medical imaging, it is common to get functional data acquired on a subpart of a surface only. Even with a rigid registration during the preprocessing, subparts with a well defined signal may be highly non-overlapping on two different subjects. This may yield to oscillating solutions if the gradient descent is used naively. We note that in practice, regularizing the gradient with respect to f may fix this issue in some cases. The method is similar to the one described at the end of Section 9.2 for the geometric part of the gradient.

Boundary problems

In this section, we discuss some problems arising when the "free" mean template method of Section 7.2 is used with fshapes containing boundaries. An issue concerns the gradient (with respect to x) of the data attachment terms g i (see Algorithm 3) -note that the following discussion is relevant both for the purely geometric and functional varifolds or currents framework. The values of the signal has no influence here and we assume without loss of generality that the signals are 0. Typical examples of data with boundaries (inspired by the dataset presented in Section 8.2) are depicted in Figure 14. To illustrate the problem that one can face, let E = R 3 and (y, 0) ∈ E Py × R Py be a fixed discrete target fshape as in Figure 14b. We consider the following minimization problem,

       min x∈E P J geo (x)
where

J geo (x) = µ (x,0) -µ (y,0) 2 W ′ , (85)
Our goal is to study the behaviour of a gradient descent in the variable x = (x k ) P k=1 ∈ E P on J geo starting from an initial discrete fshape (x init , 0). In numerical experiments, the gradient descent on J geo does not perform well if the fshape (x init , 0) contains boundaries as in Figure 14a. The boundary of a discrete fshape is the set of vertices that are an end of an edge belonging to a single triangle (see Figure 15). The main issue on the gradient of J geo is the following. The norm of ∂ ∂x k J geo has different orders of magnitude depending on the location of the points x k in the fshape: gradient of boundary points may be much larger than gradient of interior points, see Figure 16a. These unbalanced values between interior points and boundary points in the gradient induce undesirable effects such as self-crossings or changes in topology during the gradient descent as in Figure 16b. Note also that the situation becomes worse when the number P of points in x increases. In order to explain the singularities of ∇ x J geo at the boundary points, notice first that this is sufficient to focus on the gradient of the fvarifold inner product as we have ∇ x µ (x,f )µ (y,g) 2

W ′ = ∇ x µ (x,f ) , µ (x,f ) W ′ -2∇ x µ (x,f ) , µ (y,g) W ′ .
Nevertheless, it is hard to clearly understand this phenomenon by looking at the discrete formulas directly. It is then convenient to use the results of Section 3.3.2 whose aim is to analyse the variations of the fvarifold inner product in the continuous setting. Thence, we may consider that (x, 0) (resp. (y, 0)) is sampled from a continuous fshape (X, 0) (resp. (Y, 0)) and that µ (x,0) , µ (y,0 W ′ approximates its continuous counterpart µ (X,0) , µ (Y,0) W ′ . We further assume that (x, 0) is regularly discretized, so that the typical size of the interior volume elements r ℓ is of order P -1 and the typical size of the boundary volume elements is of order P -d-1 d . Formula (40) gives an expression of the variations of µ (X,0) , µ (Y,0) W ′ when (X, 0) moves. The right hand side of (40) is composed by two separated integral terms whose domains of integration are X and its boundary ∂X respectively. This formula then makes clear the differences between interior points and boundary points. A discrete approximation of the integral over the interior of X should then involve terms of order P -1 (i.e. the volume of the surface element), whereas a discrete approximation of the integral over the boundary ∂X should involve terms of order P -d-1 d (corresponding to the size of the volume element of the boundary). This heuristic explains why there exists a multiplicative factor between the size of the gradient at interior points and boundary points. We also know that this factor is of order P 1 d meaning that the difference increases when the resolution of the discrete fshape (x, 0) increases.

Thus the gradient ∇ x J geo computed thanks to the Euclidean metric involves singular terms on the boundary. A convenient way to address this problem is to compute a gradient via a more regular metric. Let us introduce a (conveniently normalized) Gaussian kernel

K reg : R n ×R n → R of scale σ reg . Let K reg x,x . = K reg (x k , x k ) 1≤k,k ′ ≤P ∈ R P ×P and define ∇x J geo (x) . = K reg x,x ∇ x J geo (x)
This new gradient is, in some sense, a convolution of the singular gradient with a Gaussian kernel, see Figure 17a. As the kernel K reg is positive definite, -∇J geo still defines a descent direction for J geo . Figure 17b shows the output of a regularized gradient descent. This new gradient can be seen as the discretization of a dense vector field belonging to the RKHS V reg of kernel K reg . The displacement of the template from its initial position in the regularized gradient descent can be then considered as the integration of a sequence of vector fields living in V reg . It is shown in [START_REF] Charon | Analysis of geometric and functional shapes with extensions of currents. Application to registration and atlas estimation[END_REF] (Property 4.2.1) that the regularized gradient flow is well defined. Thence, provided that the steps along the gradient descent are sufficiently small, the evolution of the points of x may be considered as a diffeomorphic evolution (belonging to the orbit G Vreg .(x init , 0)). Note that, in general, σ reg can be chosen to be different from the scale σ V of the deformations. This allows more flexibility in the template evolution during the gradient descent and justifies the methodology used in the algorithm described in Section 7.2.

Tuning the parameters

The procedure of atlas estimation for fshapes contains numbers of parameters which have to be tuned by the user. Although there exists some heuristics to choose the values of these parameters, it may be hard to find a good range of settings for a particular dataset. One of our goals is to provide an algorithm behaving in a similar way whatever the scale and the resolution of the data are. In particular, the user should be allowed to tune the parameters at coarse resolution (when computation times are low) and perform a matching or an atlas estimation at high resolution with only few updates for the parameters values. Even if these refinements are not fundamental in theory, it makes the algorithm usable in practice and it has the great advantage of normalizing the choice of the balance parameters γ 0 , γ 0 f , γ V , γ f , γ W between the various penalization terms and the data attachment term, as well as the initial step sizes in Algorithm 5 (these quantities may vary considerably otherwise).

In this section, J denotes either the functional J tan 0 of equation (78) (Algorithms 1 and 2) or J tan f ree of equation (80) (Algorithms 3 and 4).

Scale invariance

We first describe how we normalize the functional J and its gradient in order to provide an algorithm with outputs that are invariant if a scaling is applied to the geometry or to the signal of the fshapes. To do so, we carefully examine the effect of a scaling (in the geometrical or in the signal space) on each term composing J . Without this normalization step, the balance between the different terms of J would be modified by scale which is not desirable in practice.

Let us examine the behaviour of the functional when we apply a scaling on the geometrical space x → x ′ = λ e x ∈ E P and on the signal f → f ′ = λ f f ∈ R P for some λ e , λ f > 0. We assume that the scale parameters σ e and σ f of the Gaussian kernels k e and k f are also rescaled as they measure the size at which the fshapes are compared. These new kernels are denoted k ′ e and k ′ t and they verify k ′ e (x ′ , y ′ ) = k e (x, y) and k ′ t (f ′ , g ′ ) = k t (f , g). The scaled local volume elements r ′ ℓ now satisfy r ′ ℓ = λ d e r ℓ and the scaled momenta p ′ = λ e p so that the displacement generated by p ′ is also scaled. All the terms composing J are normalized in the same way and we just detail the case of the fvarifold data attachment term g i . By formula (66), we easily see that the scaled data attachment term is proportional to λ 2d e and does not depend on the scale of the signal, namely

g ′ i ((x ′ , f ′ ), (x i ′ , f i ′ )) = λ 2d e g i ((x, f ), (x i , f i )).
Therefore, we may compute a normalized version (R e ) -2d g i of g i with R e . = max i trace(Var(x i ))

where Var(x i ) is the covariance matrix of the points cloud x i . From there on J n denotes the functional containing the normalized terms of J (i.e. each term being multiplied by an appropriate power of R e and an appropriate power of R f . = max i Var(f i )). By construction we (formally) have

J n (x, p, f ) = J ′ n (x ′ , p ′ , f ′ ) where J ′
n is computed for σ ′ e = λ e σ e and σ ′ f = λ f σ f . We now examine the gradient's behaviour under scaling. The goal is to keep a similar dynamic during the gradient descent whatever the (geometric or functional) scales of the fshapes are. As an illustration, consider the gradient with respect to x: we need J n (x -δ∇

x J n , p, f ) = J ′ n (x ′ -δ∇ x ′ J ′ n , p ′ , f ′ ) which yields to ∇ x ′ J ′ n = λ e ∇ x J n . As ∇ x ′ g ′ i = λ -1 e ∇
x g i we use the normalized version (R e ) 2-2d ∇ x g i of the gradient ∇ x g i . Each term composing the gradient ∇J n is then normalized with an appropriate power of R e and R f and we denote ∇ n J n the normalized version ∇J n .

In our implementation of algorithm 1 and 3 (resp. 2 and 4), we compute J n (resp. ∇ n J n ) rather than J (resp. ∇J ). These normalizations then guarantee that using similar but rescaled (or translated and rotated) data will provide comparable energies and energy decreases during gradient descent.

Consistency with respect to P

We now focus on the behaviour of the template estimation algorithm when the sampling of the discretized template changes. The methodology and notations will be similar to Section 9.3.1. Assume that we are working with surfaces: we wish to understand qualitatively the asymptotic behaviour of the expressions of the functional J and its gradient when the mesh on the template is refined (i.e. when P and T increase). The idea is once again to keep a similar behaviour in the optimization process at coarse and fine resolutions. We illustrate the invariance with respect to the number of points by a numerical experiment presented in Figure 18 and 19 where the same geometrico-functional matching is performed at various resolutions.

Let (x P , f P ) P be a sequence of regularly discrete fshapes sampled from a continuous fshape (X, f ) and containing respectively P points. Let also (p P ) P be a sequence of discretized momenta (sampled from a continuous vector field p) attached to the points of (x P , f P ). If the discretization is sufficiently uniform, we may consider that the typical size of a local volume element r ℓ is of order P -1 . When the number P of points is large, we then have

|f P | 2 x P . = ℓ f 2 ℓ r ℓ ≈ X f 2 (x)dH d (x)
where f is defined by formula (64). Notice that other discretization methods may be used to approximate the L 2 norm on X as discussed in Section 6.3. If we now assume that µ (x P ,f P )µ (X ,f ) W ′ → 0, then formula (66) implies that

µ (x P ,f P ) 2 W ′ ≈ µ (X,f ) 2 W ′
when P is large. Finally, if v p P -v p V → 0 then the initial velocity fields satisfies

v p P (•) = k=1 K(•, x k )p P k ≈ X K(•, x)p(x)dH d (x) = v p (•) (86)
when P is large. This means that the magnitude of the momenta should be proportional to the inverse of the density of points to generate comparable displacements at various resolution. In particular, the update of the momentums p P should be of order P -1 as equation ( 86) suggests and ∇ p P J is multiplied, in our code, by P -1 to be at the right scale. In the same spirit, it can be shown that the gradients ∇ x P J and ∇ f P J are of order P -1 . We then multiply these terms by P to ensure the homogeneity of the update and keep a similar dynamics along the optimization procedure even if the resolution of the meshes changes. We provide an example in Figure 18b where the values of the functional along the gradient descent at various resolutions are plotted.

(a) The source (x, f ) and the target (x 1 , f 1 ). Signal f = 0.5 is constant and f 1 is equals to 1 on the head and on the tail and 0 elsewhere.

Conclusion

In this article, we have proposed a well-formalized and efficient extension of the ideas of shape spaces for the treatment of functional shapes. To do so, we have introduced the structure of fshapes' vector bundle and the metamorphosis setting to model and quantify geometricofunctional transformations within these bundles. In addition, the concepts of geometric measure theory (varifolds in particular) were generalized to provide dissimilarity metrics between such objects. The combination of these settings enabled us to express atlas estimation on populations of fshapes as a variational problem, for which we were able to prove existence of solutions.

In the second part of the paper, we addressed the issue of providing practical numerical schemes to efficiently perform the functionals' optimization. In particular, we detailed the discrete expressions corresponding to the fvarifolds' data attachment terms, proposed a gradient descent algorithm for the estimation of all variables in the atlas and carefully examined some of the important numerical issues related to these algorithms. The resulting codes includes atlas estimation on populations of classical curves and surfaces and extends it to fshapes. They shall be made publicly available shortly. This was implemented under the simplified 'tangential model' : the full metamorphosis setting derived in the theoretical part of the paper is likely to be implemented as well in the near future.

Extensions of the present framework are possible in several directions worth exploring as future work. One of them is to consider fshape bundles modeled on more regular spaces than L 2 and hopefully recover existence of solutions in a more direct way than the proof presented in this paper, at the price of more involved numerical computations. Another track is to generalize such a framework to more general spaces of signals (e.g vector fields or tensor fields) and more general group actions, which has been considered only partially in [START_REF] Charon | Analysis of geometric and functional shapes with extensions of currents. Application to registration and atlas estimation[END_REF]. Finally, one important follow-up problem to atlas estimation is the one of statistical analysis and classification based on these fshape atlases : we have deliberately postponed such questions to an upcoming paper.

have div X (ωv ⊤ ) = ω div X (v ⊤ ) + ∇ω|v ⊤ . Applying the Divergence Theorem (cf [START_REF] Simon | Lecture notes on geometric measure theory[END_REF] Section 7) on the submanifold X gives:

X ω div X (v ⊤ ) = - X ∇ω|v ⊤ + ∂X ω ν, v ⊤
where ν is the unit outward normal to the boundary.

Derivative with respect to tangent spaces: We now come to the derivative term on the tangent space part in equation (87). To explicitly compute variations with respect to variables in the Grassmann varifold, the most convenient way is to use the embedding of G d (E) into L(E) that identifies any subspace V with the orthogonal projector p V on V . With this identification, one can represent the tangent space at V of G d (E) as L(V, V ⊥ ). Then, as explained with more details in [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF], if we set V t = d x φ t (T x X), one can show that:

d dt ↾ t=0 V t = p TxX ⊥ • ∇v ↾ TxX ∈ L(T x X, (T x X) ⊥ ) .
We can now introduce ∂ω ∂V as an element of L(T x X, T x X ⊥ ) * ≈ (T x X ⊥ ) * ⊗T x X and which we can write: ∂ω ∂V = n j=d+1 η * j ⊗ α j for (η d+1 , .., η n ) an orthonormal frame of T x X ⊥ and (α j ) vectors of T x X (as usual η * denotes the linear form η, . ). Then the variation we wish to compute is:

∂ω ∂V |∇v = n j=d+1 η j , ∇ α j v .
If we introduce ∂ω ∂V |v = n j=d+1 η * j (v)α j = n j=d+1 η j , v α j which is a tangent vector field on X, we have:

div X ∂ω ∂V |v = d i=1 n j=d+1
( e i , ∇ e i α j η j , v + e i , ∇ e i η j , v α j + e i , η j , ∇ e i v α j )

The last term in the sum is also n j=d+1 η j , ∇ α j v , which is nothing else than Hence, it follows that:

∂ω ∂V |∇v = div X ∂ω ∂V |v -(div X ∂ω ∂V |v) (88)
Integrating equation (88) over the submanifold X and using the Divergence Theorem as before, we find that:

X ∂ω ∂V |∇v = ∂X ν, ∂ω ∂V |v - X (div X ∂ω ∂V |v) (89)
Synthesis: Summing all the different terms from eq.(87), we eventually obtain: B Proof of Proposition 6

X (£ (v,h) ω) = X ∂ω ∂x -div X ∂ω ∂V

B.1 Perturbation

We introduce now a perturbation process on any measure ν on E × G d (E) × R that shall be useful for the following. Let a > 0 to be fixed later and consider for any t ∈ R the function ρ t : R → R such that ρ t (z) = z + t(sgn(z)a -z)1 |z|>a

where sgn(z) is the sign of z. We have ρ 0 = Id R and ρ 1 is a symmetric threshold at level a. Now for any t ∈ R, we denote ν t the new measure defined for any ω ∈ C b (E × G d (E) × R) as: ν t (ω) = ω(x, V, ρ t (f ))dν(x, V, f ) .

Obviously ν 0 = ν and ν 1 is such that ν 1 (|f | > a) = 0 so that t → ν t is an homotopy from ν to a measure under which the signal is a.e. bounded by a.

B.2 Proof of Lemma 2

We show the existence of a fvarifold minimizer in M X (cf (54)) for the extended functional J.

For any ν ∈ M X and t ∈ R, we denote J t . = J(ν t ) where ν t is the previously defined perturbation of ν (cf B.1) and we assume that J 0 < ∞ (which is equivalent to say that ν(|f | 2 ) < ∞). We recall that µ (X i ,f i ) -ν t 2 W ′ = (µ (X i ,f i ) -ν t )(ω i ) where ω i = K W (µ (X i ,f i ) -ν t ) ∈ W . Then one easily checks that J t < ∞ and, since we assume that W is continuously embedded into C 2 0 (E × G d (E) × R), with existing derivative J ′ t at any location t given by

J ′ t = ν d dt (ρ t (f )) γ f ρ t (f ) + γ W N i=1 ∂ω i ∂f (x, V, ρ t (f )) (92)
Using again the continuous embedding of W into C 2 0 (E ×G d (E)×R) we get for a constant C > 0 that

∂ω i ∂f (x, V, ρ t (f )) ≤ C ω i W ≤ C µ (X i ,f i ) W ′ + ν t W ′ . (93)
Moreover, as we mentioned after Proposition 4, µ (X i ,f i ) W ′ ≤ cte.H d (X i ). Similarly, ν t W ′ ≤ cte.ν t (E×G d (E)×R) and, since ν t ∈ M X , we have ν t (E×G d (E)×R) = H d (X) and consequently ν t W ′ ≤ cte.H d (X). Thus, there exists a constant K > 0 such that:

N i=1 ∂ω i ∂f (x, V, ρ t (f )) ≤ K N i=1 H d (X i ) + H d (X) (94)
Noticing now that d dt (ρ t (f )) ρ t (f ) ≤ 0, that | d dt (ρ t (f )) | = 0 for |f | ≤ a and that |ρ t (f )| ≥ a for |f | ≥ a and t ∈ [0, 1], we get for t ∈ [0, 1]

J ′ t ≤ ν - d dt (ρ t (f )) 1 |f |>a γ f a -γ W K N i=1 H d (X i ) + H d (X) (95)
so that

J(ν 1 ) ≤ J(ν 0 ) if a ≥ K γ W γ f N i=1 H d (X i ) + H d (X) . ( 96
)
An important consequence of (96) is that one can restrict the search of a minimum for J on fvarifolds ν such that ν(1 |f |>a ) = 0 (97)

with a = K γ W γ f N i=1 H d (X i ) + H d (X) .
In particular, since ν ∈ M X , we will have

x ∈ X and |f | ≤ a ν a.e.

Since X is bounded and G d (E) compact, we can restrict the search of a minimum on a measure supported on a compact subset C ⊂ E × G d (E) × R so that we introduce:

M X,C . = { ν ∈ M X | (x, V, f ) ∈ C ν a.e. } . (99)
An easy check shows that J is lower semi-continuous on the set M X,C for the weak convergence topology. In addition, M X,C is sequentially compact. Indeed, if ν n is a sequence in M X,C then all ν n are supported by the compact C and in particular (ν n ) is tight. Also, as already noted, there exists a constant cte independent of n such that ν n (E ×G d (E)×R) ≤ cteH d (X) and thus the sequence is uniformly bounded for the total variation norm. It results, thanks to the Prokhorov Theorem, that there exists a subsequence of (ν n ) converging for the weak topology. These compactness and lower semicontinuity properties guarantee the existence of a minimizer ν * of J with ν * ∈ M X,C and J(ν * ) ≤ inf f ∈L 2 (X)

J X (f ) . (100)

B.3 Proof of proposition

At this point, we do not yet have a minimizer of J X . The problem is that if the marginal on E ×G d (E) of ν * is the transport of H d |X under the application x → (x, T x X), we cannot guarantee that ν * does not weight multiple signal values in the fiber above a location (x, T x X). We will now show that for large enough γ f /γ W , there exists f * ∈ L 2 (X) such that ν * = ν X,f * so that we will deduce

J X (f * ) = J(ν * ) ≤ inf f ∈L 2 (X) J X (f ) (101)
and the existence of a minimizer on L 2 (X).

Let δf ∈ C b (E × G d (E) × R) and for any t ∈ R consider the perturbation ν t ∈ M X of any ν ∈ M X,C such that for any g ∈ C b (E × G d (E) × R) we have:

ν t (g) . = g(x, V, f + tδf (x, V, f ))dν(x, V, f ) . ( 102
)
Here again, the function t → J(ν t ) is differentiable everywhere and we have for

ω i . = K W (µ (X i ,f i ) - ν) d dt J(ν t ) ↾ t=0 = ν γ f f + γ W N i=1
∂ω i ∂f (x, V, f ) δf (x, V, f ) , so that when ν = ν * we get    γ f f + γ W A(x, V, f ) = 0 ν * a.e. with A(x, V, f ) . = N i=1 ∂ω i ∂f (x, V, f ) .

(103)

The partial derivative of f → γ f f + γ W A(x, V, f ) with respect to f equals γ f + γ W ∂A ∂f (x, V, f ). As before, using the continuous embedding W ֒→ C 2 0 (E × G d (E) × R), we have once again a certain constant K such that

∂A ∂f (x, V, f ) ≤ K N i=1 H d (X i ) + H d (X) , for all (x, V, f ) ∈ E × G d (E) × R. ( 104
)
It results that for γ f /γ W large enough and for all (x, V ) ∈ E × G d (E), f → γ f f + γ W A(x, V, f ) is a strictly increasing function going from -∞ at -∞ to +∞ at +∞ and thus there is a unique solution f (x, V ) to (103). Now, since the application f → γ f f + γ W A(x, V, f ) is also

C 1 on E × G d (E) × R,
Therefore, setting f * (x) = f (x, T x X) for x ∈ X, we see that |f * | ≤ a so that f ∈ L ∞ (X) and with (105), we deduce that ν * = µ (X,f * ) which shows that the solution of the optimization is a fvarifold associated to a true fshape (X, f * ). In addition, if X is a C p submanifold then x → T x X is a C p-1 function on X and, if W ֒→ C m 0 (E × G d (E) × R) with m ≥ 2 and m ≥ p, A and f are C p-1 functions so f * is also C p-1 , which concludes the proof of Proposition 6.

C Proof of Theorem 6

We shall basically follow the same steps as in the previous simpler cases. First of all, exactly as in 5.2.2, existence of a template shape X is guaranteed with the same compacity and lower semicontinuity arguments. Thus we may assume that X is fixed and we only have to show existence of minimizers to the simplified functional:

J X (f, (ζ i ), (v i )) . = γ f 2 X |f (x)| 2 dH d (x) + 1 2 N i=1 v i 2 L 2 ([0,1],V ) + γ ζ X |ζ i (x)| 2 dH d (x) + γ W µ (X i ,f i ) -µ (φ v i 1 (X),(f +ζ i )•(φ v i 1 ) -1 ) 2 W ′
Now, as for v 0 , due to the presence of the penalizations v i 2 L 2 ([0,1],V ) , one can assume that all vector fields v i belong to a fixed closed ball B of radius r > 0 in L 2 ([0, 1], V ). As in the proof of Proposition 6, we first show existence of a minimizer in a space of fvarifolds. Namely, extending the definitions of the previous subsections, we introduce the space M X of measures ν on E × G d (E) × R × R N such that for all continuous and bounded function h on E × G d (E), we have: ν(h) = h(x, V )dν(x, V, f, (ζ i )) = X h(x, T x X)dH d (x) .

For a measure ν on E × G d (E) × R × R N and a diffeomorphism φ, we denote by φ.ν the transport of ν by φ defined by: (φ.ν)(g) = |d x φ ↾ V |g(φ(x), d x φ(V ), f, (ζ i ))dν(x, V, f, (ζ i )) .

We now introduce the extended functional:

J(ν, (v i )) . = γ f 2 ν(|f | 2 ) + 1 2 N i=1 v i 2 L 2 ([0,1],V ) + γ ζ ν(|ζ i | 2 ) + γ W µ (X i ,f i ) -(φ v i 1 ).ν i 2 W ′
for ν ∈ M X , (v i ) ∈ (L 2 ([0, 1], V )) N and for all i ∈ {1, .., N }, ν i being the fvarifold defined for all ω ∈ W by:

ν i (ω) = ω(x, V, f + ζ i )dν(x, V, f, ζ i ) .
As previously, we can consider the perturbation function ρ t acting on signals and the measures ν t (g) . = g(x, V, ρ t (f ), (ρ t (ζ i )))dν(x, V, f, (ζ i )) .

Denoting J t = J(ν t , (v i )), we have, for t ∈ [0, 1],

J ′ t = ν d dt (ρ t (f )) γ f ρ t (f ) + γ W N i=1 |d x φ v i 1 ↾ V |. ∂ω i ∂f (φ v i (x), d x φ v i (V ), ρ t (f ) + ρ t (ζ i )) + ν N i=1 d dt ρ t (ζ i ) γ f ρ t (ζ i ) + γ W |d x φ v i 1 ↾ V |. ∂ω i ∂f (φ v i (x), d x φ v i (V ), ρ t (f ) + ρ t (ζ i )) (106)
where, for all i ∈ {1, .., N }, ω i = K W (µ (X i ,f i ) -(φ v i 1 ).ν i t ). On the first hand, we know that there exists a constant cte such that for all i, x ∈ E and V ∈ G d (E),

|d x φ v i ↾ V | ≤ cte|dφ v i | ∞ .
In addition, it is a classical result on flows of differential equations (cf [START_REF] Younes | Shapes and diffeomorphisms[END_REF]) that there exists a non-decreasing continuous function τ :

R + → R + independent of v ∈ L 2 ([0, 1], V ) such that |dφ v 1 | ∞ ≤ τ ( v L 2 ([0,1],V )
). Now, using the same controls as in the previous subsections, we have, on the other hand:

∂ω i ∂f (φ v i (x), d x φ v i (V ), ρ t (f ) + ρ t (ζ i )) ≤ ∂ω i ∂f ∞ ≤ cte ω i W ≤ cte( µ (X i ,f i ) W ′ + (φ v i 1 ).ν i t W ′ ) ≤ cte(H d (X i ) + (φ v i 1 ).ν i t (E × G d (E) × R)) .
It is also straightforward that (φ v i 1 ).

ν i t (E × G d (E) × R)) ≤ cte|dφ v i | ∞ ν i t (E × G d (E) × R)
and, using the fact that ν ∈ M X as already argued in 5.2.1,

ν i t (E × G d (E) × R) = H d (X). It results,

Figure 1 :

 1 Figure 1: Example of metamorphosis transformation of a fshape.

Figure 2 :

 2 Figure 2: Fshape bundle and metamorphosis.

1 ,

 1 • • • , p N and the functionals ζ 1 , • • • , ζ N . All the fshapes should be provided with their respective connectivity matrix. The momenta p 1 , • • • , p N are usually initialized to 0 and ζ 1 , • • • , ζ N are usually initialized to a constant.

Figure 6 :

 6 Figure 6: Mapping of the template on subjects 2 and 5. From left to right, we display the template and several intermediate time steps of the geometrico-functional transformations, i.e. (φ i t (x), f + ζ i t ) where ζ i t = tζ i as in formula (76).

Figure 7 :

 7 Figure 7: Estimated template and the geometrico-functional matchings to each subject.

Figure 9 :

 9 Figure 9: The hypertemplate and the estimated mean (computed with the full dataset containing 51 observations).

  .

Figure 10 :

 10 Figure 10: Results for Data 1 (control dataset).Figure 10e and 10f should be compared with Figure 10b and 10c respectively.

  Figure 10: Results for Data 1 (control dataset).Figure 10e and 10f should be compared with Figure 10b and 10c respectively.

  Mean template and the residuals (x, f + ζ 2 ) (view 1) (e) Deformed mean template and residuals φ v p 2 .(x, f + ζ 2 ) (view 1) (f) Deformed mean template and residuals φ v p 2 .(x, f + ζ 2 ) (view 2)

Figure 11 :

 11 Figure 11: Results for Data 2 (glaucoma dataset).Figure 11e and 11f should be compared with Figure 11b and 11c respectively.

  Figure 11: Results for Data 2 (glaucoma dataset).Figure 11e and 11f should be compared with Figure 11b and 11c respectively.

  84)where the source fshape (x, f ) ∈ E Px × R Px and the target fshape (y, g) ∈ E Py × R Py are two flat overlapping squares belonging to the same plane, see Figure12. The source (x, f ) contains P x = 6400 vertices and the target (y, g) contains P y = 400 vertices distinct from the vertices of x. The meaning of (84) is the following: we are trying to register (x, f ) onto the fixed target (y, g) by tuning the signal f ∈ R Px of the source only. This is not obvious to figure out what a good solution of this problem should be. Note that Proposition 6 ensures the existence of a proper solution of the continuous version of (84) if γ f is large enough.

( a )

 a The source (x, f ) (b) The target (y, g) (c) (x, f ) and (y, g)

  γ f = 20

Figure 13 :

 13 Figure 13: The fshape (x, f * ) where f * is the output of a gradient descent on J f un of equation (84) with various value of γ f .

  (a) Initialization (x init , 0) (b) Target (y, 0)

Figure 14 :

 14 Figure 14: The initial fshape (x init , 0) is a flat square with a centered opening. The target fshape (y, 0) is a (non smooth) bended version of (x init , 0). All the signals are zero.

Figure 15 :

 15 Figure 15: An example of a typical Delaunay triangulation and the characterization of a boundary point

  The arrows represent the initial gradient ∇ x J geo (x init ) on (x init , 0). The color represent the norm of this gradient (logarithmic scale).(b) Output of the gradient descent after few steps (has to be compared with 14b).

Figure 16 :

 16 Figure 16: Gradient descent on J geo of equation (85) starting from (x init , 0) of Figure 14a.

  The arrows represent the initial regularized gradient ∇x J geo (x init ) on (x init , 0). The color represent the norm of the gradient.(b) Output of the regularized gradient descent after few steps (has to be compared with 14b).

Figure 17 :

 17 Figure 17: Regularized gradient descent on J geo of equation (85) starting from (x init , 0) of Figure 14a.

  Value of J tan 0 (y-axis) as a function of the step number (x-axis) at three different resolutions.

Figure 18 :

 18 Figure 18: Geometrico-functional registrations of an ellipsoid onto a functional version of the Stanford's Bunny. The fshapes are depicted in 18a at medium resolution (around 15000 points). The experiment consists in minimizing the functional J tan 0 of equation (78) in p 1 and ζ 1 only (here N = 1). The results are presented Figure 19.

  (a) Geometrico-functional registration (at time t = 0, t = 0.5 and t = 1) (b) Target (c) Geometrico-functional registration (at time t = 0, t = 0.5 and t = 1) (d) Target

Figure 19 :

 19 Figure 19: Results of registration presented Figure 18. First row: coarse resolution (around 1200 points) computations takes 30 seconds. Second row: high resolution (around 80000 points) computations take 2.5 hours.

  for the two other terms in the sum, it is easy to see that it equals:

  we deduce from the Implicit Function Theorem that f is a C 1 function on E × G d (E). Going back to the solution ν * , we know that for ν * almost every(x, V, f ) ∈ E × G d (E) × R, we have (x, V, f ) ∈ C and f = f (x, V ), so that | f | ≤ a a.e. For any continuous and bounded function ω: ν * (ω) = ω(x, V, f )dν * = ω(x, V, f (x, V ))dν * and if we denote by ω(x, V ) . = ω(x, V, f (x, V )) which is a continuous and bounded function on E × G d (E), we have by definition of the space M X (eq.(54)): ν * (ω) = ν * (ω) = X ω(x, T x X, f (x, T x X))dH d (x) .

2 :

 2 Inputs: A momentum p 0 , a template signal f and N momenta p i and functional residuals ζ i .

3: Begin 4:

  Now, we remind that ω(x) = ω(x, T x X, f (x)) so and applying the result of equation (88) to v ⊤ :

					(∇ω|v ⊤ ) =	∂ω ∂x	|v ⊤ +	∂ω ∂V	|∇v ⊤ +	∂ω ∂m	∇f, v ⊤
						∂ω ∂V	|∇v ⊤ = div X	∂ω ∂V	|v ⊤ -(div X	∂ω ∂V	|v ⊤ )
	We notice that		∂ω ∂V	|v ⊤ = 0 by the expression of	∂ω ∂V	and using the equality v = v ⊤ + v ⊥ we
	find eventually that:
	X	£ (v,h) ω =	X	∂ω ∂x	-div X	∂ω ∂V	-ωH X |v ⊥ +	∂ω ∂m	.(h-∇f, v ⊤ )+	∂X	ν,	∂ω ∂V	|v +ωv ⊤
	which proves the result of Theorem 5.
								X	∂ω ∂m	.h +	∂X	ν,	∂ω ∂V	|v + ωv ⊤ .

|v -X ∇ω|v ⊤ + ω H X |v ⊥ +

We are using here that for X compact, the mapping v → φ v 1 .X is continuous for the weak convergence on v ∈ L

([0, 1], V ) and the convergence for the Hausdorff metric on the set of all compact subsets of E

Acknowledgment

We would like to thank Mirza Faisal Beg, Sieun Lee, Evgeniy Lebed, Marinko Sarunic and their collaborators for providing the OCT dataset and for fruitful discussions. The authors also acknowledge the support of the French Agence Nationale de la Recherche project HM-TC (number ANR-09-EMER-006).

Declare a boolean to break the following loop: breakLoop ← False. Update each variable separately: for each i = 1, . . . , K compute Ĵnew

12:

Update every variable at the same time:

13:

Keep the best configuration: let

Update points and steps sizes:

16:

breakLoop ← True. Decrease all steps sizes: until "breakLoop = True" or "steps sizes are too small".

21:

until "Maximum iteration is reached" or "descent is to small" or "steps sizes are too small". 22: End 23: Output:

Some tricks may be also used to improve the optimization strategy. For instance, it may be convenient to regularize some part of the gradient as discussed in Section 9. Another efficient method is to change the fvarifold kernel widths σ e , σ f , σ t along the optimization process. In that case, the optimization involves several runs: a first gradient descent is performed at coarse scale by choosing large σ e , σ f , σ t and then a second gradient descent is performed with smaller σ e , σ f , σ t (with starting point equals to the end point of the first run) and so on.

8 Numerical experiments

Synthetic dataset

We now present a first set of results of the previous atlas estimation algorithms on a synthetic dataset of six textured statues constructed using Sculptris software. The set of subjects is shown on the last row of Figure 7. The atlas is estimated using the free mean fshape tangential A Variation formula for fvarifolds: proof of Theorem 5

The proof follows the same steps as the corresponding result for usual varifolds (cf [START_REF] Charon | The varifold representation of non-oriented shapes for diffeomorphic registration[END_REF]). Given a C 1 vector field v on E with compact support, we can consider the 1-parameter group of diffeomorphisms φ t with φ 0 = Id and ∂ t↾ t=0 φ t = v. Then, it follows that:

As we see, the previous leads to several terms in the derivative: differentiate the volume change term J t . = |d x φ t↾ TxX |, the function ω with respect to the position variable, with respect to the tangent space direction and to the signal part. The derivative with respect to point positions and signal values are easy to obtain and equals respectively, since ω is assumed to be C 1 , ∂ω ∂x |v and ∂ω ∂m h. The two other terms require more attention.

Derivative of the volume change: For any vector field u defined on X, we shall denote by u ⊤ and u ⊥ the tangential and normal components of u with respect to the tangent space of X at each point. We also introduce the connection ∇ • • on the ambient space and an orthonormal frame of tangent vector fields (e i ) i=1,..,d on X. Now J t = det([ d x φ t (e i ), d x φ t (e j )] i,j ) so a simple calculation shows that:

which is the tangential divergence of the vector field v ⊤ denoted usually div X (v ⊤ ). The second term becomes d i=1 e i , ∇ e i v ⊥ . For all i = 1, .., d, we have e i , v ⊥ = 0 so that after differentiation we find that e i , ∇ e i v ⊥ = -∇ e i e i , v ⊥ . Therefore:

In this last expression, we recognize the mean curvature vector to the submanifold X, which is the trace of the Weingarten map and is denoted H X . As a result, we find that:

Now, the first term can be rewritten as a boundary integral by applying the Divergence Theorem. Indeed, if we denote by ω the function defined on X by ω(x) = ω(x, T x X, f (x)) which is C 1 , we from all the previous inequalities, the existence of a non-decreasing continuous function that we will still call τ such that for all i, x, V, f, ζ i :

Following the same path that previously lead to (95)

Just as in 5.2.1, this implies that J(ν 1 , (v i )) ≤ J(ν 0 , (v i )) as soon as:

Therefore, one may restrict the search of a minimum on a set of measures ν that are supported on a compact subset C of E × G d × R × R N , which space we shall denote again M X,C . The rest of the proof is now very close to the one of 5.2.1. Due to lower semi-continuity of the functional and the compacity of M X,C and B for the weak convergence topologies (respectively on the space of measures and on L 2 ([0, 1], V )), we obtain the existence of a minimizer (ν * , (v i ) * ) for the functional J. The last step is to prove that ν * , which belongs a priori to the measure space M X,C , can be written under the form ν * = ν X,f * ,(ζ i * ) , i.e that there exists functions f * and ζ i * on X such that, for all continuous and bounded function

We then consider variations of the signals (δf, (δζ i )) all belonging to the space

) and the path t → ν t defined by:

, expressing that J ′ t↾ t=0 = 0 for all δf and (δζ i ) gives, similarly to 5.2.1, the following set of equations:

. The derivatives ∂ f A and ∂ ζ i A can be shown again to be uniformly bounded in x, V, f, ζ i , and a previous argument provides the existence of unique solutions f = f (x, V ) and ζ i = ζi (x, V ) to (110). The rest of the proof is exactly the same as in the end of appendix B.3: we set f * (x) = f (x, T x X) and ζ i * (x) = ζi (x, T x X), which are again L ∞ functions on X. In addition, one shows easily that the minimizing measure ν * equals ν X,f * ,(ζ i * ) in the sense of (109). Finally, the regularity of f * and ζ * when X is a C p submanifold is obtained again by applying the Implicit Function Theorem to (110).

D Proof of Theorem 7

We shall only sketch the essential steps to adapt the content of appendix C. We start by writing (61) in an extensive way. This gives:

Now, with Lemma 3, we know that the optimal functions h 0 * and h i * are given by (60) and thus the variational problem of (111) can be replaced by the optimization with respect to residual functions ζ 0 and ζ i living in L 2 (X) of the functional:

where C 0 (x) . = (

But we note that the previous, up to the weights in the L 2 metrics given by functions C i , becomes now extremely close to the problem examined in Theorem 6. In fact, the proof of appendix C can be adapted almost straightforwardly to this situation. As previously, the essential step is to reformulate the optimization problem in a space of measures. Defining the functions:

for i ∈ {1, ..., N } and (x, H) ∈ E × G d (E), we can set, with the same definitions as in appendix C:

(113) for ν ∈ M X . The rest of the proof follows the same path, relying on the fact that we can assume the vector fields v 0 and v i to be bounded in L 2 ([0, 1], V 0 ) and L 2 ([0, 1], V ) as we explained in the beginning of appendix C. This implies, as already argued in the same section, that we have uniform lower and upper bounds for |d x φ v i s ↾ H |, s ∈ [0, 1] and i ∈ {1, ..., N } and for the quantities |d x [φ v 0 s • (φ v 0 1 ) -1 ] ↾ H |. Consequently, we can assume that we have α, β > 0 such that for all i ∈ {0, ..., N }, α ≤ Ci ∞ ≤ β. Using these inequalities, one can check that we get equivalent controls as in the proof of Theorem 6 which allows us to conclude the existence of a measure minimizer for the extended functional J and then go back to a fshape solution for J with a similar implicit functions' argument.