Metrics, quantization and registration in varifold spaces
Hsi-Wei Hsieh Nicolas Charon







Department of Applied Mathematics and Statistics, Johns Hopkins University



hhsieh9@jhu.edu





Department of Applied Mathematics and Statistics, Johns Hopkins University
charon@cis.jhu.edu











This paper is concerned with the theory and applications of varifolds to the representation, approximation and diffeomorphic registration of shapes. One of its purpose is to synthesize and extend several prior works which, so far, have made use of this framework mainly in the context of submanifold comparison and matching. In this work, we instead consider deformation models acting on general varifold spaces, which allows to formulate and tackle diffeomorphic registration problems for a much wider class of geometric objects and lead to a more versatile algorithmic pipeline. We study in detail the construction of kernel metrics on varifold spaces and the resulting topological properties of those metrics, then propose a mathematical model for diffeomorphic registration of varifolds under a specific group action which we formulate in the framework of optimal control theory. A second important part of the paper focuses on the discrete aspects. Specifically, we address the problem of optimal finite approximations (quantization) for those metrics and show a Γ-convergence property for the corresponding registration functionals. Finally, we develop numerical pipelines for quantization and registration before showing a few preliminary results for one and two-dimensional varifolds.    

49Q20 49M25 58E50 68T10




Metrics, quantization and registration in varifold spaces
    Nicolas Charon
    
=========================================================


















.



§ INTRODUCTION

Shape is a bewildering notion: while simultaneously intuitive and ubiquitous to many scientific areas from pure mathematics to biomedicine, it remains very challenging to pin down and analyze in a systematic way. The goal of the research field known as shape/pattern analysis is precisely to provide solid mathematical and algorithmic frameworks for tasks such as automatic comparison or statistical analysis in ensembles of shapes, which is key to many applications in computer vision, speech and motion recognition or computational anatomy, among many others.

What makes shape analysis such a difficult and still largely open problem is, on the one hand, the numerous modalities and types of objects that can fall under this generic notion of shape but also the fundamental nonlinearity that is an almost invariable trait to most of the shape spaces encountered in applications. As a result, the seemingly simple issue of defining and computing distances or means on shapes is arguably a research topic of its own, which has generated countless works spanning several decades and involving concepts from various subdisciplines of mathematics. Among many important works, the model of shape space laid out by Grenander in <cit.> is especially relevant to the present paper. The underlying principle is to build distances between shapes which are induced by metrics on some deformation groups acting on those shapes. This approach has the advantage (at a theoretical level at least) of shifting the problem of metric construction from the many different cases of shape spaces to the single setting of deformation groups. One of the fundamental requirement is the right-invariance of the metrics on those groups; finding the induced distance between two given shapes then reduces to determining a deformation of minimal cost in the group, in other words to solving a registration problem.


Besides usual finite-dimensional groups like rigid of affine transformations, there is in fact a lot of practical interest in applying such an approach with groups of "large deformations", specifically groups of diffeomorphisms. This has triggered the exploration of right-invariant metrics over diffeomorphism groups. The Large Deformation Diffeomorphic Metric Mapping (LDDMM) model pioneered in <cit.> is one of such framework that defines Riemannian metrics for diffeomorphic mappings obtained as flows of time-dependent vector fields (c.f. the brief presentation of Section <ref>). In this setting, registering two shapes can be generically formulated as an optimal control problem, the functionals to optimize being typically a combination of a deformation regularization term given by the LDDMM metric on the group and a fidelity term that enforces (approximate) matching between the two shape objects. Applications of this model have been widespread in particular within the field of computational anatomy, due to the ability to adapt it to various data structures including landmarks, 2D and 3D images, tensor fields... see e.g. <cit.> for recent reviews.


Interestingly, this line of work has also been drawing many useful concepts from the seemingly distant area of mathematics known as geometric measure theory <cit.>. The key idea of representing shapes (submanifolds) as measures or distributions has been instrumental in the theoretical study of Plateau's problem on minimal surfaces and more generally in calculus of variations. It can also prove effective for computational purposes, in problems such as discrete curvature approximations <cit.> or estimation of shape medians <cit.>. With regard to the aforementioned deformation analysis problems, the potential interest of geometric measure theory has been identified early on in the works of <cit.>. Indeed, LDDMM registration of objects like geometric curves or surfaces requires fidelity terms independent of the parametrization of either of the two shapes. On the practical side, this means that one cannot usually rely on predefined pointwise correspondences between the vertices of two triangulated surfaces for instance, which makes the registration problem significantly harder than in the case of labelled objects such as landmarks or images. 

The embedding of unparametrized shapes into measure spaces provides one possible way to address the issue, by constructing parametrization-invariant fidelity metrics as restrictions of metrics on those measure spaces themselves. Several competing approaches have been introduced, each relying on embeddings into different spaces of generalized measures: <cit.> are based on the representation of oriented curves and surfaces as currents, <cit.> and <cit.> extended this model to the setting of unoriented and oriented varifolds, while <cit.> considers the higher-order representation of normal cycles, see also the recent survey <cit.>. One common feature to all those works, however, is that they are focused primarily on registration of curves or surfaces. In other words, the use of current, varifolds or normal cycles confines to the computation of a fidelity metric to guide registration algorithms but the deformation model itself remains tied to the curve/surface setting or equivalently, in the discrete situation, to objects described by point set meshes. 

The guiding theme and main objective of this paper is to investigate an alternative framework that, in contrast with those prior works, would formulate the deformation model as well as tackle the registration problem directly in these generalized measure spaces: we focus specifically on the (oriented) varifold setting of <cit.>. There are several arguments for the interest of such an approach 
but in our point of view, the primary motivation lies in the fact that, varifolds being more general than submanifolds, the proposed framework allows to extend large deformation analysis methods to a range of new geometric objects while giving more flexibility to deal with some of the flaws which are commonplace in shapes segmented from raw data. As a proof of concept, our recent work <cit.> considered the simple case of registration of discrete one-dimensional varifolds. Building on these preliminary results, the present paper intends to provide a thorough and general study of the framework.


The specific contributions and organization of this paper are the following. First, we propose a comprehensive study of the class of kernel metrics on varifold spaces initiated in <cit.>, in particular by examining the required conditions to recover true distances between all varifolds (as opposed to the subset of rectifiable varifolds) and comparing the resulting topologies with some standard metrics on measures. This is presented in Section <ref> after the brief introduction to the notion of oriented varifold of Section <ref>. In Section <ref>, we discuss the action of diffeomorphisms and from there derive a formulation of LDDMM registration of general varifolds, for which we show the existence of solutions and derive the Hamiltonian equations associated to the corresponding optimal control problem. Section <ref> addresses the issue of quantization in varifold space, namely of approximating any varifold as a finite sum of Dirac masses. We consider a novel approach in this context, that consists in computing projections onto particular cones of discrete varifolds. We then prove the Γ-convergence of the corresponding approximate registration functionals. In Section <ref>, we derive the discrete version of the optimal control problem and optimality equations, from which we deduce a geodesic shooting algorithm for the diffeomorphic registration of discrete varifolds. Finally, results on 1- and 2- varifolds are presented in Section <ref>, emphasizing the potentiality of the approach to tackle data structures which are typically challenging for previous algorithms that are designed for point sets and meshes. The MATLAB code used in this paper is made available under the GNU general public license through the Github repository <https://github.com/charoncode/Var_LDDMM>.




§ THE SPACE OF ORIENTED VARIFOLDS


The concept of varifold was originally developed in the context of geometric measure theory by <cit.>, <cit.> and <cit.> for the study of Plateau's problem on minimal surfaces. The interest in registration and shape analysis was evidenced in <cit.>. In those works, varifolds provide a convenient representation of geometric shapes such as rectifiable curves and surfaces and an efficient approach to define and compute fidelity terms for registration, or to perform clustering, classification in those shape spaces. The main purpose of this section is to introduce varifolds in this latter context. The case of non-oriented shapes was thoroughly investigated in <cit.>. Later on, the generalized framework of oriented varifold was proposed in <cit.> but only for objects of dimension or co-dimension one. In the following, we provide a fully general presentation of oriented varifolds and their properties, that also does not specifically focus on the case of rectifiable varifolds as these previous works did. Although we assume here that all the considered shapes are oriented, we emphasize that the non-oriented framework of <cit.> can be recovered almost straightforwardly through adequate choices of orientation-invariant kernels as we shall briefly point out later on.



 §.§ Definition
 

The underlying principle of varifolds is to extend measures of ^n by incorporating an additional tangent space component. In this work, we will consider such spaces to be oriented. Thus, for a given dimension 0≤ d ≤ n, we first need to introduce the set of all possible d-dimensional oriented tangent spaces in ℝ^n:

 The d-dimensional oriented Grassmannian G_d^n is the set of all oriented d-dimensional linear subspaces of ℝ^n.
 
 
 
 

The oriented Grassmannian is a compact manifold of dimension d(n-d) which can be identified to the quotient SO(n)/(SO(d) × SO(n-d)). It is also a double cover of the (non-oriented) Grassmannian G_d^n of d-dimensional subspaces of ^n. For practical purposes, a more convenient representation of G_d^n is the one detailed in the following remark.



Given T ∈G_d^n, there exists a basis {u_i }_i=1,…,d∈^n× d of T such that [u_1,⋯,u_d] has consistent orientation with T. Then the following map, called the oriented Plücker embedding, is well defined and injective,

    i_P:  G_d^n    ↦{ξ∈Λ^d(ℝ^n) : |ξ| =1}
    
        T    ↦u_1 ∧⋯∧ u_d/|u_1 ∧⋯∧ u_d |.

This allows to identify G^n_d as a subset of the unit sphere of Λ^d(ℝ^n) which inherits the topology of the inner product on Λ^d(ℝ^n). We remind that this inner product is defined for any ξ=ξ_1 ∧…∧ξ_d, η=η_1 ∧…∧η_d in Λ^d(ℝ^n) by the determinant of the Gram matrix: 

    ⟨ξ,η⟩ = (ξ_i ·η_j)_i,j=1,…,d

Through this identification, one can also define the action of linear transformations on G^n_d as follows 

    A · T := A u_1 ∧⋯∧ Au_d/|A u_1 ∧⋯∧ Au_d |

for any T ∈G^n_d and A: ℝ^n ↦ℝ^n a linear invertible map.


Similar to the definition of classical varifolds in <cit.>, we define oriented varifolds as measures on ℝ^n ×G^n_d. 


 An oriented d-varifold μ on ℝ^n is a nonnegative finite Radon measure on the space ℝ^n ×G^n_d. Its weight measure |μ| is defined by |μ|(A) := μ(A×G^n_d) for all Borel subset A of ^n. We denote by 𝒱_d the space of all oriented d-varifolds.

In the rest of the paper, with a slight abuse of vocabulary, we will often use the word varifold instead of oriented varifold for the sake of concision. Recall that from the Riesz representation theorem, we can alternatively view any varifold μ as a distribution, i.e. an element of the dual space C_0(ℝ^d ×G^n_d)^*, where C_0(ℝ^d ×G^n_d) denotes the set of continuous functions vanishing at infinity on ℝ^d ×G^n_d. It is defined for any test function ω∈ C_0(ℝ^d ×G^n_d) by:

    (μ |ω) ≐∫_ℝ^n ×G^n_dω(x,T) dμ(x,T).

As an additional note, another useful representation of a general varifold in 𝒱_d can be obtained by the disintegration theorem (see <cit.> Chap. 2). Namely, if μ∈𝒱_d, for |μ|-almost every x in ^n, there exists a probability measure ν_x on G^n_d such that x ↦ν_x is |μ|-measurable and we can write

    (μ |ω) = ∫_^n∫_G^n_dω(x,T) dν_x(T) d|μ|(x).

In other words, the varifold μ can be decomposed as its weight measure on ^n together with a family of tangent space probability measures on the Grassmannian at the different points in the support of |μ|. This is usually referred to as the Young measure representation of μ.



 §.§ Diracs and rectifiable varifolds

There are a few important families of varifolds which will be relevant for the following. First of those are the Diracs. For x ∈ℝ^n and T ∈G^n_d, the associated Dirac varifold δ_(x,T) acts on functions of C_0(ℝ^n ×G^n_d) by the relation

    (δ_(x,T)|ω) = ω(x,T),  ∀ω∈ C_0(ℝ^n ×G^n_d).

δ_(x,T) can be viewed as a singular particle at position x that carries the oriented d-plane T.


A second particular class is the one of rectifiable varifolds, which are in essence the varifolds representing an oriented shape of dimension d. More precisely, given an oriented d-dimensional submanifold X of ℝ^n of finite total d-volume, denoting by T_X(x) ∈G^n_d the oriented tangent space at x ∈ X, one can associate to X the varifold μ_X, which is defined for all Borel subset B⊂^n ×G^n_d by μ_X(B) = ℋ^d({x ∈ X | (x,T_X(x)) ∈ B }). Here, ℋ^d is the d-dimensional Hausdorff measure on ^n, i.e. the measure of d-volume of subsets of ^n (we refer the reader to <cit.> for the precise construction and properties of Hausdorff measures). It is then not hard to see that, as an element of C_0(ℝ^n ×G^n_d)^*,  

    (μ_X|ω)    = ∫_ℝ^d ×G^n_dω(x,T) d μ_X(x,T) 
       = ∫_X ω(x,T_X(x)) d ℋ^d(x).

Such a representation X ↦μ_X can be extended to slightly more general objects known as oriented rectifiable sets. A subset X of ℝ^n is said to be a countably ℋ^d-rectifiable set if ℋ^d(X ∖∪_j=1^∞F_j(ℝ^d)) = 0, where F_j: ℝ^d ↦ℝ^n are Lipschitz function for all j (c.f. <cit.>). We say that (X,T_X) is an oriented rectifiable set if X is a countably d-rectifiable set and T_X:X ↦G^n_d is a ℋ^d-measurable function such that for ℋ^d-a.e. x ∈ X, T_X(x) is the approximate tangent space of X at x with specified orientation. Rectifiable subsets include both usual submanifolds but also piecewise smooth objects like polyhedra. Given any oriented rectifiable set (X, T_X), we can associate a varifold that we also write μ_X given again by (<ref>). The set of those μ_X will be referred to as the rectifiable oriented varifolds in this paper (note that this is actually more restrictive than the standard definition of rectifiable varifold in the literature which also incorporates an additional multiplicity function). 

Rectifiable varifolds still make a very "small" subset of 𝒱_d: indeed, in the Young measure representation of (<ref>), we have in this case the very particular constraint that probability measures ν_x are Dirac masses, specifically ν_x = δ_T_X(x). 





§ METRICS ON VARIFOLDS


In this section, we address the issue of defining adequate metrics on the space 𝒱_d. After reviewing some classical metrics and their limitations for the specific applications of this work, we turn to metrics defined through positive definite kernels, for which we extend previous constructions introduced in e.g. <cit.> and derive the most relevant properties of this class of distances.   




 §.§ Standard topologies and metrics on 𝒱_d


As a measure/distribution space, 𝒱_d can be equipped with various topologies and metrics, several of which have been regularly used in various contexts. We discuss a few of those below. 

   
  ∙   
mass norm: with the previous identification of measures in 𝒱_d with elements of the dual C_0(ℝ^n ×G^n_d)^*, one can define the following dual metric on 𝒱_d:

    d_op(μ,ν) ≐sup_|ω|_∞≤ 1 (μ -ν|ω),  ∀μ∈𝒱_d.

where |ω|_∞≐sup_ℝ^n ×G^n_d |ω|. This metric is generally too strong for applications in shape analysis and leads to a discontinuous behavior. Indeed, one can easily verify that for any two Dirac masses δ_(x,T) and δ_(x',T'), d_op(δ_(x,T),δ_(x',T')) = 2 whenever (x,T) ≠ (x',T').   


  ∙
weak-* topology: a sequence of d-varifolds {μ_i}_i converges to μ∈𝒱_d in the weak-* topology (denoted by μ_i ∗⇀μ) if and only if for all ω∈ C_c(ℝ^d ×G^n_d) (continuous compactly supported function) 

    lim_i →∞(μ_i|ω) = (μ|ω).

In fact, the weak-* topology on 𝒱_d can be metrized by the following distance:

    d_*(μ,ν) = ∑_k ∈ℕ 2^-k |(μ-ν|ω_k)|,

where {ω_k}_k∈ℕ is a dense sequence in C_c(ℝ^n ×G^n_d).


  ∙
Wasserstein metric: the Wasserstein-1 distance of optimal transport can be expressed in its Kantorovitch dual formulation <cit.> as

    d_Wass^1(μ,ν) ≐sup_Lip(ω) ≤ 1 |(μ-ν|ω)|.

where the sup is taken over all Lipschitz regular functions on ℝ^n ×G^n_d with Lipschitz constant smaller than one. This metric is however well-suited for measures with the same total mass. Several recent works <cit.> have instead proposed generalized Wasserstein distances derived from unbalanced optimal transport. 


  ∙
Bounded Lipschitz metric: similar to the previous, the bounded Lipschitz distance (sometimes referred to as the flat metric) on 𝒱_d is defined by

    d_BL(μ,ν) ≐sup_ω_∞, Lip(ω) ≤ 1 |(μ-ν|ω)|.

It can be shown (cf. Ch 8 in <cit.>) that d_BL metrizes the narrow topology on 𝒱_d, namely the topology for which a sequence (μ_i) converges to μ if and only if lim_i →∞(μ_i|ω) = (μ|ω) for all bounded continuous functions ω.


Clearly, the narrow topology is stronger than the weak-* topology. Furthermore, it is also well known that d_BL locally metrizes the weak-* topology on 𝒱_d, namely:

Let μ and {μ_i}_i be varifolds such that the sequence {μ_i}_i is tight. Then μ_i ∗⇀μ if and only if d_BL(μ_i,μ) → 0.


Since d_BL metrizes the narrow topology, it suffices to show that μ_i converges to μ in the narrow topology. Let ω be a bounded continuous function defined on ℝ^n ×G^n_d and ε>0. By the tightness property, we may choose a compact set K⊂ℝ^n ×G^n_d such that μ(K^c)+ sup_i μ_i(K^c)< ε/2ω_∞. Let B be an open ball that contains K. Define

    η(x,T)  ≐{[        ω(x,T),   if (x,T) ∈ K;            0 , if (x,T) ∈ B^c ].

From Tietz extension theorem, there exists a continuous extension ω̃ of η on ℝ^n ×G^n_d such that ω̃|_K =ω|_K and ω̃∈ C_c(ℝ^n ×G^n_d). This implies that

    | ∫_ℝ^n ×G^n_dω d(μ_i-μ) | 
         ≤∫_ℝ^n ×G^n_d|ω- ω̃| d (μ_i+μ)
        + | ∫_ℝ^n ×G^n_dω̃ d (μ_i-μ) |
       ≤| ∫_ℝ^n ×G^n_dω̃ d (μ_i-μ) | + ε.

Taking lim sup on both sides, we see that 

    lim sup_i →∞| ∫_ℝ^n ×G^n_dω d(μ_i-μ) | < ε.

Since ε is arbitrary, we obtain that μ_i converges to μ in the narrow topology. 

As a direct consequence of Proposition <ref>, we have in particular that weak-* convergence and convergence in d_BL are equivalent if one restricts to varifolds that are supported in a fixed compact subset of ℝ^n ×G^n_d. Note also that a very similar result to Proposition <ref> holds when replacing the bounded Lipschitz distance by generalized Wasserstein metrics, as proved in <cit.>. 

The above metrics on varifolds all originate from classical ones in standard measure theory. Unlike the mass norm, Wasserstein and bounded Lipschitz metrics have nice theoretical properties in terms of shape comparison. However, for the purpose of diffeomorphic registration that we shall tackle below, one needs metrics that are easy to evaluate numerically. This is typically not the case of d_Wass^1 and d_BL expressed above as there is no straightforward way to compute the corresponding suprema over the respective sets of test functions. One line of work has been considering approximations of optimal transport distances with e.g. entropic regularizers for which Sinkhorn-based algorithms can be derived, see for instance the recent work <cit.>. In this paper, we focus on the alternative approach previously developed for currents in <cit.> and unoriented varifolds in <cit.> which instead relies on particular Hilbert spaces of test functions, as we detail in the next section.




 §.§ Kernel metrics


In this section, we start by defining a general class of pseudo-metrics on 𝒱_d based on positive definite kernels and their corresponding reproducing kernel Hilbert space (RKHS). We will then study sufficient conditions on such kernels to recover true metrics before examining the relationship between those kernel metrics and the ones of Section <ref>.



  §.§.§ Kernels for varifolds

We refer the reader to <cit.> for a presentation of the construction and main properties of positive kernels and Reproducing Kernel Hilbert Spaces which we do not recall in detail here for the sake of concision. In the context of varifolds, we are interested in defining positive definite kernels on the product ℝ^n ×G^n_d. Along the lines of previous works like <cit.>, we restrict to separable kernels for which we have:

 
Let k^pos and k^G be continuous positive definite kernels on ℝ^n and G^n_d respectively. Assume in addition that for any x ∈ℝ^n, k^pos(x,·) ∈ C_0(ℝ^n). Then k:= k^pos⊗ k^G is a positive definite kernel on ℝ^n ×G^n_d and the RKHS W associated to k is continuously embedded in C_0(ℝ^n ×G^n_d) i.e. there exists c_W>0 such that for any ω∈ W, we have ω_∞≤ c_W ω_W.

We recall that the tensor product kernel has the exact expression k((x,T),(x',T')) = k^pos(x,x') k^G(T,T'). The proof of Proposition <ref> is a straightforward adaptation of the same result for unoriented varifolds (cf. <cit.> Proposition 4.1).  


To simplify the rest of the presentation and in the perspective of later numerical considerations, we will also assume specific forms for k^pos and k^G, namely that k^pos is a translation/rotation invariant radial kernel k^pos(x,y) = ρ(|x-y|^2),  ∀ x,y ∈ℝ^n, with ρ(0)>0, and k^G(S,T) = γ(⟨ S,T ⟩),  ∀ S,T ∈G^n_d where ⟨·, ·⟩ is the inner product on G^n_d inherited from Λ^d(ℝ^n) introduced in remark <ref>. These assumptions are quite natural as they will eventually induce metrics on varifolds invariant to the action of rigid motion, as we shall explain later. Note that the unoriented framework of <cit.> can be also recovered in this setting by simply restricting to orientation-invariant kernels k^G i.e. such that γ(-t)=γ(t) for all t.    


Now, if we let ι_W:W ↪ C_0(ℝ^d ×G^n_d) be the continuous embedding given by Proposition <ref> and ι_W^* its adjoint, for any μ∈ C_0(ℝ^n ×G^n_d)^*, we have 

    (ι^* μ| ω) = ∫_ℝ^d ×G^n_dω(x,T) dμ(x,T),  ∀ω∈ W.

With (<ref>), we may identify μ as an element of the dual RKHS W^*. Note that ι_W^* is not injective in general, in other words one can have μ = μ' in W^* but μ≠μ' in C_0(ℝ^n ×G^n_d)^*.

In any case, one can compare any two varifolds μ,μ' ∈𝒱_d through the Hilbert norm of W^* by defining: 

    d_W^*(μ,μ')^2 = μ-μ'_W^*^2 = μ_W^*^2 -2⟨μ,μ' ⟩_W^* + μ'_W^*^2

where we use the small abuse of notation of writing μ and μ' instead of ι_W^* μ and ι_W^* μ' on the two right hand sides. Due to the potential non-injectivity of ι_W^*, in general d_W^* only induces a pseudo-metric on 𝒱_d. 

The main advantage of this construction is that d_W^* can be expressed more explicitly based on the reproducing kernel property of W. Indeed, given any μ and ν in 𝒱_d, the inner product between them is given by 

    ⟨μ, μ' ⟩_W^*   = ∫_(ℝ^d ×G^n_d)^2 k^pos(x,x') k^G(T,T') dμ(x,T) d μ'(x',T') 
       = ∫_(ℝ^d ×G^n_d)^2ρ(|x-x'|^2) γ(⟨ T,T' ⟩) dμ(x,T) d μ'(x',T')

for kernels selected as in Remark <ref>.




  §.§.§ Characterization of distances

As mentioned above, d_W^* is a priori a pseudo-distance between varifolds. It's a natural question to ask under which conditions it leads to an actual distance. 

Most past works have addressed this question focusing on the case of varifolds representing submanifolds and reunion of submanifolds <cit.>. We can first provide an extension of these results to the general case of oriented rectifiable varifolds. A key notion for the rest of this section is the one of C_0-universality of kernels:

A positive definite kernel k on a metric space ℳ is called C_0-universal when its RKHS is dense in C_0(ℳ) for the uniform convergence topology.  

C_0-universality has been studied in great length in such works as <cit.>. In particular, one can provide characterizations of C_0-universality for certain classes of kernels and spaces ℳ. In the case of translation-invariant kernels on ℳ=^n for instance, it has been established that C_0-universal kernels are the ones which can be expressed through the Fourier transform of finite Borel measures with full support on ^n, which includes: compactly-supported kernels, Gaussian kernels, Laplacian kernels... With the previous definition, we have the following sufficient condition: 


Suppose k^pos is a C_0-universal kernel on ^n, γ(1)>0 and γ(t) ≠γ(-t),  ∀ t ∈ [-1,1]. Let (X,T(·)) and (Y,S(·)) be two oriented ℋ^d-rectifiable sets with ℋ^d(X), ℋ^d(Y)< ∞. If μ_X-μ_Y _W'=0, then ℋ^d(X△ Y) =0 and T = S  ℋ^d-a.e. 

The full proof can be found in the Appendix. Note that the first part of the proof directly gives an equivalent statement for unoriented rectifiable varifolds (if one instead assumes γ(t) = γ(-t) for all t), generalizing the result of <cit.>.

However, the previous proposition does not necessarily lead to a distance on the full space 𝒱_d. Counter-examples in the case d=1 are discussed for example in <cit.>. To recover a true distance on 𝒱_d, one needs the previous map ι_W^* or equivalently the map

    μ↦∫_ℝ^d ×G^n_d k(·,(y,T)) d μ(y,T),  μ∈ C_0(ℝ^d ×G^n_d)^*

to be injective. As follows from Theorem 6 in <cit.>, this is in fact guaranteed when the kernel k on the product space ℝ^n ×G^n_d is C_0-universal, specifically


The pseudo-distance d_W^* induces a distance between signed measures of ℝ^n ×G^n_d if and only if k is C_0-universal on ℝ^n ×G^n_d. In particular, a sufficient condition for d_W^* to be a distance on 𝒱_d is that k^pos and k^G are C_0-universal kernels on ^n and G^n_d respectively.

Note that these conditions are more restrictive than in Theorem <ref>. To our knowledge, there is no simple characterization for general C_0-universal kernels on the Grassmannian. However, within the setting of Remark <ref>, one easily constructs C_0-universal kernels by restriction (based on the Plücker embedding) of C_0-universal kernels defined on the vector space Λ^d(^n).




  §.§.§ Comparison with classical metrics

We now study more precisely the topology induced by the (pseudo) distance d_W^* on 𝒱_d in comparison with the ones defined in Section <ref>. First of all, we observe that, for any ω∈ W with ω_W≤ 1, one must have ω_∞≤ c_W, where c_W is the embedding constant of Proposition <ref>. Thus, for any μ and μ' in 𝒱_d, we have

    μ-μ'_W^* = sup_ω∈ W, ω_W ≤ 1∫_ℝ^d ×G^n_d ω d(μ-μ') ≤ c_W d_op(μ,μ').

From the above inequalities we see that convergence in d_op implies convergence in d_W^*.



With more assumptions on the regularity of the kernel k, namely if W is continuously embedded in C_0^1(ℝ^d ×G^n_d), following a similar reasoning as above, one obtains the bound μ-μ'_W^*≤ c_W d_BL(μ,μ').  


Suppose μ_i  converges to μ in narrow topology. Since the map (ν_1,ν_2) ↦ν_1 ⊗ν_2 is continuous with respect to the narrow topology, we have

    μ_i_W^*^2    = ∫_(ℝ^d ×G^n_d)^2 k((x,S),(y,T)) dμ_i(x,S) d μ_i(y,T) 
       →∫_(ℝ^d ×G^n_d)^2 k((x,S),(y,T)) dμ(x,S) d μ(y,T) 
       = μ_W^*^2,

as i →∞.
Also, it's clear that lim_i →∞⟨μ_i, μ⟩_W^*→μ_W^*^2 and hence μ_i →μ with respect to d_W^*. To summarize the discussion above:


Let {μ_i}_i and μ be varifolds in 𝒱_d and assume that μ_i →μ with respect to the operator norm or the narrow topology, then μ_i →μ in W^*.


We emphasize that the result of Proposition <ref> only requires the assumptions of Proposition <ref> and thus holds whether ι is injective or not. 


As for the weak-* topology, with the C_0-universality assumption of Theorem <ref> and restricting to varifolds with bounded total mass, we show that d_W^* induces a topology stronger than weak-* convergence:


If k is C_0-universal, then the topology induced by d_W^* is finer than the weak-* topology on 𝒱_d,M≐{μ∈𝒱_d  s.t |μ|(ℝ^n) ≤ M } for any fixed M>0.




 Let {μ_i}_i and μ be varifolds in 𝒱_d,M and assume that lim_i →∞d_W^*(μ_i,μ) = 0. For any f ∈ C_0(ℝ^d ×G^n_d) and ε>0, there exists a g ∈ W such that g-f < ε/2M. Then we obtain that μ_i ∗⇀μ from the following inequalities:

    |(μ_i-μ|f)| ≤ |(μ_i| f-g)| + |(μ| g-f)| + |(μ_i - μ|g)| ≤ε + μ_i -μ_W^*g_W.




Note that the topology induced by d_W^* may be strictly finer on 𝒱_d,M. Indeed, if ρ(0),γ(1)>0, consider μ_i = δ_(x_i,S), where lim_i →∞|x_i| = ∞ and S∈G^n_d fixed. Then μ_i ∗⇀ 0 while μ_i_W^*^2 =  ρ(0) γ(1) > 0  for all i. Yet, by combining Propositions <ref>, <ref> and <ref>, we have the following 


Let M>0 and K ⊂^n ×G^n_d be a compact subset. If k is C_0-universal, then d_W^* metrizes the weak-* convergence of varifolds on 𝒱_d,M,K≐{μ∈𝒱_d  s.t |μ|(ℝ^n) ≤ M,  supp(μ)⊂ K }.


In summary, C_0-universality provides a sufficient condition to obtain actual distances between varifolds that can be expressed based on the kernel function. Furthermore, the resulting topology is locally equivalent to the weak-* topology as well as the topology induced by the bounded Lipschitz distance. This equivalence will be of importance in Section <ref>. 




§ DEFORMATION AND REGISTRATION OF VARIFOLDS
 

Having defined a way of comparing general varifolds through the above kernel metrics d_W^*, our goal is now to focus on deformation models for those objects in order to formulate and study the diffeomorphic registration problem on 𝒱_d. 



 §.§ Deformation models


In this section, we discuss different models for how varifolds can be transported by a diffeomorphism of ^n, in other words what are possible group actions of the diffeomorphism group Diff(^n) on 𝒱_d. 

Let us start by considering the case of an oriented rectifiable subset (X,T_X). A diffeomorphism ϕ∈Diff(ℝ^n) transports (X,T_X) as

    ϕ· (X,T_X) ≐ (ϕ(X), T_ϕ(X)),

where the transported orientation map writes

    T_ϕ(X)(y) ≐ d_ϕ^-1(y)ϕ· T_X(ϕ^-1(y))

the above term being well-defined from (<ref>). This suggests introducing the following pushforward action on 𝒱_d, which is defined for all μ∈𝒱_d and ϕ∈Diff(ℝ^n) by:

    (ϕ_#μ | ω) ≐∫_ℝ^d ×G^n_dω(ϕ(x),d_xϕ· T) J_T ϕ(x) d μ(x,T)

in which J_T ϕ(x) denotes the determinant of the Jacobian of ϕ along T (i.e. the change of d-volume induced by ϕ along T at x) which is given by 

    J_T ϕ(x) = ((d_xϕ(e_i) · d_xϕ(e_j))_i,j=1,…,d)

for (e_1,…,e_d) an orthonormal basis of T. One easily verifies that (ϕ,μ) ↦ϕ_#μ defines a group action which commutes with the action on oriented rectifiable sets, namely


 For any oriented rectifiable set (X,T_X) and diffeomorphism ϕ∈Diff(ℝ^n), ϕ_#μ_X = μ_ϕ(X).

This follows from the area formula for integrals over rectifiable sets, c.f. <cit.> Chapter 2. 

This pushforward action also extends the diffeomorphic transport of measures with densities on ^n. Indeed if μ=θ(x). ℒ^n with θ a measurable density function on ^n and ℒ^n the Lebesgue measure, we can extend μ to a n-varifold in 𝒱_n by taking a constant global orientation in G_n^n = {± 1}. Then, for any orientation-preserving diffeomorphism ϕ, (<ref>) writes in this case: ϕ_#μ = |Jϕ(x)| θ(x) . ℒ^n with Jϕ(x) is the full Jacobian determinant of ϕ, leading to the usual action on densities ϕ·θ(x) = |Jϕ(x)| θ(x).  



However, in contrast with past works on submanifold registration, this is not the only possible group action that could be considered on the space 𝒱_d. For instance, one can define another action by removing the above volume change term, taking instead

    (ϕ_*μ | ω) := ∫_ℝ^d ×G^n_dω(ϕ(x),d_xϕ· T) d μ(x,T).


This normalized action has the property of preserving the total mass of the varifold, i.e.,

    | ϕ_*μ | (ℝ^n) = | μ | (ℝ^n),  ∀μ∈𝒱_d  and ϕ∈Diff(ℝ^n).

Although this action is not consistent with the action on rectifiable sets as in Proposition <ref>, this model may be more adequate in applications to certain types of data in which mass preservation is natural. 

We refer the interested reader to <cit.> for a more in depth discussion on the properties (orbits, isotropy subgroups...) of these group actions in the simpler case of 1-varifolds. In the rest of the paper, we will restrict ourselves to the pushforward action model of (<ref>), although we expect the following derivations to adapt to other cases as well, which precise study is for now left as future work.




 §.§ The diffeomorphic registration problem



With the group action defined above, we are now ready to introduce the mathematical formulation of the diffeomorphic registration problem for general varifolds in 𝒱_d. As deformation model, we will rely on the Large Deformation Diffeomorphic Metric Mapping (LDDMM) setting mentioned in the introduction. 

Let us briefly sum up the basic construction of LDDMM, which details can be found in <cit.>. In this framework, deformations consist of diffeomorphisms generated by flowing time-dependent vector fields. Let V be a fixed RKHS of vector fields on ^n and L^2([0,1],V) be space of time dependent velocity fields v such that for all t∈ [0,1], v_t belongs to V. The flow map t ↦φ_t^v is defined for all t∈[0,1] by φ_0^v = id and the ODE φ̇^v_t = v_t ∘φ^v_t. If V is continuously embedded in C_0^1(^n,^n), one can show that for all t, φ_t^v is a C^1-diffeomorphism of ^n. Moreover, on the subgroup Diff_V ={φ_1^v  |  v ∈ L^2([0,1],V)} of Diff(^n), one can define the following right-invariant Riemannian metric:

    d_G_V(id,ϕ) = inf{∫_0^1 v_t_V^2 dt  |  φ_1^v = ϕ}



Let us now consider a source (or template) varifold μ_0 ∈𝒱_d as well as a target μ_tar∈𝒱_d. With the above deformation model and metric, registering μ_0 to μ_tar consists in finding a deformation ϕ that minimizes d_G_V(id,ϕ) with the constraint that ϕ_#μ_0 is close to μ_1 in the sense of a kernel metric ·_W^* defined in Section <ref>. This can be reformulated as the following optimal control problem:

    _v ∈ L^2([0,1],V){E(v) = 1/2∫_0^1v_t_V^2 dt + λμ(1) - μ_tar_W^*^2 }

with v being the control, E the total cost and the state equation is given by μ(t) ≐ (ϕ_t^v)_#μ_0 for the pushforward model. The first term in (<ref>) is the regularization term that constrains the regularity of the estimated deformation paths. The second term measures the similarity between the deformed varifold μ(1) and the target varifold μ_tar. λ is a weight parameter between the regularization and fidelity terms. Note that this is consistent with the generic inexact registration problem formulation in LDDMM that was proposed for objects like images, landmarks, submanifolds...   


The well-posedness of the optimal control problem (<ref>) holds under the following assumptions:


If V is continuously embedded in C_0^2(ℝ^n,ℝ^n), W is continuously embedded in C_0^1(ℝ^n ×G^n_d)  and supp(μ_0) ⊂ K, for some compact subset K of ℝ^n ×G^n_d, then there exists a global minimizer to the problem (<ref>).

The proof is similar to previous results of the same type on rectifiable currents and varifolds. We give it in Appendix for the sake of completeness.  

One can derive necessary and sufficient conditions on the kernels of W and V for the two embedding assumptions of Theorem <ref> to hold (see for instance Theorem 2.11 in <cit.>). In our context, in order to get W ↪ C_0^1(ℝ^n ×G^n_d) for instance, it is enough to assume that ρ and γ are C^2 functions such that all derivatives of ρ up to order 2 vanish as x → +∞. 


As an important note, the formulation of (<ref>) extends registration of submanifolds or rectifiable subsets in the sense that if μ_0 = μ_X_0 and μ_tar=μ_X_tar for two oriented d-rectifiable subsets of ^n then (<ref>) becomes equivalent, thanks to Proposition <ref>, to registering rectifiable subsets, i.e. to the problem 

    _v ∈ L^2([0,1],V){1/2∫_0^1v_t_V^2 dt + λμ_X(1) - μ_X_tar_W^*^2 }

with X(t) = φ_t^v · X_0, which is the setting of many past works as for instance <cit.>.




 §.§ General optimality conditions


A last important question we address in this section is the derivation of necessary optimality conditions for the solutions of (<ref>). In standard finite-dimensional optimal control problems, these are provided by the Pontryagin Maximum Principle (PMP) introduced originally in <cit.>. The approach generalizes, with a certain number of technicalities, to a broad class of infinite-dimensional shape matching problems, as developed in <cit.>. 

We follow the same setting as well as related works such as <cit.> by first rewriting the above problem as an optimal control problem on diffeomorphisms, i.e.

    _v ∈ L^2([0,1],V){1/2∫_0^1v_t_V^2 dt + g(φ_1^v)  |  s.t. φ̇_t^v = v_t ∘φ_t^v }

with g(φ_1^v) ≐λ(φ_1^v)_#μ_0 - μ_tar_W^*^2. The state variables are now given by the deformations φ_t^v which we view as elements of the Banach space ℬ≐id + C_0^1(^n,^n). Let us denote, for ϕ∈Diff(^n), ξ_ϕ : V → C_0^1(^n,^n) the mapping v ↦ v ∘ϕ. We then introduce the Hamiltonian functional H:  C_0^1(^n,^n)^* ×ℬ× V → defined by:

    H(p,ϕ,v) = (p|v∘ϕ) - 1/2v_V^2

where p is the costate variable which is a vector distribution of C_0^1(^n,^n)^* and (p|v∘ϕ) denotes the duality bracket in C_0^1(^n,^n)^*. With the assumptions of Theorem <ref>, it follows from the maximum principle shown in <cit.> that if (v_t,φ_t^v) is a global minimum of the optimal control problem, there exists a path of costates p ∈ H^1([0,1],C_0^1(^n,^n)^*) such that the following equations hold:

    {[                               φ̇_t^v = ∂_p H(p_t,φ_t^v,v_t);                                 ṗ_t = -∂_ϕ H(p_t,φ_t^v,v_t);                                    ∂_v H(p_t,ϕ_t^v,v_t) = 0 ].

with the end time boundary conditions p_1 = -∂_ϕ g(φ_1^v). From the last equation in (<ref>), we can attempt to deduce the form of the optimal v. Introducing the Riesz isometry operator K_V : V^*→ V and its inverse L_V=𝐊_V^-1: V → V^*, we get:

    ξ_φ_t^v^* p_t - L_V v_t = 0  ⇒ v_t = K_V ξ^*_φ_t^v p_t.


One additional consequence of (<ref>) is the following conservation of momentum again proved in <cit.>: for all u ∈ C_0^1(^n,^n) and t∈[0,1], 

    (p_t | dφ_t^v u) = (p_0 | u).

Note that (<ref>), (<ref>) and (<ref>) are generic to the LDDMM model and so far independent of the nature of the deformed objects and of the term g(φ_1^v) in the cost. This dependency is entirely encompassed by the boundary condition p_1 = -∂_ϕ g(φ_1^v) which we may describe a little more precisely based on the following:



 The end-time momentum p_1 is a vector distribution in C_0^1(^n,^n)^* of the form
 
    (p_1|u) =    ∫_^nα(x) · u(x)  d|μ_0|(x) 
     
         +   ∫_^n ×G^n_dβ(x,T)  du|_T(x)  dμ_0(x,T) 
    
         +   ∫_^n ×G^n_dγ(x,T) div_T u(x)  dμ_0(x,T)

 where α: ^n →^n, β: ^n ×G^n_d → (^n× d)^* and γ: ^n ×G^n_d → are continuous fields and for all T ∈G^n_d, div_T u and du|_T denote the divergence and differential of u restricted to T. 

A condensed proof of this proposition can be found in the Appendix, although we have left aside the technical derivations related to differential calculus on the Grassmannian (this will be discussed further in Section <ref> in the discrete setting). This result extends in a way first variation formulas for varifolds proved in <cit.> which considered variations of rectifiable varifolds resulting from variations of the underlying rectifiable sets. This corresponds to the special case in which μ_0 = μ_X_0. In that case, one can show, after some derivations, that the above expression of p_1 can be rewritten in the form of a vector distribution u ↦∫_φ_1^v(X_0) u(x) · h(x) d ℋ^d in C_0^0(^n,^n)^* with vectors h(x) normal to φ_1^v(X_0) at each x. In our more general situation, this is however not possible and p_1 is a priori a distribution that involves first order derivatives of the test function u.

Now, the conservation law of (<ref>) gives that for all t ∈ [0,1], 

    (p_t | dφ_t^v u) = (p_1 | dφ_1^v u) = (p_0 | u).

Using the expression of p_1 in Proposition <ref>, and grouping all 0-th and 1-st order terms in the resulting expressions, we may write p_t in the general form:

    (p_t | u)    = ∫_^n ×G_d^nα_t(x,T) · u(x)  dμ_0(x,T) + ∫_^n ×G_d^n B_t(x,T) d_x u|_T dμ_0(x,T)

where α_t: ^n ×G_d^n →^n and B_t: ^n ×G_d^n → (^n× d)^* are continuous fields, with α_1(x,T)=α(x) and B_1(x,T) du|_T(x) = β(x,T) du|_T(x) + γ(x,T) div_T u(x). Furthermore, optimal vector fields satisfy v_t = K_V ξ^*_φ_t^v p_t and we have

    (ξ^*_φ_t^v p_t | u)    = (p_t | u∘φ_t^v) 
       = ∫_^n ×G_d^nα_t(x,T) · u(φ_t^v(x))  dμ_0(x,T) + ∫_^n ×G_d^n B_t(x,T) d_φ_t^v(x)u|_d_xφ_t^v · T dμ_0(x,T).

Denoting K_V: ^n ×^n →^n× n the reproducing kernel of V, the reproducing kernel property implies that for all u∈ V and x,h ∈^n, u(x) · h = ⟨ K_V(x,·) h, u ⟩_V. Moreover, the similar property on the kernel first order derivatives <cit.> gives that for any h,h' ∈^n, 

    d_x u (h) · h' = ⟨∂_1 K_V(x,·)(h) · h', u ⟩_V.

Then, we rewrite the linear maps B_t as B_t(x,T) H = ∑_i=1^d b_t,i(x,T) · H_i for any H=(H_1,…,H_d) ∈^n× d and where b_i(x,T) ∈^n are the component vector fields of B_t. By the above and the linearity of K_V, we obtain the following general expression for optimal vector fields   

    v_t    = ∫_^n ×G_d^n K_V(φ_t^v(x),·) α_t(x,T)  dμ_0(x,T)    
       + ∫_^n ×G_d^n(∑_i=1^d ∂_1 K_V(φ_t^v(x),·)(d_xφ_t^v(t_i))· b_t,i(x,T)) dμ_0(x,T).

In contrast with LDDMM registration of submanifolds or point clouds, the expression of optimal deformation fields involves in general both the kernel function and its first order derivatives. We do not explicit the vector fields α and b_i at this point, it will be specified later in the discrete setting, see Section <ref>. 




§ APPROXIMATIONS BY DISCRETE VARIFOLDS



The previous derivations were so far conducted for completely general measures in the space 𝒱_d which include objects of widely different natures. In the perspective of implementing numerically the above approach, which is the subject of Section <ref>, we first need to build an adequate discretization framework in 𝒱_d with approximation guarantees, and even more importantly investigate the consistency of the discretized registration problems (Theorem <ref>), which is the main result of this section. 



 §.§ Discrete approximations


In what follows, we will consider the specific class of varifolds which can be written as finite combinations of Dirac masses:

    μ = ∑_i=1^N r_i δ_(x_i,T_i), r_i ∈ℝ_+, x_i ∈ℝ^n,  T_i ∈G̃^n_d.

for some N≥ 1. Throughout this paper, varifolds of this form will be called discrete varifolds. It is quite natural to consider this type of varifolds for the purpose of representing discrete shapes, which has been exploited in previous works on piecewise linear curves and surfaces. For example, if X = ⋃_i=1^N X_i is a triangulated surface, with X_i being the mesh triangles with specified orientations, one can write μ_X = ∑_i=1^N μ_X_i and for each i ∈{1,⋯,N} approximate μ_X_i by r_i δ_(x_i,T_i), where x_i is the center of X_i, T_i the oriented plane containing X_i and r_i = ℋ^d(X_i). This leads to the approximation μ_X := ∑_i=1^N r_i δ_(x_i,d_i). As proved in <cit.>, this approximation provides an acceptable error bound for d_W^*:

    d_W^*(μ_X,μ_X) ≤ Cte  ℋ^d(X) max_idiam(X_i).

The main interest of such discrete varifold approximations is that the expression of the metric (<ref>) becomes particularly simple to compute numerically. Indeed, given two discrete varifolds μ = ∑_i=1^N r_i δ_(x_i,S_i) and μ' = ∑_j=1^M r'_j δ_(x'_j,T'_j), we have as a particular case of (<ref>):

    ⟨μ,μ' ⟩_W^* =  ∑_i=1^N∑_j=1^M  r_i r'_j ρ(|x_i-x'_j|^2) γ(⟨ T_i,T'_j ⟩).


The above approximation scheme only applies to the case of piecewise linear shapes given by meshes such as polygonal curves or triangulated surfaces. In the more general context of this work, a key issue is to construct similar discrete varifold approximations for more general and less structured objects. Specifically, given a varifold μ with finite total weight, can it be approximated by discrete varifolds and will approximations converge as N→ +∞? This is the problem known as quantization, which has been studied intensively in the case of probability measures over Euclidean spaces <cit.> or manifolds <cit.>, under specific regularity assumptions on those measures. In the situation of varifolds, an interesting recent work on this question is <cit.>. The authors prove that any rectifiable varifold with finite mass can be approximated by a sequence of discrete varifolds for the bounded Lipschitz distance and propose a numerical approach to approximate mean curvature measures based on discrete varifolds. 

In this section, we first wish to extend approximation results to general oriented varifolds of finite mass for both d_BL and d_W^* metrics.
 

Let
 
    𝒱^N_d := {∑_i=1^N r_i δ_(x_i,T_i) | r_i ∈ℝ_+, x_i ∈ℝ^n,  T_i ∈G^n_d }

 be the (non-convex) cone of discrete varifolds with at most N Diracs. For any oriented varifold μ∈𝒱_d with |μ |(ℝ^n) < ∞, there exists a sequence μ_N ∈𝒱_d^Nsuch that lim_N →∞ d_BL(μ_N,μ) =0. Moreover, if μ has compact support, then we can assume that for all N, supp(μ_N) ⊂ K for some compact set K ⊂ℝ^n ×G^n_d and 
 
    d_BL(μ_N,μ) < C/N^1/(n+d(n-d)),

 where C is a constant that only depends on n, d and supp(μ).



 We first tackle the case of compactly supported μ. Without loss of generality, we may also assume that μ is a probability measure. Let D = n+ d(n-d) and B ⊂ℝ^n be a closed ball that contains supp|μ|. For brevity, we write M ≐ B ×G^n_d. Since we can view M as a compact D-dimensional submanifold of ^n ×Λ^d(ℝ^n) (using Plücker embedding), M is also regular of dimension D (cf. <cit.>), i.e., 0<ℋ^D(M)< ∞ and there exist c,r_0>0, such that
  
    1/cr^D ≤ℋ^D  M(B_r(a)) ≤ c r^D, ∀ a ∈ M,  r ∈ (0,r_0).

  Given ε∈ (0,5 r_0), by the 5-Times Covering Lemma (cf. <cit.>), these exists a subset ℐ⊂ M, such that M ⊂∪_x ∈ℐ B_ε(x) and B_ε/5(x) ∩ B_ε/5(y) = ∅ for all x ≠ y ∈ℐ. Therefore,
  
    ℋ^D(M) ≥∑_x ∈ℐℋ^D(M ∩ B_ε/5(x)) 
          ≥|ℐ| ε^D/c 5^D.

 We can thus obtain a partition {A_i}_i=1,⋯,|ℐ| of M from the the collection {B_ε(x) ∩ M}_x∈ℐ which satisfies sup_i diam(A_i) < ε and 
 
    |ℐ| ≤c 5^D ℋ^D(M)/ε^D.

Let r_i = μ(A_i) and (x_i,T_i) ∈ A_i and define ν=∑_i=1^|ℐ| r_i δ_(x_i,T_i). For any φ∈Lip_1(ℝ^n ×G^n_d), with φ_∞≤ 1, we have


    | ∫_ℝ^n ×G_d^nφ(x,T)d ν - ∫_ℝ^n ×G_d^nφ(x,T) dμ| 
       = | ∑_i=1^|ℐ|( μ(A_i) φ(x_i,T_i) - ∫_A_iφ(x,T) dμ) | 
       ≤∑_i=1^|ℐ|∫_A_i| φ(x_i,T_i) - φ(x,T)| dμ
       <∑_i=1^|ℐ|εμ(A_i) = ε.

Taking the supremum over all ω∈Lip_1(ℝ^n ×G^n_d) with φ_∞≤ 1, we obtain d_BL(μ,ν)< ε. Then for each N ∈ℕ, we can choose ε_N = 5(C ℋ^D(M)/N)^1/D and we obtain μ_N ∈𝒱_d^N such that

    d_BL(μ,μ_N)< 5 C^1/D (ℋ^D(M))^1/D/N^1/D

and in particular lim_N → +∞ d_BL(μ,μ_N) = 0.  

Suppose now that supp(μ) is not compact: we show that for any ε>0, there exists a discrete varifold ν such that d_BL(μ,ν)< ε. Choose a compact set K ⊂ℝ^n ×G^n_d such that μ(ℝ^n ×G^n_d ∖ K)< ε/2. From the previous case, we can find a discrete varifold ν such that d_BL(μ K,ν )<ε/2, and hence d_BL(μ,ν )<ε.

Note that the proposition clearly holds for non-oriented varifolds as well. Another direct consequence, thanks to proposition <ref> and remark <ref>, is the following corresponding statement for d_W^*:


With the assumptions from proposition <ref>, one also has lim_N →∞ d_W^*(μ_N, μ) =0. If in addition W ↪ C_0^1(^n×G̃_d^n), an equivalent upper bound as in Theorem <ref> holds for d_W^*(μ,μ_N).    

We should point out that the asymptotic convergence rate given by the previous upper bound is rather slow, especially as the dimensions d and n grow. This is however under very mild assumptions on the varifold μ. We expect much better convergence properties for certain specific classes of varifolds, for instance assuming Alfors regularity as in <cit.>, although we leave such questions for future investigation.




 §.§ Optimal approximating sequence

In addition to the asymptotic approximation results of the previous section, we now want to construct such sequences of discrete approximating varifolds. Given any μ∈𝒱_d and N ∈, a natural idea is to look for the optimal discrete varifold in 𝒱_d^N that approximates μ in terms of the metric d_W^*. Due to the intricate structure of the set 𝒱_d^N (infinite-dimensional non-convex cone), this is far from a straightforward problem. Several different approaches in some simpler contexts have been proposed to circumvent this issue, which we briefly recap. One possibility is to restrict to finite-dimensional vector spaces of 𝒱_d (e.g. generated by finite sets of Diracs). Works such as <cit.> for instance, which are focused on the model of currents, consider dictionaries of Diracs defined on a predefined grid of point positions in ^n. Then the problem can be recast as the one of finding sparse approximations of μ in such a dictionary. It remains a NP hard problem but solutions can be approached either through greedy algorithms like orthogonal matching pursuit as proposed in <cit.> or by considering the L^1 relaxation formulation leading to a standard convex LASSO program. Such ideas apply well to the specific situation of currents mainly as a result of the inherent linearity of this model: indeed, at any iteration of a matching pursuit procedure, once the optimal position of a Dirac is found, the corresponding direction vector and weight are explicitly determined. This allows to limit the search over grid of points in the spatial domain only. Unfortunately, for the general oriented varifold metrics we consider in this paper, such a property no longer holds and, as a result, these methods would involve very large dictionaries defined on grids on the product ^n ×G_d(^n). Such an increase in dimension makes the approach numerically impractical as soon as n≥ 3 and d≥2. Another downside is that the use of finite dictionaries and greedy algorithms like matching pursuit is not guaranteed to give an optimal approximation of varifolds for a given number N of Diracs. The approach we develop in this section consists instead in directly tackling the non-convex problem of computing the projection onto 𝒱_d^N for the class of kernel metrics d_W^*. It shares some connection with the recent work of <cit.> that considers a related problem for standard measures defined on the torus ^n / ℤ^n. 

Fix a varifold μ_* ∈𝒱_d. For any N ∈ℕ,  N≥ 1, we seek μ_N ∈𝒱_d^N that is closest to μ_* for d_W^*, namely

    μ_N = *argmin_μ∈𝒱^N_dμ - μ_* _W^*


By construction, if |μ|(ℝ^n)< ∞ then Corollary <ref> will imply that (μ_N) converges to μ in the metric d_W^*. We only need to ensure that such a projection is well defined, which is the object of the following proposition:

 
Suppose all assumptions in proposition <ref> and remark <ref> hold. We further assume  that the functions ρ and γ defining the kernels are non-negative.
Then for any μ∈𝒱_d and N ∈ℕ, there exists μ_N∈𝒱^N_d such that μ_N = *argmin_ν∈𝒱^N_dμ - μ_* _W^* 



  Let μ_m = ∑_i=1^N r_i^m δ_(x_i^m,T_i^m) be a minimizing sequence, K_W: W^* ↦ W be the dual operator and f := K_W(μ). Without loss of generality, we may assume that there is a N_1 ≤ N such that sup_1≤ i ≤ N_1 |x_i^m| remains bounded and inf_ M+1≤ i ≤ N |x_i^m| tends to ∞ as m →∞.

Observe that sup_1 ≤ i ≤ N{r_i^m} must be bounded. If it's not bounded, then from the assumptions that ρ, γ≥ 0, we obtain 


    μ_m -μ_* _W^*   ≥μ_m  _W^* -μ_*  _W^* 
                             = √(∑_i,j=1^N r_i^m r_j^m ρ(|x_i^m - x_j^m|^2) γ(⟨ T_i^m,T_j^m ⟩)) - μ_* _W^*
       ≥√(∑_i=1^N (r_i^m)^2 ρ(0) γ(1)) - μ_* _W^*→∞

as m →∞, which is absurd.

Since r_i^m, T_i^m, f(x_i^m,T_i^m) and 

    A_m := (ρ(|x_i^m -x_j^m|^2) γ(⟨ T_i^m, T_j^m ⟩))_1 ≤ i,j ≤ N

 are all bounded sequences of m, we may replace them by convergent subsequences, thus we could assume that 


    lim_m →∞ r_i^m = r_i,   lim_m →∞ T_i^m = T_i,  lim_m →∞ f(x_i^m,T_i^m) = f_i,   lim_m →∞ A_m = A.

Since ρ∈ C_0(ℝ), the matrix A must has the following form:


    A = ([ B_1   0;   0 B_2 ]),

where B_1 and B_2 are N_1-by-N_1 and N-N_1-by-N-N_1 semi-positive definite matrices. Combining this with the assumption f ∈ C_0(ℝ^n ×G^n_d), we obtain


    lim_m →∞μ_m -μ_* _W^*^2 
        = r'^T B_1 r' + r”^T B_2 r” - 2 f'^T  r' + μ_* _W^*^2,

where r' = (r_1,⋯,r_N_1), r” = (r_N_1+1,⋯,r_N), and f' = (f_1,⋯,f_N_1). Since sup_1 ≤ i ≤ N_1 |x_i^m| is bounded we can assume that  lim_m →∞ x_i^m = x_i,  1 ≤ i ≤ N_1. Let μ := ∑_i=1^N_1 r_i δ_(x_i,u_i), then 

    μ -μ_* _W^*^2    = r'^T B_1 r' - 2 f'^T  r' + μ_* _W^*^2 
         ≤lim_m →∞μ_m -μ_* _W^*^2.

Hence μ is a minimizer. 

However, in general this projection is not unique. We also point out that the existence is a priori not guaranteed if kernels ρ and γ take negative values. It is so far an open question to determine to what extent one could generalize the result of Proposition <ref>, one particular but important case being the one of current metrics obtained for γ(t)=t which is not covered by our result.

As written in the proof of Proposition <ref>, (<ref>) is equivalent to the optimization problem: 

    (r_i,x_i,T_i) = *argmin_(w_i,y_i,S_i)∑_i=1^N w_i δ_(y_i,S_i) - μ_* _W^*^2

Any solution must satisfy first order optimality conditions obtained by differentiating μ_N-μ_W^* with respect to the (r_k,x_k,T_k). In particular, we have

    0    = ∂μ_N-μ_* _W^*^2/∂ r_k
       = 2 ( ∑_i=1^N r_i ρ(|x_i-x_k|^2) γ(⟨ T_i,T_k ⟩)  - ∫_ℝ^n ×G^n_dρ(|x_k-x|^2) γ(⟨ T_k,T ⟩) d μ_*(x,T) ).

which gives ⟨μ_N-μ_*,μ_N ⟩_W^* = 0. It shows that for any N ∈, μ_N_W^*≤μ_*_W^*.




 §.§ Γ-convergence of registration functionals

Ultimately, our purpose is to use the previous approximating discrete varifolds μ_N to approximate the diffeomorphic registration problem (<ref>). The natural question that arises is whether replacing the source varifold μ_0 by its projections μ_N in (<ref>) still leads to reasonable approximations (at least asymptotically) of optimal deformation fields for the original problem. In this section, we address this by showing a Γ-convergence property for these variational problems. We point out that our setting and the following proof differ quite a bit from previous results of the same type that were dealing with the specific case of surface triangulations such as <cit.>.      
   
To obtain such convergence results for solutions of variational problems, one usually requires the approximating sequence to possess certain nice properties. Specifically, assuming μ∈𝒱_d with compact support and finite mass and {μ_N}⊂𝒱^N_d such that lim_N →∞μ_N - μ_W^* = 0, we will need that ⋃_Nsupp(μ_N) ⊂ K for some compact set K ⊂ℝ^n or that sup_N |μ_N|(ℝ^n) < ∞. Unfortunately, this does not hold in general since convergence in d_W^* does not allow to control the support nor the total mass of the sequence μ_N.

Yet, provided that ⋃_Nsupp(μ_N) ⊂ K, we can actually retrieve the boundedness of the total mass. We assume in what follows that the kernels are such that ρ(0)>0 and γ(1)>0. 
   
      Let {μ_N} be a sequence of discrete varifolds with finite mass such that there exists a compact K ⊂^n with supp(|μ_N|) ⊂ K for all N. We assume that {μ_N_W^*} is bounded. Then {|μ_N|(ℝ^n)} is bounded. 
   
   

We prove it by contradiction. Assume that (|μ_N|(ℝ^n))_N ≥ 1 is unbounded. Then, up to extracting a subsequence, we can assume that |μ_N|(ℝ^n) → +∞. Let's write μ_N = ∑_i=1^p_N r_i,Nδ_(x_i,N,T_i,N). Thus |μ_N|(ℝ^n) = ∑_i=1^p_N r_i,N→ +∞. Since ρ and γ are continuous and ρ(0),γ(1)>0, we can find compact subsets A ⊂ K and B ⊂G̃^n_d with diameters small enough, so that: inf_x,y ∈ Aρ(|x-y|^2)>m>0, inf_u,v ∈ Bγ(⟨ u, v ⟩) > m' >0 and lim_N →∞∑_i∈ℐ_N r_i,N = ∞, 
where ℐ_N := {i: (x_i,N,u_i,N) ∈ A × B}. It follows that, as N →∞,

    μ_N_W^*^2    = ∑_i=1^p_N∑_j=1^p_N r_i,N r_j,Nρ(|x_i,N-x_j,N|^2) γ(⟨ T_i,N,T_j,N⟩) 
       ≥ mm' ∑_i,j ∈ℐ_Nr_i,N r_j,N
       = mm' ( ∑_i ∈ℐ_Nr_i,N)^2 →∞,

which is a contradiction. 


Lemma <ref> suggests that one should enforce the uniform compactness of the supports of the μ_N. To do so in the context of the projection approach of the previous sections, we consider solving the optimization problem (<ref>) with the additional constraint that the support of μ_N stays in a compact set containing supp(|μ|). We still have to verify the convergence of the resulting sequence:  


   Let μ_0 be a varifold with finite mass and K be a compact set in ℝ^n which contains supp(|μ_0|). Construct the approximating sequence of μ_0 by solving the following constrained optimization problem:
   
    μ_K,N = *argmin_ν∈𝒱^N_dν - μ_W^*
       subject  to  supp(|ν|) ⊂ K.

   Then μ_K,N converges to μ_0 in d_W^* and, if the kernel k is C_0-universal, it also converges in d_BL. 


   Thanks to Theorem <ref> and Lemma <ref>, we immediately get that μ_K,N - μ_0_W^*→ 0 as N →∞ and sup_N(|μ_K,N|(ℝ^n))< ∞. Moreover, if k is C_0-universal, then by Proposition <ref> it implies that μ_K,N∗⇀μ_0. Since ⋃_N supp(|μ_K,N|) ⊂ K and sup_N(|μ_K,N|(ℝ^n))< ∞, weak-* convergence implies that μ_K,N converges to μ_0 in d_BL by Proposition <ref>. 


We are now able to state the main result of this section. We assume that the source/template varifold μ_0 is compactly supported and we fix K is a compact subset of ^n that contains supp(|μ_0|). Then for any N ∈,  N≥ 1, μ_K,N is defined as in Proposition <ref> and we introduce the following energy functionals E_N :  L^2([0,1],V) →_+:

    E_N(v) ≐1/2∫_0^1v_t_V^2 dt + λμ_K,N(1) - μ_tar_W^*^2 
       subj  to  {[ ∂_t φ^v_t = v_t ∘φ^v_t,  φ^v_0 = id;          μ_K,N(t) = ( φ^v_t)_#μ_K,N ].

which are the equivalent to the energy E of the original problem (<ref>) but replacing the template varifold μ_0 by its approximations μ_K,N. 


      With the above notations, we assume that the reproducing kernel k of W is C_0-universal and satisfies all the conditions of Proposition <ref>. We also assume the continuous embedding V ↪ C_0^2(ℝ^n,ℝ^n). Then, the sequence of functionals E_N Γ-converges to E for the weak topology on L^2([0,1],V). Consequently, if v_N is a global minimizer of E_N for each N≥ 1, then (v^N) is bounded in L^2([0, 1],V) and every cluster point for the weak topology of L^2([0,1],V) is a global minimum of E.  






We first show that whenever v^N converges to v̅ weakly in L^2([0,1],V), we have

    E(v̅) ≤lim inf_N →∞ E_N(v^N).

Since v ↦∫_0^1 v_V^2 dt is lower semicontinuous with respect to the weak topology, we only need to prove the following,

    lim_N →∞(φ_1^v^N)_#μ_K,N - φ_1^v̅·μ_0 _W^* = 0.







For all ω∈ W with ω_W≤ 1, we have


    |((φ_1^v^N)_#μ_K,N - (φ_1^v^N)_#μ_0 | ω) | 
       = |∫_K ×G̃^n_d J_S φ_1^v^N(x) 
     ω (φ_1^v^N(x),d_x φ_1^v^N· S ) d (μ_K,N - μ_0) |
       ≤ C_1  ∫_K ×G̃^n_dsup_N ≥ 1 J_S φ_1^v^N d (μ_K,N - μ_0) 
       ≤ C_1 ∫_ℝ^n ×G̃^n_d g(x,T) d (μ_K,N - μ_0),

where g ∈ C_c(ℝ^n ×G̃^n_d) and sup_N ≥ 1 J_S φ_1^v^N≤ g(x,S), for all (x,S) ∈ K ×G̃^n_d. Similar to the computation done in the proof of Theorem <ref>, we see that


    |((φ_1^v^N)_#μ_0 - (φ_1^v̅)_#μ_0 | ω) | 
    ≤ C_2  (φ_1^v^N - φ_1^v̅)|_K _1,∞.

Taking supremum over all ω∈ W with ω_W ≤ 1, we obtain the following inequality,

    (φ_1^v^N)_#μ_K,N - (φ_1^v̅)_#μ_0 _W^*≤  C_1  (μ_K,N-μ_0| g) 
    +C_2  (φ_1^v^N - φ_1^v̅)|_K _1,∞.

From Proposition <ref>, μ_K,N converges to μ_0 in the narrow topology. Hence the right hand side in the equation above converges to 0 as N →∞. This proves (<ref>). 

Second, we need to show that for each v̅∈ L^2([0,1],V), there exists a sequence v^N converging to v̅ weakly such that

    E(v̅) ≥lim sup_N →∞ E_N(v^N).

In fact, it suffices here to take v^N to be the constant sequence v^N = v̅ since, by a similar argument to the proof of (<ref>), it leads to  
 
    lim_N →∞(φ_1^v̅)_#μ^N - (φ_1^v̅)_#·μ_W^* = 0

and thus implies that
  
    lim sup_N →∞ E_N(v^N) = lim_N →∞ E_N(v̅) = E(v̅).



Note that we stated the result of Theorem <ref> in the situation where only the source varifold μ_0 is approximated by the projection approach that we presented in the previous sections but one can easily extend it to the scenario in which both source and target are replaced by discrete approximating sequences, the conclusion being the same in that case.  




§ NUMERICAL CONSIDERATIONS


Having introduced a variational formulation for the varifold registration problem together with an approach for projecting onto the space of discrete varifolds with fixed number of Diracs, we now turn more specifically to the numerical implementation of methods for solving those problems. The first hurdle, which we start by addressing in Section <ref>, is to define an adequate framework for representing and computing with elements of the oriented Grassmannian.



 §.§ Frame representation for metric computation and quantization


In order to come up with a computationally effective representation of G_d(^n) and by extension of discrete oriented varifolds, we consider a slightly different setting than the Plücker embedding idea of Remark <ref>, primarily because the dimension of the embedding vector space Λ^d(^n) may become prohibitively large in practice. We may instead choose to represent an element T ∈G_d(^n) by an oriented frame (u^(1),…,u^(d)) ∈^n × d of independent vectors for which T=Span(u^(1),…,u^(d)). Such a representation is of course not unique since elements of G_d(^n) are equivalence classes of oriented frames but we leave to the next section the more thorough analysis of the additional invariances that this representation will imply.

We will in fact go one step further by also incorporating the weight of Dirac varifolds in this frame representation itself, which is done as follows. Let μ be a discrete varifold of the form μ = ∑_i=1^N r_i δ_(x_i,T_i). For each i, we consider a frame {u_i^(1),⋯,u_i^(d)} such that

    T_i = u_i^(1)∧⋯∧ u_i^(d)/|u_i^(1)∧⋯∧ u_i^(d)| and  r_i = |u_i^(1)∧⋯∧ u_i^(d)|.

In other words, the oriented space spanned by the frame {u_i^(1),⋯,u_i^(d)} corresponds to T_i while its d-volume matches the weight r_i. Given such a choice of frame for each i, we can then identify μ with the (non-unique) state variable q = (x_i,u^(1)_i,⋯,u^(d)_i)_i=1,⋯,N in the vector space ℝ^Nn(d+1). Conversely, such a frame q with (u^(1)_i,⋯,u^(d)_i) a matrix of rank k for all i, corresponds to the (unique) discrete oriented varifold defined by the relations of (<ref>); we will denote it by μ^q in what follows.   

In this representation, the kernel metrics for discrete varifolds expressed in (<ref>) can be explicitly written as

    ⟨μ, μ' ⟩_W^*^2 =∑_i=1^N∑_j=1^M  r_i r'_j ρ(|x_i-x'_j|^2) γ(1/r_i r'_j(u_i^(k)·u'_j^(l))_k,l)

where r_i = |u_i^(1)∧⋯∧ u_i^(d)| = √((u_i^(k)· u_i^(l))_k,l). Note that this expression does not depend on the choice of frames that satisfy the conditions of (<ref>) for μ (and similarly for μ'). In the case where μ' is a more general non-discrete varifold in 𝒱_d, the computation of ⟨μ, μ' ⟩_W^*^2 involves integrals over ^n ×G_d(^n) of the kernel functions, which requires introducing specific quadrature schemes for approximating them. We do not address those issues in more details in this work as it needs particular discussion depending on the nature, regularity and dimension of the varifolds under consideration. Provided such adequate quadrature schemes have been defined, the W^* metric then formally reduces to an expression equivalent to (<ref>) in which the x_j',u'_j and r'_j are now the quadrature nodes and associated weights of the scheme.

In this setting, the solution to the projection problem (<ref>) can be computed by an iterative descent strategy on the vector q=(x_i,u^(1)_i,⋯,u^(d)_i)_i=1,⋯,N. The gradient of q ↦μ^q-μ_*_W^*^2 can be computed by direct differentiation of expressions like (<ref>) with respect to the x_i and u_i^(l). In practice, computations of varifold kernel metrics for different classes of kernels and gradients of the metrics can be conveniently implemented with automatic differentiation pipelines. In our MATLAB implementation, we make use of the recent KeOps library <cit.> which allows to generate CUDA functions for the low-level kernel sum evaluations and their automatic differentiation. The optimization itself is done using a limited memory BFGS algorithm from the HANSO library <cit.> which we typically initialize by taking a random subset of N Diracs composing the varifold μ_*. Note that one of the main downside of this projection algorithm, in contrast with the previously mentioned approach of fixing a dictionary and solving a convex sparse decomposition problem, is that we can provide no general guarantees of convergence to a global minimum of (<ref>). Results of this algorithm are discussed below in Section <ref>.




 §.§ Discrete registration model


This frame representation also provides a convenient setting to express the diffeomorphism action and registration problem on discrete varifolds. Indeed, let φ be a diffeomorphism of ^n and μ∈𝒱_d^N, the pushforward action φ_#μ in (<ref>) is equivalent to the following action in the frame model:

    φ_# q := (φ(x_i),d_x φ (u^(1)_i),⋯,d_x φ (u^(d)_i))_i=1,⋯,N .

Now, this allows us to rewrite the former infinite-dimensional optimal control problem by considering instead the finite-dimensional state variable q ∈^Nn(d+1). In the next paragraphs, we give a direct derivation of the optimality conditions in this discrete setting, in order to arrive at simpler and more explicit equations than the general abstract derivations presented in Section <ref>. Note that the resulting Hamiltonian equations we obtain are eventually very similar to the ones appearing in the 1st-order jets model studied in <cit.>, although there are a few notable differences due to the specific extra invariances attached to the varifold framework (c.f.  <cit.> for a more detailed discussion in the d=1 case).      

Following once again the Pontryagin maximum principle approach, the Hamiltonian for this discrete representation is given by:

    H(q,p,v) ≐∑_i=1^N [ p_i^x· v(x_i) + ∑_k=1^d  p_i^u_k· d_x_i v(u_i^(k)) ]  - 1/2v^2_V

with p^x,p^u_k∈^n denoting respectively the costates for the position x and frame vector u^(k) variables.
The PMP then shows that optimal trajectories of the registration problem are governed by the dynamical system:

    {[                                                      ẋ_i = v_t(x_i);                                         u̇_i^(k) = d_x_i v(u_i^(k)); ṗ_i^x =  - d_x_iv^T p_i^x - ∑_k=1^d d_x_i^(2)v(·,u_i^(k))^T p_i^u_k;                                        ṗ_i^u_k = - d_x_iv^T p_i^u_k;                                                                     ].

while optimal vector fields v satisfy

    v_t(·)    = ∑_i=1^N ( K(x_i(t),·)p^x_i(t)  
        + ∑_k=1^d ∂_1 K(x_i(t),·)(u_i^(k)(t)) · p_i^u_k).

Plugging the above expression of v with respect to (p,q) in the Hamiltonian (<ref>) gives the reduced Hamiltonian H_r(p,q) ≐ H(p,q,v) which writes:

    H_r(p,q) = 1/2∑_i,j=1^N[    p_i^x · K(x_i,x_j) p_j^x  + p_i^x ·∑_k=1^d ∂_1 K(x_j,x_i)(u_j^(k)) · p_j^u_k
       + ∑_k=1^d  p_i^u_k·∂_2 K(x_j,x_i)(u_j^(k)) p_j^x + ∑_k,l=1^d p_i^u_k·∂^2_1,2 K(x_j,x_i)(u_j^(l),u_i^(k)) p_j^u_l]

and (<ref>) becomes a coupled system in the variables q and p called the reduced Hamiltonian equations. Consequently, the set of optimal paths is entirely determined by the initial values (q(0),p(0)) and the value of the reduced Hamiltonian H_r(p(t),q(t)) = 1/2v_t_V^2 is conserved along an optimal trajectory.  

There are in addition several other conserved quantities in such a system as evidenced by the following lemma: 


For any i=1,…,N, the matrix 

    D^i(t)  ≐( ⟨ u_i^(k)(t), p_i^u_ℓ(t) ⟩)_1 ≤ k, ℓ≤ d,

is constant in time.


Using the Hamiltonian equations written above, we have for all k,l=1,…,d

    d/dt(D^i(t) )_k,ℓ = ⟨ d_x_iv(u_i^(k)(t)), p_i^u_ℓ⟩ - ⟨ u_i^(k), d_x_iv^T p_i^u_ℓ⟩ = 0.

Hence D^i(t) is a constant matrix. 
Note that, at this point, all those equations are fundamentally modelling the deformation of the frames {x_i,(u_i^(k))} but are not yet taking into account the invariances that result from the representation of the discrete oriented varifolds as oriented frames. Those extra invariances can be derived from the boundary conditions of the PMP:

    p(1) = - ∂_q g(q)|_q=q(1),  with g(q) = λμ^q - μ_tar_W^*^2.

As a clear consequence of (<ref>), μ^q and thus g(q) are independent of the choices of the frame vectors (u_i^(k))_k=1,…,d that span the same oriented vector spaces T_i with the same d-volumes r_i. This in turn leads to a set of conditions satisfied by the different components of the final costate p(1) and, with Lemma <ref>, of the full path p(t). These are summed up by the following result:



Let (q(t),p(t)) be optimal trajectory, then for all i, the matrices D^i(t) as defined above are constant scalar matrices. In particular, we have p_i^u_k(t) ⊥Span({u_i^(ℓ)(t)}_ℓ≠ k) for all t ∈ [0,1], i=1,…,N and k=1,…,d.

This result, which proof can be found in Appendix, is particularly interesting from a computational point of view as it allows to partly alleviate the redundancy introduced by the frame representation of Grassmannians. Indeed, we see that the costates p(t) actually lie in affine subspaces of ^Nn(d+1) of lower dimensions N(n + d(n-d)+1), which is precisely the dimension of the 'true' state space (^n ×G_d(^n)×)^N. 



 §.§ Registration algorithm


Based on the optimality equations of the previous section, we can now easily design an algorithm to solve the discrete registration problem. As mentioned earlier, optimal trajectories are completely determined, through the Hamiltonian equations (<ref>) and (<ref>), by the initial conditions q(0)=q_0, which is known, and p(0). One of the standard class of methods in optimal control, known as shooting methods, consist in directly optimizing the cost function over p(0), which has been the approach of choice in many past works on shape registration such as <cit.>. We adopt a similar strategy for our particular problem.

The main issue is to compute the gradient of the total energy E with respect to the initial costate p(0). The regularization term being equal to H_r(p(0),q(0)) thanks to the conservation of the reduced Hamiltonian, its gradient can be obtained by direct differentiation of (<ref>). The fidelity term g(q(1)) on the other hand depends indirectly on the initial costate p(0) via the integration of the forward reduced Hamiltonian equations. As standard for this type of optimal control problems, c.f. <cit.> or <cit.>, the gradient of g(q(1)) with respect to p(0) can be computed by flowing backward in time the adjoint Hamiltonian system 

    Z(t) = -dF(q(t),p(t))^T Z(t)

where F(q,p)=(∂_p H_r(p,q), -∂_q H_r(p,q)), Z(t)=(q̃(t),p̃(t)) ∈ℝ^n × (ℝ^n)^d the adjoint variables of the system, together with the end-time conditions q̃(1) = -∂_q g(q)|_q=q(1) and p̃(1)=0. Although being a linear system of ODEs, the adjoint equations can be tedious to derive and implement, in particular given the rather intricate expression of the reduced Hamiltonian function considered here. Instead, the differential appearing on the right hand side of (<ref>) can be approximated efficiently based on the finite difference trick proposed in <cit.> (Section 4.1). Indeed, it can be rewritten as follows:

    dF(q,p)^T Z = [ ∂_p(∂_q H_r) ·q̃ - ∂_q (∂_q H_r) ·p̃;  ∂_p(∂_p H_r)·q̃ - ∂_q (∂_p H_r) ·p̃ ]

which only involves directional derivatives of the components of the function F in the directions of q̃ and p̃. We then specifically approximate the above by centered finite difference 

    dF(q,p)^T Z ≈[  α; -β ],   with [ α; β ] = F(q-ϵp̃,p+ϵq̃)-F(q+ϵp̃,p-ϵq̃)/2ϵ

for some small ϵ>0, which only requires at each time t two evaluations of the same function F that appears in the forward reduced Hamiltonian equations.   

With the above approach to compute the gradient with respect to p(0), the registration algorithm then consists of essentially the same steps as the aforementioned works:   

For the numerical ODE integration steps of lines 2 and 4, we use a standard RK4 scheme with regular time samples in [0,1], where we typically take T=15 time steps in most of the experiments that we present in the next section. Note that one can easily replace the RK4 scheme by even higher order or adaptive step methods although in practice we have found this to be unnecessary for the types of ODEs involved here. The optimization update in line 6 follows the limited memory BFGS algorithm, specifically the implementation provided by the HANSO library <cit.>. One can further take additional advantage of the dimensionality reduction provided by Proposition <ref> by restricting each of the components p_i^u_k(0) to the linear subspace Span({u_i^(ℓ)(0)}_ℓ≠ k)^⊥. Lastly, as in Section <ref>, all kernel summation and differentiation operations appearing in both the varifold fidelity terms and Hamiltonian equations are coded in CUDA using the KeOps library <cit.>. The full implementation of the varifold approximation and diffeomorphic registration approach is available at <https://github.com/charoncode/Var_LDDMM> together with the scripts and data of some of the simulations presented in the next section. 




§ RESULTS


We now present some results of the previous algorithms on discrete varifolds of dimension d=1 and d=2. In all these experiments, we choose the deformation kernel K of V to be a diagonal Gaussian kernel K(x,y)=exp(-|x-y|^2/σ_V^2) Id. The kernel function ρ is a Gaussian of scale σ_ρ. The choice of these scales is adapted to the sizes of the shapes in each of the experiment. We will not discuss these questions more in detail here, since this is not our main topic and it has been more thoroughly analyzed in previous works such as <cit.>. The function γ is chosen, depending on the situation, in the different classes of functions discussed in detail in <cit.>, the main distinction being whether the considered varifolds are rectifiable or not according to the conditions given by Theorem <ref> and Theorem <ref> and whether the shapes carry a relevant orientation or not. In particular, one can use γ(t) = t^2 to recover an orientation invariant fidelity metric, or γ(t) = e^-2/σ_γ^2(1-t) which leads to an orientation sensitive distance that satisfy the conditions of Theorem <ref>. All simulations are run on a desktop computer equipped with a NVIDIA Quadro P5000 graphics card.  






 §.§ Diffeomorphic registration




We start with results of registration obtained from the algorithm of Section <ref>. In this section, we will mostly focus on examples involving 2-varifolds, the reader may refer to <cit.> for additional examples in the case d=1. First, as a sanity check, we compare our 2-varifold registration approach applied to triangulated surfaces with the previous LDDMM mesh surface matching implementation of <cit.> using the same kernel size parameters, in which case we expect both approaches to be theoretically equivalent as pointed out in the last paragraph of Section <ref>. Shown in Fig. <ref> are triangulated surfaces of amygdala segmented from two different subjects of the BIOCARD database <cit.>, containing 563 vertices, 1122 triangles and 488 vertices, 972 triangles respectively. Following the simple procedure outlined at the beginning of Section <ref>, we obtain discrete 2-varifolds (one Dirac for each triangle). The first row in the figure shows the optimal deformation estimated with our approach through the evolution of the discrete varifold of the source shape (red) to the target varifold (blue). Discrete varifolds are here displayed in the form of tangent patches and normal vectors (instead of 2-frames) for the purpose of better visualization. Now, the estimated vector fields v_t define a path of dense deformations of the full space which we can also apply to deform the original triangulated surface, which we show on the second row of Fig. <ref>. This is very comparable to the result of the surface mesh LDDMM registration approach displayed on the third row. In terms of computation times, the varifold registration takes a total of 494s (0.99s per iteration of BFGS) against 92.5s (0.18s per iterations) for the surface LDDMM algorithm. This difference comes from mainly two factors: the fact that the numerical complexities are quadratic in the number of Diracs (i.e. triangles) for varifold matching as opposed to the number of vertices for surface LDDMM, and from the increased dimensionality of the Hamiltonian systems in our model. 


In Fig. <ref>, we consider a more challenging registration scenario which was originally studied in <cit.>. Here, one of the two shape is a triangulated surface of a heart membrane segmented from high resolution CT imaging while the second one only consists of a sparse set of cross-sectional curves of the heart contour obtained from lower resolution clinical cardiac MRI data. The varifold framework of this paper leads to an alternative registration approach to the one proposed in <cit.> that relies on a tailored closest point fidelity cost for the surface to curve set comparison. In our case, we instead represent both shapes as 2-varifolds and register them using the exact same varifold registration algorithm as in the previous example. The triangulated surface is again associated to a discrete 2-varifold in the same way as above. As for the set of cross-sectional curve set, we first obtain its 1-varifold representation {x_i,u_i^(1)} which involve the tangent vectors u_i^(1) to the curve that passes through x_i. We then complete it into a 2-varifold by adding a second "vertical" (i.e. inter-sectional) frame vector u_i^(2), which can be estimated in this case by simply finding the projection of x_i onto the corresponding curve in the section immediately above (note that this does involve any attempt to estimate an actual surface mesh of the data). We show the 2-varifolds associated to each shape in the first row of Fig. <ref> and as well as the result of the 2-varifold registration both from curve set to surface and surface to curve set. In each case, we have again applied the estimated deformation between varifolds on the original shapes for visualization.




Along the same lines, we finally look into the case of even less structured data objects. Specifically, as displayed on the first row of Fig. <ref>, we consider two noisy point clouds which are obtained by first randomly selecting vertices from the groundtruth surfaces (with replacement) and then adding some Gaussian noises (σ=0.028) to the position of each sampled point. A first possible registration approach could be to treat such point clouds as standard measures of ^3 (i.e. 0-varifolds) and follow the simple point distribution LDDMM algorithm for unlabelled point sets proposed in <cit.>. The result shown on the third row of Fig. <ref> illustrates the shortcomings of such a model for this type of data. Indeed, one can see that, in the absence of any tangential information, many details of the target shape are not well-recovered. Furthermore, this point set model is not robust to sampling changes and imbalances which results in the mismatches observed below the ear region. An arguably more adequate method would be to exploit the fact that these point clouds are close to their underlying surfaces. However, due to noise and the presence of outliers, estimating triangulations of the point clouds with standard meshing algorithms can prove particularly challenging and inefficient. Instead, our approach consists in directly learning the 2-varifold structure from the point clouds based on the geometric multi-resolution analysis (GMRA) framework developed in <cit.>. Here, we fix a specific scale and GMRA then provides local partitions with estimates of tangent planes to the point clouds which eventually gives us an approximate representation as a 2-varifold illustrated on the second row of Fig. <ref>. Besides its robustness and numerical efficiency, such manifold learning algorithm is also particularly well suited for our proposed registration framework since it naturally leads to approximations in the form of 2-varifolds (and generally not meshes). In the last row of Fig. <ref>, we show the deformed point cloud resulting from the deformation estimated by the 2-varifold registration algorithm. It obviously outperforms the direct point cloud registration described above both in terms of quality of matching but also computation time (10 mins vs 39 mins in total).





 §.§ Approximation and registration


In this second part, we examine some results of the varifold quantization procedure proposed in Section <ref>, and in particular its interplay with the registration algorithm. Specifically, we wish to numerically validate the statements of Corollary <ref> and Theorem <ref>. We shall consider the following protocol. Starting from a highly sampled shape (that we treat as the groundtruth) for which the associated varifold μ_0 is composed of a very high number of Diracs, we compute the compressed varifolds given by the μ_N of (<ref>) for increasing values of N and evaluate the resulting quantization error in terms of the d_W^* metric. Then we solve the registration problems to a fixed target μ_tar from the source varifolds given by the μ_N in lieu of μ_0, and compare the estimated solutions to the registration of the groundtruth. For comparison, we will evaluate the total energy E(v^N) of the estimated deformation fields v^N for the original problem, i.e. 

    E(v^N) = ∫_0^1 v^N_t_V^2 dt + λ(φ_1^v^N)_#μ_0 - μ_tar_W^*^2.

We shall also compare this overall approach against the alternative idea of directly subsampling the original meshes and registering those subsampled shapes with point set mesh LDDMM.



We begin with a 1-varifold toy example given by the curves shown in Fig. <ref> from the Kimia database. These very simple curves segmented from binary images have a relatively high number of points to start with (368 vertices and edges). We look first at how well they can be approximated with smaller number of Diracs through the quantization approach described above. The upper row shows the plot of the relative approximation error μ_N - μ_0_W^*/μ_0_W^* of the source curve as a function of N (blue) as well as the same error in varifold norm when instead the curve is uniformly subsampled (green). Consistent with the fact that varifold quantization should provide the optimal error rate at a given N, we observe that the error is indeed smaller than with the subsampling approach. We also display a few of the quantized μ_N for several values of N. As a second step, we compute the optimal deformations from the reduced shapes to the fixed target and compare their registration energies to the "groundtruth" E(v^*) estimated from the full resolution source shape. The corresponding plots for the quantization versus subsampling methods are shown on the lower row in blue and green respectively. It suggests again a faster convergence to the optimal energy E(v^*) with the quantization strategy, although the difference between the two methods is rather tenuous in this example.     

Those effects can be much more significant in the two-dimensional case. We emphasize it with the triangulated heart surfaces of Fig. <ref> (data courtesy of C. Chnafa, S. Mendez and F. Nicoud, University of Montpellier). The source surface has a total of 42448 triangles leading to the same number of Diracs for the source 2-varifold μ_0 and thus compressing the representation may be in that case quite critical from a computational standpoint. Indeed, computing the groundtruth matching at full resolution takes more than 7 hours (68s per iteration) in this case. We again compare two approaches: our quantization algorithm applied to μ_0 versus directly subsampling the triangulated surface itself (we use here the reducepatch function in MATLAB to reduce the initial mesh to a given number of triangles). For both methods, we compute the relative approximation error μ_N - μ_0_W^*/μ_0_W^* with different values of N, the number of Diracs (resp. triangles) of the compressed varifold (resp. mesh). This is shown on the left second row in Fig. <ref>. Unsurprisingly, we see that the quantization approach leads to a much faster decrease in the error as a function of N but that in addition we obtain a very good approximation of μ_0 with only a small fraction of the initial number of Diracs. Some of the quantized varifolds μ_N are displayed in the figure. We also evaluate how well the solution of the registration problem to the target varifold or surface can be approximated based on the quantized source shapes. With v^* being a numerical solution for the groundtruth and v^N the solutions based on the quantized source shapes, the third row of Fig. <ref> shows the difference of the energies E(v^N)-E(v^*). We observe again a faster convergence towards the groundtruth optimal energy with the varifold quantization than with mesh subsampling.    




§ DISCUSSION

In this paper, we proposed a registration framework between varifolds that goes beyond the previous restrictions of such models to the registration of discrete or smooth submanifolds of ^n. To achieve so, we studied a general class of distances between oriented varifolds based on reproducing kernels and derived a deformation model on the space 𝒱_d, which are combined into an optimal control formulation of the registration problem between any two varifolds. We also examined the possibility to couple this approach with a quantization/compression methodology in order to eventually tackle the registration problem, in practice, on discrete varifolds with a relatively low number of Dirac masses.     

We showed that first of all this setting leads to an equivalent yet alternative formulation to the diffeomorphic registration of rectifiable sets such as continuous or discrete curves and surfaces; the resulting higher-order Hamiltonian systems in our model provides richer local patterns for the deformations but at the price of a higher numerical cost. From an application standpoint, however, the main advantage we expect from this framework is that it applies very naturally to more general geometric objects, in particular to typical situations where well-defined and reliable meshes are not available. We gave a taste of it through some of the examples of Section <ref>, although future work on a larger scale will be needed in order to evaluate such benefits more thoroughly. Besides the cases mentioned here, there are also several types of data that could constitute interesting test applications for this setting. This includes for instance high-angular resolution diffusion MRI in which the data is effectively modeled as spatially distributed orientation probability distribution functions consistent with the Young measure representation of varifolds in (<ref>), or the case of contrast-invariant image registration c.f. <cit.>. 

At the theoretical level, there are several questions left open by this work which we believe can constitute interesting tracks for future work. One is to study the possibility of extending all or part of the results of Section <ref> to more general kernel metrics (in particular currents) and determining tighter quantization error bounds. Moreover, the registration model at play in this paper is based on the pushforward group action of Diff(^n) on 𝒱_d. Yet, other group actions could be have been considered, as briefly evoked in Section <ref>, that involve different choices of reweighing factor, for which we could expect very different properties of the solutions to the registration problem. 

Lastly, some additional work on the numerical side is likely needed for potential future applications to large scale databases, most notably to generalize this work to the estimation of means and atlases over populations of many high resolution shapes. Indeed, as we pointed out, even with the ability to compress the size of varifolds in the registration pipeline using the quantization approach, the higher complexity of the dynamical equations involved in the registration model has a non-negligible numerical toll. This could be improved in the future by using more efficient computational schemes for the repeated evaluations of sums of kernels and derivative of kernels appearing in the Hamiltonian equations, possibly along the lines of fast multiple methods.         



§ ACKNOWLEDGEMENTS

The authors would like to thank Benjamin Charlier, Siamak Ardekani, Laurent Younes and the BIOCARD team for sharing the data used in some of the examples of this paper. This work was supported by NSF grant No 1819131.






§ APPENDIX

Proof of Theorem <ref>
1ex
We first prove that ℋ^d(X△ Y) =0. Let us denote by W^pos and W^G the RKHS associated to kernels k^pos and k^G respectively. Suppose that X and Y are rectifiable sets as above such that μ_X-μ_Y _W^*=0 and ℋ^d(X△ Y) >0. Without loss of generality, we may assume that ℋ^d(X ∖ Y) >0. From Lusin's theorem, there exists a subset U of X such that T|_U is continuous and ℋ^d(X∖ U) <  ℋ^d(X∖ Y). Let us denote by E := U ∩ (X ∖ Y), we see that ℋ^d(E)>0. Since for ℋ^d  a.e.  x∈ E, 


    lim sup_r → 0ℋ^d(B_r(x)∩ E )/π^d/2/Γ(d/2+1)r^d≥1/2^d,

(cf <cit.>), there exists x_0 ∈ E,  ℋ^d(B_r(x_0) ∩ E)>0 for any r >0.  

Let g:G^n_d → be defined by g(·) = k^G(T(x_0),·). Since x ⟼ g(T(x)) is continuous on E and g(T(x_0))>0, there exists r_0>0 such that ∀ x ∈ B_r_0(x_0) ∩ E,  g(T(x))>0. Let A ≐ B_r_0(x_0) ∩ E and h(x) := 1_A(x), then ℋ^d(A)>0 and g(T(x))>0,  ∀ x ∈ A. Using the density of C_c(ℝ^n) in L^1(ℝ^n,ℋ^d  (X ∪ Y)) together with the fact that k^pos is C_0-universal, there exist {f_j}_j=1^∞⊂ C_c(ℝ^n) and {h_j}_j=1^∞⊂ W^pos 
such that lim_j →∞ f_j =h in L^1(ℝ^n,ℋ^d  (X ∪ Y)) and f_j-h_j _∞< 1/j. Now, since h_j ⊗ g ∈ W and μ_X = μ_Y in W^*, we have


    0 = (μ_X-μ_Y)(h_j ⊗ g) = ∫_X h_j(x) g(T(x)) d ℋ^d(x) - ∫_Y h_j(y) g(S(y)) d ℋ^d(y) →∫_A g(T(x)) d ℋ^d(x) >0,

which is a contradiction. Hence we have ℋ^d(X△ Y) = 0
 
Next, we show that T(x) = S(x) ℋ^d-a.e.. Let F := {x ∈ X | T(x) = -S(x)} and assume that ℋ^d(F)>0. From Lusin's theorem, there exists subset F' ⊂ F such that T|_F' is continuous and ℋ^d(F')>0. Using the upper density argument as above, we can find z_0 ∈ F' such that ℋ^d(B_r(z_0) ∩ F')>0 for all r>0. Since the map x ↦⟨ T(x),T(z_0) ⟩ restricted to F' is continuous, there exists a δ_0>0 satisfying:

    ⟨ T(x),T(z_0) ⟩ >0,  ∀ x ∈ B_δ_0(z_0) ∩ F'.

Define B:= B_δ_0(z_0) ∩ F', η(·):= γ(⟨·,T(z_0) ⟩) and u(x) := η(T(x)) - η(S(x)). Observe that, from the assumption γ(t) ≠γ(-t),  ∀ t ∈ [-1,1],

    u(x) = η(T(x)) - η(-T(x)) ≠ 0,  ∀ x ∈ F'.

From this, we may assume that u(x)>0,  ∀ x ∈ F'. Let {f_j'}_j and {h_j'}_j be sequences in C_c(ℝ^n) and W_pos such that f_j' converges to 1_B in L^1(ℝ^n,ℋ^d  F) and f_j'-h_j'_∞ < 1/j. We obtain
 
 
    0 = (μ_X - μ_Y| h_j' ⊗η) = ∫_X h_j'(x) u(x) dℋ^d(x) →∫_B u(x) dℋ^d(x) > 0,

which is impossible. 3ex


Proof of Theorem <ref>
1ex
 Thanks to the first term in E, any minimizing sequence of E is bounded in L^2([0,1],V). Let {v_j} be a subsequence of such minimizing sequence which converges weakly to some v̅ in L^2([0,1],V). Using the results of <cit.> Chapter 7.2, we know that
 
    lim_j →∞(φ_1^v_j - φ_1^v̅)|_K _1,∞ = 0.

Furthermore, for any ω∈ W, we have

    |( (φ_1^v_j)_#μ_0 - (φ_1^v̅)_#μ_0 | ω) | 
       =|∫_K J_S φ_1^v_j(x) ω(φ_1^v_j(x),d_x φ_1^v_j· S)  
    - J_Sφ_1^v̅(x) ω(φ_1^v̅(x),d_x φ_1^v̅· S) d μ_0 |
       ≤∫_K |J_S φ_1^v_j(x)| 
     |ω(φ_1^v_j(x),d_x φ_1^v_j· S) -  ω(φ_1^v̅(x),d_x φ_1^v̅· S) |  d μ_0 
       + ∫_K|J_S φ_1^v_j(x) - J_Sφ_1^v̅(x) |  |ω(φ_1^v̅(x),d_x φ_1^v̅· S)  | d μ_0

Now, using the embedding W ↪ C_0^1(ℝ^n ×G^n_d)

    |( (φ_1^v_j)_#μ_0 - (φ_1^v̅)_#μ_0 | ω) | 
       ≤( ∫_K |J_S φ_1^v_j(x)| d μ_0 ) ω_1,∞ (φ_1^v^N - φ_1^v̅)|_K _1,∞ 
    + C  (φ_1^v_j - φ_1^v̅)|_K _1,∞
       ≤ C'  (φ_1^v_j - φ_1^v̅)|_K _1,∞.

Taking supremum over all ω∈ W with ω_W ≤ 1, we obtain that

    (φ_1^v_j)_#μ_0 - (φ_1^v̅)_#μ_0 _W^*≤  C'  (φ_1^v_j - φ_1^v̅)|_K _1,∞→ 0

as j →∞. Combining this with lower semicontinuity of v ↦v_L^2([0,1],V)^2, we finally obtain that

    E(v̅) ≤lim inf_j →∞ E(v_j)

and hence v̅ is a global minimizer. 


3ex

Proof of Proposition <ref>
1ex
Recall that for all ϕ∈Diff(^n), g(ϕ) = λϕ_#μ_0 - μ_tar_W^*^2 which we may rewrite as

    g(ϕ) = λ (ϕ_#μ_0 | K_W(ϕ_#μ_0 -2μ_tar)) + λμ_tar_W^*^2.

Thus, the variation with respect to ϕ in the Banach space ℬ writes

    ∂_ϕ g(ϕ)  = ∂_ϕ (ϕ_#μ_0 | ω_0)

where ω_0 ≐ 2λ  K_W(ϕ_#μ_0 -μ_tar) ∈ W. Moreover

    (ϕ_#μ_0 | ω_0) = ∫_^n ×G_d^nω_0(ϕ(x),d_x ϕ· T) J_T ϕ(x) d μ_0(x,T).

Taking the variation with respect to ϕ along any u ∈ C_0^1(^n,^n), we obtain:

    (∂_ϕ g(ϕ) | u)    = ∫_^n ×G_d^n∂_x ω_0(ϕ(x),d_x ϕ· T) · u(x) J_T ϕ(x) d μ_0(x,T) 
       +∫_^n ×G_d^n∂_T ω_0(ϕ(x),d_x ϕ· T) · (d_x u |_T) J_T ϕ(x) d μ_0(x,T) 
       +∫_^n ×G_d^nω_0(ϕ(x),d_x ϕ· T). div_T u(x). J_T ϕ(x) d μ_0(x,T)

where the last term follows from the differentiation of Gram determinant matrices while the notation ∂_T in the second term is a shortcut notation for differentiation on the Grassmannian which we do not explicit further here, we however refer to the similar computations done in <cit.> and to the developments in Section <ref> for more details. For the first term, we can rely on the Young measure decomposition μ_0 = |μ_0| ⊗ν_x introduced at the end of Section <ref> which gives:

    (1) = ∫_^nα̃(ϕ,x) · u(x)  d |μ_0|(x),   where  α̃(ϕ,x) = ∫_G_d^n∂_x ω_0(ϕ(x),d_x ϕ· T) J_T ϕ(x) dν_x(T).

We can also rewrite the third term as:

    (3) = ∫_^n ×G_d^nγ̃(ϕ,x,T)   div_T u(x)  d μ_0(x,T),   with  γ̃(ϕ,x,T) =  ω_0(ϕ(x),d_x ϕ· T)  J_T ϕ(x).

As for the second term in (<ref>), for each (x,T) the integrand involves a linear combination (depending on ϕ) of the partial derivatives of u along the subspace T i.e. of the elements of the matrix d_x u|_T ∈^n× d. Thus, without attempting to specify this term explicitly, we can in general write it as β̃(ϕ,x,T) d_x u|_T where B̃ is a continuous map from ℬ×^n ×G̃_d^n into ℒ(^n× d,) giving us 

    (2) = ∫_^n ×G_d^nB̃(ϕ,x,T) d_x u|_T  d μ_0(x,T).

The result of the theorem then follows by setting α(x) ≐α̃(φ_1^v,x), β(x,T) ≐β̃(φ_1^v,x,T) and γ(x,T)=γ̃(φ_1^v,x,T). 


3ex

Proof of Proposition <ref>
1ex
We can treat the case of each particle i separately and thus, without loss of generality, we may directly assume that N=1. We write q(t)=(x(t),u^(1)(t),⋯,u^(d)(t)), p(t)=(p^x(t),p^u_1(t),…,p^u_d(t)) for the state and costate variables along an optimal trajectory and 

    U ≐Span{u^(1)(1),⋯,u^(d)(1)}.

Consider the group of linear transformations, G ≐ SL (U) ⊕ GL (U^⊥), i.e., for any g∈ G,

    g(x) = g_(x_U) + g_⊥(x_U^⊥),

where x_U and x_U^⊥ are the orthogonal projections of x on U and U^⊥, with   g_∈ SL(U) and g_⊥∈ GL(U^⊥). The Lie algebra of G is 𝔤 = 𝔰𝔩(U) ×ℒ(U^⊥) and 𝔰𝔩(U) is the set of all zero trace linear transformations of U. Now, consider the action of G on ℝ^(d+1)n defined as:

    g· q := (q_0,g(q_1),⋯,g(q_d)).

for any q=(q_0,…,q_d) ∈ℝ^(d+1)n. We see that μ^g· q(1) = μ^q(1) for all g∈ G and therefore g(g· q(1)) = g(q(1)). 

Now, if we let {g_t} be a smooth curve in G that satisfies g_0 = id and d/dτ|_τ=0g_τ = h ∈𝔤, differentiating the equality g(g_τ· q(1)) = g(q(1)) shows that for any h ∈𝔤, we have

    0 = (p(1)|h · q(1)) = ∑_k=1^d⟨ p^u_k(1) , h(u^(k)(1)) ⟩

Since h ∈𝔤, we must have that h|_U is a zero trace linear map. For any 1 ≤ i < j ≤ d, we may choose h such that h(u^(i)(1)) = -h(u^(j)(1)) and h(u^(k)(1)) = 0 ,  ∀ k ∉{i,j}, which leads to ⟨ u^(i)(1), p^u_i(1) ⟩ = ⟨ u^(j)(1), p^u_j(1) ⟩. Consequently,

    ⟨ u^(1)(1), p^u_1(1) ⟩ = ⋯ = ⟨ u^(d)(1), p^u_d(1) ⟩ =α

for some constant α. In addition, for any i ≠ j, we can also choose h such that h(u^(i)(1)) = u^(j)(1) and h(u^(k)(1))= 0 ,  ∀ k ∉{i,j}, which gives ⟨ u^(i)(1),p^u_j(1) ⟩ = 0. It results that D(1) = α.I_d × d.  


Finally, since D(t) is constant by Lemma <ref>, we obtain that

    D(t) = ( [ α   0;   ⋱  ; 0   α ]),
 
for all t ∈ [0,1]. 


apacite



