Metrics, quantization and registration in varifold spaces

TL;DR

This paper advances the theory and algorithms for shape analysis using varifold spaces, enabling registration of diverse geometric objects through new metrics, models, and numerical methods with proven convergence properties.

Contribution

It introduces a generalized framework for diffeomorphic registration of varifolds, including new metrics, a mathematical model, and quantization techniques with convergence guarantees.

Findings

Development of kernel metrics on varifold spaces

A mathematical model for diffeomorphic registration under group actions

Numerical pipelines demonstrating initial results in 1D and 2D

Abstract

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 $\Gamma$-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.

Figures

Equations

$$i_P:{\tilde{G}}_d^n$$

$$T$$

$$⟨\xi ,\eta ⟩=\det ({\xi }_i\cdot {\eta }_j{)}_{i,j=1,\dots ,d}$$

$$A\cdot T:=\frac{A{u}_1\wedge \cdots \wedge A{u}_d}{|A{u}_1\wedge \cdots \wedge A{u}_d|}$$

$$(\mu |\omega )\doteq {\int }_{{\mathbb{R}}^n\times {\tilde{G}}_d^n}\omega (x,T)d\mu (x,T).$$

$$(\mu |\omega )={\int }_{{\mathbb{R}}^n}{\int }_{{\tilde{G}}_d^n}\omega (x,T)d{\nu }_x(T)d|\mu |(x).$$

$$({\delta }_{(x,T)}|\omega )=\omega (x,T),\forall \omega \in C_0({\mathbb{R}}^n\times {\tilde{G}}_d^n).$$

$$({\mu }_X|\omega )$$

$$={\int }_X\omega (x,T_X(x))d{\mathcal{H}}^d(x).$$

$$d_{op}(\mu ,\nu )\doteq \underset{|\omega {|}_{\infty }\le 1}{\sup }(\mu -\nu |\omega ),\forall \mu \in {\mathcal{V}}_d.$$

$$\underset{i\to \infty }{\lim }({\mu }_i|\omega )=(\mu |\omega ).$$

$$d_{\ast }(\mu ,\nu )=\sum _{k\in \mathbb{N}}{2}^{-k}|(\mu -\nu |{\omega }_k)|,$$

$$d_{Was{s}^1}(\mu ,\nu )\doteq \underset{\text{Lip}(\omega )\le 1}{\sup }|(\mu -\nu |\omega )|.$$

$$d_{BL}(\mu ,\nu )\doteq \underset{\Vert \omega {\Vert }_{\infty },\text{Lip}(\omega )\le 1}{\sup }|(\mu -\nu |\omega )|.$$

$$\displaystyle\eta(x,T)\doteq\left\{\begin{array}[]{lc}\omega(x,T),&\ \text{if}\ (x,T)\in K\\ 0,&\ \text{if}\ (x,T)\in B^{c}\end{array}\right.$$

$$\left|{\int }_{{\mathbb{R}}^n\times {\tilde{G}}_d^n}\omega d({\mu }_i-\mu )\right|$$

$$\le \left|{\int }_{{\mathbb{R}}^n\times {\tilde{G}}_d^n}\tilde{\omega }d({\mu }_i-\mu )\right|+\epsilon .$$

$$\underset{i\to \infty }{\mathrm{limsup}}\left|{\int }_{{\mathbb{R}}^n\times {\tilde{G}}_d^n}\omega d({\mu }_i-\mu )\right|<\epsilon .$$

$$({\iota }^{\ast }\mu |\omega )={\int }_{{\mathbb{R}}^d\times {\tilde{G}}_d^n}\omega (x,T)d\mu (x,T),\forall \omega \in W.$$

$$d_{W^{\ast }}(\mu ,{\mu }^{\prime }{)}^2=\Vert \mu -{\mu }^{\prime }{\Vert }_{W^{\ast }}^2=\Vert \mu {\Vert }_{W^{\ast }}^2-2⟨\mu ,{\mu }^{\prime }{⟩}_{W^{\ast }}+\Vert {\mu }^{\prime }{\Vert }_{W^{\ast }}^2$$

$$⟨\mu ,{\mu }^{\prime }{⟩}_{W^{\ast }}$$

$$={\int }_{({\mathbb{R}}^d\times {\tilde{G}}_d^n{)}^2}\rho (|x-x^{\prime }{|}^2)\gamma (⟨T,T^{\prime }⟩)d\mu (x,T)d{\mu }^{\prime }(x^{\prime },T^{\prime })$$

$$\mu \mapsto {\int }_{{\mathbb{R}}^d\times {\tilde{G}}_d^n}k(\cdot ,(y,T))d\mu (y,T),\mu \in C_0({\mathbb{R}}^d\times {\tilde{G}}_d^n{)}^{\ast }$$

$$\Vert \mu -{\mu }^{\prime }{\Vert }_{W^{\ast }}=\underset{\omega \in W,\Vert \omega {\Vert }_{W\le 1}}{\sup }{\int }_{{\mathbb{R}}^d\times {\tilde{G}}_d^n}\omega d(\mu -{\mu }^{\prime })\le c_W{d}_{op}(\mu ,{\mu }^{\prime }).$$

$$\Vert {\mu }_i{\Vert }_{W^{\ast }}^2$$

$$\to {\int }_{({\mathbb{R}}^d\times {\tilde{G}}_d^n{)}^2}k((x,S),(y,T))d\mu (x,S)d\mu (y,T)$$

$$=\Vert \mu {\Vert }_{W^{\ast }}^2,$$

$$|({\mu }_i-\mu |f)|\le |({\mu }_i|f-g)|+|(\mu |g-f)|+|({\mu }_i-\mu |g)|\le \epsilon +\Vert {\mu }_i-\mu {\Vert }_{W^{\ast }}\Vert g{\Vert }_W.$$

$$\varphi \cdot (X,T_X)\doteq (\varphi (X),T_{\varphi (X)}),$$

$$T_{\varphi (X)}(y)\doteq d_{{\varphi }^{-1}(y)}\varphi \cdot T_X({\varphi }^{-1}(y))$$

$$({\varphi }_{\#}\mu |\omega )\doteq {\int }_{{\mathbb{R}}^d\times {\tilde{G}}_d^n}\omega (\varphi (x),d_x\varphi \cdot T)J_T\varphi (x)d\mu (x,T)$$

$$J_T\varphi (x)=\det \left((d_x\varphi (e_i)\cdot d_x\varphi (e_j){)}_{i,j=1,\dots ,d}\right)$$

$$({\varphi }_{\ast }\mu |\omega ):={\int }_{{\mathbb{R}}^d\times {\tilde{G}}_d^n}\omega (\varphi (x),d_x\varphi \cdot T)d\mu (x,T).$$

Full text

Mathematics Subject Classification (2010) 49Q20  $\cdot$ 49M25  $\cdot$ 58E50  $\cdot$ 68T10

.

1. 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 Grenander (\APACyear1993) 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 Beg \BOthers. (\APACyear2005); Younes (\APACyear2019) 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 4.2). 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. Miller \BBA Qiu (\APACyear2009); Miller, Trouvé\BCBL \BBA Younes (\APACyear2015) 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 Federer (\APACyear1969). 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 Cohen-Steiner \BBA Morvan (\APACyear2003); Buet \BOthers. (\APACyear2018) or estimation of shape medians Hu \BOthers. (\APACyear2018). With regard to the aforementioned deformation analysis problems, the potential interest of geometric measure theory has been identified early on in the works of Glaunès \BOthers. (\APACyear2004); Glaunès \BBA Vaillant (\APACyear2006). 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: Glaunès \BBA Vaillant (\APACyear2006); Glaunès \BOthers. (\APACyear2008); Durrleman \BOthers. (\APACyear2010) are based on the representation of oriented curves and surfaces as currents, Charon \BBA Trouvé (\APACyear2013) and Kaltenmark \BOthers. (\APACyear2017); Bauer \BOthers. (\APACyear2019) extended this model to the setting of unoriented and oriented varifolds, while Roussillon \BBA Glaunès (\APACyear2016, \APACyear2017) considers the higher-order representation of normal cycles, see also the recent survey Charon \BOthers. (\APACyear2020). 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 Kaltenmark \BOthers. (\APACyear2017). 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 Hsieh \BBA Charon (\APACyear2019) 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 Charon \BBA Trouvé (\APACyear2013); Kaltenmark \BOthers. (\APACyear2017), 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 3 after the brief introduction to the notion of oriented varifold of Section 2. In Section 4, 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 5 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 $\Gamma$-convergence of the corresponding approximate registration functionals. In Section 6, 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 7, 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.

2. The space of oriented varifolds

The concept of varifold was originally developed in the context of geometric measure theory by Young (\APACyear1942), Almgren (\APACyear1966) and Allard (\APACyear1972) for the study of Plateau’s problem on minimal surfaces. The interest in registration and shape analysis was evidenced in Charon \BBA Trouvé (\APACyear2013); Kaltenmark \BOthers. (\APACyear2017). 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 Charon \BBA Trouvé (\APACyear2013). Later on, the generalized framework of oriented varifold was proposed in Kaltenmark \BOthers. (\APACyear2017) 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 Charon \BBA Trouvé (\APACyear2013) can be recovered almost straightforwardly through adequate choices of orientation-invariant kernels as we shall briefly point out later on.

2.1. Definition

The underlying principle of varifolds is to extend measures of $\mathbb{R}^{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\leq d\leq n$, we first need to introduce the set of all possible $d$-dimensional oriented tangent spaces in $\mathbb{R}^{n}$:

Definition 2.1**.**

The $d$-dimensional oriented Grassmannian $\widetilde{G}_{d}^{n}$ is the set of all oriented $d$-dimensional linear subspaces of $\mathbb{R}^{n}$.

The oriented Grassmannian is a compact manifold of dimension $d(n-d)$ which can be identified to the quotient $SO(n)/(SO(d)\times SO(n-d))$. It is also a double cover of the (non-oriented) Grassmannian $G_{d}^{n}$ of $d$-dimensional subspaces of $\mathbb{R}^{n}$. For practical purposes, a more convenient representation of $\widetilde{G}_{d}^{n}$ is the one detailed in the following remark.

Remark 2.2*.*

Given $T\in\widetilde{G}_{d}^{n}$, there exists a basis $\{u_{i}\}_{i=1,\ldots,d}\in\mathbb{R}^{n\times d}$ of $T$ such that $[u_{1},\cdots,u_{d}]$ has consistent orientation with $T$. Then the following map, called the oriented Plücker embedding, is well defined and injective,

[TABLE]

This allows to identify $\widetilde{G}^{n}_{d}$ as a subset of the unit sphere of $\Lambda^{d}(\mathbb{R}^{n})$ which inherits the topology of the inner product on $\Lambda^{d}(\mathbb{R}^{n})$. We remind that this inner product is defined for any $\xi=\xi_{1}\wedge\ldots\wedge\xi_{d}$, $\eta=\eta_{1}\wedge\ldots\wedge\eta_{d}$ in $\Lambda^{d}(\mathbb{R}^{n})$ by the determinant of the Gram matrix:

[TABLE]

Through this identification, one can also define the action of linear transformations on $\widetilde{G}^{n}_{d}$ as follows

[TABLE]

for any $T\in\widetilde{G}^{n}_{d}$ and $A:\mathbb{R}^{n}\mapsto\mathbb{R}^{n}$ a linear invertible map.

Similar to the definition of classical varifolds in Simon (\APACyear1983), we define oriented varifolds as measures on $\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$.

Definition 2.3**.**

An oriented $d$-varifold $\mu$ on $\mathbb{R}^{n}$ is a nonnegative finite Radon measure on the space $\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$. Its weight measure $|\mu|$ is defined by $|\mu|(A):=\mu(A\times\widetilde{G}^{n}_{d})$ for all Borel subset $A$ of $\mathbb{R}^{n}$. We denote by $\mathcal{V}_{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 $\mu$ as a distribution, i.e. an element of the dual space $C_{0}(\mathbb{R}^{d}\times\widetilde{G}^{n}_{d})^{*}$, where $C_{0}(\mathbb{R}^{d}\times\widetilde{G}^{n}_{d})$ denotes the set of continuous functions vanishing at infinity on $\mathbb{R}^{d}\times\widetilde{G}^{n}_{d}$. It is defined for any test function $\omega\in C_{0}(\mathbb{R}^{d}\times\widetilde{G}^{n}_{d})$ by:

[TABLE]

As an additional note, another useful representation of a general varifold in $\mathcal{V}_{d}$ can be obtained by the disintegration theorem (see Ambrosio \BOthers. (\APACyear2000) Chap. 2). Namely, if $\mu\in\mathcal{V}_{d}$, for $|\mu|$-almost every $x$ in $\mathbb{R}^{n}$, there exists a probability measure $\nu_{x}$ on $\widetilde{G}^{n}_{d}$ such that $x\mapsto\nu_{x}$ is $|\mu|$-measurable and we can write

[TABLE]

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

2.2. 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\in\mathbb{R}^{n}$ and $T\in\widetilde{G}^{n}_{d}$, the associated Dirac varifold $\delta_{(x,T)}$ acts on functions of $C_{0}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})$ by the relation

[TABLE]

$\delta_{(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 $\mathbb{R}^{n}$ of finite total $d$-volume, denoting by $T_{X}(x)\in\widetilde{G}^{n}_{d}$ the oriented tangent space at $x\in X$, one can associate to $X$ the varifold $\mu_{X}$, which is defined for all Borel subset $B\subset\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$ by $\mu_{X}(B)=\mathcal{H}^{d}(\{x\in X|(x,T_{X}(x))\in B\})$. Here, $\mathcal{H}^{d}$ is the $d$-dimensional Hausdorff measure on $\mathbb{R}^{n}$, i.e. the measure of $d$-volume of subsets of $\mathbb{R}^{n}$ (we refer the reader to Simon (\APACyear1983) for the precise construction and properties of Hausdorff measures). It is then not hard to see that, as an element of $C_{0}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})^{*}$,

[TABLE]

Such a representation $X\mapsto\mu_{X}$ can be extended to slightly more general objects known as oriented rectifiable sets. A subset $X$ of $\mathbb{R}^{n}$ is said to be a countably $\mathcal{H}^{d}$-rectifiable set if $\mathcal{H}^{d}(X\setminus\cup_{j=1}^{\infty}F_{j}(\mathbb{R}^{d}))=0$, where $F_{j}:\mathbb{R}^{d}\mapsto\mathbb{R}^{n}$ are Lipschitz function for all $j$ (c.f. Simon (\APACyear1983)). We say that $(X,T_{X})$ is an oriented rectifiable set if $X$ is a countably $d$-rectifiable set and $T_{X}:X\mapsto\widetilde{G}^{n}_{d}$ is a $\mathcal{H}^{d}$-measurable function such that for $\mathcal{H}^{d}-a.e.\ x\in 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 $\mu_{X}$ given again by (2.2). The set of those $\mu_{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).

Remark 2.4*.*

Rectifiable varifolds still make a very ”small” subset of $\mathcal{V}_{d}$: indeed, in the Young measure representation of (4), we have in this case the very particular constraint that probability measures $\nu_{x}$ are Dirac masses, specifically $\nu_{x}=\delta_{T_{X}(x)}$.

3. Metrics on varifolds

In this section, we address the issue of defining adequate metrics on the space $\mathcal{V}_{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. Charon \BBA Trouvé (\APACyear2013); Kaltenmark \BOthers. (\APACyear2017) and derive the most relevant properties of this class of distances.

3.1. Standard topologies and metrics on $\mathcal{V}_{d}$

As a measure/distribution space, $\mathcal{V}_{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.

• $\bullet$

mass norm: with the previous identification of measures in $\mathcal{V}_{d}$ with elements of the dual $C_{0}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})^{*}$, one can define the following dual metric on $\mathcal{V}_{d}$:

[TABLE]

where $|\omega|_{\infty}\doteq\sup_{\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}}|\omega|$. 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 $\delta_{(x,T)}$ and $\delta_{(x^{\prime},T^{\prime})}$, $d_{op}(\delta_{(x,T)},\delta_{(x^{\prime},T^{\prime})})=2$ whenever $(x,T)\neq(x^{\prime},T^{\prime})$.

• $\bullet$

weak-* topology: a sequence of $d$-varifolds $\{\mu_{i}\}_{i}$ converges to $\mu\in\mathcal{V}_{d}$ in the weak-* topology (denoted by $\mu_{i}\overset{\ast}{\rightharpoonup}\mu$) if and only if for all $\omega\in C_{c}(\mathbb{R}^{d}\times\widetilde{G}^{n}_{d})$ (continuous compactly supported function)

[TABLE]

In fact, the weak-* topology on $\mathcal{V}_{d}$ can be metrized by the following distance:

[TABLE]

where $\{\omega_{k}\}_{k\in\mathbb{N}}$ is a dense sequence in $C_{c}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})$.

• $\bullet$

Wasserstein metric: the Wasserstein-1 distance of optimal transport can be expressed in its Kantorovitch dual formulation Villani (\APACyear2008) as

[TABLE]

where the sup is taken over all Lipschitz regular functions on $\mathbb{R}^{n}\times\widetilde{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 Piccoli \BBA Rossi (\APACyear2014); Chizat \BOthers. (\APACyear2018) have instead proposed generalized Wasserstein distances derived from unbalanced optimal transport.

• $\bullet$

Bounded Lipschitz metric: similar to the previous, the bounded Lipschitz distance (sometimes referred to as the flat metric) on $\mathcal{V}_{d}$ is defined by

[TABLE]

It can be shown (cf. Ch 8 in Bogachev (\APACyear2007)) that $d_{BL}$ metrizes the narrow topology on $\mathcal{V}_{d}$, namely the topology for which a sequence $(\mu_{i})$ converges to $\mu$ if and only if $\lim_{i\rightarrow\infty}(\mu_{i}|\omega)=(\mu|\omega)$ for all bounded continuous functions $\omega$.

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 $\mathcal{V}_{d}$, namely:

Proposition 3.1**.**

Let $\mu$ and $\{\mu_{i}\}_{i}$ be varifolds such that the sequence $\{\mu_{i}\}_{i}$ is tight. Then $\mu_{i}\overset{\ast}{\rightharpoonup}\mu$ if and only if $d_{BL}(\mu_{i},\mu)\rightarrow 0$.

Proof.

Since $d_{BL}$ metrizes the narrow topology, it suffices to show that $\mu_{i}$ converges to $\mu$ in the narrow topology. Let $\omega$ be a bounded continuous function defined on $\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$ and $\varepsilon>0$. By the tightness property, we may choose a compact set $K\subset\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$ such that $\mu(K^{c})+\sup_{i}\mu_{i}(K^{c})<\varepsilon/2\|\omega\|_{\infty}$. Let $B$ be an open ball that contains $K$. Define

[TABLE]

From Tietz extension theorem, there exists a continuous extension $\tilde{\omega}$ of $\eta$ on $\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$ such that $\tilde{\omega}|_{K}=\omega|_{K}$ and $\tilde{\omega}\in C_{c}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})$. This implies that

[TABLE]

Taking $\limsup$ on both sides, we see that

[TABLE]

Since $\varepsilon$ is arbitrary, we obtain that $\mu_{i}$ converges to $\mu$ in the narrow topology. ∎

As a direct consequence of Proposition 3.1, 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 $\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$. Note also that a very similar result to Proposition 3.1 holds when replacing the bounded Lipschitz distance by generalized Wasserstein metrics, as proved in Piccoli \BBA Rossi (\APACyear2014).

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 Feydy \BOthers. (\APACyear2017). In this paper, we focus on the alternative approach previously developed for currents in Glaunès \BOthers. (\APACyear2008) and unoriented varifolds in Charon \BBA Trouvé (\APACyear2013) which instead relies on particular Hilbert spaces of test functions, as we detail in the next section.

3.2. Kernel metrics

In this section, we start by defining a general class of pseudo-metrics on $\mathcal{V}_{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 3.1.

3.2.1. Kernels for varifolds

We refer the reader to Aronszajn (\APACyear1950); Hastie \BOthers. (\APACyear2001); Micheli \BBA Glaunès (\APACyear2014) 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 $\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$. Along the lines of previous works like Charon \BBA Trouvé (\APACyear2013); Kaltenmark \BOthers. (\APACyear2017), we restrict to separable kernels for which we have:

Proposition 3.2**.**

Let $k^{pos}$ and $k^{G}$ be continuous positive definite kernels on $\mathbb{R}^{n}$ and $\widetilde{G}^{n}_{d}$ respectively. Assume in addition that for any $x\in\mathbb{R}^{n}$, $k^{pos}(x,\cdot)\in C_{0}(\mathbb{R}^{n})$. Then $k:=k^{pos}\otimes k^{G}$ is a positive definite kernel on $\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$ and the RKHS $W$ associated to $k$ is continuously embedded in $C_{0}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})$ i.e. there exists $c_{W}>0$ such that for any $\omega\in W$, we have $\|\omega\|_{\infty}\leq c_{W}\|\omega\|_{W}$.

We recall that the tensor product kernel has the exact expression $k((x,T),(x^{\prime},T^{\prime}))=k^{pos}(x,x^{\prime})k^{G}(T,T^{\prime})$. The proof of Proposition 3.2 is a straightforward adaptation of the same result for unoriented varifolds (cf. Charon \BBA Trouvé (\APACyear2013) Proposition 4.1).

Remark 3.3*.*

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)=\rho(|x-y|^{2}),\ \forall x,y\in\mathbb{R}^{n}$, with $\rho(0)>0$, and $k^{G}(S,T)=\gamma(\langle S,T\rangle),\ \forall S,T\in\widetilde{G}^{n}_{d}$ where $\langle\cdot,\cdot\rangle$ is the inner product on $\widetilde{G}^{n}_{d}$ inherited from $\Lambda^{d}(\mathbb{R}^{n})$ introduced in remark 2.2. 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 Charon \BBA Trouvé (\APACyear2013) can be also recovered in this setting by simply restricting to orientation-invariant kernels $k^{G}$ i.e. such that $\gamma(-t)=\gamma(t)$ for all $t$.

Now, if we let $\iota_{W}:W\hookrightarrow C_{0}(\mathbb{R}^{d}\times\widetilde{G}^{n}_{d})$ be the continuous embedding given by Proposition 3.2 and $\iota_{W}^{*}$ its adjoint, for any $\mu\in C_{0}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})^{*}$, we have

[TABLE]

With (10), we may identify $\mu$ as an element of the dual RKHS $W^{*}$. Note that $\iota_{W}^{*}$ is not injective in general, in other words one can have $\mu=\mu^{\prime}$ in $W^{*}$ but $\mu\neq\mu^{\prime}$ in $C_{0}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})^{*}$.

In any case, one can compare any two varifolds $\mu,\mu^{\prime}\in\mathcal{V}_{d}$ through the Hilbert norm of $W^{*}$ by defining:

[TABLE]

where we use the small abuse of notation of writing $\mu$ and $\mu^{\prime}$ instead of $\iota_{W}^{*}\mu$ and $\iota_{W}^{*}\mu^{\prime}$ on the two right hand sides. Due to the potential non-injectivity of $\iota_{W}^{*}$, in general $d_{W^{*}}$ only induces a pseudo-metric on $\mathcal{V}_{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 $\mu$ and $\nu$ in $\mathcal{V}_{d}$, the inner product between them is given by

[TABLE]

for kernels selected as in Remark 3.3.

3.2.2. 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 Charon \BBA Trouvé (\APACyear2013); Kaltenmark \BOthers. (\APACyear2017). 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:

Definition 3.4**.**

A positive definite kernel $k$ on a metric space $\mathcal{M}$ is called $C_{0}$-universal when its RKHS is dense in $C_{0}(\mathcal{M})$ for the uniform convergence topology.

$C_{0}$-universality has been studied in great length in such works as Carmeli \BOthers. (\APACyear2010); Sriperumbudur \BOthers. (\APACyear2011). In particular, one can provide characterizations of $C_{0}$-universality for certain classes of kernels and spaces $\mathcal{M}$. In the case of translation-invariant kernels on $\mathcal{M}=\mathbb{R}^{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 $\mathbb{R}^{n}$, which includes: compactly-supported kernels, Gaussian kernels, Laplacian kernels… With the previous definition, we have the following sufficient condition:

Theorem 3.5**.**

Suppose $k^{pos}$ is a $C_{0}$-universal kernel on $\mathbb{R}^{n}$, $\gamma(1)>0$ and $\gamma(t)\neq\gamma(-t),\ \forall t\in[-1,1]$. Let $(X,T(\cdot))$ and $(Y,S(\cdot))$ be two oriented $\mathcal{H}^{d}$-rectifiable sets with $\mathcal{H}^{d}(X)$, $\mathcal{H}^{d}(Y)<\infty$. If $\left\|\mu_{X}-\mu_{Y}\right\|_{W^{\prime}}=0$, then $\mathcal{H}^{d}(X\bigtriangleup Y)=0$ and $T=S\ \mathcal{H}^{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 $\gamma(t)=\gamma(-t)$ for all $t$), generalizing the result of Charon \BBA Trouvé (\APACyear2013).

However, the previous proposition does not necessarily lead to a distance on the full space $\mathcal{V}_{d}$. Counter-examples in the case $d=1$ are discussed for example in Hsieh \BBA Charon (\APACyear2019). To recover a true distance on $\mathcal{V}_{d}$, one needs the previous map $\iota_{W}^{*}$ or equivalently the map

[TABLE]

to be injective. As follows from Theorem 6 in Sriperumbudur \BOthers. (\APACyear2011), this is in fact guaranteed when the kernel $k$ on the product space $\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$ is $C_{0}$-universal, specifically

Theorem 3.6**.**

The pseudo-distance $d_{W^{*}}$ induces a distance between signed measures of $\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$ if and only if $k$ is $C_{0}$-universal on $\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$. In particular, a sufficient condition for $d_{W^{*}}$ to be a distance on $\mathcal{V}_{d}$ is that $k^{pos}$ and $k^{G}$ are $C_{0}$-universal kernels on $\mathbb{R}^{n}$ and $\widetilde{G}^{n}_{d}$ respectively.

Note that these conditions are more restrictive than in Theorem 3.5. To our knowledge, there is no simple characterization for general $C_{0}$-universal kernels on the Grassmannian. However, within the setting of Remark 3.3, 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 $\Lambda^{d}(\mathbb{R}^{n})$.

3.2.3. Comparison with classical metrics

We now study more precisely the topology induced by the (pseudo) distance $d_{W^{*}}$ on $\mathcal{V}_{d}$ in comparison with the ones defined in Section 3.1. First of all, we observe that, for any $\omega\in W$ with $\|\omega\|_{W}\leq 1$, one must have $\|\omega\|_{\infty}\leq c_{W}$, where $c_{W}$ is the embedding constant of Proposition 3.2. Thus, for any $\mu$ and $\mu^{\prime}$ in $\mathcal{V}_{d}$, we have

[TABLE]

From the above inequalities we see that convergence in $d_{op}$ implies convergence in $d_{W^{*}}$.

Remark 3.7*.*

With more assumptions on the regularity of the kernel $k$, namely if $W$ is continuously embedded in $C_{0}^{1}(\mathbb{R}^{d}\times\widetilde{G}^{n}_{d})$, following a similar reasoning as above, one obtains the bound $\|\mu-\mu^{\prime}\|_{W^{*}}\leq c_{W}d_{BL}(\mu,\mu^{\prime})$.

Suppose $\mu_{i}$ converges to $\mu$ in narrow topology. Since the map $(\nu_{1},\nu_{2})\mapsto\nu_{1}\otimes\nu_{2}$ is continuous with respect to the narrow topology, we have

[TABLE]

as $i\rightarrow\infty$. Also, it’s clear that $\lim_{i\rightarrow\infty}\langle\mu_{i},\mu\rangle_{W^{*}}\rightarrow\|\mu\|_{W^{*}}^{2}$ and hence $\mu_{i}\rightarrow\mu$ with respect to $d_{W^{*}}$. To summarize the discussion above:

Proposition 3.8**.**

Let $\{\mu_{i}\}_{i}$ and $\mu$ be varifolds in $\mathcal{V}_{d}$ and assume that $\mu_{i}\rightarrow\mu$ with respect to the operator norm or the narrow topology, then $\mu_{i}\rightarrow\mu$ in $W^{*}$.

Remark 3.9*.*

We emphasize that the result of Proposition 3.8 only requires the assumptions of Proposition 3.2 and thus holds whether $\iota$ is injective or not.

As for the weak-* topology, with the $C_{0}$-universality assumption of Theorem 3.6 and restricting to varifolds with bounded total mass, we show that $d_{W^{*}}$ induces a topology stronger than weak-* convergence:

Proposition 3.10**.**

If $k$ is $C_{0}$-universal, then the topology induced by $d_{W^{*}}$ is finer than the weak- topology on $\mathcal{V}_{d,M}\doteq\{\mu\in\mathcal{V}_{d}\ \text{s.t}\ |\mu|(\mathbb{R}^{n})\leq M\}$ for any fixed $M>0$.*

Proof.

Let $\{\mu_{i}\}_{i}$ and $\mu$ be varifolds in $\mathcal{V}_{d,M}$ and assume that $\lim_{i\rightarrow\infty}d_{W^{*}}(\mu_{i},\mu)=0$. For any $f\in C_{0}(\mathbb{R}^{d}\times\widetilde{G}^{n}_{d})$ and $\varepsilon>0$, there exists a $g\in W$ such that $\|g-f\|<\varepsilon/2M$. Then we obtain that $\mu_{i}\overset{\ast}{\rightharpoonup}\mu$ from the following inequalities:

[TABLE]

∎

Note that the topology induced by $d_{W^{*}}$ may be strictly finer on $\mathcal{V}_{d,M}$. Indeed, if $\rho(0),\gamma(1)>0$, consider $\mu_{i}=\delta_{(x_{i},S)}$, where $\lim_{i\rightarrow\infty}|x_{i}|=\infty$ and $S\in\widetilde{G}^{n}_{d}$ fixed. Then $\mu_{i}\overset{\ast}{\rightharpoonup}0$ while $\|\mu_{i}\|_{W^{*}}^{2}=\rho(0)\gamma(1)>0$ for all $i$. Yet, by combining Propositions 3.1, 3.8 and 3.10, we have the following

Corollary 3.11**.**

Let $M>0$ and $K\subset\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$ be a compact subset. If $k$ is $C_{0}$-universal, then $d_{W^{*}}$ metrizes the weak- convergence of varifolds on $\mathcal{V}_{d,M,K}\doteq\{\mu\in\mathcal{V}_{d}\ \text{s.t}\ |\mu|(\mathbb{R}^{n})\leq M,\ \text{supp}(\mu)\subset 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 5.

4. 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 $\mathcal{V}_{d}$.

4.1. Deformation models

In this section, we discuss different models for how varifolds can be transported by a diffeomorphism of $\mathbb{R}^{n}$, in other words what are possible group actions of the diffeomorphism group $\textrm{Diff}(\mathbb{R}^{n})$ on $\mathcal{V}_{d}$.

Let us start by considering the case of an oriented rectifiable subset $(X,T_{X})$. A diffeomorphism $\phi\in\textrm{Diff}(\mathbb{R}^{n})$ transports $(X,T_{X})$ as

[TABLE]

where the transported orientation map writes

[TABLE]

the above term being well-defined from (2). This suggests introducing the following pushforward action on $\mathcal{V}_{d}$, which is defined for all $\mu\in\mathcal{V}_{d}$ and $\phi\in\textrm{Diff}(\mathbb{R}^{n})$ by:

[TABLE]

in which $J_{T}\phi(x)$ denotes the determinant of the Jacobian of $\phi$ along $T$ (i.e. the change of d-volume induced by $\phi$ along $T$ at $x$) which is given by

[TABLE]

for $(e_{1},\ldots,e_{d})$ an orthonormal basis of $T$. One easily verifies that $(\phi,\mu)\mapsto\phi_{\#}\mu$ defines a group action which commutes with the action on oriented rectifiable sets, namely

Proposition 4.1**.**

For any oriented rectifiable set $(X,T_{X})$ and diffeomorphism $\phi\in\textrm{Diff}(\mathbb{R}^{n})$, $\phi_{\#}\mu_{X}=\mu_{\phi(X)}$.

This follows from the area formula for integrals over rectifiable sets, c.f. Simon (\APACyear1983) Chapter 2.

Remark 4.2*.*

This pushforward action also extends the diffeomorphic transport of measures with densities on $\mathbb{R}^{n}$. Indeed if $\mu=\theta(x).\mathcal{L}^{n}$ with $\theta$ a measurable density function on $\mathbb{R}^{n}$ and $\mathcal{L}^{n}$ the Lebesgue measure, we can extend $\mu$ to a n-varifold in $\mathcal{V}_{n}$ by taking a constant global orientation in $\widetilde{G}_{n}^{n}=\{\pm 1\}$. Then, for any orientation-preserving diffeomorphism $\phi$, (15) writes in this case: $\phi_{\#}\mu=|J\phi(x)|\theta(x).\mathcal{L}^{n}$ with $J\phi(x)$ is the full Jacobian determinant of $\phi$, leading to the usual action on densities $\phi\cdot\theta(x)=|J\phi(x)|\theta(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 $\mathcal{V}_{d}$. For instance, one can define another action by removing the above volume change term, taking instead

[TABLE]

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

[TABLE]

Although this action is not consistent with the action on rectifiable sets as in Proposition 4.1, 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 Hsieh \BBA Charon (\APACyear2019) 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 (15), although we expect the following derivations to adapt to other cases as well, which precise study is for now left as future work.

4.2. 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 $\mathcal{V}_{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 Beg \BOthers. (\APACyear2005); Younes (\APACyear2019). In this framework, deformations consist of diffeomorphisms generated by flowing time-dependent vector fields. Let $V$ be a fixed RKHS of vector fields on $\mathbb{R}^{n}$ and $L^{2}([0,1],V)$ be space of time dependent velocity fields $v$ such that for all $t\in[0,1]$, $v_{t}$ belongs to $V$. The flow map $t\mapsto\varphi_{t}^{v}$ is defined for all $t\in[0,1]$ by $\varphi_{0}^{v}=\text{id}$ and the ODE $\dot{\varphi}^{v}_{t}=v_{t}\circ\varphi^{v}_{t}$. If $V$ is continuously embedded in $C_{0}^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$, one can show that for all $t$, $\varphi_{t}^{v}$ is a $C^{1}$-diffeomorphism of $\mathbb{R}^{n}$. Moreover, on the subgroup $\text{Diff}_{V}=\{\varphi_{1}^{v}\ |\ v\in L^{2}([0,1],V)\}$ of $\text{Diff}(\mathbb{R}^{n})$, one can define the following right-invariant Riemannian metric:

[TABLE]

Let us now consider a source (or template) varifold $\mu_{0}\in\mathcal{V}_{d}$ as well as a target $\mu_{tar}\in\mathcal{V}_{d}$. With the above deformation model and metric, registering $\mu_{0}$ to $\mu_{tar}$ consists in finding a deformation $\phi$ that minimizes $d_{G_{V}}(\text{id},\phi)$ with the constraint that $\phi_{\#}\mu_{0}$ is close to $\mu_{1}$ in the sense of a kernel metric $\|\cdot\|_{W^{*}}$ defined in Section 3.2. This can be reformulated as the following optimal control problem:

[TABLE]

with $v$ being the control, $E$ the total cost and the state equation is given by $\mu(t)\doteq(\phi_{t}^{v})_{\#}\mu_{0}$ for the pushforward model. The first term in (16) is the regularization term that constrains the regularity of the estimated deformation paths. The second term measures the similarity between the deformed varifold $\mu(1)$ and the target varifold $\mu_{tar}$. $\lambda$ 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 (16) holds under the following assumptions:

Theorem 4.3**.**

If $V$ is continuously embedded in $C_{0}^{2}(\mathbb{R}^{n},\mathbb{R}^{n})$, $W$ is continuously embedded in $C_{0}^{1}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})$ and $\rm{supp}(\mu_{0})\subset K$, for some compact subset $K$ of $\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$, then there exists a global minimizer to the problem (16).

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.

Remark 4.4*.*

One can derive necessary and sufficient conditions on the kernels of $W$ and $V$ for the two embedding assumptions of Theorem 4.3 to hold (see for instance Theorem 2.11 in Micheli \BBA Glaunès (\APACyear2014)). In our context, in order to get $W\hookrightarrow C_{0}^{1}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})$ for instance, it is enough to assume that $\rho$ and $\gamma$ are $C^{2}$ functions such that all derivatives of $\rho$ up to order 2 vanish as $x\rightarrow+\infty$.

As an important note, the formulation of (16) extends registration of submanifolds or rectifiable subsets in the sense that if $\mu_{0}=\mu_{X_{0}}$ and $\mu_{tar}=\mu_{X_{tar}}$ for two oriented d-rectifiable subsets of $\mathbb{R}^{n}$ then (16) becomes equivalent, thanks to Proposition 4.1, to registering rectifiable subsets, i.e. to the problem

[TABLE]

with $X(t)=\varphi_{t}^{v}\cdot X_{0}$, which is the setting of many past works as for instance Glaunès \BOthers. (\APACyear2008); Charon \BBA Trouvé (\APACyear2013); Kaltenmark \BOthers. (\APACyear2017).

4.3. General optimality conditions

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

We follow the same setting as well as related works such as Sommer \BOthers. (\APACyear2013) by first rewriting the above problem as an optimal control problem on diffeomorphisms, i.e.

[TABLE]

with $g(\varphi_{1}^{v})\doteq\lambda\|(\varphi_{1}^{v})_{\#}\mu_{0}-\mu_{tar}\|_{W^{*}}^{2}$. The state variables are now given by the deformations $\varphi_{t}^{v}$ which we view as elements of the Banach space $\mathcal{B}\doteq\text{id}+C_{0}^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$. Let us denote, for $\phi\in\text{Diff}(\mathbb{R}^{n})$, $\xi_{\phi}:V\rightarrow C_{0}^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$ the mapping $v\mapsto v\circ\phi$. We then introduce the Hamiltonian functional $H:\ C_{0}^{1}(\mathbb{R}^{n},\mathbb{R}^{n})^{*}\times\mathcal{B}\times V\rightarrow\mathbb{R}$ defined by:

[TABLE]

where $p$ is the costate variable which is a vector distribution of $C_{0}^{1}(\mathbb{R}^{n},\mathbb{R}^{n})^{*}$ and $(p|v\circ\phi)$ denotes the duality bracket in $C_{0}^{1}(\mathbb{R}^{n},\mathbb{R}^{n})^{*}$. With the assumptions of Theorem 4.3, it follows from the maximum principle shown in Arguillère \BOthers. (\APACyear2015) that if $(v_{t},\varphi_{t}^{v})$ is a global minimum of the optimal control problem, there exists a path of costates $p\in H^{1}([0,1],C_{0}^{1}(\mathbb{R}^{n},\mathbb{R}^{n})^{*})$ such that the following equations hold:

[TABLE]

with the end time boundary conditions $p_{1}=-\partial_{\phi}g(\varphi_{1}^{v})$. From the last equation in (18), we can attempt to deduce the form of the optimal $v$. Introducing the Riesz isometry operator $\mathbb{K}_{V}:V^{*}\rightarrow V$ and its inverse $\mathbb{L}_{V}=\mathbf{K}_{V}^{-1}:V\rightarrow V^{*}$, we get:

[TABLE]

One additional consequence of (18) is the following conservation of momentum again proved in Arguillère \BOthers. (\APACyear2015): for all $u\in C_{0}^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$ and $t\in[0,1]$,

[TABLE]

Note that (18), (19) and (20) are generic to the LDDMM model and so far independent of the nature of the deformed objects and of the term $g(\varphi_{1}^{v})$ in the cost. This dependency is entirely encompassed by the boundary condition $p_{1}=-\partial_{\phi}g(\varphi_{1}^{v})$ which we may describe a little more precisely based on the following:

Proposition 4.5**.**

The end-time momentum $p_{1}$ is a vector distribution in $C_{0}^{1}(\mathbb{R}^{n},\mathbb{R}^{n})^{*}$ of the form

[TABLE]

where $\alpha:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, $\beta:\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}\rightarrow(\mathbb{R}^{n\times d})^{*}$ and $\gamma:\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}\rightarrow\mathbb{R}$ are continuous fields and for all $T\in\widetilde{G}^{n}_{d}$, $\text{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 6 in the discrete setting). This result extends in a way first variation formulas for varifolds proved in Charon \BBA Trouvé (\APACyear2013); Charlier \BOthers. (\APACyear2017) which considered variations of rectifiable varifolds resulting from variations of the underlying rectifiable sets. This corresponds to the special case in which $\mu_{0}=\mu_{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\mapsto\int_{\varphi_{1}^{v}(X_{0})}u(x)\cdot h(x)d\mathcal{H}^{d}$ in $C_{0}^{0}(\mathbb{R}^{n},\mathbb{R}^{n})^{*}$ with vectors $h(x)$ normal to $\varphi_{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 (20) gives that for all $t\in[0,1]$,

[TABLE]

Using the expression of $p_{1}$ in Proposition 4.5, and grouping all [math]-th and $1$-st order terms in the resulting expressions, we may write $p_{t}$ in the general form:

[TABLE]

where $\alpha_{t}:\mathbb{R}^{n}\times\widetilde{G}_{d}^{n}\rightarrow\mathbb{R}^{n}$ and $B_{t}:\mathbb{R}^{n}\times\widetilde{G}_{d}^{n}\rightarrow(\mathbb{R}^{n\times d})^{*}$ are continuous fields, with $\alpha_{1}(x,T)=\alpha(x)$ and $B_{1}(x,T)du|_{T}(x)=\beta(x,T)du|_{T}(x)+\gamma(x,T)\text{div}_{T}u(x)$. Furthermore, optimal vector fields satisfy $v_{t}=K_{V}\xi^{*}_{\varphi_{t}^{v}}p_{t}$ and we have

[TABLE]

Denoting $K_{V}:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n\times n}$ the reproducing kernel of $V$, the reproducing kernel property implies that for all $u\in V$ and $x,h\in\mathbb{R}^{n}$, $u(x)\cdot h=\langle K_{V}(x,\cdot)h,u\rangle_{V}$. Moreover, the similar property on the kernel first order derivatives Micheli \BBA Glaunès (\APACyear2014) gives that for any $h,h^{\prime}\in\mathbb{R}^{n}$,

[TABLE]

Then, we rewrite the linear maps $B_{t}$ as $B_{t}(x,T)H=\sum_{i=1}^{d}b_{t,i}(x,T)\cdot H_{i}$ for any $H=(H_{1},\ldots,H_{d})\in\mathbb{R}^{n\times d}$ and where $b_{i}(x,T)\in\mathbb{R}^{n}$ are the component vector fields of $B_{t}$. By the above and the linearity of $\mathbb{K}_{V}$, we obtain the following general expression for optimal vector fields

[TABLE]

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 $\alpha$ and $b_{i}$ at this point, it will be specified later in the discrete setting, see Section 6.2.

5. Approximations by discrete varifolds

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

5.1. Discrete approximations

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

[TABLE]

for some $N\geq 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=\bigcup_{i=1}^{N}X_{i}$ is a triangulated surface, with $X_{i}$ being the mesh triangles with specified orientations, one can write $\mu_{X}=\sum_{i=1}^{N}\mu_{X_{i}}$ and for each $i\in\{1,\cdots,N\}$ approximate $\mu_{X_{i}}$ by $r_{i}\delta_{(x_{i},T_{i})}$, where $x_{i}$ is the center of $X_{i}$, $T_{i}$ the oriented plane containing $X_{i}$ and $r_{i}=\mathcal{H}^{d}(X_{i})$. This leads to the approximation $\widetilde{\mu}_{X}:=\sum_{i=1}^{N}r_{i}\delta_{(x_{i},d_{i})}$. As proved in Kaltenmark \BOthers. (\APACyear2017), this approximation provides an acceptable error bound for $d_{W^{*}}$:

[TABLE]

The main interest of such discrete varifold approximations is that the expression of the metric (3.2.1) becomes particularly simple to compute numerically. Indeed, given two discrete varifolds $\mu=\sum_{i=1}^{N}r_{i}\delta_{(x_{i},S_{i})}$ and $\mu^{\prime}=\sum_{j=1}^{M}r^{\prime}_{j}\delta_{(x^{\prime}_{j},T^{\prime}_{j})}$, we have as a particular case of (3.2.1):

[TABLE]

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 $\mu$ with finite total weight, can it be approximated by discrete varifolds and will approximations converge as $N\rightarrow+\infty$? This is the problem known as quantization, which has been studied intensively in the case of probability measures over Euclidean spaces Graf \BBA Luschgy (\APACyear2007) or manifolds Kloeckner (\APACyear2012), under specific regularity assumptions on those measures. In the situation of varifolds, an interesting recent work on this question is Buet \BOthers. (\APACyear2018). 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.

Theorem 5.1**.**

Let

[TABLE]

be the (non-convex) cone of discrete varifolds with at most $N$ Diracs. For any oriented varifold $\mu\in\mathcal{V}_{d}$ with $|\mu|(\mathbb{R}^{n})<\infty$, there exists a sequence $\mu_{N}\in\mathcal{V}_{d}^{N}$such that $\lim_{N\rightarrow\infty}d_{BL}(\mu_{N},\mu)=0$. Moreover, if $\mu$ has compact support, then we can assume that for all $N$, $\rm{supp}(\mu_{N})\subset K$ for some compact set $K\subset\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$ and

[TABLE]

where $C$ is a constant that only depends on $n$, $d$ and $\rm{supp}(\mu)$.

Proof.

We first tackle the case of compactly supported $\mu$. Without loss of generality, we may also assume that $\mu$ is a probability measure. Let $D=n+d(n-d)$ and $B\subset\mathbb{R}^{n}$ be a closed ball that contains $\rm{supp}|\mu|$. For brevity, we write $M\doteq B\times\widetilde{G}^{n}_{d}$. Since we can view $M$ as a compact $D$-dimensional submanifold of $\mathbb{R}^{n}\times\Lambda^{d}(\mathbb{R}^{n})$ (using Plücker embedding), M is also regular of dimension $D$ (cf. Graf \BBA Luschgy (\APACyear2007)), i.e., $0<\mathcal{H}^{D}(M)<\infty$ and there exist $c,r_{0}>0$, such that

[TABLE]

Given $\varepsilon\in(0,5r_{0})$, by the 5-Times Covering Lemma (cf. Simon (\APACyear1983)), these exists a subset $\mathcal{I}\subset M$, such that $M\subset\cup_{x\in\mathcal{I}}B_{\varepsilon}(x)$ and $B_{\varepsilon/5}(x)\cap B_{\varepsilon/5}(y)=\varnothing$ for all $x\neq y\in\mathcal{I}$. Therefore,

[TABLE]

We can thus obtain a partition $\{A_{i}\}_{i=1,\cdots,|\mathcal{I}|}$ of $M$ from the the collection $\{B_{\varepsilon}(x)\cap M\}_{x\in\mathcal{I}}$ which satisfies $\sup_{i}\textrm{diam}(A_{i})<\varepsilon$ and

[TABLE]

Let $r_{i}=\mu(A_{i})$ and $(x_{i},T_{i})\in A_{i}$ and define $\nu=\sum_{i=1}^{|\mathcal{I}|}r_{i}\delta_{(x_{i},T_{i})}$. For any $\varphi\in\textrm{Lip}_{1}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})$, with $\|\varphi\|_{\infty}\leq 1$, we have

[TABLE]

Taking the supremum over all $\omega\in\textrm{Lip}_{1}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})$ with $\|\varphi\|_{\infty}\leq 1$, we obtain $d_{BL}(\mu,\nu)<\varepsilon$. Then for each $N\in\mathbb{N}$, we can choose $\varepsilon_{N}=5(C\mathcal{H}^{D}(M)/N)^{1/D}$ and we obtain $\mu_{N}\in\mathcal{V}_{d}^{N}$ such that

[TABLE]

and in particular $\lim_{N\rightarrow+\infty}d_{BL}(\mu,\mu_{N})=0$.

Suppose now that $\rm{supp}(\mu)$ is not compact: we show that for any $\varepsilon>0$, there exists a discrete varifold $\nu$ such that $d_{BL}(\mu,\nu)<\varepsilon$. Choose a compact set $K\subset\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}$ such that $\mu(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d}\setminus K)<\varepsilon/2$. From the previous case, we can find a discrete varifold $\nu$ such that $d_{BL}(\mu\mathbin{\vrule height=6.88889pt,depth=0.0pt,width=0.55974pt\vrule height=0.55974pt,depth=0.0pt,width=5.59721pt}K,\nu)<\varepsilon/2$, and hence $d_{BL}(\mu,\nu)<\varepsilon$. ∎

Note that the proposition clearly holds for non-oriented varifolds as well. Another direct consequence, thanks to proposition 3.8 and remark 3.7, is the following corresponding statement for $d_{W^{*}}$:

Corollary 5.2**.**

With the assumptions from proposition 3.2, one also has $\lim_{N\rightarrow\infty}d_{W^{*}}(\mu_{N},\mu)=0$. If in addition $W\hookrightarrow C_{0}^{1}(\mathbb{R}^{n}\times\tilde{G}_{d}^{n})$, an equivalent upper bound as in Theorem 5.1 holds for $d_{W^{*}}(\mu,\mu_{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 $\mu$. We expect much better convergence properties for certain specific classes of varifolds, for instance assuming Alfors regularity as in Kloeckner (\APACyear2012), although we leave such questions for future investigation.

5.2. 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 $\mu\in\mathcal{V}_{d}$ and $N\in\mathbb{N}$, a natural idea is to look for the optimal discrete varifold in $\mathcal{V}_{d}^{N}$ that approximates $\mu$ in terms of the metric $d_{W^{*}}$. Due to the intricate structure of the set $\mathcal{V}_{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 $\mathcal{V}_{d}$ (e.g. generated by finite sets of Diracs). Works such as Durrleman \BOthers. (\APACyear2009); Gori \BOthers. (\APACyear2016) for instance, which are focused on the model of currents, consider dictionaries of Diracs defined on a predefined grid of point positions in $\mathbb{R}^{n}$. Then the problem can be recast as the one of finding sparse approximations of $\mu$ 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 Durrleman \BOthers. (\APACyear2009) 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 $\mathbb{R}^{n}\times\widetilde{G}_{d}(\mathbb{R}^{n})$. Such an increase in dimension makes the approach numerically impractical as soon as $n\geq 3$ and $d\geq 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 $\mathcal{V}_{d}^{N}$ for the class of kernel metrics $d_{W^{*}}$. It shares some connection with the recent work of Chauffert \BOthers. (\APACyear2017) that considers a related problem for standard measures defined on the torus $\mathbb{R}^{n}/\mathbb{Z}^{n}$.

Fix a varifold $\mu_{*}\in\mathcal{V}_{d}$. For any $N\in\mathbb{N},\ N\geq 1$, we seek $\mu_{N}\in\mathcal{V}_{d}^{N}$ that is closest to $\mu_{*}$ for $d_{W^{*}}$, namely

[TABLE]

By construction, if $|\mu|(\mathbb{R}^{n})<\infty$ then Corollary 5.2 will imply that $(\mu_{N})$ converges to $\mu$ 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:

Proposition 5.3**.**

Suppose all assumptions in proposition 3.2 and remark 3.3 hold. We further assume that the functions $\rho$ and $\gamma$ defining the kernels are non-negative. Then for any $\mu\in\mathcal{V}_{d}$ and $N\in\mathbb{N}$, there exists $\mu_{N}\in\mathcal{V}^{N}_{d}$ such that $\mu_{N}=\operatorname*{argmin}_{\nu\in\mathcal{V}^{N}_{d}}\|\mu-\mu_{*}\|_{W^{*}}$

Proof.

Let $\mu_{m}=\sum_{i=1}^{N}r_{i}^{m}\delta_{(x_{i}^{m},T_{i}^{m})}$ be a minimizing sequence, $K_{W}:W^{*}\mapsto W$ be the dual operator and $f:=K_{W}(\mu)$. Without loss of generality, we may assume that there is a $N_{1}\leq N$ such that $\sup_{1\leq i\leq N_{1}}|x_{i}^{m}|$ remains bounded and $\inf_{M+1\leq i\leq N}|x_{i}^{m}|$ tends to $\infty$ as $m\rightarrow\infty$.

Observe that $\sup_{1\leq i\leq N}\{r_{i}^{m}\}$ must be bounded. If it’s not bounded, then from the assumptions that $\rho,\gamma\geq 0$, we obtain

[TABLE]

as $m\rightarrow\infty$, which is absurd.

Since $r_{i}^{m},$ $T_{i}^{m}$, $f(x_{i}^{m},T_{i}^{m})$ and

[TABLE]

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

[TABLE]

Since $\rho\in C_{0}(\mathbb{R})$, the matrix $A$ must has the following form:

[TABLE]

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\in C_{0}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})$, we obtain

[TABLE]

where $\boldsymbol{r^{\prime}}=(r_{1},\cdots,r_{N_{1}})$, $\boldsymbol{r^{\prime\prime}}=(r_{N_{1}+1},\cdots,r_{N})$, and $\boldsymbol{f^{\prime}}=(f_{1},\cdots,f_{N_{1}})$. Since $\sup_{1\leq i\leq N_{1}}|x_{i}^{m}|$ is bounded we can assume that $\lim_{m\rightarrow\infty}x_{i}^{m}=x_{i},\ 1\leq i\leq N_{1}$. Let $\mu:=\sum_{i=1}^{N_{1}}r_{i}\delta_{(x_{i},u_{i})}$, then

[TABLE]

Hence $\mu$ 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 $\rho$ and $\gamma$ take negative values. It is so far an open question to determine to what extent one could generalize the result of Proposition 5.3, one particular but important case being the one of current metrics obtained for $\gamma(t)=t$ which is not covered by our result.

As written in the proof of Proposition 5.3, (24) is equivalent to the optimization problem:

[TABLE]

Any solution must satisfy first order optimality conditions obtained by differentiating $\|\mu_{N}-\mu\|_{W^{*}}$ with respect to the $(r_{k},x_{k},T_{k})$. In particular, we have

[TABLE]

which gives $\langle\mu_{N}-\mu_{*},\mu_{N}\rangle_{W^{*}}=0$. It shows that for any $N\in\mathbb{N}$, $\|\mu_{N}\|_{W^{*}}\leq\|\mu_{*}\|_{W^{*}}$.

5.3. $\Gamma$-convergence of registration functionals

Ultimately, our purpose is to use the previous approximating discrete varifolds $\mu_{N}$ to approximate the diffeomorphic registration problem (16). The natural question that arises is whether replacing the source varifold $\mu_{0}$ by its projections $\mu_{N}$ in (16) 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 $\Gamma$-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 Arguillère \BOthers. (\APACyear2016).

To obtain such convergence results for solutions of variational problems, one usually requires the approximating sequence to possess certain nice properties. Specifically, assuming $\mu\in\mathcal{V}_{d}$ with compact support and finite mass and $\{\mu_{N}\}\subset\mathcal{V}^{N}_{d}$ such that $\lim_{N\rightarrow\infty}\|\mu_{N}-\mu\|_{W^{*}}=0$, we will need that $\bigcup_{N}\textrm{supp}(\mu_{N})\subset K$ for some compact set $K\subset\mathbb{R}^{n}$ or that $\sup_{N}|\mu_{N}|(\mathbb{R}^{n})<\infty$. 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 $\mu_{N}$.

Yet, provided that $\bigcup_{N}\textrm{supp}(\mu_{N})\subset K$, we can actually retrieve the boundedness of the total mass. We assume in what follows that the kernels are such that $\rho(0)>0$ and $\gamma(1)>0$.

Lemma 5.4**.**

Let $\{\mu_{N}\}$ be a sequence of discrete varifolds with finite mass such that there exists a compact $K\subset\mathbb{R}^{n}$ with $\rm{supp}(|\mu_{N}|)\subset K$ for all $N$. We assume that $\{\|\mu_{N}\|_{W^{*}}\}$ is bounded. Then $\{|\mu_{N}|(\mathbb{R}^{n})\}$ is bounded.

Proof.

We prove it by contradiction. Assume that $(|\mu_{N}|(\mathbb{R}^{n}))_{N\geq 1}$ is unbounded. Then, up to extracting a subsequence, we can assume that $|\mu_{N}|(\mathbb{R}^{n})\rightarrow+\infty$. Let’s write $\mu_{N}=\sum_{i=1}^{p_{N}}r_{i,N}\delta_{(x_{i,N},T_{i,N})}$. Thus $|\mu_{N}|(\mathbb{R}^{n})=\sum_{i=1}^{p_{N}}r_{i,N}\rightarrow+\infty$. Since $\rho$ and $\gamma$ are continuous and $\rho(0),\gamma(1)>0$, we can find compact subsets $A\subset K$ and $B\subset\tilde{G}^{n}_{d}$ with diameters small enough, so that: $\inf_{x,y\in A}\rho(|x-y|^{2})>m>0$, $\inf_{u,v\in B}\gamma(\langle u,v\rangle)>m^{\prime}>0$ and $\lim_{N\rightarrow\infty}\sum_{i\in\mathcal{I}_{N}}r_{i,N}=\infty$, where $\mathcal{I}_{N}:=\{i:(x_{i,N},u_{i,N})\in A\times B\}$. It follows that, as $N\rightarrow\infty$,

[TABLE]

which is a contradiction. ∎

Lemma 5.4 suggests that one should enforce the uniform compactness of the supports of the $\mu_{N}$. To do so in the context of the projection approach of the previous sections, we consider solving the optimization problem (24) with the additional constraint that the support of $\mu_{N}$ stays in a compact set containing $\rm{supp}(|\mu|)$. We still have to verify the convergence of the resulting sequence:

Proposition 5.5**.**

Let $\mu_{0}$ be a varifold with finite mass and $K$ be a compact set in $\mathbb{R}^{n}$ which contains $\rm{supp}(|\mu_{0}|)$. Construct the approximating sequence of $\mu_{0}$ by solving the following constrained optimization problem:

[TABLE]

Then $\mu_{K,N}$ converges to $\mu_{0}$ in $d_{W^{*}}$ and, if the kernel $k$ is $C_{0}$-universal, it also converges in $d_{BL}$.

Proof.

Thanks to Theorem 5.1 and Lemma 5.4, we immediately get that $\|\mu_{K,N}-\mu_{0}\|_{W^{*}}\rightarrow 0$ as $N\rightarrow\infty$ and $\sup_{N}(|\mu_{K,N}|(\mathbb{R}^{n}))<\infty$. Moreover, if $k$ is $C_{0}$-universal, then by Proposition 3.10 it implies that $\mu_{K,N}\overset{\ast}{\rightharpoonup}\mu_{0}$. Since $\bigcup_{N}\rm{supp}(|\mu_{K,N}|)\subset K$ and $\sup_{N}(|\mu_{K,N}|(\mathbb{R}^{n}))<\infty$, weak-* convergence implies that $\mu_{K,N}$ converges to $\mu_{0}$ in $d_{BL}$ by Proposition 3.1. ∎

We are now able to state the main result of this section. We assume that the source/template varifold $\mu_{0}$ is compactly supported and we fix $K$ is a compact subset of $\mathbb{R}^{n}$ that contains $\rm{supp}(|\mu_{0}|)$. Then for any $N\in\mathbb{N},\ N\geq 1$, $\mu_{K,N}$ is defined as in Proposition 5.5 and we introduce the following energy functionals $E_{N}:\ L^{2}([0,1],V)\rightarrow\mathbb{R}_{+}$:

[TABLE]

which are the equivalent to the energy $E$ of the original problem (16) but replacing the template varifold $\mu_{0}$ by its approximations $\mu_{K,N}$.

Theorem 5.6**.**

With the above notations, we assume that the reproducing kernel $k$ of $W$ is $C_{0}$-universal and satisfies all the conditions of Proposition 5.3. We also assume the continuous embedding $V\hookrightarrow C_{0}^{2}(\mathbb{R}^{n},\mathbb{R}^{n})$. Then, the sequence of functionals $E_{N}$ $\Gamma$-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\geq 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$.

Proof.

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

[TABLE]

Since $v\mapsto\int_{0}^{1}\|v\|_{V}^{2}dt$ is lower semicontinuous with respect to the weak topology, we only need to prove the following,

[TABLE]

For all $\omega\in W$ with $\|\omega\|_{W}\leq 1$, we have

[TABLE]

where $g\in C_{c}(\mathbb{R}^{n}\times\tilde{G}^{n}_{d})$ and $\sup_{N\geq 1}J_{S}\varphi_{1}^{v^{N}}\leq g(x,S)$, for all $(x,S)\in K\times\tilde{G}^{n}_{d}$. Similar to the computation done in the proof of Theorem 4.3, we see that

[TABLE]

Taking supremum over all $\omega\in W$ with $\|\omega\|_{W}\leq 1$, we obtain the following inequality,

[TABLE]

From Proposition 5.5, $\mu_{K,N}$ converges to $\mu_{0}$ in the narrow topology. Hence the right hand side in the equation above converges to [math] as $N\rightarrow\infty$. This proves (29).

Second, we need to show that for each $\bar{v}\in L^{2}([0,1],V)$, there exists a sequence $v^{N}$ converging to $\bar{v}$ weakly such that

[TABLE]

In fact, it suffices here to take $v^{N}$ to be the constant sequence $v^{N}=\bar{v}$ since, by a similar argument to the proof of (29), it leads to

[TABLE]

and thus implies that

[TABLE]

∎

Note that we stated the result of Theorem 5.6 in the situation where only the source varifold $\mu_{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.

6. 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 6.1, is to define an adequate framework for representing and computing with elements of the oriented Grassmannian.

6.1. Frame representation for metric computation and quantization

In order to come up with a computationally effective representation of $\widetilde{G}_{d}(\mathbb{R}^{n})$ and by extension of discrete oriented varifolds, we consider a slightly different setting than the Plücker embedding idea of Remark 2.2, primarily because the dimension of the embedding vector space $\Lambda^{d}(\mathbb{R}^{n})$ may become prohibitively large in practice. We may instead choose to represent an element $T\in\widetilde{G}_{d}(\mathbb{R}^{n})$ by an oriented frame $(u^{(1)},\ldots,u^{(d)})\in\mathbb{R}^{n\times d}$ of independent vectors for which $T=\text{Span}(u^{(1)},\ldots,u^{(d)})$. Such a representation is of course not unique since elements of $\widetilde{G}_{d}(\mathbb{R}^{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 $\mu$ be a discrete varifold of the form $\mu=\sum_{i=1}^{N}r_{i}\delta_{(x_{i},T_{i})}$. For each $i$, we consider a frame $\{u_{i}^{(1)},\cdots,u_{i}^{(d)}\}$ such that

[TABLE]

In other words, the oriented space spanned by the frame $\{u_{i}^{(1)},\cdots,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 $\mu$ with the (non-unique) state variable $q=(x_{i},u^{(1)}_{i},\cdots,u^{(d)}_{i})_{i=1,\cdots,N}$ in the vector space $\mathbb{R}^{Nn(d+1)}$. Conversely, such a frame $q$ with $(u^{(1)}_{i},\cdots,u^{(d)}_{i})$ a matrix of rank $k$ for all $i$, corresponds to the (unique) discrete oriented varifold defined by the relations of (31); we will denote it by $\mu^{q}$ in what follows.

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

[TABLE]

where $r_{i}=|u_{i}^{(1)}\wedge\cdots\wedge u_{i}^{(d)}|=\sqrt{\det(u_{i}^{(k)}\cdot u_{i}^{(l)})_{k,l}}$. Note that this expression does not depend on the choice of frames that satisfy the conditions of (31) for $\mu$ (and similarly for $\mu^{\prime}$). In the case where $\mu^{\prime}$ is a more general non-discrete varifold in $\mathcal{V}_{d}$, the computation of $\langle\mu,\mu^{\prime}\rangle_{W^{*}}^{2}$ involves integrals over $\mathbb{R}^{n}\times\widetilde{G}_{d}(\mathbb{R}^{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 (32) in which the $x_{j}^{\prime},u^{\prime}_{j}$ and $r^{\prime}_{j}$ are now the quadrature nodes and associated weights of the scheme.

In this setting, the solution to the projection problem (24) can be computed by an iterative descent strategy on the vector $q=(x_{i},u^{(1)}_{i},\cdots,u^{(d)}_{i})_{i=1,\cdots,N}$. The gradient of $q\mapsto\|\mu^{q}-\mu_{*}\|_{W^{*}}^{2}$ can be computed by direct differentiation of expressions like (32) 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 Charlier \BOthers. (\APACyear2018) 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 Overton (\APACyear2016) which we typically initialize by taking a random subset of $N$ Diracs composing the varifold $\mu_{*}$. 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 (24). Results of this algorithm are discussed below in Section 7.2.

6.2. Discrete registration model

This frame representation also provides a convenient setting to express the diffeomorphism action and registration problem on discrete varifolds. Indeed, let $\varphi$ be a diffeomorphism of $\mathbb{R}^{n}$ and $\mu\in\mathcal{V}_{d}^{N}$, the pushforward action $\varphi_{\#}\mu$ in (15) is equivalent to the following action in the frame model:

[TABLE]

Now, this allows us to rewrite the former infinite-dimensional optimal control problem by considering instead the finite-dimensional state variable $q\in\mathbb{R}^{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 4.3. Note that the resulting Hamiltonian equations we obtain are eventually very similar to the ones appearing in the 1st-order jets model studied in Sommer \BOthers. (\APACyear2013); Jacobs (\APACyear2013), although there are a few notable differences due to the specific extra invariances attached to the varifold framework (c.f. Hsieh \BBA Charon (\APACyear2019) 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:

[TABLE]

with $p^{x},p^{u_{k}}\in\mathbb{R}^{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:

[TABLE]

while optimal vector fields $v$ satisfy

[TABLE]

Plugging the above expression of $v$ with respect to $(p,q)$ in the Hamiltonian (33) gives the reduced Hamiltonian $H_{r}(p,q)\doteq H(p,q,v)$ which writes:

[TABLE]

and (38) 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))=\frac{1}{2}\|v_{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:

Lemma 6.1**.**

For any $i=1,\ldots,N$, the matrix

[TABLE]

is constant in time.

Proof.

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

[TABLE]

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:

[TABLE]

As a clear consequence of (31), $\mu^{q}$ and thus $g(q)$ are independent of the choices of the frame vectors $(u_{i}^{(k)})_{k=1,\dots,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 6.1, of the full path $p(t)$. These are summed up by the following result:

Proposition 6.2**.**

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)\perp\text{Span}(\{u_{i}^{(\ell)}(t)\}_{\ell\neq k})$ for all $t\in[0,1]$, $i=1,\ldots,N$ and $k=1,\ldots,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 $\mathbb{R}^{Nn(d+1)}$ of lower dimensions $N(n+d(n-d)+1)$, which is precisely the dimension of the ’true’ state space $(\mathbb{R}^{n}\times\widetilde{G}_{d}(\mathbb{R}^{n})\times\mathbb{R})^{N}$.

6.3. 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 (38) and (39), 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 Vialard \BOthers. (\APACyear2012); Sommer \BOthers. (\APACyear2013); Charon \BBA Trouvé (\APACyear2013). 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 (6.2). 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. Vialard \BOthers. (\APACyear2012) or Arguillère \BOthers. (\APACyear2015), the gradient of $g(q(1))$ with respect to $p(0)$ can be computed by flowing backward in time the adjoint Hamiltonian system

[TABLE]

where $F(q,p)=(\partial_{p}H_{r}(p,q),-\partial_{q}H_{r}(p,q))$, $Z(t)=(\tilde{q}(t),\tilde{p}(t))\in\mathbb{R}^{n}\times(\mathbb{R}^{n})^{d}$ the adjoint variables of the system, together with the end-time conditions $\tilde{q}(1)=-\partial_{q}g(q)|_{q=q(1)}$ and $\tilde{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 (42) can be approximated efficiently based on the finite difference trick proposed in Arguillère \BOthers. (\APACyear2015) (Section 4.1). Indeed, it can be rewritten as follows:

[TABLE]

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

[TABLE]

for some small $\epsilon>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:

1:repeat

2:     From $(q(0),p(0))$ compute $(q(t),p(t))$ by forward integration of the reduced Hamiltonian system given by (38) and (39).

3:     Compute $g(q(1))$ and $-\partial_{q}g(q)|_{q=q(1)}$.

4:     Integrate backward the adjoint Hamiltonian equations (42) to obtain $\partial_{p(0)}g(q(1))$.

5:     Deduce the gradient of the full cost function with respect to $p(0)$.

6:     Update $p(0)$.

7:until convergence

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 Overton (\APACyear2016). One can further take additional advantage of the dimensionality reduction provided by Proposition 6.2 by restricting each of the components $p_{i}^{u_{k}}(0)$ to the linear subspace $\text{Span}(\{u_{i}^{(\ell)}(0)\}_{\ell\neq k})^{\perp}$. Lastly, as in Section 6.1, all kernel summation and differentiation operations appearing in both the varifold fidelity terms and Hamiltonian equations are coded in CUDA using the KeOps library Charlier \BOthers. (\APACyear2018). 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.

7. 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(-\frac{|x-y|^{2}}{\sigma_{V}^{2}})Id$. The kernel function $\rho$ is a Gaussian of scale $\sigma_{\rho}$. 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 Bruveris \BOthers. (\APACyear2012); Kaltenmark \BOthers. (\APACyear2017); Hsieh \BBA Charon (\APACyear2019). The function $\gamma$ is chosen, depending on the situation, in the different classes of functions discussed in detail in Kaltenmark \BOthers. (\APACyear2017), the main distinction being whether the considered varifolds are rectifiable or not according to the conditions given by Theorem 3.5 and Theorem 3.6 and whether the shapes carry a relevant orientation or not. In particular, one can use $\gamma(t)=t^{2}$ to recover an orientation invariant fidelity metric, or $\gamma(t)=e^{-\frac{2}{\sigma_{\gamma}^{2}}(1-t)}$ which leads to an orientation sensitive distance that satisfy the conditions of Theorem 3.6. All simulations are run on a desktop computer equipped with a NVIDIA Quadro P5000 graphics card.

7.1. Diffeomorphic registration

We start with results of registration obtained from the algorithm of Section 6.3. In this section, we will mostly focus on examples involving 2-varifolds, the reader may refer to Hsieh \BBA Charon (\APACyear2019) 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 Charon \BBA Trouvé (\APACyear2013); Kaltenmark \BOthers. (\APACyear2017) 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 4.2. Shown in Fig. 1 are triangulated surfaces of amygdala segmented from two different subjects of the BIOCARD database Miller, Ratnanather\BCBL \BOthers. (\APACyear2015), containing 563 vertices, 1122 triangles and 488 vertices, 972 triangles respectively. Following the simple procedure outlined at the beginning of Section 5.1, 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. 1. 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. 2, we consider a more challenging registration scenario which was originally studied in Ardekani \BOthers. (\APACyear2011). 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 Ardekani \BOthers. (\APACyear2011) 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. 2 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. 3, 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 ($\sigma=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 $\mathbb{R}^{3}$ (i.e. 0-varifolds) and follow the simple point distribution LDDMM algorithm for unlabelled point sets proposed in Glaunès \BOthers. (\APACyear2004). The result shown on the third row of Fig. 3 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 Allard \BOthers. (\APACyear2012). 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. 3. 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. 3, 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).

7.2. Approximation and registration

In this second part, we examine some results of the varifold quantization procedure proposed in Section 5, and in particular its interplay with the registration algorithm. Specifically, we wish to numerically validate the statements of Corollary 5.2 and Theorem 5.6. We shall consider the following protocol. Starting from a highly sampled shape (that we treat as the groundtruth) for which the associated varifold $\mu_{0}$ is composed of a very high number of Diracs, we compute the compressed varifolds given by the $\mu_{N}$ of (24) 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 $\mu_{tar}$ from the source varifolds given by the $\mu_{N}$ in lieu of $\mu_{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.

[TABLE]

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. 4 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 $\|\mu_{N}-\mu_{0}\|_{W^{*}}/\|\mu_{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 $\mu_{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. 5 (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 $\mu_{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 $\mu_{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 $\|\mu_{N}-\mu_{0}\|_{W^{*}}/\|\mu_{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. 5. 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 $\mu_{0}$ with only a small fraction of the initial number of Diracs. Some of the quantized varifolds $\mu_{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. 5 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.

8. 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 $\mathbb{R}^{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 $\mathcal{V}_{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 7, 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 (4), or the case of contrast-invariant image registration c.f. Hsieh \BBA Charon (\APACyear2019).

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 5 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 $\text{Diff}(\mathbb{R}^{n})$ on $\mathcal{V}_{d}$. Yet, other group actions could be have been considered, as briefly evoked in Section 4.1, 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 3.5 We first prove that $\mathcal{H}^{d}(X\bigtriangleup 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 $\left\|\mu_{X}-\mu_{Y}\right\|_{W^{*}}=0$ and $\mathcal{H}^{d}(X\bigtriangleup Y)>0$. Without loss of generality, we may assume that $\mathcal{H}^{d}(X\setminus Y)>0$. From Lusin’s theorem, there exists a subset $U$ of $X$ such that $T|_{U}$ is continuous and $\mathcal{H}^{d}(X\setminus U)<\mathcal{H}^{d}(X\setminus Y)$. Let us denote by $E:=U\cap(X\setminus Y)$, we see that $\mathcal{H}^{d}(E)>0$. Since for $\mathcal{H}^{d}\ a.e.\ x\in E$,

[TABLE]

(cf Evans (\APACyear2018)), there exists $x_{0}\in E,\ \mathcal{H}^{d}(B_{r}(x_{0})\cap E)>0$ for any $r>0$.

Let $g:\widetilde{G}^{n}_{d}\rightarrow\mathbb{R}$ be defined by $g(\cdot)=k^{G}(T(x_{0}),\cdot)$. Since $x\longmapsto g(T(x))$ is continuous on $E$ and $g(T(x_{0}))>0$, there exists $r_{0}>0$ such that $\forall\ x\in B_{r_{0}}(x_{0})\cap E,\ g(T(x))>0$. Let $A\doteq B_{r_{0}}(x_{0})\cap E$ and $h(x):=\mathbf{1}_{A}(x)$, then $\mathcal{H}^{d}(A)>0$ and $g(T(x))>0,\ \forall\ x\in A$. Using the density of $C_{c}(\mathbb{R}^{n})$ in $L^{1}(\mathbb{R}^{n},\mathcal{H}^{d}\mathbin{\vrule height=6.88889pt,depth=0.0pt,width=0.55974pt\vrule height=0.55974pt,depth=0.0pt,width=5.59721pt}(X\cup Y))$ together with the fact that $k^{pos}$ is $C_{0}$-universal, there exist $\{f_{j}\}_{j=1}^{\infty}\subset C_{c}(\mathbb{R}^{n})$ and $\{h_{j}\}_{j=1}^{\infty}\subset W^{pos}$ such that $\lim\limits_{j\rightarrow\infty}f_{j}=h$ in $L^{1}(\mathbb{R}^{n},\mathcal{H}^{d}\mathbin{\vrule height=6.88889pt,depth=0.0pt,width=0.55974pt\vrule height=0.55974pt,depth=0.0pt,width=5.59721pt}(X\cup Y))$ and $\|f_{j}-h_{j}\|_{\infty}<\frac{1}{j}$. Now, since $h_{j}\otimes g\in W$ and $\mu_{X}=\mu_{Y}$ in $W^{*}$, we have

[TABLE]

which is a contradiction. Hence we have $\mathcal{H}^{d}(X\bigtriangleup Y)=0$

Next, we show that $T(x)=S(x)$ $\mathcal{H}^{d}$-$a.e.$. Let $F:=\{x\in X|T(x)=-S(x)\}$ and assume that $\mathcal{H}^{d}(F)>0$. From Lusin’s theorem, there exists subset $F^{\prime}\subset F$ such that $T|_{F^{\prime}}$ is continuous and $\mathcal{H}^{d}(F^{\prime})>0$. Using the upper density argument as above, we can find $z_{0}\in F^{\prime}$ such that $\mathcal{H}^{d}(B_{r}(z_{0})\cap F^{\prime})>0$ for all $r>0$. Since the map $x\mapsto\langle T(x),T(z_{0})\rangle$ restricted to $F^{\prime}$ is continuous, there exists a $\delta_{0}>0$ satisfying:

[TABLE]

Define $B:=B_{\delta_{0}}(z_{0})\cap F^{\prime}$, $\eta(\cdot):=\gamma(\langle\cdot,T(z_{0})\rangle)$ and $u(x):=\eta(T(x))-\eta(S(x))$. Observe that, from the assumption $\gamma(t)\neq\gamma(-t),\ \forall t\in[-1,1]$,

[TABLE]

From this, we may assume that $u(x)>0,\ \forall x\in F^{\prime}$. Let $\{f_{j}^{\prime}\}_{j}$ and $\{h_{j}^{\prime}\}_{j}$ be sequences in $C_{c}(\mathbb{R}^{n})$ and $W_{pos}$ such that $f_{j}^{\prime}$ converges to $\mathbf{1}_{B}$ in $L^{1}(\mathbb{R}^{n},\mathcal{H}^{d}\mathbin{\vrule height=6.88889pt,depth=0.0pt,width=0.55974pt\vrule height=0.55974pt,depth=0.0pt,width=5.59721pt}F)$ and $\|f_{j}^{\prime}-h_{j}^{\prime}\|_{\infty}<1/j$. We obtain

[TABLE]

which is impossible. ∎

Proof of Theorem 4.3 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 $\bar{v}$ in $L^{2}([0,1],V)$. Using the results of Younes (\APACyear2019) Chapter 7.2, we know that

[TABLE]

Furthermore, for any $\omega\in W$, we have

[TABLE]

Now, using the embedding $W\hookrightarrow C_{0}^{1}(\mathbb{R}^{n}\times\widetilde{G}^{n}_{d})$

[TABLE]

Taking supremum over all $\omega\in W$ with $\|\omega\|_{W}\leq 1$, we obtain that

[TABLE]

as $j\rightarrow\infty$. Combining this with lower semicontinuity of $v\mapsto\|v\|_{L^{2}([0,1],V)}^{2}$, we finally obtain that

[TABLE]

and hence $\bar{v}$ is a global minimizer. ∎

Proof of Proposition 4.5 Recall that for all $\phi\in\text{Diff}(\mathbb{R}^{n})$, $g(\phi)=\lambda\|\phi_{\#}\mu_{0}-\mu_{tar}\|_{W^{*}}^{2}$ which we may rewrite as

[TABLE]

Thus, the variation with respect to $\phi$ in the Banach space $\mathcal{B}$ writes

[TABLE]

where $\omega_{0}\doteq 2\lambda K_{W}(\phi_{\#}\mu_{0}-\mu_{tar})\in W$. Moreover

[TABLE]

Taking the variation with respect to $\phi$ along any $u\in C_{0}^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$, we obtain:

[TABLE]

where the last term follows from the differentiation of Gram determinant matrices while the notation $\partial_{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 Charon \BBA Trouvé (\APACyear2013) and to the developments in Section 6 for more details. For the first term, we can rely on the Young measure decomposition $\mu_{0}=|\mu_{0}|\otimes\nu_{x}$ introduced at the end of Section 2.1 which gives:

[TABLE]

We can also rewrite the third term as:

[TABLE]

As for the second term in (Appendix), for each $(x,T)$ the integrand involves a linear combination (depending on $\phi$) of the partial derivatives of $u$ along the subspace $T$ i.e. of the elements of the matrix $d_{x}u|_{T}\in\mathbb{R}^{n\times d}$. Thus, without attempting to specify this term explicitly, we can in general write it as $\tilde{\beta}(\phi,x,T)d_{x}u|_{T}$ where $\tilde{B}$ is a continuous map from $\mathcal{B}\times\mathbb{R}^{n}\times\tilde{G}_{d}^{n}$ into $\mathcal{L}(\mathbb{R}^{n\times d},\mathbb{R})$ giving us

[TABLE]

The result of the theorem then follows by setting $\alpha(x)\doteq\tilde{\alpha}(\varphi_{1}^{v},x)$, $\beta(x,T)\doteq\tilde{\beta}(\varphi_{1}^{v},x,T)$ and $\gamma(x,T)=\tilde{\gamma}(\varphi_{1}^{v},x,T)$. ∎

Proof of Proposition 6.2 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),\cdots,u^{(d)}(t))$, $p(t)=(p^{x}(t),p^{u_{1}}(t),\ldots,p^{u_{d}}(t))$ for the state and costate variables along an optimal trajectory and

[TABLE]

Consider the group of linear transformations, $G\doteq{\rm SL}(U)\oplus{\rm GL}(U^{\perp})$, i.e., for any $\mathrm{g}\in G$,

[TABLE]

where $x_{U}$ and $x_{U^{\perp}}$ are the orthogonal projections of $x$ on $U$ and $U^{\perp}$, with $\mathrm{g}_{\mathbin{\!/\mkern-5.0mu/\!}}\in{\rm SL}(U)$ and $\mathrm{g}_{\perp}\in{\rm GL}(U^{\perp})$. The Lie algebra of $G$ is $\mathfrak{g}=\mathfrak{sl}(U)\times\mathcal{L}(U^{\perp})$ and $\mathfrak{sl}(U)$ is the set of all zero trace linear transformations of $U$. Now, consider the action of $G$ on $\mathbb{R}^{(d+1)n}$ defined as:

[TABLE]

for any $q=(q_{0},\ldots,q_{d})\in\mathbb{R}^{(d+1)n}$. We see that $\mu^{\mathrm{g}\cdot q(1)}=\mu^{q(1)}$ for all $\mathrm{g}\in G$ and therefore $g(\mathrm{g}\cdot q(1))=g(q(1))$.

Now, if we let $\{\mathrm{g}_{t}\}$ be a smooth curve in $G$ that satisfies $\mathrm{g}_{0}=id$ and $\frac{d}{d\tau}|_{\tau=0}\mathrm{g}_{\tau}=h\in\mathfrak{g}$, differentiating the equality $g(\mathrm{g}_{\tau}\cdot q(1))=g(q(1))$ shows that for any $h\in\mathfrak{g}$, we have

[TABLE]

Since $h\in\mathfrak{g}$, we must have that $h|_{U}$ is a zero trace linear map. For any $1\leq i<j\leq d$, we may choose $h$ such that $h(u^{(i)}(1))=-h(u^{(j)}(1))$ and $h(u^{(k)}(1))=0,\ \forall k\notin\{i,j\}$, which leads to $\langle u^{(i)}(1),p^{u_{i}}(1)\rangle=\langle u^{(j)}(1),p^{u_{j}}(1)\rangle$. Consequently,

[TABLE]

for some constant $\alpha$. In addition, for any $i\neq j$, we can also choose $h$ such that $h(u^{(i)}(1))=u^{(j)}(1)$ and $h(u^{(k)}(1))=0,\ \forall k\notin\{i,j\}$, which gives $\langle u^{(i)}(1),p^{u_{j}}(1)\rangle=0$. It results that $D(1)=\alpha.I_{d\times d}$.

Finally, since $D(t)$ is constant by Lemma 6.1, we obtain that

[TABLE]

for all $t\in[0,1]$. ∎