Arguill Sylvain

Emmanuel Tr Ère
email: emmanuel.trelat@upmc.fr

Alain Élat

Laurent Trouv É
email: laurent.younes@jhu.edu

Younes

E Tr

A Trouv
email: trouve@cmla.ens-cachan.fr

Ecole Normale

MULTIPLE SHAPE REGISTRATION USING CONSTRAINED OPTIMAL CONTROL

Keywords: Subject Classification. 58D05 49N90 49Q10 68E10 Shape analysis, optimal control, deformations, groups of diffeomorphisms

Lagrangian particle formulations of the large deformation diffeomorphic metric mapping algorithm (LDDMM) can be numerically challenging when the number of particles is large. In this paper, we introduce and discuss numerical schemes useful for surface and image matching, which are based on representing the Eulerian velocity over a finite-dimensional basis that deforms in time. The method is described within the optimal control formalism, and optimality conditions are derived, together with the equations that are needed to implement gradient-based methods. Experimental results are shown, both with surfaces and images.

2000 Mathematics

INTRODUCTION

The large deformation diffeomorphic metric mapping (LDDMM) approach to shape matching is a powerful topology-preserving registration method with an increasing record of successful applications in medical imaging. It was first described in [START_REF] Sarang | Landmark matching via large deformation diffeomorphisms[END_REF] for point sets and in [START_REF] Dupuis | Variation Problems on Flows of Diffeomorphisms for Image Matching[END_REF][START_REF] Trouvé | Diffeomorphisms groups and pattern matching in image analysis[END_REF][START_REF] Michael I 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] for images and has become widely used in the medical imaging literature and other applications. While deeper understanding and extensions of the underlying theoretical framework was pursued [START_REF] Michael I Miller | The metric spaces, euler equations, and normal geodesic image motions of computational anatomy[END_REF][START_REF] Glaunès | Landmark Matching via Large Deformation Diffeomorphisms on the Sphere[END_REF][START_REF] Cao | Large deformation diffeomorphic metric mapping of vector fields[END_REF][START_REF] Michael I Miller | Geodesic shooting for computational anatomy[END_REF][START_REF] Younes | Jacobi fields in groups of diffeomorphisms and applications[END_REF][START_REF] Glaunès | Large deformation diffeomorphic metric curve mapping[END_REF][START_REF] Younes | Shapes and diffeomorphisms[END_REF][START_REF] Bruveris | The momentum map representation of images[END_REF][START_REF] Sylvain Arguillère | Shape deformation analysis from the optimal control viewpoint[END_REF] and alternative numerical methods were designed [START_REF] Brian B Avants | Lagrangian frame diffeomorphic image registration: Morphometric comparison of human and chimpanzee cortex[END_REF][START_REF] Cotter | Singular solutions, momentum maps and computational anatomy[END_REF][START_REF] Trouvé | Sparse approximation of currents for statistics on curves and surfaces[END_REF][START_REF] Niethammer | Geodesic regression for image time-series[END_REF][START_REF] Günther | Direct lddmm of discrete currents with adaptive finite elements[END_REF][START_REF] Franc ¸ois-Xavier | Diffeomorphic 3D Image Registration via Geodesic Shooting Using an Efficient Adjoint Calculation[END_REF][START_REF] Ashburner | Diffeomorphic registration using geodesic shooting and Gauss-Newton optimisation[END_REF][START_REF] Günther | Flexible shape matching with finite element based lddmm[END_REF][START_REF] Stanley Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF], LDDMM has been applied to medical imaging data including brain [START_REF] Michael I Miller | Increasing the power of functional maps of the medial temporal lobe by using large deformation diffeomorphic metric mapping[END_REF][START_REF] Zhang | Characterization of mouse brain and its development using diffusion tensor imaging and computational techniques[END_REF][START_REF] Qiu | Cortical hemisphere registration via large deformation diffeomorphic metric curve mapping[END_REF]11,[START_REF] Risser | Simultaneous multi-scale registration using large deformation diffeomorphic metric mapping[END_REF], heart [START_REF] Ardekani | Cardiac motion analysis in ischemic and non-ischemic cardiomyopathy using parallel transport[END_REF][START_REF] Azencott | Diffeomorphic matching and dynamic deformable surfaces in 3d medical imaging[END_REF] and lung [START_REF] Vidal | Template registration with missing parts: Application to the segmentation of m. tuberculosis infected lungs[END_REF] images. This algorithm provides a non-rigid registration method between various types of objects (point sets, curves surfaces, functions, vector fields. . . ) within a unified framework driven by Grenander's concept of deformable templates [START_REF] Grenander | General pattern theory: A mathematical study of regular structures[END_REF]. It optimizes a flow of diffeomorphisms that transform an initial object (shape) into a target one.

The practical importance of shape registration is underlined by the increasing amount of work that has flourished in the literature over the past few years. LDDMM is one among many methods that have been proposed to perform this task. Several such methods are based on elastic matching energies [START_REF] Bajcsy | A computerized system for the elastic matching of deformed radiographic images to idealized atlas images[END_REF][START_REF] Droske | A variational approach to nonrigid morphological image registration[END_REF], and other, like LDDMM, inspired by viscous fluid dynamics [12,[START_REF] Thirion | Image matching as a diffusion process: an analogy with maxwell's demons[END_REF][START_REF] Vercauteren | Non-parametric diffeomorphic image registration with the demons algorithm[END_REF][START_REF] Ashburner | A fast diffeomorphic image registration algorithm[END_REF][START_REF] Vercauteren | Diffeomorphic demons: Efficient non-parametric image registration[END_REF]. For surfaces, which will be our main focus, several authors have developed approaches to find approximate conformal parametrizations with respect to the unit disk or sphere [START_REF] Hurdal | Quasi-conformally flat mapping the human cerebellum[END_REF][START_REF] Haker | Conformal surface parameterization for texture mapping[END_REF][START_REF] Hurdal | Coordinate systems for conformal cerebellar flat maps[END_REF][START_REF] Gu | Genus zero surface conformal mapping and its application to brain surface mapping[END_REF][START_REF] Miao | Optimal global conformal surface parameterization[END_REF][START_REF] Wang | Intrinsic brain surface conformal mapping using a variational method[END_REF][START_REF] Monica | Cortical cartography using the discrete conformal approach of circle packings[END_REF][START_REF] Gu | Genus zero surface conformal mapping and its application to brain surface mapping[END_REF][START_REF] Monica | Discrete conformal methods for cortical brain flattening[END_REF]. More recently, quasi-conformal parametrizations based on the minimization of the Beltrami coefficient have been designed [START_REF] Zeng | Registration for 3d surfaces with large deformations using quasi-conformal curvature flow[END_REF][START_REF] Lok | Shape-based diffeomorphic registration on hippocampal surfaces using Beltrami holomorphic flow[END_REF][START_REF] Wang | Teichmüller shape space theory and its application to brain morphometry[END_REF]. Another class of non-rigid registration methods include those based on optimal mass transportation [START_REF] Haber | 3d nonrigid registration via optimal mass transport on the gpu[END_REF][START_REF] Haker | Optimal mass transport for registration and warping[END_REF][START_REF] Lipman | Surface comparison with mass transportation[END_REF][START_REF] Lipman | Conformal wasserstein distances: Comparing surfaces in polynomial time[END_REF], while [START_REF] Mémoli | A theoretical and computational framework for isometry invariant recognition of point cloud data[END_REF][START_REF] Mémoli | On the use of gromov-hausdorff distances for shape comparison[END_REF][START_REF] Mémoli | Gromov-wasserstein distances and the metric approach to object matching[END_REF] introduce comparison methods based on Gromov-Hausdorff or Gromov-Wasserstein distances.

• M = {1, . . . , m} is finite, and q can be identified with a collection q 1 , . . . , q m of distinct points in R d . • M = [0, 1] and q is a curve in R d .

• M = S d-1 (the unit sphere in R d ) and q is a hypersurface.

Our goal is to discuss models in which several shapes can deform, while being subject to contact constraints. The deformation process will be similar to the one designed for large deformation diffeomorphic metric mapping (LDDMM), which can be formulated as an optimal control problem. Before introducing our general framework, it will be easier to start with a description of the now well explored single-shape problem upon which we will build. For this, we let (V, • V ) be a Hilbert space of vector fields on R d , assumed to be continuously imbedded in the space B := C p 0 (R d , R d ), the completion of the space of smooth compactly supported vector fields for the norm • p,∞ , which denotes the sum of supremum norms of derivatives of order p or less, with p ≥ 1. Then V possesses a reproducing kernel, that is a mapping K : (x, y) → K(x, y), defined over R d × R d , with values in the space of d × d matrices, such that all partial derivatives with order less than p with respect to each variable exist and K(•, y)a ∈ V with K(•, y)a , w V = a T w(y), for all (a, y) ∈ (R d ) 2 . The LDDMM algorithm uses flows of time-dependent vector fields v(•) ∈ L 2 (0, 1; V ).

2.2.

Registering Two Shapes Using LDDMM.

General Problem.

The general LDDMM problem is formulated as the infinite-dimensional optimal control problem consisting of minimizing the cost functional (1)

F (v) = 1 2 1 0 u(t) 2 V dt + U (q(1)),
subject to the constraint ∂ t q(t) = u(t) • q(t) for a.e. t ∈ [0, 1],

q(0) = q init . (2)
This differential constraint is a control system, where the control is the time-dependent vector field u(•) ∈ L 2 (0, 1; V ), the solution of which is q(t, •) = ϕ(t, q init (•)) where ϕ is the flow of diffeomorphisms generated by u(•), defined as the unique solution of the Cauchy problem ∂ t ϕ(t) = u(t) • ϕ(t), ϕ(0) = id R d . For every time t, we have ϕ(t, •) ∈ Diff p , the set of p-times differentiable diffeomorphisms in R d .

The function U is a matching cost function, that is, a penalization that pushes the solution of (1)-(2) towards a target. It will be assumed to be Fréchet differentiable from Q to R. To simplify the discussion, and because this covers most of the interesting cases in practice, we will assume that there exists some fixed measure ν M on M such that its derivative, denoted dU (q) or dU q when evaluated at q ∈ Q, can be expressed in the form dU q = z q ν M for some (ν M -measurable)

z q : M → R d , meaning that ∀h ∈ C p (M, R d ), (dU q | h) = M h(m) • z q (x) dν M (x).
Throughout the paper, for any Banach space X, the notation (µ | v) will be used to designate the application µ(v) of a linear form µ ∈ X * to a vector v ∈ X.

Under these assumptions, one can prove that the gradient of the objective function F defined by (1) (which is a mapping on L 2 ([0, 1], V )) is given by

∇ V F (u)(t, ξ) = u(t, ξ) - M K(ξ, q(t, x))α(t, x) dν M (x),
where K is the reproducing kernel of V , and α :

[0, 1] × M → R d is a time-dependent function defined by α(1, •) = -z q(1) and (3) ∂ t α = -(du • q) T α.
This result implies, in particular, that the solutions of ( 1)-( 2) must satisfy the Pontryagin maximum principle (see [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF][START_REF] Trélat | Contrôle optimal[END_REF]), which is the following first-order necessary condition for optimality. Let H u be the Hamiltonian defined by

H u (ρ, q) = (ρ | u • q) - 1 2 u 2 V ,
for every u ∈ V , every q ∈ Q and every ρ ∈ Q * . If u(•) is an optimal control, solution of the optimal control problem (1)-( 2), then it must be such that (4) u(t) = argmax w H w (ρ(t), q(t)),

where (ρ, q) are solutions of

∂ t q(t) = ∂ ρ H u (ρ(t), q(t)), ∂ t ρ(t) = -∂ q H u (ρ(t), q(t)),
(ρ is the so-called co-state, or adjoint state) and ρ(1) = -dU q(1) . Indeed, it suffices to take ρ(•) = α(•)ν M and to use properties of the reproducing kernel to check that all the conditions are satisfied. Moreover, (4) then implies that u = M K(•, q(x))α(x) dν M at every time.

Examples of shapes and matching cost functions.

Example 1. To start with a simple example, let M = {1, . . . , m} so that a shape q = (q(1), . . . , q(m)) is a collection of landmarks, and consider the landmark matching cost function defined by

U (q) = m k=1 |q(k) -y k | 2 , for fixed y = (y 1 , . . . , y m ) ∈ R N , in which | • | is the Euclidean norm on R d . We then have (dU q | h) = 2 m k=1 (q(k) -y k ) T h(k) ,
or, to interpret this result in the general form provided above, dU q = zν M with z(k) = 2(q(k)-y k ) and ν M the counting measure on {1, . . . , m}.

Example 2. If µ is a scalar measure on R d and z a µ-integrable R d -valued function defined on support(µ), we will denote by zµ the vector measure such that

(zµ | w) = support(µ) z • w dµ,
where z • w denotes the standard euclidean dot product.

Vector measures of the form zµ are continuous linear forms over any space that is continuously imbedded in C 0 0 (R d , R d ), and in particular over any reproducing kernel Hilbert space W . For such a space, equipped with a reproducing kernel χ, the operator norm of zµ is given by

zµ 2 χ = z(x) • (χ(x, y)z(y)) dµ(x) dµ(y),
and more generally, the norm of the difference between two such measures is

zµ -z μ 2 χ = z(x) • (χ(x, y)z(y)) dµ(x) dµ(y) -2 z(x) • (χ(x, y)z(y)) dµ(x) dμ(y) + z(x) • (χ(x, y)z(y)) dμ(x) dμ(y).
Note that W and V have no relationship one to each other, except that both have a continuous inclusion in C 0 0 (R d , R d ), so that χ is different from K. One can deduce from this the surface-matching cost function introduced in [START_REF] Vaillant | Surface Matching via Currents[END_REF], in which an oriented surface S is represented as a geometric current and a dual-RKHS norm between currents is used. Identifying surface currents with vector measures, this leads to the representation of S given by µ S = N S σ S , where N S is the unit normal to S and σ S its volume form. Assume that M (the parameter space) is an oriented 2D manifold so that S = q(M ) is a surface, and let η be a positively oriented volume form on M . For m ∈ M let N q (x) ∈ R 3 denote the "area-weighted normal" to S = q(M ) at q(x), defined by N q (x) = dq x e 1 × dq x e 2 where (e 1 , e 2 ) is an arbitrary basis of T x M such that η x (e 1 , e 2 ) = 1. Then

(µ S | w) = M (w • q • N q ) dη, for every w ∈ C 0 0 (R 3 , R 3
). Now, given a reproducing kernel χ and a target surface S = q(M ), we define the surfacematching cost by

U (q) = µ q(M ) -µ q(M ) 2 χ .
Example 3. This cost function is actually a special case of the most general framework in which one compares compact k-dimensional oriented submanifolds of R d , which we briefly discuss hereafter. Given such a manifold, S, with a global parametrisation q : M → S, one can associate to any ω ∈ C p 0 (R d , (Λ k R d ) * ) (the set of C p differential k-forms on R d that vanish at infinity), its integral

(C S | ω) = S ω = M q * ω,
where q * ω denotes the pull-back of ω on M . An RKHS, W , of such forms, is a Hilbert space continuously embedded in C p 0 (R d , (Λ k R d ) * ), with kernel χ(x, y) taking values in the space of bilinear functions on

Λ k R d × Λ k R d .
The linear form C S then belongs in W * , and is a special form of a geometric current, as defined in [START_REF] Federer | Geometric measure theory[END_REF]. If S = q(M ) and S is a target manifold, they can be compared using the operator norm ( 5)

U (q) = C q(M ) -C S 2 W * .
Now, if we consider q : M . = (-1, 1) × M → R d a smooth perturbation of q such that q ε = q + εδq + o(ε) where q ε (x) = q(ε, x) for ε ∈ (-1, 1) and x ∈ M , we have, for

M ε . = {ε} × M , (dU q | δq) = d d M q * ω | =0 = d d M q * ω | =0 = M 0 L ∂/∂ε (q * ω),
where L ∂/∂ε is the Lie derivative along the vector field ∂ ∂ε on M (which is equal to (1, 0) ∈ R × T m M at any location (ε, x) ∈ M ) and ( 6)

ω = 2K W (C q(M ) -C S ),
with K W the isometry from W * to W . We next show that

(dU q | δq) = M α q • δq vol M + ∂M β q • δq vol ∂M ,
where vol M and vol ∂M are the positive Riemannian volume forms on M and ∂M , and α q : M → R d (resp. β q : ∂M → R d ) is such that α q (x) is normal to S = q(M ) (resp. β q (x) is normal to ∂S = q(∂M )) at q(x). Using the Cartan magic formula we get

L ∂/∂ε (q * ω) = i ∂/∂ε d(q * ω) + d(i ∂/∂ε (q * ω)) ,
so that, applying the Stokes theorem,

(dU q | δq) = M 0 i ∂/∂ε (q * dω) + ∂M 0 i ∂/∂ε (q * ω),
where i ∂/∂ε denotes the contraction operator. Note that, for ξ 1 , . . . ,

ξ k ∈ T x M , i ∂/∂ε (q * dω) (0,x) (ξ 1 , • • • , ξ k ) = dω q(x) (δq(x), dq x ξ 1 , . . . , dq x ξ k ) = dω q(x) (δq ⊥ (x), dq x ξ 1 , . . . , dq x ξ k ),
where δq ⊥ (x) denotes the projection of δq(x) on T q(x) S ⊥ (since the form vanishes if δq(x) ∈ T q(x) S = T q(x) q(M )). Since the set of k-forms is a one-dimensional space on M , this means that we can write, for some function α q such that α q (x) ∈ (T q(x) S) ⊥ , with x ∈ M , i ∂/∂ε (q * dω) (0,x) = ((δq

• α q ) vol M ) x .
Similarly, i ∂/∂ε (q * ω q ) (0,x) = ((δq

• β q )vol ∂M ) x , for β q (x) ∈ (T q(x) ∂S) ⊥ .
The data term defined in ( 5) is derived from this general construction, using the fact that twoforms in R 3 (or (d -1)-forms in R d ) can be identified with vector fields via ω x (e 1 , e 2 ) = w(x) • (e 1 × e 2 ). The current C S is then identified with the vector measure µ S . The form ω introduced in (6) becomes, introducing a global parametrizaion q : M → S of the target S, the vector field

w(•) = 2 R d χ(•, y) d(µ q(M ) -µ q(M ) )(y).
With this identification, we have α q = div(w) • q N q where N q is the oriented "area-weighted" normal to S at q(x) defined previously, and β q = τ q × w • q, where τ q is the oriented "lengthweighted" tangent to q(∂M ) given as τ q (x) = dq x e 1 , where e 1 is the unit positively oriented tangent vector at x along ∂M .

Example 4. Returning to surfaces, the discrete case, in which triangulated surfaces are compared, is, for practical purposes, even more important. We here also consider the case M = {1, . . . , m}, with an additional family F of facets, which are ordered triples (i, j, k) with i, j, k ∈ M . (We assume that F is a consistent with a manifold structure: The set V i of indices that share a facet with i must form a chain and no pair of indices can be included in more than two facets.)

Given a one-to-one mapping q : M → R 3 , define S q as the collection of triangles S q = {(q(i), q(j), q(k)), (i, j, k) ∈ F }. If f = (i, j, k), let q(f ) = (q(i), q(j), q(k)), N q (f ) = (q(j)q(i)) × (q(k) -q(i)) and c q (f ) = (q(i) + q(j) + q(k))/3 respectively denote the triangle, areaweighted normal and center associated to the facet f . Following [START_REF] Vaillant | Surface Matching via Currents[END_REF], we define the vector measure associated to q by

µ q = f ∈F N q (f )δ cq(f ) .
Here, δ x (with x ∈ R 3 ) denotes the atomic measure of mass 1 with support {x}. The (discrete) surface matching cost associated to a target q is then defined by U (q) = µ q -µ q 2 χ . Note that q does not need to be consistent with q, and can be defined on a different set of indices, M = {1, . . . , m} and triangle structure F . One then has dU q = α q ν M , where, as above, ν M is the counting measure on M , and

α q (i) = f ∈F,i∈f (dZ T cq(f ) N q (f )/3 + e q (f, i) × Z(c q (f ))),
with e q (f, i) = q(k) -q(j) the oriented edge opposed to q(i) in q(f ), and

Z(•) = f ∈F χ(•, c q (f ))N q (f ) - f ∈ F χ(•, c q( f )N q( f ).
2.2.3. Reduced Problem. Since the optimal control must satisfy

(7) v = M K(•, q(t, x))α(t, x) dν M
for some function α defined on M , it is natural to parametrize v by α and use this function as a new control. We define the inner product

α , β q = M ×M α(x) T K(q(x), q(x))β(x) dν M (x) dν M (x)
between two measurable functions α and β defined on M . If v is given by ( 7), the reproducing property of the kernel implies that v 2 V = α 2 q . The optimal control problem ( 1)-( 2) is then equivalent to the reduced problem consisting of minimizing the cost functional ( 8)

F (α) = 1 2 1 0 α(t) 2 q(t) dt + U (q(1)), subject to the constraint (control system) (9) ∂ t q(t, x) = M K(q(t, x), q(t, x))α(t, x)dν M (x),
almost everywhere over the time interval [0, 1].

According to [START_REF] Azencott | Diffeomorphic matching and dynamic deformable surfaces in 3d medical imaging[END_REF]13,[START_REF] Vaillant | Surface Matching via Currents[END_REF][START_REF] Younes | Evolutions equations in computational anatomy[END_REF], we have ∇F (α) = α-p, where p (the co-state) is a time-dependent vector-valued measurable function on M such that p(1)ν M = -dU q(1) and

∂ t p(t) = -∂ q p(t) , α(t) q -α(t) 2 q(t) /2 ,
where q is defined by [START_REF] Bruveris | The momentum map representation of images[END_REF]. Here, the gradient is computed with respect to the inner product • , • q .

3. MULTIPLE SHAPE PROBLEMS 3.1. Motivating Examples. In the previous formulation, the shape evolution was controlled by a single, smooth vector field v, inducing a single diffeomorphism of R d restricted to the considered shape. This approach has been successfully used to model variations of single, homogeneous shapes, and led to important applications in computational anatomy, including, among many other examples, the impact of pathologies like Huntington disease [START_REF] Younes | Regionally selective atrophy of subcortical structures in prodromal hd as revealed by statistical shape analysis[END_REF], schizophrenia [START_REF] Qiu | Neuroanatomical asymmetry patterns in individuals with schizophrenia and their non-psychotic siblings[END_REF], and Alzheimer's disease [START_REF] Wang | Large deformation diffeomorphism and momentum based hippocampal shape discrimination in dementia of the alzheimer type[END_REF][START_REF] Durrleman | Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets[END_REF][START_REF] Risser | Simultaneous multi-scale registration using large deformation diffeomorphic metric mapping[END_REF] on brain structures. This deformation model, however, is not well adapted in situations in which several shapes interact, or situations in which shapes have heterogeneous parts. Let us review some motivating examples.

(1) Consider a schematic representation of a kite, or a manta ray, composed with a twodimensional surface, representing the body, and an open curve attached to it representing the tail. When comparing two such objects, the body is assumed to only show small differences in shape, while the tail can vary widely. (2) Consider a two-dimensional representation of a mouth, with two curves representing the upper and lower lip. Because the mouth can be wide open or closed, it is not possible to consider its deformations as resulting from the restriction of a smooth diffeomorphism of R 2 . (3) Finally, it is natural, when analyzing multiple organs in the human body, to consider multiple shapes, each of them being relatively stable (only subject to small deformations) while their position with respect to each other is subject to larger variations, so that the background (the intersection of their complements) is subject to very large deformations. Here again, modeling the whole process with a single diffeomorphism is not adequate.

These examples suggest using multiple deformations applied to each component of the considered model. Generalizing ( 1)-( 2), consider parameter spaces M 1 , . . . , M n for an n-component model. Each shape, or component, is a mapping

q (k) ∈ Q k : M k → R d . The shape space will then be Q = Q 1 × • • • × Q n . To each shape, associate a control u k ∈ V k , where V k is an RKHS embedded in C p 0 (R d , R d ) with the state evolution equation ∂ t q (k) = u k • q (k)
. We can then choose each V k according to how "wildly" we want to allow the k-th shape to deform. These evolutions, however, must be consistent with each other, implying contact constraints that we will consider in two forms:

• Identity constraints: These are constraints that make a subset of the k-th shape stay stitched to a subset of the l-th shape, so that these subsets coincide in R d and move identically along the deformation. Given some pair (k, l) ∈ {1, . . . , n} 2 , and given a one-to-one mapping g kl :

A kl ⊂ M k → A lk = g kl (A kl ) ⊂ M l , one has q (k) (x) = q (l) (g kl (x)) for every x ∈ A kl .
• Sliding constraints: These are constraints that force a closed submanifold of the k-th shape to slide on a corresponding submanifold of the l-th shape along the deformation, for all (k, l). Here, we assume that all parameter spaces are orientable differential manifolds, and that all q (k) 's are immersions. Given some pair (k, l) ∈ {1, . . . , n}, and a closed submanifold without boundary A kl ⊂ M k , there exists a diffeomorphism g kl :

A kl → A lk ⊂ M l onto a fixed closed submanifold A lk of M l such that q (k) (x) = q (l) (g kl (x)) for every x ∈ A kl .
Let us turn back to the examples mentioned at the beginning of the section. For Example (1), we can take M 1 = S 2 and M 2 = [0, 1], and, letting x 0 represent the north pole in S 2 , impose q (1) (x 0 ) = q (2) (0). We can then assign different deformation models to q (1) and q (2) via the metrics on V 1 and V 2 .

Example (2) requires a slightly more complex construction, that only imperfectly addresses the issue. Let (1) represent the upper lip, q (2) the lower one and q (3) their union. We use the identity constraints q (1) (1) = q (2) (0) and q (2) (1) = q (1) (0) for the extremities of each lip, and q

M 1 = M 2 = [0, 1] and M 3 = {1, 2} × [0, 1]. Let q
(3) (1, •) = q (1) (•), q (3) (2, •) = q (2) (•). We take V 1 = V 2
and choose V 3 such that the latter allows for large deformations at a small cost. With this model, it becomes easier to "almost" close the mouth, although the deformation inside the mouth must remain diffeomorphic, so that the closing cannot go all the way.

For Example (3), there are n shapes, n -1 of which are associated with the organs, and the last of which represents the background. For example, we can take M k = S 2 for k = 1, . . . , n -1, and M n = {1, . . . , n -1} × S 2 . Assuming that the shapes do not intersect, we can define identity or sliding constraints for the background, enforcing

q (n) (k, •) = q (k) (•) or q (n) ({k}×S 2 ) = q (k) (S 2 ) = S k for k ∈ {1, . . . , n -1} during the deformation.
3.2. Induced Constraints. The previous constraints can be reformulated as equality constraints involving the state and control. Identity constraints q (k) (x) = q (l) (g(x)) are equivalent (taking time derivatives) to u k (t, q (k) (t, x)) = u l (t, q (l) (t, g(x))) as soon as the constraints are satisfied at time t = 0, which we obviously assume.

Making the same assumption, sliding constraints can be expressed as ( 10)

N (k) (t, q (k) (t, x)) T (u k (t, q (k) (t, x)) -u l (t, q (k) (t, x))) = 0,
where

N (k) (t, q (k) ) is a d × (d -dim(A kl )
) matrix consisting of independent vectors perpendicular to T q (k) B kl (t) (e.g, a normal frame to B kl ), with B kl = q (k) (A kl ) for every (k, l). Let us briefly justify this statement.

We express the sliding constraint as q (k) (t, x) = q (l) (t, g(t, x)) for some diffeomorphism g(t, •) : A kl → A lk , assuming a differentiable dependency on time. Taking time derivatives, we get

u k (t, q (k) (t, x)) = u l (t, q (l) (t, g(t, x))) + dq (l) (t, g(t, x))∂ t g(t, x), x ∈ M k .
Since q (l) (t, g(t, x)) = q (k) (t, x), we obtain

u k (t, q (k) (t, x)) -u l (t, q (k) (t, x)) = dq (l) (t, g(t, x))∂ t g(t, x) = dq (k) (t, x)dg(t, x) -1 ∂ t g(t, x),
which is tangent to B kl at x. Note that, since the image of g(t, •) is A lk for every time t, we do have

∂ t g(t, x) ∈ T g(t,x) A lk = dg(t, x)(T x A kl ), so dg(t, x) -1 ∂ t g(t, x) is well-defined.
Conversely, assume that (10) holds for every x ∈ A kl , with q (k) (0, x) = q (l) (0, g 0 (x)) for some diffeomorphism g 0 : A kl → A lk ⊂ M l . Then for every time t, the mapping

w : x ∈ A kl → dq (k) (t, x) -1 (u k (t, q (k) (t, x)) -u l (t, q (k) (t, x)) ∈ T q (k) (x) B kl ) ∈ T x A kl
defines a time-dependent vector field on A kl . Since A kl is a closed manifold, this vector field is complete, and we denote its flow by h(t, •) :

A kl → A kl . Then ∂ t q (k) (t, h(t, x)) = u k (t, q (k) (t, h(t, x))) + dq (k) (t, h(t, x))∂ t h(t, x) = u l (q (k) (t, h(t, x))),
where the last identity holds for every x ∈ A kl , so that q (k) (t, h(t, x)) and q (l) (t, g 0 (x)) satisfy the same differential equation with the same initial condition and therefore coincide. Hence, letting g(t, x) = g 0 (h -1 (t, x)), we obtain q (l) (t, g(t, x)) = q (k) (t, x) for every x ∈ A kl .

It is possible to extend this construction to the case where A kl is a compact manifold with boundary. In this case, the matrix N (k) (t) must consist of a normal frame along ∂B kl (t) and possesses therefore an extra column, and we get an additional constraint along the boundary of A kl .

We will need the constraints to depend smoothly on q, and therefore we will need a smooth representation of the normal space to q (a smooth map q → N (q)) in order to be able to use [START_REF] Cao | Large deformation diffeomorphic metric mapping of vector fields[END_REF]. When this is not possible (or convenient), one can also use the alternative approach of introducing a new state, say N (k) (t, x), evolving according to

(11) ∂ t N (k) = -du k (q (k) ) T N (k) ,
which ensures that N (k) (t, x) remains perpendicular to T q (k) (t,x) S k as soon as this holds true at t = 0. The constraint

N (k) (t, x) T (u k (t, q (k) (t, x)) -u l (t, q (k) (t, x))
) is now a smooth function of the extended state.

The two problems that we consider are therefore special cases of the general problem considered in [START_REF] Sylvain Arguillère | Shape deformation analysis from the optimal control viewpoint[END_REF], which is the problem of minimizing the cost functional

(12) 1 2 n k=1 1 0 u k (t) 2 V k dt + n k=1 U k (q (k) (1)),
subject to the constraints (13) ∂ t q (k) (t) = u k (t, q (k) (t)), and C(q(t))u(t) = 0, almost everywhere over the time interval [0, 1], where C : M → L(V, Y) takes values in the space of bounded linear operators from V to a Banach space Y. Here, we have

V = V 1 × • • • × V n and q = (q (1) , . . . , q (n) ).
The study of this constrained optimal control problem, and in particular, the derivation of its first-order optimality conditions (of the type of Pontryagin maximum principle), is challenging in this infinite-dimensional setting. In [START_REF] Sylvain Arguillère | Shape deformation analysis from the optimal control viewpoint[END_REF], it is proved that, under some differentiability conditions, and under the important assumption that C(q) is surjective for every q ∈ M, optimal solutions must be such that there exist p = (p (1) , . . . , p

(n) ) ∈ H 1 ([0, 1], Q * ) and λ ∈ L 2 ([0, 1], Y * ) that satisfy (14)                  ∂ t q (k) = u k (q (k) ), ∂ t p (k) = -∂ q (k) p (k) | u k (q (k) ) -∂ q (k) (λ | C(q)u), u k , v V k = -p (k) | v • q (k) -(λ | C k (q)v), v ∈ V k , n k=1 C k (q)u k = 0,
where C k (q)u k = C(q)(0, . . . , 0, u k , 0, . . . , 0).

Unfortunately, the constraints C(q) that correspond to our identity or contact constraints are, in general, not surjective, and the results of [START_REF] Sylvain Arguillère | Shape deformation analysis from the optimal control viewpoint[END_REF] cannot be applied in a fully general infinitedimensional context. However, surjectivity becomes almost straightforward when these constraints are discretized to a finite number. They are true as soon as the points involved in the constraints are all distinct, which is a mild assumption. We now proceed to the description of a discrete version of this approach.

DISCRETE APPROXIMATIONS

4.1. Augmented Lagrangian. As an example, and to simplify the presentation, we detail our implementation for multi-shape problems in which shapes interact (through constraints) with a background, but not directly with each other. Direct interactions between shapes can be handled in a similar way. Our constrained optimization method uses the augmented Lagrangian method (see, e.g., [START_REF] Nocedal | Numerical Optimization[END_REF]). In a nutshell, in order to minimize a function u → F (u) subject to multi-dimensional equality constraints C(u) = 0, the augmented Lagrangian method consists of considering the functional

L(u) = F (u) -λ • C(u) + µ 2 |C(u)| 2 ,
in which λ lives in the dual space of the space of constraints Y , and µ is a positive real number. Each iteration of the algorithm consists in minimizing L with fixed λ and µ (our implementation using nonlinear conjugate gradient) until the gradient norm passes below some upper bound, and then in updating λ according to the rule

λ ← λ -µC(u),
before running a new minimization of L. The constant µ is increased only if needed, i.e., if the norm of the constraint did not decrease enough during the minimization. More details can be found in [START_REF] Nocedal | Numerical Optimization[END_REF]. We first apply this to identity constraints, which only require the shapes to be discretized into a sets of points. We will then discuss sliding constraints, which will require more structure in order to define normal frames to the boundary. 4.2. Identity Constraints. We consider n -1 objects, discretized into point sets, so that M k is a finite set of indices for each k. Let x

(k) j = q (k) (j) and x (k) = (x (k) 1 , . . . , x (k) m k ), for k = 1, . . . , n -1, with m k = |M k |.
We add as n-th object the background, defined on

M n = ({1} × M 1 ) ∪ • • • ∪ ({n -1} × M n-1 ). We let z (k) j = q (n) (k, j), z (k) = (z (k) 1 , . . . , z (k) m k ) and z = (z (1) , . . . , z (n-1) ) (a collection of m = m n = m 1 + . . . + m n-1 points).
Assume that end-point cost functions U 1 (x (1) ), . . . , U n-1 (x (n-1) ) are defined, typically measuring the discrepancy between each collection of points and an associated target. We assume similar functions Ũ1 (z (1) ), . . . , Ũn-1 (z (n-1) ) for the background, typically using U j = Ũj . The associated constrained optimal control problem consists in minimizing the cost functional

1 2 n k=1 1 0 u k (t) 2 V k dt + n-1 k=1 U k (x (k) (1)) + n-1 k=1 Ũk (z (k) (1)),
subject to the constraints (almost everywhere along [0, 1])

     ∂ t x (k) j = u k (x (k) j ) j = 1, . . . , m k , k = 1, . . . , n -1, ∂ t z (k) j = u n (z (k) j ), j = 1, . . . , m k , k = 1, . . . , n -1, z (k) = x (k) , k = 1, . . . , n -1.
For y and y ordered families of points in R d , let K (k) (y, y ) be the matrix formed with all d × d blocks K V k (y i , y j ), and let K (k) (y) = K (k) (y, y), for k = 1, . . . , n, where K V k is the kernel of V k . Since the problems only depend on the values taken by u 1 , . . . , u n on their corresponding point set trajectories x (1) , . . . , x (n-1) , z, the optimal vector fields take the form

u k (•) = K (k) (•, x (k) )α (k) , k = 1, . . . , n -1, u n (•) = K (n) (•, z)β,
for some families of d-dimensional vectors α (1) , . . . , α (n-1) , β. The problem can therefore be reduced to the finite-dimensional optimal control problem consisting in minimizing the cost functional

E(α, β, x, z) = 1 2 n-1 k=1 1 0 α (k) • (K (k) (x (k) )α (k) ) dt + 1 2 1 0 β • (K (n) (z)β) dt + n-1 k=1 U k (x (k) (1)) + n-1 k=1 Ũk (z (k) (1))
subject to the constraints (almost everywhere along [0, 1])

     ∂ t x (k) = K (k) (x (k) )α (k) , ∂ t z = K (n) (z)β, z (k) = x (k) , k = 1, . . . , n -1.
Extending E with the augmented Lagrangian method, we introduce coefficients λ (k) , k = 1, . . . , n -1 (where λ (k) has the same dimension as x (k) ) and µ > 0, defining

L(α, β, x, z) = 1 2 n-1 k=1 1 0 α (k) • (K (k) (x (k) )α (k) ) dt + 1 2 1 0 β • (K (n) (z)β) dt + n-1 k=1 U k (x (k) (1)) + n-1 k=1 Ũk (z (k) (1)) - n-1 k=1 1 0 λ (k) • (x (k) -z (k) ) dt + µ 2 n-1 k=1 1 0 |x (k) -z (k) | 2 dt,
which will be minimized subject to the constraints (almost everywhere along [0, 1])

∂ t x (k) = K (k) (x (k) )α (k) , ∂ t z = K (n) (z)β.
From the constraints, L can be considered as a function of α and β only, and its differential with respect to these variables can be computed via the adjoint method as follows. Denoting the co-states by p x,k , k = 1, . . . , n -1, and p z , the associated Hamiltonian is

H = n-1 k=1 1 0 p x,k • K (k) (x (k) )(α (k) ) dt + 1 0 p z • K (n) (z)β dt -L.
The computation of the gradient of L follows the same general scheme as the one described in Section 2.2 for the basic LDDMM algorithm. Given α and β and the associated trajectories x and z, one has solve the adjoint equations

∂ t p x,k = -∂ x (k) H, p x,k (1) = -∇U k (x (k) (1)), k = 1, . . . , n -1, ∂ t p z = -∂ z H, p z,k (1) = -∇ Ũk (z (k) (1)), k = 1, . . . , n -1.
The computation of the differential system gives

∂ t p x,k i = - m k j=1 ∇ 1 (p x,k i • K (k) (x (k) i , x (k) j )α (k) j ) - m k j=1 ∇ 1 (α (k) i • K (k) (x (k) i , x (k) j )p x,k j ) + 2 m k j=1 ∇ 1 (α (k) i • K (k) (x (k) i , x (k) j )α (k) j ) -(λ (k) -µ(x (k) -z (k) )), and
∂ t p z,k i = - n-1 l=1 m l j=1 ∇ 1 (p z,k i • K (n) (z (k) i , z (l) j )β (l) j ) - n-1 l=1 m l j=1 ∇ 1 (β (k) i • K (n) (z (k) i , z (l) j )p z,l j ) + 2 n-1 l=1 m l j=1 ∇ 1 (β (k) i • K (n) (z (k) i , z (l) j )β (l) j ) + n l=1 (λ (l) -µ(x (l) -z (l)
))

The gradient of L with respect to (α, β) is then deduced from the partial differentials of H with respect to these variables, yielding

∇ α (k) L = K (k) (x (k) )(α (k) -p x,k ), ∇ β L = K (n) (z)(β -p z ).
Alternatively, on may choose to use the gradient relative to the dot product on V 1 × . . . × V n , which is simply given by

∇ α (k) L = α (k) -p x,k , ∇ β L = β -p z .
The latter choice is simpler, and generally more efficient numerically.

4.3. Sliding Interface. Assume that the parameter sets M k are vertices of pure oriented, simplicial complexes T k of dimension r k < d (we will however only provide implementation details for codimension d -r k = 1). We let F k denote the set of facets of the k-th complex. We also assume that T 1 , . . . , T n-1 are disjoint and that T n is their union, T n = n-1 k=1 T k . We also let

F = n-1 k=1 F k (disjoint union).
The associated shape space is formed by functions q k : M k → R d such that q k (f ) is not degenerate (i.e., has maximal dimension) for all f ∈ F k . Each object is allowed to slide against the background. We will write x (k) = q (k) (M k ), k = 1, . . . , n -1, and z (k) = q (n) (M k ), z = (z (1) , . . . , z (n-1) ), in accordance with our previous notation. If f ∈ F k is a facet in T k ⊂ T n , we discretize ( 10) into ( 15)

N (n) (f ) • j∈f (u k (z (k) j ) -u n (z (k) j )) = 0,
where N (n) (f ) is a d×(d-r k ) matrix spanning the normal space to q (n) (f ), assumed to be defined as a smooth function of q (n) . If r k = d -1, this is always possible, since N (n) is a vector that can be taken as the cross product of z f,2 -z f,1 , . . . , z f,d -z f,1 where z f,1 , . . . , z f,d is any labeling of the vertices of q (n) (f ) ordered consistently with the orientation.

We now restrict to this case, with d = 3, so that shapes are triangulated surfaces in R 3 , as discussed in Section 2.2. For f ∈ F k and j ∈ f , we denote by e j,f the edge (z

(k) j -z (k) j )
where j and j are the other two vertices of f such that (j, j , j ) is positively oriented. Similarly, let e j,f = (z

(k) j -z (k) j ) and e j,f = (z (k) j -z (k) j ) be the two edges stemming from z (k) j so that e j,f × e j,f = 2 area(q (n) (f ))N (n) (f ) =: Ñ (n) (f )
is the area-weighted positively oriented normal to f in q (n) (M k ). Note that e j,f = e j,f -e j,f .

With this notation, we can rewrite the constraint in the form j∈f det(e j,f , e j,f , u k (z

(k) j ) -u n (z (k) j )) = 0.
holding for all f ∈ F k and k = 1, . . . , n -1.

Introducing a Lagrange multiplier λ f for each of these constraints, after reduction of the vector fields, which proceeds similarly to the identity constraints case, the augmented Lagrangian takes the form

L(α, β, x, z) = 1 2 n-1 k=1 1 0 α (k) • (K (k) (x (k) )α (k) ) dt + 1 2 1 0 β • (K (n) (z)β) dt + n-1 k=1 U k (x (k) (1)) + n-1 k=1 Ũk (z (k) (1)) - n-1 k=1 f ∈F k 1 0 (λ f Γ (k) f (x (k) , z) - µ 2 Γ (k) f (x (k) , z) 2 ) dt, with Γ (k) f (x (k) , z) := j∈f det e j,f , e j,f , K (k) (z (k) j , x (k) )α (k) -K (n) (z (k) j , z)β .
We now compute the evolution equations for the co-states, as done with identity constraints. For f ∈ F k and i ∈ M k , we have ( 16)

∂ x (k) i Γ (k) f = j∈f ∇ 1 (α (k) i • K (k) (x (k) i , z (k) j ) Ñ (n) (f )).
Denoting

δ (k) (f ) := j∈f (u k (z (k) j ) -u n (z (k) j )), if i ∈ f ∈ F k , then (17) ∂ z (k) i Γ (k) f = -e i,f × δ (k) (f ) - n-1 l=1 m l j=1 ∇ 1 ( Ñ (n) (f ) • K (n) (z (k) i , z (l) j )β (l) j ) + m k j=1 ∇ 1 ( Ñ (n) (f ) • K (k) (z (k) i , x (k) j )α (k) j ) - n-1 l=1 m l j=1 ∇ 1 (β (k) i • K (n) (z (k) i , z (l) j ) Ñ (n) (f )).
Let p x,1 , . . . , p x,n-1 and p z = (p z,1 , . . . , p z,n-1 ) be the co-states. Let γ

(k) f = λ f -µΓ (k) f . For i ∈ M k , let F k (i) = {f ∈ M k : i ∈ f }. Then ∂ t p x,k i = - N k j=1 ∇ 1 (p x,k i • K (k) (x (k) i , x (k) j )α (k) j ) - N k j=1 ∇ 1 (α (k) i • K (k) (x (k) i , x (k) j )p x,k j ) + 2 N k j=1 ∇ 1 (α (k) i • K (k) (x (k) i , x (k) j )α (k) j ) - f ∈F k γ (k) f ∂ x (k) i Γ (k) f , and
∂ t p z,k i = - N j=1 ∇ 1 (p z,i i • K (n) (z (k) i , z j )β j ) - N j=1 ∇ 1 (β (k) i • K (n) (z (k) i , z j )p z j ) + 2 N j=1 ∇ 1 (β (k) i • K (n) (z (k) i , z j )β j ) - f ∈F k γ (k) f ∂ z (k) i Γ (k) f ,
where

∂ x (k) i Γ (k) f and ∂ z (k) i Γ (k)
f are given by ( 16) and [START_REF] Dupuis | Variation Problems on Flows of Diffeomorphisms for Image Matching[END_REF].

For f ∈ F k (k = 1, . . . , n -1), we have

∂ α (k) Γ (k) f = j∈f K (k) (x (k) , z (k) j ) Ñ (n) j (f ), ∂ β Γ (k) f = - j∈f K (n) (z, z (k) j ) Ñ (n) j (f ), Letting θ (k)
j = f ∈F k :j∈f γ (k) f N (n) j (f )
, the gradient of L in α and in β is then given by

∇ α (k) L = K (k) (x (k) )(α (k) -p x,k ) -K (k) (x (k) , z)θ, ∇ β L = K (n) (z)(β -p z ) + K (n) (z, z)θ
or, taking the Hilbert gradient,

∇ α (k) L = α (k) -p x,k -K (k) (x (k) ) -1 K (k) (x (k) , z)θ, ∇ β L = β -p z + K (n) (z, z)θ.
In spite of it requiring the inversion of a linear system in the first equation, we found the latter version preferable to the L 2 gradient in our experiments.

Remarks.

Existence of constrained solutions. It is important to note that, according to [2, Theorem 1], there always exists at least one solution of ( 12)-( 13) satisfying the constraints.

Convergence to surfaces. A question naturally arising is whether our discrete approximation using triangulations converges to the smooth setting as triangles get smaller and smaller. More precisely, assume that smooth initial surfaces

S k init = q (k)
init (M k ) are triangulated, with increasingly fine triangulations q (k, ) init : M k, → R 3 , = 1, 2, . . ., where M k, labels the vertices of a simplicial complex T k, whose faces are F k, . We discuss whether minimizers (u k, , k = 1, . . . , n) of the discrete problems have a subsequence that converges to a minimizer (u k , k = 1, . . . , n) of the limit problem.

Assume that the following condition holds for the sequence of triangulations:

(i) We assume that for all k and , and for every f ∈ F k, , there exists an embedding ψ f k, :

t f k, → S (k) init = q (k) init (M k ) (where t f k, is the interior of the triangle q k, (f )) such that (ψ f k, (t f k, ), f ∈ F k, ) partitions S (k)
init up to a negligible set and max f ψ f k, -id t k, 1,∞ → 0 when → ∞. This conditions ensure that data attachment terms like those described in section 2.2.2 computed at diffeomorphic transformations ϕ k, • q (k, ) init converge to the same term computed at ϕ k • q (k) init as soon as ϕ (k, ) converges to ϕ (k) in C 1 (R 3 ). Given this, we sketch the argument leading to the consistency of the discrete approximations.

For identity constraints, one can use [2, Proposition 5], which proves that, if the triangulations are nested (every vertex at step lies on the limit surface and is also a vertex at step + 1), then one can extract, from a corresponding sequence of identity-constrained optimal vector fields, a subsequence that converges towards an identity-constrained solution of ( 12)-(13).

For sliding constraints, one cannot directly apply [2, Proposition 5], because the constraints are not nested, even when the triangulations are. To obtain a consistent approximation, we need to relax the discrete problems. More precisely, let t → (u

1 (t), . . . , u n (t)) ∈ V 1 × • • • × V n
be a minimizer of the continuous problem with sliding constraints, and let (ϕ 1 , . . . , ϕ n ) denote the corresponding flow with q

(k) (t) = ϕ k (t) • q (k)
init the corresponding deformation of q (k) init , and

N (k) (t, x) = dϕ k (t) T N (k) init (x) a normal to S k (t) = q (k) (t)(M k ) at q (k) (t, x)
. In particular we have, for every k = 1, . . . , n, every x ∈ M k , and almost every time t,

N (k) (t, x) • (u k (t, q (k) (t, x)) -u n (t, q (k) (t, x)) = 0.
Moreover, as n k=1 u k (t) 2 is constant, both u k (t) and du k (t) are α-Lipschitz for some positive constant α that does not depend on t or k.

Now let q (k, ) (t) = ϕ k (t) • q (k, )
init be the corresponding deformation of the discretization at step . Recall that N (k, ) (t, f ) denotes the unit normal to the triangle q (k, ) (t, f ). We will prove:

(ii) The discretized deformations at step satisfy the following relaxed sliding constraints

N (k, ) (f ) • j∈f (u k (z (k, ) j ) -u n (z (k, ) j
) ≤ ε at almost every t and for every face f in F k, , for a suitably chosen sequence ε > 0 going to 0 as goes to ∞. Indeed, fix a face f and an integer . Define

N (k, ) (t, f ) = dϕ k (t) T N (k, ) (0, f ) |dϕ k (t) T N (k, ) (0, f )|
for every time t. Note that assumption (i) implies that

|N (k, ) (t, f ) -N (k, ) (t, f )| ≤ γ ,
for some sequence γ , independent of f and t, and going to 0 as goes to infinity. Assumption (i) also implies that there exists a sequence η , independent of f and going to 0 as → ∞ such that for y ∈ ψ f k, (f ), and q

(k) init (x) = y, we have |N (k) init (x) -N (k, ) (0, f )| + j∈f |q (k) init (x) -z (k, ) j (0)| ≤ η .
Some triangle inequalities, Gronwall's lemma, and the fact that u k and du k are α-Lipschitz then imply

N (k, ) (t, f ) • j∈f (u k (t, z (k, ) j ) -u n (t, z (k, ) j ) ≤ ε ,
with ε = αη e α + γ going to 0 as goes to infinity.

Consequently, if (u 1, , . . . , u k, ) is a sequence of minimizers of the discretized problem at step with relaxed sliding constraints

N (k, ) (t, f ) • j∈f (u k, (t, z (k, ) j (t)) -u n, (t, z (k, ) j (t)) ≤ ε ,
we see that the infimum limit over of the respective discretized costs of (u 1, , . . . , u n, ) is smaller than or equal to the cost of a minimizer of the continuous problem with sliding constraints. So to prove that a limit point of that sequence is a minimizer of the cost for the continuous problem with sliding constraints, all we need is to check that any such limit point does satisfy the constraints.

So let (u 1, , . . . , u n, ) be a sequence of minimizers of the relaxed discrete problem at step 1 that weakly converges to (u 1 , . . . , u n ) in V 1 ו • •×V n (which is true for at least one subsequence of any minimizing sequence), then the associated flows ϕ k, and their first two space derivatives converge uniformly in time and space to the flows ϕ k and their first two derivatives.

From this it is easy to see that any sequence approximating q (k) init (x) as in assumption (i) is such that q (k, ) (t, j ) → q (k) (t, x) and N (k, ) (t, j ) → N (k) (t, x) at all times. Moreover, such a minimizing sequence must be, like all geodesics, such that n k=1 u k, (t) 2 V k is constant in time, and smaller than the cost function associated to, say, u k = 0 for all k. This implies that the vector fields u k, (t) are continuous uniformly in k, and t, which, combined with the continuity of the evaluation functionals in an RKHS implies that u k, (q (k, ) (t, j )) → u k (q (k) (t, x)) for all times. Consequently, one easily checks that each (u 1, , . . . , u n, ) satisfy a relaxed version of the continuous version of the constraints, with a precision that goes to 0 as goes to infinity. This finishes proving that the constraints are satisfied with exactitude at the limit.

Kernel derivatives. Expressions similar to ∇ 1 (n • K(x, y)α) appear at multiple times in the previous computation (for some vectors n and α). For radial kernels (K(x,

y) = G(|x -y| 2 )Id R d ), we have ∇ 1 (n • K(x, y)α) = 2G (|x -y| 2 )(n • α)(x -y), 1
It is easy to prove that such minimizers exist using the same method as that of [START_REF] Sylvain Arguillère | Shape deformation analysis from the optimal control viewpoint[END_REF], and replacing equality constraints with inequality constraints. which (slightly) simplifies the expressions.

Sliding Interface -Alternate Version. As discussed in Section 3, the sliding constraint can also be handled by introducing a new state variable N that tracks a vector (or frame) normal to the interface via (11). In the discrete case, one can discretize this equation by introducing states N (f ), f ∈ F , indexed by the facets of M , and evolving according to

∂ t N (f ) = - 1 |f | i∈f du (n) (z i ) T N (f ),
where |f | is the number of vertices in f . The sliding constraints are now expressed in terms of the state variables in a more direct way, but with a new co-state variable for the normals, bringing in an extra degree of complexity and increasing the computational cost. Note that the finite-dimensional reduction is still possible in this case, so that u (n) (•) = K (n) (•, z)β, and the evolution of the normals can be expressed in a form involving the differential of the kernel. This yields an adjoint system involving second derivatives of the kernel. We will not detail the computations in this paper, since they follow the same pattern as the other two that were already discussed (see [START_REF] Sommer | Higher-order momentum distributions and locally affine lddmm registration[END_REF] for more examples on how higher-order variables can be handled in similar contexts). Note that this alternate version of the sliding constraints is slightly more general than the one discussed in the previous section, since it does not require a definition of a normal field that smoothly depends on the manifolds.

EXPERIMENTAL RESULTS

5.1. Synthetic Example. The first example is described in Figure 1. In this synthetic example, the template has two identical balls initially close to each other. In the target, the first ball (referred to as "Ball A") gets bigger, and "impacts" the other one (referred to as "Ball B"), which assumes an oblong, non-convex shape (the target shapes slightly overlap, so that an exact homeomorphic match cannot be achieved). (i) the tangent Jacobian, which is the Jacobian determinant of the surface-to-surface transformations, and which measures the ratio between the areas of elementary surface patches at each point before and after deformation; (ii) the normal Jacobian, which is the ratio of the Jacobian determinant (of the 3D diffeomorphism) to the tangent Jacobian, and which measures the ratio between the length of an infinitesimal line element normal to the surface after and before transformation. These indexes are mapped on the deformed template image, which is close to the target.

Figure 2 compares the normal Jacobian of the shape and background deformations when using identity constraints. While shape diffeomorphisms characterize each shape transformation (uniform variation for Ball A, expansion at the top and compression otherwise for Ball B), the effect of compressing the space is clearly visible in the background deformation, when the two shapes get close to each other.

Figure 3 provides the corresponding tangent Jacobian, which is identical for shape and background transformation since we are using identity constraints.

Figures 4 and5 compare the normal and tangent Jacobians for the synthetic experiment with sliding constraints. Regarding the former, the most notable difference is with Ball B, which shows an expansion pattern at the tips in its shape diffeomorphism which is inverse of the one observed with identity constraints. One plausible explanation is that sliding constraints allow the two shapes to use translation-like motion to position themselves differently, without the need for limiting the amount of shear in the background that would have resulted from identity constraints. The second notable difference can be noted in the background diffeomorphism, in which compression is mostly observed with Ball B. In contrast with the identity constraints, the tangent Jacobians are very different between shape and background diffeomorphisms. Note that Figure 5 uses two different color scales for the left and right panels because of the strong difference between the ranges of the Jacobians in each case. The background deformation, in particular, has a huge tangent expansion around the impact location, which cannot be observed in the shape deformations. Note that both patterns in the sliding case are very different from the one that was observed in the identity case.

For comparison purposes, Figure 6 provides the result of the LDDMM algorithm using a single diffeomorphism. One observes a very strong compression effect for the normal jacobian resulting in an expansion observed on the tangent jacobian on Ball B, that was not observed in any of the constrained examples. The nice uniform expansion in Ball A that could be observed in the sliding constraint case is not observed either. 5.2. Subcortical Structures. We now describe an example mapping a group of three subcortical structures: hippocampus, amygdala and entorhinal cortex (ERC). The template and target sets are represented in Figure 7. One can observe shape changes in each structure, combined with a significant displacement of the ERC relative to the other two structures when comparing template to target. Because the structures were segmented independently, there is some overlap between the target hippocampus and amygdala.

Figures 8 and9 provide the normal and tangent Jacobian obtained with identity constraints, while Figures 10 and 11 The hippocampus is the central shape, with the amygdala on its left and the ERC on its right.

DISCUSSION

The previous approach provides a solution, using constrained optimal control, of the important issue of dealing with multiple objects with varying deformation properties for registration. We have focused on surface matching, numerically dealing with constraints using an augmented Lagrangian method. Note that a similar approach was introduced for plane curves in [START_REF] Sylvain Arguillère | Shape deformation analysis from the optimal control viewpoint[END_REF].

The formulation is quite general and can accommodate constraints in various forms, including the examples discussed in Section 3. The investigation of these additional applications will be the subject of future work. One of the limitations of the present implementation is the slow convergence of the augmented Lagrangian procedure, for which each minimization step is, in addition, high dimensional and computationally demanding. One possible alternative can be based on solving the optimality conditions ( 14) (which hold in the discrete case) by means of a numerical shooting method. This approach has, however, its own numerical challenges, because solving [START_REF] Cotter | Singular solutions, momentum maps and computational anatomy[END_REF] requires the determination of λ such that the last equation (constraint) is satisfied, and this leads to a possibly ill-posed problem for systems in large dimension (see [START_REF] Sylvain Arguillère | Shape deformation analysis from the optimal control viewpoint[END_REF] for additional details).

We have illustrated our examples using deformation markers derived from the jacobian determinant. This markers are routinely used in shape analysis studies and led to important conclusion in computational anatomy. When dealing with multiple shapes, however, figures 6 and 12 show that, when using the classical LDDMM method with multiple shapes, these markers becomes as much, if not more, influenced by interactions between the shapes as by the changes in the shapes themselves. For this reason, multi-shape computational anatomy studies have applied registration methods separately to each shape, without ensuring that the obtained diffeomorphisms are consistent with each other. This limitation is addressed in the present paper, in which we exhibit deformation markers that are meaningful in describing tangential and normal surface stretching, while being consistently associated to a global transformation of the space.

FIGURE 1 .

 1 FIGURE 1. Template and target shapes for synthetic example

FIGURE 2 .

 2 FIGURE 2. Three views of the normal Jacobian with identity constraints: shape diffeomorphisms (left) and background diffeomorphism (right).

FIGURE 3 .

 3 FIGURE 3. Tangential Jacobian: shape and background diffeomorphisms (identity constraints).

  provide this information for sliding constraints. The two types of constraints provide similar deformation indices, especially for the normal jacobians (Figures8 and

FIGURE

  FIGURE 4. Three views of the normal Jacobian with sliding constraints: shape diffeomorphisms (left) and background diffeomorphism (right).

4 .

 4 FIGURE 4. Three views of the normal Jacobian with sliding constraints: shape diffeomorphisms (left) and background diffeomorphism (right).

FIGURE 5 .

 5 FIGURE 5. Three views of the tangent Jacobian with sliding constraints: shape diffeomorphisms (left) and background diffeomorphism (right).

FIGURE 6 .

 6 FIGURE 6. Three views of the Normal (left) and tangential Jacobians (right) when using a single diffeomorphism.

  ). The deformation patterns associated to using a single diffeomorphism (Figure12) are significantly different, though, exhibiting very strong compression, for example, where shapes are close to each other.

FIGURE 7 .

 7 FIGURE 7. Template (blue) and target (red) shapes for subcortical structures. The hippocampus is the central shape, with the amygdala on its left and the ERC on its right.

FIGURE 8 .

 8 FIGURE 8. Three views of the normal Jacobian with identity constraints: First Row: Shape diffeomorphisms; Second Row: Background diffeomorphism.

FIGURE 9 .

 9 FIGURE 9. Three views of the tangent Jacobian with identity constraints: Shape and background diffeomorphisms.

FIGURE 10 .

 10 FIGURE 10. Three views of the normal Jacobian with sliding constraints: First Row: Shape diffeomorphisms; Second Row: Background diffeomorphism.

[ 11 ]

 11 Can Ceritoglu, Kenichi Oishi, Xin Li, Ming-Chung Chou, Laurent Younes, Marilyn Albert, Constantine Lyketsos, Peter van Zijl, Michael I Miller, and Susumu Mori. Multi-contrast large deformation diffeomorphic metric mapping for diffusion tensor imaging. Neuroimage, 47(2):618-627, 2009. [12] Gary E. Christensen, Richard D. Rabbitt, and Michael I. Miller. Deformable templates using large deformation kinematics. IEEE Transactions on Image Processing, 5(10):1435-1447, 1996. [13] CJ Cotter and DD Holm. Discrete momentum maps for lattice epdiff. In Temam and Tribbia, editors, Handbook of Numerical Analysis, pages 247-278. North-Holland, 2009.

FIGURE 11 .

 11 FIGURE 11. Three views of the tangent Jacobian with sliding constraints: First Row: Shape diffeomorphisms; Second Row: Background diffeomorphism.

This work was partially supported by the ONR award N000140810606.

S. ARGUILL ÈRE, E. TR ÉLAT, A. TROUV É, AND L. YOUNES