I R Ène

Alain Trouv

ESTIMATION OF A GROWTH DEVELOPMENT WITH PARTIAL DIFFEOMORPHIC MATCHINGS

In the field of computational anatomy, the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework has proved to be highly efficient for addressing the problem of modeling and analyzing of the variability of populations of shapes, allowing for the direct comparison and quantization of diffeomorphic morphometric changes. However, with the progress achieved in medical imaging analysis, the interest for longitudinal data set has substantially increased in the last years and requires the processing of more complex changes, which especially appear during growth or aging phenomena. The observed organisms are subject to transformations over time that are no longer diffeomorphic, at least in a biological sense. One reason might be a gradual creation of new material. The evolution of the shape can then be described by the joint action of a deformation process and a creation process.

In this paper, we extend the LDDMM framework to address the problem of non diffeomorphic structural variations in longitudinal data. We keep the geometric central concept of a group of deformations acting on embedded shapes. The need for partial mappings leads to a timevarying dynamic that modifies the action of the group of deformations. We develop a theoretical framework and two algorithms to estimate realistic individual growth scenarios from a set of observations sparsely distributed in time. We present few numerical experiments on animal horns where the shapes are modeled by oriented varifolds. Each computed scenario is parametrized by low-dimensional variables providing the support for statistical analysis.

1. Issues about growth modeling 1.1. Context. In the field of anatomy, the massive investment in the acquisition of medical imaging calls for the development of new numerical techniques to model and analyze the variability of large databases. Already a few decades ago, the willingness to help neuroscientists and clinicians in the analysis of the substructures of the human brain led to a new discipline named Computational Anatomy by U. Grenander and M. Miller [27]. Various mathematical frameworks are at the foundations of this new field [6,58]. The developed theories and methods for registering and comparing shapes have been successfully applied to, among many others examples, the study of the shape of Hippocampus in relation to the evolution of Alzheimer disease, similar works on the planum temporale for schizophrenia, Down syndrome, the analysis of brain connectivity based on DTI imaging, studies of heart shapes and malformations. Miller, Trouvé, and Younes recently presented a review [START_REF] Miller | Hamiltonian Systems and Optimal Control in Computational Anatomy: 100 Years Since D'Arcy Thompson[END_REF] of the approach, called Diffeomorphometry, that consists in the definition of shape spaces as homogeneous spaces under the action of a group of diffeomorphisms. This construction has provided theoretically sound and numerically efficient tools, like the Large Deformation Diffeomorphic Metric Mapping (LDDMM), allowing to consider a wide variety of databases of shapes, as images, landmarks, curves, surfaces, fiber sets, or more recently, functional shapes that are shapes equipped with a signal [START_REF] Charon | Functional Currents: A New Mathematical Tool to Model and Analyse Functional Shapes[END_REF][START_REF] Charlier | The fshape framework for the variability analysis of functional shapes[END_REF].

Besides the cross-sectional variability analysis emerges the study of longitudinal data sets. Each subject of a population is represented by a sequence of measurements at different times. Among many other examples, the interest for these more complex data is motivated by the clinical studies of diseases or treatments that have a progressive impact over time and therefore entail changes on these evolution scenarios [START_REF] Tang | Miller Shape abnormalities of subcortical and ventricular structures in mild cognitive impairment and Alzheimer's disease: detecting, quantifying, and predicting[END_REF]. Although modeling evolution scenarios and analyzing their variations appear as two different processes, in a lot of situations they can both be achieved by diffeomorphic registration. Shape spaces as Riemannian manifolds have provided various methods ranging from parallel transport [46], Riemannian splines [START_REF] Trouvé | Shape Splines and Stochastic Shape Evolutions: A Second Order Point of View[END_REF], geodesic regression [43, [START_REF] Vialard | 3d image registration via geodesic shooting using and efficient adjoint calculation[END_REF]22] including the inference from a population of a prototype scenario of evolution and its spatiotemporal variability [START_REF] Durrleman | Toward a comprehensive framework for the spatiotemporal statistical analysis of longitudinal shape data[END_REF].

In the past decades, growth modeling has been mainly addressed with biophysical models based on partial differential equations (see for example the numerous studies on tumor growth [START_REF] Benzekry | Classical Mathematical Models for Description and Prediction of Experimental Tumor Growth[END_REF] or plant growth). These models provides highly accurate description for each patient as they are able to take into account a large set of patient-specific parameters. The difficulty comes then from multi-subject comparisons and the lack of inter-subjects correspondences. Interestingly, among few others, Zacharaki et al. [START_REF] Zacharaki | ORBIT: A Multiresolution Framework for Deformable Registration of Brain Tumor Images[END_REF] developed a registration method driven by a biophysical model.

In the field of optimal transport, a lot of new metrics have emerged to capture growth phenomena and relax the conservation of mass of the generic framework. For this approach, called unbalanced optimal transport, Piccoli and Rossi [44], and Lombardi and Maitre [START_REF] Lombardi | Eulerian models and algorithms for unbalanced optimal transport[END_REF] consider a source term in the transportation equation of the density. In order to model growth tumor, Lombardi and Maitre control the mass increase by the source term that can either be optimized or used to integrate prior biological information. Another method introduced by Figalli and Gigli [START_REF] Figalli | A new transportation distance between non-negative measures, with applications to gradients flows with Dirichlet boundary conditions[END_REF] is to consider the boundary as a source or a sink of mass. However, this reserve of mass can only be used at the ends of a geodesic path. In a recent work, Feydy et al. [START_REF] Feydy | Optimal Transport for Diffeomorphic Registration[END_REF] use unbalanced regularized optimal transport methods as a fidelity term for diffeomorphic registration purpose.

These new metrics allow to adjust the contributions of the transport part and of the mass creation/destruction in the optimal transformation. Up to now, the longitudinal analysis has been limited to the study of data sets with homologous observations. Yet, in some situations this assumption seems inappropriate. During the growth or the degeneration of an organism, the changes occurring over time cannot always be modeled by diffeomorphic transformations, at least in a biological sense.

1.2. Two main types of growth process. Growth refers to a positive change in size of an organism by adding material. On the macroscopic scale, one can globally identify two main types of growth process:

• Type I: a growth homogeneously distributed.

• Type II: a growth process that involves new material on specific areas, usually on the organism's boundary. Although the first case seems to be the most common, many examples illustrate the second type of growth as crystal growth or mineralized tissues as bone, horn, mollusc shells, tendon, cartilage, tooth enamel. Plant growth offers also examples of the two processes as illustrated in Figure 1 (see also [10]). On the first row, besides various local growth rates, the growth of the leaf mainly involves a scaling process (growth of type I). One could consider that we have a creation of new material stricto sensus but the homology structure remains stable at a macroscopic and anatomic viewpoint. The growth can thus be modeled by a one-to-one deformation process. The situation is more complex on the second row. Although the leaf evolves through a deformation globally similar to a scaling, one can observe at the bottom of the leaf the emergence of new material highlighted by some geometrical features: new veins and new indentations on the boundary of the leaf (growth of types I & II). Note that it is not the emergence of new veins that justifies the growth of type II but the creation of new areas providing the place for the new veins. Once again, one could argue that theses areas contains new cells and old cells but for registration purpose at macroscopic scale, it is more coherent to say that theses new lands admit no equivalent in the prior ages of the leaf. Note at last that a creation process is not always highlighted by geometrical features as illustrated in Figure 1 by the development of bone. The red doted areas are called epiphyseal plates. Although the ossification is a complex diffused process, the main part of the elongation comes from these areas where the long bones grow longer.

If one wants to analyze the evolution of a small region of the leaf or the bone, a biologically coherent registration is mandatory, i.e. the precise identification of corresponding points between multiple observations. When a creation process occurs, this registration should not be achieved with one-to-one deformations. Even if the region of interest exists in the youngest shape, e.g. around the yellow point in the bone, a deformation that would not consider the creation of new areas, e.g. similar to a scaling, would unlikely retrieve the accurate image of this point in the second bone. Remark at last that the second leaf grows through both a deformation process and a creation process. The nonlinear combination of these two processes is a challenging problem.

In [START_REF] Gonzalez | Leaf size control: complex coordination of cell division and expansion[END_REF], authors express the importance of dissecting the genetic networks underlying plant growth and organ size. They highlight parameters as the number of the source cells at the beginning of the growth, the rate and duration of both cell division process and cell expansion process. Modeling the evolution of a shape during a growth process requires thus to be able to integrate either the emergence of new geometrical features or more generally the main biological priors driving the growth process. The model introduced in this paper aims at finding a trade-off between the reliable replication of the growth dynamic and the simplicity of diffeomorphisms allowing subject comparison. It will be applied to the development of animal horns which are considered as rigid shapes. Their growth process, illustrated hereafter in Figure 2, is strictly of type II and although it induces rigid displacements, it is not affected by distortions due to other causes. At last, horn growth still illustrates a general situation by the absence of geometrical features, providing thus a perfect case study.

The examples in Figure 1 raise the question of expliciting the homology relations between two shapes.

1.3. Coordinate systems in growth modeling. In his pioneer work, D'Arcy Thompson illustrated through several examples that the differences between related species could be explained by simple geometrical deformations highlighting a common coordinate system to all the related shapes of a given population. Before proceeding further, we present below a brief survey of the literature on the topic of coordinate system for growth modeling.

Seashells grow through the local creation of mass along the aperture. Considering that a seashell is a rigid shape, its development is then completely described by the longitudinal evolution of its aperture. Moulton et al. [START_REF] Moulton | Mechanical growth and morphogenesis of seashells[END_REF] describe this growth in a local basis attached to each point of the aperture. We can say that the authors are introducing a dynamic coordinate system to study this growth process. More generally, Hejnowicz [START_REF] Hejnowicz | Trajectories of principal directions of growth, natural coordinate system in growing plant organ[END_REF] lists different coordinate systems to express the principal directions of an organ's growth and to identify different patterns of nonrandom alignments in the cell wall network. With a more local approach, Grenander et al. [START_REF] Grenander | Characterization of biological growth using iterated diffeomorphisms[END_REF] seek to decompose the growth process by typical biologically-meaningful local deformations to capture and characterize different patterns. They apply their model on growth data of a rat's brain. Aiming for future extensions of the model, they introduce the notion of a Darcyan coordinate system as a natural curvilinear system evolving jointly with the shape.

A coordinate system also emerges in brain surface analysis where modeling the development of the cortical surface from gestation to adulthood remains a challenge. The large inter-subject variability prevents to readily extract anatomical similarities between brains and the registration of cortical surfaces is still limited by the ill-defined problem of homologies between different cortices. However, a large set of geometrical features of different scales and dimensions, like sulci (folds of the surface), sulcal lines, sulcal basins or sulcal pits, are mostly common to every individual and describe a global cortical organization. They allow to exhibit a natural coordinate system to parametrize the cortical surface [START_REF] Auzias | Model-driven harmonic parameterization of the cortical surface: HIP -HOP[END_REF][START_REF] Auzias | Model-driven parameterization of fetal cortical surfaces[END_REF][START_REF] Clouchoux | Model-driven parameterization of the cortical surface for localization and inter-subject matching[END_REF].

An interesting point in this last situation is that a unique coordinate system allows to describe the complete formation of the final organization of the cortical surface with sulcal lines. Returning to the examination of leaf growth, a unique coordinate system could also describe the vein pattern of the oak's leaf. However, for the basil, it would shrink the last veins to appear in a small region of the coordinate system. In this case, we would rather consider a coordinate system where new coordinates appear over time with the primary purpose to prevent correspondences between non homologous points. In addition to integrating prior information to the parametrized shapes, e.g. an organization with geometric features as vein patterns or sulcal lines, it should contains the time of creation of its new coordinates.

1.4. Growth evolution by foliation. At this stage, we need to specify a creation process for which we can anticipate the creation of new coordinates and their organization in a structure common to the whole population of interest. In Figure 2, we extrapolate the drawings of D'Arcy Thompson to examine creation or degeneration processes on some examples. On the left on top, we compare two ages of an animal horn. Below, we represent an organism similar to an onion or the cross section of a tree trunk. Finally on the right, an other representation of a membrane, such as the human skin where the double arrow indicates that we can read this example in both directions, as the growth or the degeneration of an organic tissue. In each case, we assume that new material is progressively created during the growth at the boundary of the shape. With the regularity of this growth process, the subsets of the new points created at each time define a natural decomposition of the complete set of all the coordinates. These subsets all have the same shape: similar to lines for the horns or circles for the onion. They induce a collection of foliated leaves meaning that the shape is locally similar to a connected disjoint union of parallel lines. The creation of new material is linked to the appearance of new coordinates in order to provide a biological coordinate system that ensures the pointwise homology between two shapes of different ages. Consequently, the grids are either extended or shortened with the red or blue lines. Moreover, one can associate to each leave its time of creation and likewise each coordinate can inherit the associated time marker. At last, note that the three deformations induced by the coordinates in Figure 2 are partial mappings delimited by the restriction to the black grids.

Such a type of growth process by foliation is a strong biological prior that compensates the lack of homology between individuals and the lack of geometrical features to drive the evolution of a shape. Although this paper is focused on application to animal horns, we will discuss at the end how to exploit it for plant growth.

Remark 1.1. A foliation can be simply described as a decomposition of a manifold into pathconnected submanifolds, called leaves, such that the manifold looks locally like a parallel union of these leaves. We refer to [47,[START_REF] Lee | Introduction to Smooth Manifolds[END_REF] for more details.

1.5. Contributions. As summarized in [START_REF] Durrleman | Toward a comprehensive framework for the spatiotemporal statistical analysis of longitudinal shape data[END_REF], a statistical framework for longitudinal data analysis can be divided in three components:

(1) For each subject of the data set, the estimation of a continuous shape evolution that fits the few observations at different times of the subject. (2) A temporal re-alignment of all individual trajectories.

(3) The estimation of a mean scenario representative of the population called the template allowing to establish the variance of the population.

This paper investigates the first item for growth modeling. A core hypothesis is that all the scenarios of a studied population share a common growth pattern. We introduce a new type of coordinate systems, called biological coordinate systems, to explicit this growth process. We then extend the LDDMM methods to integrate scenarios modeled on a given biological coordinate system. This system plays a key role to restrict the diffeomorphic registrations and to anticipate the creation of new material. At the end of the modeling, one should get a representation of the population as illustrated in Figure 3 where the biological coordinate system is some kind of canonical growth scenario. Note that its role is different from the template defined as a mean scenario. A practical goal of this paper is, given a temporal sequence of shapes (S tar i ) i of a horn at times (t i ) i , to provide numerical algorithms able to retrieve its continuous development from its youngest state to its oldest one. This means to generate a scenario t → S t such that S ti ≈ S tar i for all i. The theoretical approach will be developed in Section 3 and 4. Section 6 will present some numerical results for which the shapes are represented by oriented varifolds introduced in Section 5.

We will see that our model can also produce a path modeling the complete development of a horn from only one final observation. If we imagine the horn at its birth as reduced to a single point, we can construct a continuous path from this point to a nontrivial shape at the final time matching the given observation. In each application, the complete development of a horn produced by the algorithm is encoded in a low-dimensional forecast initial condition, providing the support to a statistical analysis.

Throughout this paper, the growth scenarios are all defined on a time interval fixed to [0, 1].

Modeling growth evolutions

2.1. Biological coordinate system. In the sketch of Figure 3, all individuals follow the same development pattern. The body is created first, head comes next and at last, hands and feet appear. The biological coordinate system codes this information. We will model the evolution of a shape during a growth process by a curve of mappings t → (q t : X → R d ) defined on a compact space X called the coordinate space. The coordinates of X should describe the homologous correspondences between two ages of an individual's shape and we will thus call them biological coordinates. The space X contains all the biological coordinates needed to model the complete evolution, but each coordinate will not necessary be used during the whole evolution. For any time t ∈ [0, 1], we denote X t the set of coordinates that actually parametrized the shape at time t. An individual scenario (t → S t ) t∈[0,1] is then modeled by:

(1) S t = q t (X t ) .

X t is called the set of active coordinates at time t. If new material is created in a time interval [s, t], then new coordinates are activated, i.e. X s = X t . Since we aim to model growth evolutions with no degeneration process, we assume that for any s ≤ t ∈ [0, 1],

(2)

X s ⊂ X t .
More precisely, this nested sequence is described by a function τ : X → [0, 1], called the birth tag, that indicates when a coordinate x ∈ X starts to be active, i.e. at what time the material associated to this coordinate is created. For any time t ∈ [0, 1], we define

(3) X t . = {x ∈ X | τ (x) ≤ t} .
A pair (X, τ ) as defined above is called a biological coordinate system. Equation [START_REF] Arguillère | Shape Deformation Analysis From the Optimal Control Viewpoint[END_REF] implies that a population parametrized by a given biological coordinate system share a common growth process. Intuitively, the curve of nested shapes t → X t forms a canonical scenario that describes this growth process.

Example 2.1. Figure 4 explicits the growth process associated to the biological coordinate system

(4)    X . = [0, 1] × S 1 , τ (x) . = t for any x = (t, x 0 ) ∈ X .
At time 0, the active coordinates form a circle. They grow into a cylinder by creation at each time t ∈]0, 1] of a new circle that is added at the bottom of the shape. The regularity of τ gives to X its foliated structure. The coordinate space X is a connected disjoint reunion of leaves that are the level lines of τ . We will denote

(5) X {t} . = {x ∈ X | τ (x) = t}
the new leaf at time t in the general case. Here, we have X {t} = {t} × S 1 .

To induce this growth process to the whole population, we assume that X is a subset of R d and that any individual scenario is given by a curve of continuous and injective mappings. An exception is yet allowed to model a population of horns with the biological coordinate system given by equation (4): the image of the first set of coordinates X {0} is reduced to a point to form the tip of the horn. Figure 5 illustrates two examples of scenarios built on this coordinate system. We retrieve for both scenarios the growth process described in Example 2.1. At each time t > 0, a new leaf, homeomorphic to a circle, appears at the bottom of the shape and the points of this new leaf have no biological correspondence with the points of the younger shapes.

Figure 5. Two examples of individual scenarios built on a given biological coordinate system (X, τ ). Each image scenario inherits the foliation of the biological coordinate system induced by the birth tag τ and enlightened by the color gradient.

2.2. Partial diffeomorphic registration under constraint. We now investigate how to generate the curve of mappings (q t : X → R d ) t , modeling an individual scenario, by diffeomorphic registration of each time step. The central idea of this paper is to produce a biologically coherent model of the evolution of a shape during a growth process. The action of any diffeomorphism on the shape should thus exhibit the homologous points between the shape at ages s and t. Consider a flow (a curve) of diffeomorphisms t → φ t on R d , the action of these diffeomorphisms during a time interval [s, t] is given by the deformation φ s,t . = φ t • φ -1 s as illustrated in Figure 6. Any coordinate x ∈ X models a point of the scenario that is created at time s = τ (x) at the position q s (x) ∈ S s . Its position at any time t > s is then given by

q t (x) = φ s,t (q s (x)) ∈ S t . (6)
Since the points q s (x) and q t (x) inherit the birth tag of the coordinate x, the action of φ s,t does preserve the birth tag. This parametrization of the shapes by the biological coordinate system via the mappings (q t ) t is the key ingredient that ensures the preservation of the birth tag under the action of the flow. Once a coordinate x ∈ X is activated, its image is thus completely determined Figure 6. Partial mappings and preservation of the birth tag. In this sketch with three discrete times, we simplify and display the values of the birth tag with only three colors to highlight the preservation of this birth tag under the action of the flow. by the flow. The subtlety of this model is to define the position of all the coordinates when they are activated. For this purpose, we introduce the birth place function q : X → R d mapping the coordinates to their respective place of birth. This function is an individual parameter of each subject of the population. We will see that given any horn at its final state, we will always have a good estimation of this parameter.

Any coordinate x ∈ X is activated at time τ (x) and its position is then fully determined by

q τ (x) (x) = q(x) , for any t > τ (x) , q t (x) = φ τ (x),t (q(x)) . (7)
For a unified framework, we want to extend the mappings q t : X t → R d into a single space of mappings on X. We thus anticipate the position of coordinates that are not activated yet. We choose to say that for any x ∈ X, (8)

q t (x) =    q t (x) if τ (x) ≤ t q(x) otherwise
.

If the point associated to a coordinate x ∈ X does not exist at time t, q t (x) returns its future place of birth. It follows that ( 9) q 0 = q .

Remark 2.1. The development of a horn presents a specificity in the location of the new points. They all emerge in a plane area that is implicitly the boundary of the animal's head. This planar constraint helps us to estimate the birth place function from the final age of a horn as illustrated in Figure 7 and justify the choice of equation (8).

Figure 7. First row : two ages of a horn given at an intermediate time t 0 ∈ [0, 1] and at final time 1. The small horn is the image q t0 (X t0 ) and the large horn q 1 (X 1 ) = q 1 (X). Second row: on the left, we additionally see the initial position of all the inactive coordinates that will gradually appear to form the large horn. On the right, we display the virtual horn at time 0. It is a flat disc, given by q 0 (X) = q(X), that contains the (anticipated) initial positions of every coordinates. Note that this flat horn can be roughly estimated from the bottom boundary of the full horn and the rings, or leaves, are linearly deduced. The sizes and points of view have been adjusted for better visibility but the colors give the pointwise correspondences between the figures.

Remark 2.2. The concept of growth mapped evolution (GME) has been introduced in [START_REF] Kaltenmark | Trouvé Partial Matchings and Growth Mapped Evolutions in Shape Spaces[END_REF] to formalize the idea of partial diffeomorphic registrations under the constraint of a tagging function (Figure 6). This work led to the introduction of a coordinate system that parametrizes the shapes as presented above. Given a nested sequence (X t ) t , we exhibited the topological conditions required to ensure the existence of a tagging function τ that satisfies equation (3). See also [START_REF]Kaltenmark Geometrical Growth Models for Computational Anatomy[END_REF] for more details.

2.3.

Illustration of the new dynamic for growth scenarios. In order to estimate later a cost for the action of any flow of diffeomorphisms and as diffeomorphometry methods usually proceed, a flow of diffeomorphisms (φ t ) t∈[0,1] is produced by integration of a time-varying vector field v : [0, 1] → V where V is a space of vector fields on the fixed ambient R d :

(10)

φ t = Id + t 0 v s • φ s ds .
Given a biological coordinate system (X, τ ), the growth scenarios modeled on this system can then be parametrized by an initial position q 0 : X → R d and a time-varying vector field v : [0, 1] → V . By derivation of ( 7) and ( 8), the scenario represented by a pair (q 0 , v) is given by the solution q of qt (

x) =    v t (q t (x)) if τ (x) ≤ t 0 otherwise = 1 1 τ (x)≤t v t (q t (x)) . ( 11
)
The shapes of the individual's scenario are then retrieved by the restriction of the mappings (q t ) t∈[0,1] to the active coordinates: [START_REF] Charon | Functional Currents: A New Mathematical Tool to Model and Analyse Functional Shapes[END_REF] S t = q t|X t (X t ) .

Definition 2.1. We will refer throughout this paper to equation (11) as the growth dynamic.

Remark 2.3. When τ ≡ 1, the shape evolves through a pure deformation process and we retrieve the standard dynamic through the flow.

We present now few examples to illustrate how to build scenarios from birth place functions and time-varying vector fields. We consider again the biological coordinate system given by X = [0, 1] × S 1 ⊂ R 3 and τ the projection on the first coordinate. We recall that Figure 4 highlights the trivial scenario induced by this system and Figure 5 offers two examples of image scenario.

Example 2.2. The initial position q 0 is here given by embedding X in R 3 then projecting it on the horizontal plane so that each leaf is sent on the unit circle of R 2 × {0}. This localizes the area where occurs all the creation of the growth process. The time-varying vector field acting on the ambient space is a simple constant upward translation. The resulting scenario is displayed in Figure 8. One can see in this example that the initial position of the shape is not an embedding in

0 -0.2 0 0.2 -0.2 0 0.2 0.2 0.4 0.6 0.8 1 1.2 0 -0.2 0 0.2 -0.2 0 0.2 0.2 0.4 0.6 0.8 1 1.2 0 -0.2 0 0.2 -0.2 0 0.2 0.2 0.4 0.6 0.8 1 1.2 Figure 8
. Generation of a cylinder. From left to right, the initial position q 0 (X), an intermediate position q t0 (X) at time t 0 ∈]0, 1[, the final position q 1 (X) and the biological coordinate system (X, τ ). This example shows a situation where q 0 is not an embedding and yet q 1 is one. The cylinder at time 0 is completely folded flat on itself and it unfolds gradually until time 1 when it is fully grown. All the creation process occurs at the base of the cylinder. Each newly created leaf pushes upwards the rest of the cylinder.

the ambient space. This initial shape can be seen as a compressed accordion that will progressively be unfolded. Even if the flow is diffeomorphic, it only sees gradually the shape. These mingled leaves can thus be separated over time. The flow will yet never be able to separate two points that appear at the same time at the same position.

Example 2.3. In order to produce a population of horn's evolutions, we consider now a family of initial positions of parameter ρ whose images fit the unit disk:

q 0 : [0, 1] × [0, 2π[ → {0} × R 2 (t, θ) → 0, ρ(t) cos(θ), ρ(t) sin(θ) , (13)
where ρ : [0, 1] → [0, 1] is an increasing homeomorphism. We consider three specific v : [0, 1] → R to generate vertical translations in R 3 . One is constant, one is increasing, and the last one is decreasing. Figure 9 shows these three vector fields and the associated scenarios that are generated with also three different initial positions of parameter ρ where ρ is either constant, increasing or decreasing. This figure only displays the final state of the scenarios. However, since their flow The deformations are restrained to vertical upward translations whose amplitudes are displayed in the first row. The first column shows the initial positions q 0 (X). We display in the center of the table the final cones q 1 (X) resulting from this nine configurations. One can compare, regarding to the growth process, the variations of the solutions with respect to the cone on the top left corner.

are built with rigid deformations, the foliations allow to retrieve each scenario (these foliations are enlightened by the meshes of the shapes and the color gradient). For a given initial position, when the amplitude of v increases, the growth accelerates and the final horn has a smooth tip. Conversely, when the amplitude decreases, the growth decreases and the final horn is sharp. The middle column illustrates that the choice of the initial position can have the effect of a time warping on a scenario.

2.4.

Theoretical study of the generative model. In this section, we formalize the previous construction in a general mathematical framework where a shape is defined as a mapping q from X to R d and the set of shapes, denoted B, is assumed to be a Banach space but also a shape space in the sense formalized by Arguillière [START_REF] Arguillère | The abstract setting for shape deformation analysis and lddmm methods[END_REF][START_REF] Arguillère | The general setting for shape analysis[END_REF]. This means that the group Diff(R d ) of smooth diffeomorphisms of R d continuously acts on B and that this action on any shape q is differentiable at Id ∈ Diff(R d ). It implies the existence of an operator ξ q : V → B modeling the infinitesimal action of smooth vector fields of R d . However, to model our growth dynamic, this infinitesimal action needs to depend on the time. Consequently, we cannot directly apply the general results of Arguillière's framework to generate new scenarios. We introduce a general time-dependent function

ξ : B × [0, 1] → L(V, B) ,
induced by the differentiation of the action of a group of diffeomorphisms, where V is a Hilbert space of vector fields on the ambient space and L(V, B) is the space of continuous linear operators from V to B.

Remark 2.4. An element ξ (q,t) (v) models the small variation of q induced by v. It belongs thus to the tangent space at q of B. Since B here is a linear space, we have ξ (q,t) (v) ∈ B.

The previous introduction (see equation ( 11)) invites to study general integral equations of the type [START_REF] Charlier | The fshape framework for the variability analysis of functional shapes[END_REF] ∀t ∈ [0, 1],

q t = q 0 + t 0 ξ (qs,s) (v s ) ds .
V can be more generally seen as a space of control. The introduction of the operator ξ highlights the linearity with respect to the control v. The specificity of this model lies on the fact that ξ depends on both time and spatial control.

Remark 2.5. Equation ( 14) is the integral version of the standard equation

qt = ξ (qt,t) (v t ) , (15)
with a given initial condition q 0 at time 0. However, we will see that the function t → ξ (qt,t) (v t ) is not regular enough to imply the existence of a derivative at all time. Theorem 2.1 below says that a solution q of the integral equation ( 14) admits a time derivative qt at almost all time and that this derivative is integrable. This property is shared by any absolutely continuous function (definition recalled below) .

Definition 2.2. A function

F : [0, 1] → B with values in a Banach space B is said absolutely continuous if there exists a function f ∈ L 1 ([0, 1], B) such that for any t ∈ [0, 1], F t = F 0 + t 0 f s ds.
The space of absolutely continuous functions with values in B will be denoted AC([0, 1], B).

We denote (H ξ 1 ) and (H V 1 ) few sets of regularity conditions on ξ and V . The (H ξ 1 ) conditions ensure the existence of time integrable controls on ξ and its derivative. The (H V 1 ) conditions allow to model both diffeomorphisms that tends to the identity at the infinity and linear deformations like rotations and translations.

(H ξ 1 ) (i) ξ t ∈ C 1 (B, L(V, B) for any t ∈ [0, 1] .
(ii) There exists c > 0 such that for any (q, t)

∈ B × [0, 1] , we have |ξ(q, t)| L(V,B) ≤ c(|q| B + 1) and |∂ q ξ(q, t)| L(B,L(V,B)) ≤ c .
where we use the shorthand notation

∂ q . = ∂ ∂q . (H V 1 ) (i) V ⊂ C 2 (R d , R d ) .
(ii) There exists c > 0 such that for any (

x, v) ∈ R d × V , we have |v(x)| R d ≤ c|v| V (|x| R d + 1) , |dv(x)| L(R d ,R d ) + |d 2 v(x)| L(R d ⊗R d ,R d ) ≤ c|v| V .
The main result of this section guarantees that any square integrable control v, i.e. v ∈ L 2 ([0, 1], V ) (also denoted L 2 V ), generates a scenario q on the complete time interval [0, 1].

Theorem 2.1. Under the (H ξ 1 ) and (H V 1 ) conditions, one can define the function Φ ξ that returns, for any initial condition q 0 ∈ B and any control v ∈ L 2 ([0, 1], V ), the unique solution q ∈ AC([0, 1], B) of the integral equation [START_REF] Charlier | The fshape framework for the variability analysis of functional shapes[END_REF].

Φ ξ : B × L 2 ([0, 1], V ) -→ AC([0, 1], B) (q 0 , v) -→ q : t → q 0 + t 0 ξ (qs,s) (v s )ds . Proof.
The definition of this integral is grounded in the theory of the Bochner integral (integration of functions with values in a Banach space) and the absolute continuity of the solution results from the integrability of v and the notion of Bochner-Lebesgue point (straightforward extension of Lebesgue point). See Theorem 4.2 and Proposition 4.2 in [START_REF]Kaltenmark Geometrical Growth Models for Computational Anatomy[END_REF] for more details.

We can now apply this general result for the growth dynamic that we introduced in Definition 2.1. When a space X is equipped with a temporal marker τ : X → [0, 1] and B is a space of mappings q : X → R d , the growth dynamic is written

ξ (q,t) (v) = (x → 1 1 τ (x)≤t v(q(x))) . (16) Proposition 2.1. If B = L ∞ (X, R d ), the (H V 1 )
conditions ensure that the operator ξ induced by the growth dynamic, given by equation ( 16), takes values in L(V, B) and satisfies the (H ξ 1 ) conditions.

Proof. See Section 4.4 in [START_REF]Kaltenmark Geometrical Growth Models for Computational Anatomy[END_REF].

Although the vector fields are spatially smooth, the indicator given by the temporal marker implies that ξ (q,t) (v) is not always continuous. The reference space to model the shapes is thus fixed to L ∞ (X, R d ) and the spatial regularity of q t demands a special attention. Nevertheless, Proposition 2.2 says that for any initial position q 0 ∈ C(X, R d ), the continuity is preserved. That is to say that the development q belongs to C([0, 1], C(X, R d )). Moreover, unlike the standard dynamic, a new characteristic appears here. The spatial regularity of q t depends on the temporal regularity of the vector field. Indeed, without more assumption on v, Proposition 2.3 exhibits a null set N ⊂ X where q t might not be differentiable. This set corresponds to the level lines of τ for times that are not Bochner-Lebesgue points of t → v t . Proposition 2.2 (Spatial continuity of the solution). If q 0 ∈ C b (X, R d ) (continuous and bounded) and τ ∈ C(X, [0, 1]) then q belongs to C([0, 1], C b (X, R d )).

Proposition 2.3 (Differentiability of the solution). Assume that q

0 ∈ C 1 b (X, R d ) and that τ ∈ C 1 (X, [0, 1]) is a submersion. Assume the (H V 1 )
conditions. (i) There exists a null set N ⊂ X of level lines of τ such that for any t ∈ [0, 1],

• the restriction of q t to the subset X c t = {x ∈ X | τ (x) > t} is of class C 1 and we have there dq t (x) = dq 0 (x) ,

• for any x ∈ X \ N such that τ (x) < t, q t : X → R d is differentiable at x and dq t (x) is the solution at time t of the integral equation

L t (x) = dq 0 (x) -v τ (x) (q τ (x) (x))dτ (x) + t τ (x) dv s (q s (x)) • L s (x) ds (17) defined on [0, 1]. (ii) Moreover, if v is continuous, i.e. v ∈ C([0, 1], V ), then for any t ∈ [0, 1[, q t is of class C 1 on the two level sets {x ∈ X | τ (x) > t} and {x ∈ X | τ (x) < t}. At last, q 1 belongs to C 1 (X, R d ).
Proof. See Section 4.4.2 in [START_REF]Kaltenmark Geometrical Growth Models for Computational Anatomy[END_REF].

Remark 2.6. To conclude on the regularity of the solution, one should remember that when v is continuous, the shape is of class C 1 on its two main components of active and inactive points as illustrated in Figure 10.

Figure 10. Horn in the middle of its development. The colors correspond to the level lines of the birth tag τ . The shape q t (X) is divided in two parts. The active part (in blue) of real points and the inactive part (from green to red) of fictional points that will progressively appear. At the boundary between these two parts, the shape admits two half tangent planes.

At last, to anticipate on the optimization problem presented in the next section, we need to study how a solution q = Φ ξ (q 0 , v) reacts to some variations of the initial condition and of the vector field. More precisely, consider two small variations δq 0 ∈ B and δv ∈ L 2 ([0, 1], V ) that define for > 0 a new set of parameters q 0 = q 0 + δq 0 and v = v + δv. We want to establish the link between the two solutions q = Φ ξ (q 0 , v) and q = Φ ξ (q 0 , v ).

Let us recall the definition of directional derivatives in Banach spaces.

Definition 2.3 (Gâteaux-derivative). Let f : E → F be an application between two Banach spaces E and F . Let be x 0 , δx 0 ∈ E. Define the application g : R → F given by g(h) = f (x 0 + hδx 0 ). If g is derivable at 0, we say that f is Gâteaux-differentiable at x 0 in the direction δx 0 and in this case, we denote and define the Gâteaux-derivative f (x 0 ; δx 0 ) . = g (0).

This definition leads to consider the function g : R → B defined with the previous notation by g( ) = q . Then Proposition 2.4 says that g is derivable at 0 and g (0) is the Gâteaux-derivative of Φ ξ at point (q 0 , v) in the direction (δq 0 , δv). Proposition 2.4. Consider ξ and Φ ξ as defined in Theorem 2.1 and such that ξ satisfies the (H ξ 1 ) conditions. For any (q 0 , v) ∈ B × L 2 V and any (δq 0 , δv) ∈ B × L 2 V , the Gâteaux-derivative Φ ξ (q 0 , v; δq 0 , δv) ∈ AC([0, 1], B) of Φ ξ is defined by the unique solution of the linearized equation

∀t ∈ [0, 1], δq t = δq 0 + t 0 ∂ξ ∂q (q s , s)(v s ) • δq s + ξ (qs,s) • δv s ds . (18)
Proof. See Theorem 4.4 in [START_REF]Kaltenmark Geometrical Growth Models for Computational Anatomy[END_REF].

In conclusion, if we denote δq = Φ ξ (q 0 , v; δq 0 , δv), we have for any small and at any time t ∈ [0, 1], q t ≈ q t + δq t .

Note at last again that δq is not derivable at all time but absolutely continuous.

3.

Optimal matching with a time dependent dynamic 3.1. Reconstitution of a growth scenario. Consider a longitudinal data set where each individual is represented by a sample of its evolution at a finite number of times. The main problem addressed in this paper is to estimate, for any individual, its complete evolution on the time interval [0, 1]. Additionally, we assume that all the individuals share a common growth pattern.

Each evolution can thus be represented by a growth scenario parametrized by a common biological coordinate system (X, τ ) . Consider a target scenario given by a collection of shapes (S tar i ) i at a finite number of intermediate times (t i ) i ⊂ [0, 1] (with max{t i , i} = 1) . With the notation of Theorem 2.1, we aim thus to find q ∈ C([0, 1], B) in the image of the generating function Φ ξ such that q ti (X ti ) ≈ S tar i for all (t i ) i . The discrepancy between the data and a solution q ∈ C([0, 1], B) is estimated at the different times t i with a data attachment term A of the form

n i=1 d(S ti , S tar i ) 2 ,
where the shape S ti is induced by q ti and d is a distance that depends on the type of the data. To simplify the problem, we will assume throughout this section that n = 1 so that A only depends on q 1 . We consider thus a functional A : B → R + , of class C 1 . In our experiments, A will measure the difference between the image of q 1 and S tar with the square norm of an Reproducing Kernel Hilbert Space (RKHS) modeling a current space or a varifold space (see Section 5).

We can finally precise our inexact registration problem where we use the most common regularization term. The sought-after scenario is represented by a pair (q 0 , v) that minimizes the energy ( 19)

E(q 0 , v) = 1 2 1 0 |v t | 2 V dt + A(q 1 ) ,
where v is a time-varying vector field that belongs to L 2 ([0, 1], V ), q 0 ∈ B is the initial position and q = Φ ξ (q 0 , v) ∈ C([0, 1], B) is the development generated by (q 0 , v) as previously introduced (see Theorem 2.1). We recall that given ξ :

B × [0, 1] → L(V, B), we have at almost any time t ∈ [0, 1], (20) qt
= ξ(q t , t)(v t ) .
Remark 3.1. Let us recall that the target S tar will usually induce a good estimation of the initial position q 0 thanks to a biological prior that restricts the area where the new points of the scenario can appear (see Figure 7) . We consider yet q 0 as a variable to optimize this estimation, if required. See Section 3.2.3 and ?? for more details.

3.2.

Expression of the gradient via the momentum map. The minimization of E is achieved by a gradient descent. A prerequisite to establish the gradient of E is the introduction of a new variable p called the momentum.

Momentum and momentum map.

Proposition 3.1 (Existence of the Momentum). Under the (H ξ 1 ) conditions (defined for Theorem 2.1), the momentum curve p ∈ AC([0, 1], B * ) associated to q = Φ ξ (q 0 , v) is the unique solution of

ṗt = - ∂ξ ∂q (q t , t)(v t ) * • p t , for almost all t ∈ [0, 1] , (21)
with the final condition

p 1 = -dA(q 1 ) ∈ B * . ( 22
)
The uniqueness of the solutions q and p allows us to denote hereafter (q, p) = Φ ξ (q 0 , v).

Proof. Under the (H ξ 1 ) conditions, t → ∂ξ ∂q (q t , t)(v t ) belongs to L 2 ([0, 1], L(B)). The adjoint operator curve t → ∂ξ ∂q (q t , t)(v t ) * belongs thus to L 2 ([0, 1], L(B * )). The existence and uniqueness of p are given by the linear Cauchy-Lipschitz formulation given in Corollary 4.2 in [START_REF]Kaltenmark Geometrical Growth Models for Computational Anatomy[END_REF].

As in the LDDMM framework (with the standard dynamic), the momentum is a key variable to solve the registration problem presented here. The optimal vector field will be parametrized by the trajectory q, its momentum p and the time variable through a function, J ξ , usually called the momentum map. Its definition is based on the linearity of ξ with respect to v. Proposition 3.2 (Definition of the Momentum Map). The momentum map associated to ξ is the application J ξ defined under the (H ξ 1 ) conditions by

J ξ : B × B * × [0, 1] -→ V * (q, p, t)
-→ ξ * (q,t) • p . Explicitly, we have for any (q, p, t) ∈ B × B * × [0, 1] and any v ∈ V (J ξ (q, p, t) | v) = p ξ (q,t) (v) .

Proof. Let us verify that J ξ is well defined. Under the (H ξ ) conditions, there exists c > 0 such that for any (q, p, t)

∈ B × B * × [0, 1], |ξ (q,t) | L(V,B) ≤ c(|q| B + 1). Hence, |J ξ (q, p, t)| V * ≤ c(|q| B + 1)|p| B * (23)
and J ξ takes indeed its values in V * . Remark 3.2. Equation (23) implies that for any continuous curve (q, p)

∈ C([0, 1], B × B * ), t → J ξ (q t , p t , t) is square integrable on [0, 1], i.e. t → J ξ (q t , p t , t) ∈ L 2 ([0, 1], V * ) . ( 24
)
The regularity with respect to time of this momentum map will play an important role in the next sections. At this stage, we start with the following proposition. Proposition 3.3. Under the (H ξ 1 ) conditions, the momentum map is C 1 with respect to its first two variables.

Proof. Under (H ξ 1 ), q → ξ (q,t) is of class C 1 for any t ∈ [0, 1], and p → (p | δq) is smooth for any δq ∈ B.

Expression of the gradient.

Theorem 3.1 (Expression of the gradient). Assume the (H ξ 1 ) and (H V 1 ) conditions and consider A : B → R of class C 1 . Let be (q 0 , v) ∈ B × L 2 ([0, 1], V ) and (q, p) = Φ ξ (q 0 , v) the trajectory defined by (q 0 , v). Consider the energy

E(q 0 , v) = 1 2 1 0 |v t | 2 V dt + A(q 1 ) .
Then the Gâteaux-derivative of the energy at any (q 0 , v) is given in any direction (δq 0 , δv) ∈ B×L 2

V by

E (q 0 , v); (δq 0 , δv) = (p 0 | δq 0 ) + 1 0 ( v t -J ξ (q t , p t , t)| δv t ) dt .
Hence, given q 0 ∈ B, any minimizer

v * ∈ L 2 V of E satisfies at almost any time t ∈ [0, 1] (25) v * t = K V J ξ (q t , p t , t) ,
where K V : V * → V is the canonical isomorphism (V being a Hilbert space).

Proof. Consider and denote δq = Φ ξ (q 0 , v; δq 0 , δv) the Gateaux-derivative of Φ ξ given in Proposition 2.4. Since δq and p are absolutely continuous, t → (p t | δq t ) is also absolutely continuous. We deduce from the expressions of qt , ṗt , and δ qt given by equations ( 20), [START_REF] Figalli | A new transportation distance between non-negative measures, with applications to gradients flows with Dirichlet boundary conditions[END_REF], and ( 18) that

(p 1 | δq 1 ) = (p 0 | δq 0 ) + 1 0 d dt (p t | δq t ) dt = (p 0 | δq 0 ) + 1 0 (p t | δ qt ) + ( ṗt | δq t ) dt = (p 0 | δq 0 ) + 1 0 p t ∂ξ ∂q (q t , t)(v t ) • δq t + ξ (qt,t) • δv t - ∂ξ ∂q (q t , t)(v t ) * • p t δq t dt = (p 0 | δq 0 ) + 1 0 p t ∂ξ ∂q (q t , t)(v t ) • δq t + ξ (qt,t) • δv t -p t ∂ξ ∂q (q t , t)(v t ) • δq t dt = (p 0 | δq 0 ) + 1 0 ξ * (qt,t) • p t δv t dt = (p 0 | δq 0 ) + 1 0 ( J ξ (q t , p t , t) | δv t ) dt ,
where J ξ is the momentum map introduced in Proposition 3.2. As noted in Remark 3.2, this integral is finite under the (H ξ 1 ) conditions. Define then

δE = 1 0 v t , δv t V dt + (dA(q 1 ) | δq 1 ) . (26) It follows that δE = 1 0 v t , δv t V dt -(p 1 | δq 1 ) = (p 0 | δq 0 ) + 1 0 v t , δv t V -( J ξ (q t , p t , t) | δv t ) dt = (p 0 | δq 0 ) + 1 0 v t -K V J ξ (q t , p t , t), δv t V dt . ( 27
)
δE is finite and if e( ) = E(q 0 , v ), the Gâteaux-derivative of E is then equal to E (q 0 , v); (δq 0 , δv) = e (0) = δE .

Note that this theorem implies neither the existence of local minimizers of the energy E nor the uniqueness. Remark 3.3 (Time regularity of an optimal vector field). An important issue with time-varying dynamics is that t → ξ t has no reason to be continuous. Consider for example X = [0, 1], τ = Id and v ≡ y, with y ∈ (R d ) * , a constant vector field. Then we have with the growth dynamic for any

q in B = L ∞ ([0, 1], R d ) and any t < t in [0, 1] |ξ (q,t) (v) -ξ (q,t ) (v)| ∞ = | x → 1 1 t<x≤t y | ∞ = |y| R d ,
so that ξ (q,t ) does not tend to ξ (q,t) when t tends to t. Consequently and as mentioned before, we have a priori no control on the time regularity of the momentum map and thus on the continuity of an optimal vector field.

3.2.3.

Algorithm for the gradient descent. Algorithm 1 explicits a gradient descent very similar to the LDDMM model to construct a minimizer v * . The main difference in practice is to trace at each discrete time t i ∈ [0, 1] the set of active points.

Additionally, Theorem 3.1 allows to optimize the initial condition q 0 if necessary (step 6 in the algorithm). Typically, if q 0 is partially known and a reconstruction has been guessed, we can Algorithm 1 Gradient descent on v 1 -Initiate v 0 ∈ L 2 ([0, 1], V ) at zero. Then for any n ∈ N, given q 0 and v n , 2 -Compute q n ∈ C([0, 1], (R d ) k ) the path generated by v n ∈ L 2 ([0, 1], V ). 3 -Compute p n 1 = -∇A(q n 1 ) and integrate it backward to construct p n ∈ C([0, 1], (R d ) k ). 4 -Compute at any time the gradient at v n t : δv n t = v n t -K V J (q n t , p n t , t). 5 -Update the vector field by v n+1 = v n + δv n for a small > 0. 6 -(Optional) Update q n 0 by q n+1 0 = q n 0 + p n 0 for a small > 0.

optimize it under some constraints (for example, inside a subset of the ambient space) . This optimization should of course be controlled, otherwise the initial condition would just straightforwardly tend to the target. See an example in Section ??.

Hamiltonian framework.

The Hamiltonian approach leads in the standard LDDMM framework to an elegant characterization of the solutions to register two shapes. The problem can then be rewritten as an optimization over the initial momentum p 0 that what we refer to as the shooting method. Another direct result of this approach is the norm conservation of the optimal vector field. With time-varying dynamics like the growth dynamic, we need to rewrite the complete Hamiltonian framework with an important subtlety: the Hamiltonian functions depend on the time variable.

3.3.1.

Hamiltonian formulation. The central Theorem 3.1 that gives the gradient of the energy E, says that given an initial condition q 0 in B, any local minimizer v ∈ L 2

V of E(q 0 , •) must satisfy at almost any time t ∈ [0, 1] the equation

0 = J ξ (q t , p t , t) -v t (28) = ξ * (qt,t) • p t -v t
, where (q, p) = Φ ξ (q 0 , v) are the spatial mapping and the momentum curve associated to q 0 and v. This leads to introduce the following Hamiltonian function [START_REF] Gris | Modular approach on shape spaces, sub-Riemannian geometry and computational anatomy[END_REF] H :

B × B * × V × [0, 1] -→ R (q, p, v, t) -→ (J ξ (q, p, t) | v) -1 2 |v| 2
V , for which we have for any (q, p, v, t)

∈ B × B * × V × [0, 1] ∂H ∂v (q, p, v, t) = J ξ (q, p, t) -v .
An optimal trajectory (q, p) = Φ ξ (q 0 , v) associated to a local minimizer v ∈ L 2 V is thus a local extrema of this Hamiltonian at almost any time. Indeed, equation ( 28) is equivalent to ∂H ∂v (q t , p t , v t , t) = 0 which is nothing but a weak form of the Pontryagin Maximum Principle [45].

This allows to define the reduced Hamiltonian as follows:

(30)

H r : B × B * × [0, 1] -→ R (q, p, t)
-→ max v∈V H(q, p, v, t) .

and we have explicitly (31) H r (q, p, t) = H(q, p, J ξ (q, p, t), t) .

We can now rewrite Theorem 3.1 and characterize the optimal solutions with the reduced Hamiltonian.

Theorem 3.2. Assume the (H ξ 1 ) and (H V 1 ) conditions and consider A : B → R of class C 1 . Let be v ∈ L 2 ([0, 1], V ), q 0 ∈ B and denote (q, p) = Φ ξ (q 0 , v) the associated trajectory.

Then the Gâteaux-derivative of the energy

E(q 0 , v) = 1 2 1 0 |v t | 2 V dt + A(q 1 ) is given in any direction δv ∈ L 2 ([0, 1], V ) by (32) ∂E ∂v (q 0 , v) • δv = 1 0 - ∂H ∂v (q t , p t , v t , t) | δv t dt ,
where for any t ∈ [0, 1],

H(q t , p t , v t , t) = (J ξ (q t , p t , t) | v t ) - 1 2 |v t | 2 V .
Moreover, if v locally minimizes E, then the trajectory (q, p) ∈ AC([0, 1], B × B * ) satisfies at almost any time the following Hamiltonian system [START_REF]Kaltenmark Geometrical Growth Models for Computational Anatomy[END_REF] (H)

   qt = ∂Hr ∂p (q t , p t , t) ṗt = -∂Hr ∂q (q t , p t , t) ,
where H r (q t , p t , t) = max u∈V H(q t , p t , u, t) = H(q t , p t , v t , t) .

Proof. Equation ( 32) is a direct consequence of Theorem 3.1. The partial derivatives of H r are given for any (q, p, t)

∈ B × B * × [0, 1] by ∂H r ∂q (q, p, t) = ∂H ∂q (q, p, J ξ (q, p, t), t) + ∂H ∂v (q, p, J ξ (q, p, t), t) • ∂J ξ ∂q (q, p, t) = ∂H ∂q (q, p, J ξ (q, p, t), t) , ∂H r ∂p (q, p, t) = ∂H ∂p (q, p, J ξ (q, p, t), t) + ∂H ∂v (q, p, J ξ (q, p, t), t) • ∂J ξ ∂p (q, p, t)
= ∂H ∂p (q, p, J ξ (q, p, t), t) . Now, if v locally minimizes E and (q, p) ∈ AC([0, 1], B × B * ) denotes the associated optimal trajectory, we have for almost any t ∈ [0, 1], v t = J ξ (q t , p t , t) and

∂H r ∂q (q t , p t , t) = - ∂ξ ∂q (q t , t)(v t ) * • p t = -ṗt , ∂H r ∂p (q t , p t , t) = ξ qt,t (v t ) = qt .
Remark 3.4. From equation (31), we retrieve a similar result from the standard LDDMM framework that is for any optimal trajectory (q, p)

∈ AC([0, 1], B × B * ) of control v ∈ L 2 V and at almost any time t ∈ [0, 1] (34) H r (q t , p t , t) = 1 2 |J ξ (q t , p t , t)| 2 V * = 1 2 |v t | 2 V .
The important novelty here is that this Hamiltonian is not constant in time and the norm of the optimal control v follows the time variations of H r . From ∂H r ∂t (q, p, t) = ∂H ∂t (q, p, J ξ (q, p, t), t) + ∂H ∂v (q, p, J ξ (q, p, t), t) • ∂J ξ ∂t (q, p, t) = ∂H ∂t (q, p, J ξ (q, p, t), t) , we deduce that the derivative of H r with respect to t equals dH r dt (q t , p t , t) = ∂H r ∂q (q t , p t , t)

• qt + ∂H r ∂p (q t , p t , t) • ṗt + ∂H r ∂t (q t , p t , t) = -ṗt • qt + ṗt • qt + ∂H ∂t (q t , p t , v t , t) = ∂J ξ ∂t (q t , p t , t) v t .
We will see in Section 4 that with the growth dynamic this partial derivative of the momentum map is link to the mass of the new points. Hence, when there is no creation at time t, we retrieve the classic LDDMM case and the norm of the vector field is constant. Otherwise, when new points appear, this norm might vary and usually increase. Moreover, since the norm of v varies, this new model on time-varying dynamics does no longer generate geodesics of the group of deformations. This results from the fact that with the growth dynamic for example, the final shape q 1 does not depend anymore only on the final deformation φ 1 generated by v but on the complete curve t → φ t (or at least at every time when creation occurs).

3.3.2.

Existence and uniqueness of the solutions of the reduced Hamiltonian system. Among all the solutions of the Hamiltonian system (H), defined by equation [START_REF]Kaltenmark Geometrical Growth Models for Computational Anatomy[END_REF], generated from any initial condition (q 0 , p 0 ) ∈ B × B * is the sought-after trajectory (q, p) ∈ C([0, 1], B) × C([0, 1], B * ) associated to an optimal time-varying vector field that minimizes E. Therefore, the optimization problem on v can be replaced by an optimization problem on p 0 (Shooting Method ). Before presenting this new viewpoint, we need to guarantee the existence and uniqueness of a solution of (H) for any initial condition (q 0 , p 0 ). It requires to extend the regularity conditions on ξ and V as follows

(H ξ 2 ) (i) For any t ∈ [0, 1], ξ t ∈ C 2 (B, L(V, B)) .
(ii) There exists c > 0 such that for any (q, t) ∈ B × [0, 1], we have

∂ 2 ξ ∂q 2 (q, t) op ≤ c . (H V 2 ) (i) V ⊂ C 3 (R d , R d ) .
(ii) There exists c > 0 such that for any (x, v) ∈ R d × V , we have

|v(x)| R d ≤ c|v| V (|x| R d + 1) , |dv(x)| ∞ + |d 2 v(x)| ∞ + |d 3 v(x)| ∞ ≤ c|v| V .
The local existence of the solutions lies on the spatial regularity of the momentum map, i.e. with respect to q and p. The most challenging part is to show the global existence. It requires to control the reduced Hamiltonian with respect to time and likewise the momentum map, see Remark 3.4. We introduce therefore this last condition:

(H J ξ ) (i) For any (q, p) ∈ B × B * , t → J ξ (q, p, t) is piecewise constant on [0, 1].
(ii) For any local solution (q, p) ∈ C(I, B × B * ) of the reduced Hamiltonian system, there exists c > 0 such that for any t ∈ I,

lim s→t + J ξ (q s , p s , s) -lim s→t -J ξ (q s , p s , s) V * ≤ c|p 0 | B * (|q 0 | B + 1) .
Remark 3.5. We will see in Section 4.2 that this condition is satisfied with the growth dynamic in cases where the shapes are given by a discrete mesh. When the shapes are submanifolds, the most natural condition should write t → J ξ (q, p, t) is absolutely continuous and t 0 ∂J ∂t (q s , p s , s)

V * ds ≤ c|p 0 | B * (|q 0 | B + 1) .
With the growth dynamic, these conditions require to control the amount of creation, i.e. the total mass of the coordinate space X, which is thus assumed compact.

Theorem 3.3 (Global Solutions of the Reduced Hamiltonian System ). Assume the (H ξ 1 ), (H ξ 2 ), (H J ξ ) and (H V 2 ) conditions. Then for any initial condition (q 0 , p 0 ) ∈ B × B * , the reduced Hamiltonian system (H) admits a unique solution (q, p) ∈ AC([0, 1], B × B * ).

Moreover, there exists an increasing function ϕ : R + → R + such that for any (q 0 , p 0 ) ∈ B × B * and any t ∈ [0, 1], we have

|q t | B + |p t | B * ≤ ϕ(|q 0 | B + |p 0 | B * ) . (35)
Proof. See Section 5.4 and 5.2 in [START_REF]Kaltenmark Geometrical Growth Models for Computational Anatomy[END_REF].

Remark 3.6. The theorem gives an interesting property of the solutions of the reduced Hamiltonian. They are locally bounded with respect to the initial condition.

3.3.3. Shooting method. This section presents a formal approach of the shooting method. The following results are proved in Section 5.5 and 5.6 of [START_REF]Kaltenmark Geometrical Growth Models for Computational Anatomy[END_REF].

We saw in the previous section that the curve (q, p) generated by an optimal vector field can be retrieved as a solution of the reduced Hamiltonian system. These solutions are parametrized by an initial position q 0 and an initial momentum p 0 . Instead of playing with the vector field as the parameter of the set of curves, we can thus consider the initial momentum as the new parameter that can be optimized.

Denote y = Ψ ξ (q 0 , p 0 ) the unique solution of the reduced Hamiltonian system associated to the initial condition (q 0 , p 0 ) ∈ B × B * . We have thus y ∈ AC([0, 1], B × B * ) and at any time t ∈ [0, 1], y t = (q t , p t ) and

y t = y 0 + t 0 h(y s , s) ds ,
where h is the symplectic gradient of H r with respect to (q, p) defined as follows [START_REF] Lee | Introduction to Smooth Manifolds[END_REF]

h : (B × B * ) × [0, 1] -→ B × B * ((q, p), t) -→   ∂Hr ∂p (q, p, t)
-∂Hr ∂q (q, p, t)

  .
With this notation, we introduce a new expression of the energy Ê(y 0 ) = 1 0 Ĉ(y t , t) dt + Â(y 1 ) , where Ĉ(y t , t) = 1 2 |J ξ (q t , p t , t)| 2 V is equivalent to the old cost function and Â(y 1 ) = A(q 1 ) is equivalent to the old data attachment term. Therefore, if an initial momentum p 0 and a timevarying vector field v generate the same solution (q, p) ∈ AC([0, 1], B×B * ), i.e. v t = K V J ξ (q t , p t , t) a.e., then the two respective energies are equal

Ê(y 0 ) = Ê(q 0 , p 0 ) = E(q 0 , v) . (37)
The method to explicit the gradient of the energy is the same as before, reduced to two main steps and its conclusion as follows:

• As the momentum p is the covariable of q in the previous framework, the first step here is to define the covariable of y. We introduce thus z 1 = -d Â(y 1 ) ∈ (B × B * ) * and we integrate it backward through the equation

(38) żt = ∂ Ĉ ∂y (y t , t) - ∂h ∂y (y t , t) * • z t .
• The second step is to establish that the Gâteaux-derivative Ψ ξ (y 0 ; δy 0 ) is given by δy t = δy 0 + t 0 ∂h ∂y (y s , s) • δy s ds .

• Then we can write

E (y 0 ; δy 0 ) = 1 0 ∂ Ĉ ∂y (y t , t) δy t dt + (d Â(y 1 ) | δy 1 ) = 1 0 ∂ Ĉ ∂y (y t , t) δy t dt -(z 1 | δy 1 ) = 1 0 ∂ Ĉ ∂y (y t , t) δy t dt -(z 0 | δy 0 ) + 1 0 ( żt | δy t ) -(z t | δ y t ) dt = -(z 0 | δy 0 ) + 1 0 ∂ Ĉ ∂y (y t , t) -żt - ∂h ∂y (y t , t) * • z t δy t dt = -(z 0 | δy 0 ) .
At last, if we write at any time t the covariable z t as (Q t , P t ) ∈ B * × B * * , we get more explicitly

E (q 0 , p 0 ); (δq 0 , δq 0 ) = -(Q 0 | δq 0 ) -(P 0 | δp 0 ) . (39)
The Gâteaux-derivative of the energy has thus a particularly simple expression leading to a new algorithm of gradient descent. An interest of this approach is to parametrize the solution with variables of smaller dimension paving the way for a statistical analysis.

Algorithm 2 Gradient descent on p 0 1 -Given q 0 0 ∈ B, initialize p 0 0 ∈ B at zero. Then for any n ∈ N, given q n 0 and p n 0 : 2 -Integrate forward with the Hamiltonian system [START_REF] Lee | Introduction to Smooth Manifolds[END_REF]

to get (q n , p n ) ∈ C([0, 1], B 2 ). 3 -Compute Q n 1 = -dA(q n 1 )
, defined P n 1 = 0. 4 -Integrate backward with the second order Hamiltonian system [START_REF] Michor | An overview of the Riemannian metrics on spaces of curves using the Hamiltonian approach[END_REF] to get (Q n , P n ) ∈ C([0, 1], B 2 ). 4 -Update p n 0 by p n+1 0 = p n 0 + P n 0 for a small > 0. 5 -(Optional) Update q n 0 by q n+1 0 = q n 0 + Q n 0 for a small > 0.

Applications with the growth dynamic

We will now apply the previous results in the setting of mappings from a biological coordinate system (X, τ ) and with the growth dynamic that defines the operator ξ :

B × [0, 1] → L(V, B) by ξ (q,t) (v) = x → 1 1 τ (x)≤t v(q(x)
) . In order to remember that ξ is fixed throughout this section, the momentum map will be denoted J instead of J ξ , likewise with Φ. 4.1. Momentum variables with a discrete coordinate space. We assume that the coordinate space X is given as a finite set of k points with a mesh. At any time t ∈ [0, 1], q t is an element of B = (R d ) k with a mesh. One can easily show that under the (H V 2 ) conditions, the (H ξ 1 ) and (H ξ 2 ) conditions are satisfied (see Proposition 4.3 and 5.13 in [START_REF]Kaltenmark Geometrical Growth Models for Computational Anatomy[END_REF]). 4.1.1. The momentum. The general definition of the momentum curve p and its evolution are given in Proposition 3.1. Here, with the Riez representation theorem, p can be seen as an element of C([0, 1], (R d ) k ). At time t = 1, p 1 is given by definition as the gradient of the data attachment term A dA(q

1 ) • δq 1 = x∈X ∇ q1(x) A(q 1 ), δq 1 (x) R d
and p 1 can thus be parametrized by X as follows

p 1 (x) = ∇ q1(x) A(q 1 ) , so that (p 1 | δq 1 ) = p 1 , δq 1 (R d ) k .
This pointwise expression is conserved by the backward integration and we have at any time t ∈ [0, 1] for any

y ∈ B = (R d ) k (p t | y) = x∈X p t (x), y(x) R d .
Remark 4.1. In this configuration, the behaviors of q and p over time share a common pattern. For any x in X, the point q t (x) and its momentum p t (x) are both static when the coordinate x is inactive, i.e. t < τ (x), and jointly moving once x has been activated. Indeed, their dynamics are explicitly given for any x ∈ X and at almost any time t ∈ [0, 1] by

qt (x) = 1 1 τ (x)≤t v t (q t (x)) ṗt (x) = -1 1 τ (x)≤t dv t (q t (x)) T • p t (x) . (40)
4.1.2. Expression of the momentum map. As we saw in Theorem 3.1, optimal vector fields are given by the momentum map J . If (q, p) ∈ C([0, 1], B × B * ) is an optimal trajectory associated to a vector field v ∈ L 2 ([0, 1], V ) then v satisfies at almost all time:

(41)

v t = K V J (q t , p t , t) ,
where K V : V * → V the canonical isomorphism. For any (q, p, t)

∈ (R d ) k × (R d ) k × [0, 1] and any v ∈ V , we have J (q, p, t) | v . = p, ξ (q,t) (v) (R d ) k = p, 1 1 τ ≤t v • q (R d ) k = x∈X,τ (x)≤t p(x), v(q(x
)) R d .

Equivalently, J can be written

J (q, p, t) = x∈X,τ (x)≤t δ p(x) q(x)
,

where for any (x, y) ∈ (R d ) 2 and any v ∈ V , the Dirac delta δ y x ∈ V * is defined by

δ y x (v) = y, v(x) R d .
When V is a RKHS with a kernel denoted k V , we have an explicit expression of v = K V J (q, p, t). For any y ∈ R d , v(y) = x∈X,τ (x)≤t k V (q(x), y)p(x) .

In our numerical experiments, we will use a Gaussian kernel of scale σ > 0 for which we get

v(y) = x∈X,τ (x)≤t exp - |q(x) -y| 2 2σ p(x) .
Remark 4.2. The reader familiar with the LDDMM framework will recognize all these equations. The difference with the classical model only resides in the addition of the indicator function. At any time t, the coordinate space X and likewise the shape q t (X) are divided in two parts:

(1) the active coordinates, resp. the associated points of the shape that actually exist

X t = {x ∈ X | τ (x) ≤ t} , resp. q t (X t ) , (42)
(2) the inactive coordinates (the complementary set of X t ), resp. the anticipated points. At time t, the vector field v t carries only the active points of the shape. Therefore, it is natural to obtain an optimal vector field localized on the positions q t (x) and driven by the momenta p t (x) but only for the active coordinates x ∈ X t . See Figure 11.

Figure 11. The active part of the shape modeling the horn is blue. The inactive part of anticipated points goes from green to red. The arrows are the respective momenta p t (x) at points q t (x). 4.2. Specific behavior of the momentum map with the growth dynamic. The sets X t and X {t} of active coordinates and new coordinates at time t, recalled by ( 42) and (43) (or see Section 2.1), play an important role to understand the construction of an optimal vector field. Since X is a finite set, the image of τ is also a finite subset of [0, 1] that we denote {t 0 , t 1 , . . . , t n } . In other words, the times at which new points in the mesh are activated is a finite subset of [0, 1]. It follows that t → X t is a piecewise constant function on [0, 1] and the jumps are due to the sets

(X {ti} = {x ∈ X | τ (x) = t i }) i (43)
of new coordinates that are activated at times (t i ) i .

Remark 4.3 (Jumps of the momentum map). Given q, p ∈ (R d ) k , the function

t → J (q, p, t) = Xt δ p(x) q(x)
is constant on any interval [t i , t i+1 [. The jump at time t i+1 is given by

x∈X {t i+1 } δ p(x) q(x) .
This jump is thus due to the contribution of the new coordinates that are activated at time t i+1 . Remark 4.4 (Continuity of an optimal vector field). Proposition 3.3 states that the momentum map is of class C 1 with respect to its first two variables. Then we deduce from the previous remark that given a trajectory (q, p) ∈ C([0, 1], (R d ) k × (R d ) k ) solution of the Hamiltonian system, the time-varying vector field t → K V J (q t , p t , t) is continuous on any interval [t i , t i+1 [ where t i and t i+1 are two consecutive values of τ (X). The jump at time t i+1 is given by

x∈X {t i+1 } δ pt i+1 (x) qt i+1 (x) .
Since the position and momentum of any inactive coordinate remain static (see Remark 4.1), we have at any time t and for any new coordinate x ∈ X {t} , q t (x) = q 0 (x) and p t (x) = p 0 (x). The jump at time t i+1 is thus equal to

x∈X {t i+1 } δ p0(x) q0(x) .
Therefore, all jumps depend only on q 0 , p 0 and τ which guarantees the control on the evolution with respect to time of the momentum map as required by the (H J ξ ) conditions (see Remark 3.5) .

Inside these intervals [t i , t i+1 [, the evolution of K V J is the same as in the LDDMM framework. Otherwise, the jumps result from the extension of the support of K V J with the set of points that progressively appear.

The previous analysis leads to suppose that the norm of an optimal vector field increases with the apparition of new points. Although we cannot explicitly show it, we can yet propose an upper bound of this norm by an increasing function of time. With a discrete coordinate space, we saw that the momentum map is given for any (q, p, t) ∈ B × B * × [0, 1] and any v ∈ V by

J (q, p, t) | v = x∈Xt p(x), v(q(x)) R d (44) (note that B * = B = (R d ) k ) .
It results that for any continuous curve (q, p) ∈ C([0, 1], B × B * ), there exists, under the (H V 1 ) conditions, a constant c > 0 such that for any t ∈ [0, 1],

|K V J (q t , p t , t)| V ≤ c|p| ∞ (1 + |q| ∞ ) x∈Xt 1 ≤ c|p| ∞ (1 + |q| ∞ )|X t | , (45)
where |X t | is the cardinality of X t .

Data attachment term

5.1. RKHS of currents and varifolds. The inexact matching setting with the presence of an attachment term is justified by the fact that the curve q t (X) is not intended to accurately describe the real data as we do not want to capture too small details that could result from very specific characteristics of an individual or from noise. The flexibility given by the group of deformations and the precision of the attachment term (mostly the typical scale of these two elements) will determine the level of detail of the model and the independence with respect to local noises. For shapes like curves or surfaces A can be chosen as the distance on currents presented in [23,[START_REF] Glaunès | Diffeomorphic matching of distributions: A new approach for unlabelled point-sets and sub-manifolds matching[END_REF] or the distance on unoriented varifolds introduced in [13]. In our experiments, the shapes will be considered as oriented varifolds, more recently presented in [START_REF] Kaltenmark | Charon A general framework for curve and surface comparison and registration with oriented varifolds[END_REF]. Currents and varifolds enable to measure the discrepancy between shapes regardless of the parameterization.

Let us recall how to identify shapes to linear forms on some spaces of test functions. Reproducing Kernel Hilbert Spaces (RKHSs) allow then to compute very efficiently distances between these shapes. The focus is set here on triangulated surfaces embedding in R 3 .

Locally, a surface S can be encoded by a position and a unit normal vector, i.e. a pair (x,

-→ N x ) ∈ R 3 × S 2 .
A test function is then a function ω : R 3 × S 2 → R and the pair (x,

-→ N x ) is associated to the Dirac δ -→ Nx x defined by δ -→ Nx x (ω) = ω(x, -→ N x ).
A small triangle is approximated by a linear form of the type µ x = x δ -→ Nx x where x is the center of the triangle, -→ N x its normal unit vector and x its area. Given a space of test functions W , this representative µ x : W → R is applied to any ω ∈ W by µ x (ω) = x ω(x, -→ N x ) . The union of two triangles is then represented by a sum µ x + µ y . More generally, a triangulated surface can be approximated by a finite sum (46) S ≈ µ = The construction of metrics via RKHSs becomes from there rather simple in practice. It can be induced by the choice of two real positive kernels on the ambient space E = R n and on the set of tangential data T (general situation). These two kernels generate the RKHS W , i.e. the set of test functions. A kernel k E measures the distance between the positions of two infinitesimal shapes (triangles) and a kernel k T measures the distance between their respective tangential data (here, a unit normal vector with eventually the orientation).

We can now explicit a scalar product between two shapes represented in W by two sums of Diracs as in equation (47). Thanks to the RKHS properties, we have

δ -→ Nx x , δ -→ Ny y W = δ -→ Nx x k E (•, y)k T (•, -→ N y ) = k E (x, y)k T ( -→ N x , -→ N y ) . (48)
The scalar product between two shapes S ≈ µ S = x x δ -→ Nx x and S ≈ µ S = y y δ -→ Ny y is then deduced by linearity:

µ S , µ S W = x y x y k E (x, y)k T ( -→ N x , -→ N y ) . (49)
Finally, we return now to the definition of a data attachment term to compare a target shape to a deformed source shape. Once we have fixed a RKHS W , we denote µ tar ∈ W the representative of the target shape and µ v ∈ W the representative of the final state of the solution generated by a vector field v ∈ L 2 ([0, 1], V ). The data attachment term is then given by

A(q 0 , v) = 1 2 µ v -µ tar 2 W . (50)
Note that this brief overview is common to currents and varifolds. The only differences between these models lie in the choice of the RKHS's kernels. The metric in our experiments is mainly the product of two Gaussian kernels:

δ -→ Nx x , δ -→ Ny y W . = exp - |x -y| 2 R 3 2σ 2 exp - | -→ N x - -→ N y | 2 R 3 2σ 2 N -→ N x , -→ N y R 3 . (51)
where σ > 0 and σ N > 0 are two scalar parameters. This kernel is sensitive to the orientation of -→ N x and -→ N y and it defines thus a metric on oriented varifolds. Note that when σ N tends to +∞ it retrieves the scalar product on currents.

The advantage of varifolds in our situation is their non linearity with respect to the tangent space that enables them to capture sharp tails as the tip of horns. Moreover, since shapes like horns are easily oriented, there is no reason to call for non-oriented varifolds.

Additional landmark.

For the specific case of the horn, we add a penalization on the tip of the horns. This tip is the point image by q of the first leaf X {0} . Consider any x ∈ X {0} , the new attachment term can then be defined by

A(q 0 , v) = c head 2 |q 1 (x) -ŷtar | 2 R d + c data 2 µ v -µ tar 2 W ,
where ŷtar is the tip of the target horn. Without this additional landmark, the tips of both horns have no particular reason to be matched together. This single landmark gives a strong input toward the natural structure of the temporal foliation.

5.3.

Intermediate times in the input data. When intermediate states are available for real data, their integration to the model should improve the matching of the whole evolution. Assume that a target scenario is given by a collection (S tar i ) i of shapes at a finite number of intermediate times (t i ) i ⊂ [0, 1] (with max{t i , i} = 1). Let v ∈ L 2 ([0, 1], V ) be a vector field that generates a scenario t → S v t . Given a distance d on shapes, the discrepancy between the two scenarios can be estimated at the different times t i by an attachment term A of the form

A(q 0 , v) = n i=1 d(S v
ti , S tar i ) .

6.1.2. Deformation spaces. The deformations involved in the model are determined by the choice of the space of vector fields V . This space is usually a Reproducing Kernel Hilbert Space (RKHS) and the representations of the data rely thus on the choice of the kernel.

The choice to study animal horns intended to avoid the ill-posed distinction between diffeomorphic and intrinsic changes, i.e. the unknown balance between the emergence of new coordinates and the stretching of a shape. In order to exclude this issue, the most natural choice would be to use affine transformations. These transformations respect the idea that once a portion of the horn has appeared, it behaves as a solid. It is not deformed, only displaced. However, this setting requires an additional adaptation of the cost function that would complicate the presentation of the global model introduced in this paper. Morever, the case study of horns with rigid deformations is a first step to initiate the study of growth models. For clarity and in order to pave the way for more general applications, we will thus present experiments using a RKHS with a scalar Gaussian kernel: [START_REF] Trélat | Contrôle optimal : théorie et applications[END_REF] k

V : R d × R d -→ R (x, y) -→ exp -|x-y| 2 2σ 2 .
To retrieve the rigidity of the rotations and translations, we used a large σ. can be extracted at any time t ∈ [0, 1]. To compare the solution with the input data, Figure 13 displays them both, with their respective foliation, at the specific time steps of the input and the shapes are, as intended, restricted to the active coordinates (see equation ( 12)). Although the target is known with its foliation, this foliation is not seen by the optimization problem.

The aim of the model is to retrieve this hidden variable. We recall that this variable is described by the color gradient in all the figures of this paper and it intends to specify the correspondences between shapes. One can observe that the pattern of correspondences illustrated in Figure 13 could not be retrieved by a standard registration method. When new coordinates are activated, the norm of the optimal vector field tends to increase as illustrated in Figure 14. Note yet that, if necessary, the new coordinates also allows this norm to decrease: one can observe in Figure 15 that the momenta of the fifth leaf have opposite directions with the momenta of the younger leaves. This flexibility is new as compared to the standard LDDMM method and it illustrates that the deformation can be disturbed from its natural geodesic evolution at each creation of a new leaf. To display the norm of the optimal vector field, we compute the associated trajectory one more time at the end of the algorithm with a finer time discretization to show that between the appearance of two successive leaves, the dynamic is classic (no creation occurs during these intervals) and the norm is thus constant. Note at last that if we The arrangement of the momentum (concentrated around the tip) is strongly linked to the scale of the kernel. Here, we have : σ V = 0.9 constrain the initial momentum p 0 to the vertical axis, the foliation is perfectly retrieved but the result presented in Figure 13 is more interesting to understand the behavior of the model with more degrees of freedom.

6.3. Example 2 -Horn. Here, we apply the model to estimate the development of a horn. We study on this second example the impact of intermediate times. Figure 16, 17 and 18 display the results of the experiment, with and without intermediate times. As for the previous example, the model is able to globally reproduce the growth process of the shape with smooth large deformations. The orientation of the leaves is better retrieved with intermediate times. Note that this toy example only aims to illustrate the bases of the model. In application to real data, the user should, when possible, optimize the model according to his data, e.g. with additional landmarks. Another option is to find the most appropriate kernel to generate the deformation. We can indeed observe that the Gaussian kernel favors local translations over rotations which prevents to better match the orientation of the horn's leaves. As we will discuss in the next section, the group of rotations and translations allows to perform much more complex matching of horn's development.

Figure 19 displays the initial momenta p 0 on the initial position q 0 . In both cases, the momenta have been smoothed during their optimization. It led to a significant improvement especially with intermediate input data. Let us report at last that the norm of the optimal vector field increases in both experiments. 7. Extensions 7.1. Adapted deformation model. The control we have on the deformation model built on a RKHS with a Gaussian kernel is limited. We saw in our experiments that it is not completely adapted to reproduce the deformations involved in the growth of a horn. The most natural way would be to consider the group of rotations and translations. However, its successful integration to the variational problem of matching with the growth dynamic requires subtle modulations that go beyond the scope of this paper. Figure 20 illustrates yet the large improvement resulting from these changes. This type of long and highly curved horns could hardly be produced with a Gaussian kernel but more importantly, the foliation of the horn is significantly better retrieved. Another promising approach should be to use the modular framework for diffeomorphisms recently introduced in [START_REF] Gris | Modular approach on shape spaces, sub-Riemannian geometry and computational anatomy[END_REF] that builds combinations of local deformation modules with few degrees of freedom allowing a strong control and restriction of the optimal deformations. 7.2. Plant growth. We propose the heuristics for a straightforward application of our model to plant growth. In some of the few examples we observed, the bottom of the leaves' boundary (petiole excluded) has a shape either flat or similar to the letter v. In this case, this area can be consider as the birth place for the new coordinates. The growth is then modeled exactly as for horns and illustrated in Figure 21. However, in some cases, this base of the leaf can be too small meaning that the deformation required to model the evolution of the new coordinates is too important. We then suggest to systematically cut the bottom of the leaf and consider this cut as the new place of creation. 7.3. Other evolution types. We presented a theoretical model more general than the growth dynamic. An important class to explore in a future work consists in degenerative processes. However, other situations can be direct applications of the model with the growth dynamic even if they do not involve creation of material. The next example aims to sketch the opening of a pin cone. The discontinuity of the growth dynamic allows to separate the scales of the cone.

With intermediate times

Without intermediate times

Figure 22 presents the matching between two sets of curves. The source, in blue, represents 5 scales of a close pin cone and consists in 5 vertical segments aligned on a single axis but slightly shifted upwards. The target, in red, represents the open pin cone. The biological coordinate system is a union of 5 segments with for each of them a specific time of activation given by a function τ : X → {t 1 , . . . , t 5 }. Although all the scales actually exist at initial time t = 0, they are progressively activated (t 1 < . . . < t 5 ).

Conclusion

The general notion of biological coordinate systems that respect the biological homology in a population of related shapes allowed us to present in this paper a new concept to generate curves of shape that evolves under a growth process by foliation. These curves are produced with time-varying vector fields as usually achieved with diffeomorphic registration methods. The model required a specific theoretical analysis since unlike the standard approach, the infinitesimal action of the vector field depends on time. The main issue raised by this novelty lies in the spatial regularity of the generated mappings that represent a new scenario.

We successfully integrated this generative model in the problem to estimate the continuous evolution of a time-varying shape from its final state. It led to a new optimal control problem where the time dependency played again an important role. We introduced a general framework extending the LDDMM methods where the Hamiltonian approach is maintained despite the lost of the conservation property of an optimal vector field. This approach parametrizes the shape evolutions with a forecast initial position and initial momentum (q 0 , p 0 ), providing the support to a statistical analysis. We implemented this new registration method and illustrated how it provides an original and biologically coherent modeling of growth phenomena.

At last, note that the general construction of this work is paving the way for new models in computational anatomy able to integrate priors on the self-organization of living organisms.

Figure 1 .

 1 Figure 1. (a) Top row: growth of an oak's leaf. Bottom row: growth of a basil's leaf. Secondary veins are enumerated to highlight their increasing number in the second case. (b)The red doted areas are called epiphyseal plates. Although the ossification is a complex diffused process, the main part of the elongation comes from these areas where the long bones grow longer.

Figure 2 .

 2 Figure 2. The creation of new material during a growth process is linked to the appearance of new coordinates.

Figure 3 .

 3 Figure 3. Sketch of the modeling of a longitudinal shape data set where individuals share a common growth pattern.

  0

Figure 4 .

 4 Figure 4. The nested shapes (X t ) t are individually displayed on a time-line to highlight the growth process induced by the biological coordinate system.

Figure 9 .

 9 Figure 9. "Horns" generated with different initial positions and vector fields. The deformations are restrained to vertical upward translations whose amplitudes are displayed in the first row. The first column shows the initial positions q 0 (X). We display in the center of the table the final cones q 1 (X) resulting from this nine configurations. One can compare, regarding to the growth process, the variations of the solutions with respect to the cone on the top left corner.

6. 2 .Figure 13 .

 213 Figure13. Development of a cone. The solution is displayed against the target scenario (faded shapes). At each step of the gradient descent, the data attachment term measures the discrepancy of the six pairs of homologous horns.

Figure 15 .

 15 Figure 15. Initial momenta p 0 on the initial position q 0 . The arrangement of the momentum (concentrated around the tip) is strongly linked to the scale of the kernel. Here, we have : σ V = 0.9

Figure 16 .

 16 Figure 16. Development of the horn with and without intermediate times. The faded evolution is the target. On the first row, the matching only included the final state of the target. On the second row, all the intermediate states displayed have been used in the gradient descent.

Figure 17 .

 17 Figure 17. Final state of the horn without intermediate times. Left: result at time 1. Right: target at time 1. Middle: overlay of both.

Figure 18 .Figure 19 .

 1819 Figure 18. Final state of the horn with intermediate times. Left: result at time 1. Right: target at time 1. Middle: overlay of both.

Figure 20 .

 20 Figure 20. We display the development of the solution overlayed with the target's development. In addition on the top row, the final state from two different points of view.

Figure 21 .

 21 Figure 21. The base of the leaf, pointed by blue arrows, defines the birth place of the new coordinates. Doted lines represent the successive sets of new coordinates.

Figure 22 .

 22 Figure 22. Pin cone opening.

If these shapes are modeled by currents or varifolds in a RKHS W , respectively denoted µ v i ∈ W for the solution generated by v and µ tar i ∈ W for the target, the attachment term reads

Once again, we can add the L 2 distance on each pair of tips of the horns. With the previous notation, the final attachment term is given by

.

Numerical experiments

This section applies the framework presented in this paper to few examples of synthetic horns. The input of these experiments consist in a final state of a horn with eventually few prior states and the common biological coordinate system modeling the growth process of the population. We used the system introduced in Example 2.1. In each case, we performed the gradient descent algorithm described in Section 3.3.3 to build the complete growth scenario fitting the input data. The solution is initialized with the flat disc modeling the initial position as detailed in the next section.

The meshes all have the same structure : 3×l ×L points to highlight the leaves of the foliations. Shapes are thus discrete foliations of L leaves where each leaf contains l vertices. These vertices are merged for the first leaf to model the tip of the horn. The leaves always appear at regular time intervals. The time discretization of the interval [0, 1] used for the integration steps (forward and backward) follows the creation times of the leaves, which means that this time step equals 1/(L -1). Refining this discretization did not seem to improve significantly the results in the following examples. 6.1. General settings.

6.1.1. Initial position. Once the numbers of leaves and points of the meshes are fixed, the first step to create a scenario is to generate its initial position. Let us recall that the initial position can be seen as the pull backward through the flow of each leaf X {t} from its final position to its position at creation at time t. In order to generate a synthetic horn target, we define the form of the last leaf that will be the base of the target at its final state and compute a linear reduction of this leaf toward the center that will be the tip of the horn at its birth. To initialize the algorithm, the initial position q 0 is defined likewise via the base of the target. We compute a linear reduction of the base (closed polyline) to generate the L initial leaves.