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…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 20
Figure 21
Figure 22
Figure 23
Figure 24
Figure 25
Figure 26
Figure 27
Figure 28
Figure 29
Figure 30
Figure 31
Figure 32
Figure 33
Figure 34
Figure 35
Figure 36
Figure 37
Figure 38
Figure 39
Figure 40Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Metrics, quantization and registration in varifold spaces
Hsi-Wei Hsieh
Department of Applied Mathematics and Statistics, Johns Hopkins University
and
Nicolas Charon
Department of Applied Mathematics and Statistics, Johns Hopkins University
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 -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.
Mathematics Subject Classification (2010) 49Q20 49M25 58E50 68T10
Key words and phrases:
varifolds, diffeomorphic registration, reproducing kernels, quantization, optimal control, -convergence
.
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 -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 - and - 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 by incorporating an additional tangent space component. In this work, we will consider such spaces to be oriented. Thus, for a given dimension , we first need to introduce the set of all possible -dimensional oriented tangent spaces in :
Definition 2.1**.**
The -dimensional oriented Grassmannian is the set of all oriented -dimensional linear subspaces of .
The oriented Grassmannian is a compact manifold of dimension which can be identified to the quotient . It is also a double cover of the (non-oriented) Grassmannian of -dimensional subspaces of . For practical purposes, a more convenient representation of is the one detailed in the following remark.
Remark 2.2*.*
Given , there exists a basis of such that has consistent orientation with . Then the following map, called the oriented Plücker embedding, is well defined and injective,
[TABLE]
This allows to identify as a subset of the unit sphere of which inherits the topology of the inner product on . We remind that this inner product is defined for any , in by the determinant of the Gram matrix:
[TABLE]
Through this identification, one can also define the action of linear transformations on as follows
[TABLE]
for any and a linear invertible map.
Similar to the definition of classical varifolds in Simon (\APACyear1983), we define oriented varifolds as measures on .
Definition 2.3**.**
An oriented -varifold on is a nonnegative finite Radon measure on the space . Its weight measure is defined by for all Borel subset of . We denote by the space of all oriented -varifolds.
In the rest of the paper, with a slight abuse of vocabulary, we will often use the word varifold instead of oriented varifold for the sake of concision. Recall that from the Riesz representation theorem, we can alternatively view any varifold as a distribution, i.e. an element of the dual space , where denotes the set of continuous functions vanishing at infinity on . It is defined for any test function by:
[TABLE]
As an additional note, another useful representation of a general varifold in can be obtained by the disintegration theorem (see Ambrosio \BOthers. (\APACyear2000) Chap. 2). Namely, if , for -almost every in , there exists a probability measure on such that is -measurable and we can write
[TABLE]
In other words, the varifold can be decomposed as its weight measure on together with a family of tangent space probability measures on the Grassmannian at the different points in the support of . This is usually referred to as the Young measure representation of .
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 and , the associated Dirac varifold acts on functions of by the relation
[TABLE]
can be viewed as a singular particle at position that carries the oriented -plane .
A second particular class is the one of rectifiable varifolds, which are in essence the varifolds representing an oriented shape of dimension . More precisely, given an oriented -dimensional submanifold of of finite total -volume, denoting by the oriented tangent space at , one can associate to the varifold , which is defined for all Borel subset by . Here, is the -dimensional Hausdorff measure on , i.e. the measure of -volume of subsets of (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 ,
[TABLE]
Such a representation can be extended to slightly more general objects known as oriented rectifiable sets. A subset of is said to be a countably -rectifiable set if , where are Lipschitz function for all (c.f. Simon (\APACyear1983)). We say that is an oriented rectifiable set if is a countably -rectifiable set and is a -measurable function such that for , is the approximate tangent space of at with specified orientation. Rectifiable subsets include both usual submanifolds but also piecewise smooth objects like polyhedra. Given any oriented rectifiable set , we can associate a varifold that we also write given again by (2.2). The set of those 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 : indeed, in the Young measure representation of (4), we have in this case the very particular constraint that probability measures are Dirac masses, specifically .
3. Metrics on varifolds
In this section, we address the issue of defining adequate metrics on the space . 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
As a measure/distribution space, can be equipped with various topologies and metrics, several of which have been regularly used in various contexts. We discuss a few of those below.
mass norm: with the previous identification of measures in with elements of the dual , one can define the following dual metric on :
[TABLE]
where . 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 and , whenever .
weak-* topology: a sequence of -varifolds converges to in the weak-* topology (denoted by ) if and only if for all (continuous compactly supported function)
[TABLE]
In fact, the weak-* topology on can be metrized by the following distance:
[TABLE]
where is a dense sequence in .
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 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.
Bounded Lipschitz metric: similar to the previous, the bounded Lipschitz distance (sometimes referred to as the flat metric) on is defined by
[TABLE]
It can be shown (cf. Ch 8 in Bogachev (\APACyear2007)) that metrizes the narrow topology on , namely the topology for which a sequence converges to if and only if for all bounded continuous functions .
Clearly, the narrow topology is stronger than the weak-* topology. Furthermore, it is also well known that locally metrizes the weak-* topology on , namely:
Proposition 3.1**.**
Let and be varifolds such that the sequence is tight. Then if and only if .
Proof.
Since metrizes the narrow topology, it suffices to show that converges to in the narrow topology. Let be a bounded continuous function defined on and . By the tightness property, we may choose a compact set such that . Let be an open ball that contains . Define
[TABLE]
From Tietz extension theorem, there exists a continuous extension of on such that and . This implies that
[TABLE]
Taking on both sides, we see that
[TABLE]
Since is arbitrary, we obtain that converges to in the narrow topology. ∎
As a direct consequence of Proposition 3.1, we have in particular that weak-* convergence and convergence in are equivalent if one restricts to varifolds that are supported in a fixed compact subset of . 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 and 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 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 . 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 and be continuous positive definite kernels on and respectively. Assume in addition that for any , . Then is a positive definite kernel on and the RKHS associated to is continuously embedded in i.e. there exists such that for any , we have .
We recall that the tensor product kernel has the exact expression . 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 and , namely that is a translation/rotation invariant radial kernel , with , and where is the inner product on inherited from 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 i.e. such that for all .
Now, if we let be the continuous embedding given by Proposition 3.2 and its adjoint, for any , we have
[TABLE]
With (10), we may identify as an element of the dual RKHS . Note that is not injective in general, in other words one can have in but in .
In any case, one can compare any two varifolds through the Hilbert norm of by defining:
[TABLE]
where we use the small abuse of notation of writing and instead of and on the two right hand sides. Due to the potential non-injectivity of , in general only induces a pseudo-metric on .
The main advantage of this construction is that can be expressed more explicitly based on the reproducing kernel property of . Indeed, given any and in , 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, 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 -universality of kernels:
Definition 3.4**.**
A positive definite kernel on a metric space is called -universal when its RKHS is dense in for the uniform convergence topology.
-universality has been studied in great length in such works as Carmeli \BOthers. (\APACyear2010); Sriperumbudur \BOthers. (\APACyear2011). In particular, one can provide characterizations of -universality for certain classes of kernels and spaces . In the case of translation-invariant kernels on for instance, it has been established that -universal kernels are the ones which can be expressed through the Fourier transform of finite Borel measures with full support on , which includes: compactly-supported kernels, Gaussian kernels, Laplacian kernels… With the previous definition, we have the following sufficient condition:
Theorem 3.5**.**
Suppose is a -universal kernel on , and . Let and be two oriented -rectifiable sets with , . If , then and -.
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 for all ), generalizing the result of Charon \BBA Trouvé (\APACyear2013).
However, the previous proposition does not necessarily lead to a distance on the full space . Counter-examples in the case are discussed for example in Hsieh \BBA Charon (\APACyear2019). To recover a true distance on , one needs the previous map 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 on the product space is -universal, specifically
Theorem 3.6**.**
The pseudo-distance induces a distance between signed measures of if and only if is -universal on . In particular, a sufficient condition for to be a distance on is that and are -universal kernels on and respectively.
Note that these conditions are more restrictive than in Theorem 3.5. To our knowledge, there is no simple characterization for general -universal kernels on the Grassmannian. However, within the setting of Remark 3.3, one easily constructs -universal kernels by restriction (based on the Plücker embedding) of -universal kernels defined on the vector space .
3.2.3. Comparison with classical metrics
We now study more precisely the topology induced by the (pseudo) distance on in comparison with the ones defined in Section 3.1. First of all, we observe that, for any with , one must have , where is the embedding constant of Proposition 3.2. Thus, for any and in , we have
[TABLE]
From the above inequalities we see that convergence in implies convergence in .
Remark 3.7*.*
With more assumptions on the regularity of the kernel , namely if is continuously embedded in , following a similar reasoning as above, one obtains the bound .
Suppose converges to in narrow topology. Since the map is continuous with respect to the narrow topology, we have
[TABLE]
as . Also, it’s clear that and hence with respect to . To summarize the discussion above:
Proposition 3.8**.**
Let and be varifolds in and assume that with respect to the operator norm or the narrow topology, then in .
Remark 3.9*.*
We emphasize that the result of Proposition 3.8 only requires the assumptions of Proposition 3.2 and thus holds whether is injective or not.
As for the weak-* topology, with the -universality assumption of Theorem 3.6 and restricting to varifolds with bounded total mass, we show that induces a topology stronger than weak-* convergence:
Proposition 3.10**.**
If is -universal, then the topology induced by is finer than the weak- topology on for any fixed .*
Proof.
Let and be varifolds in and assume that . For any and , there exists a such that . Then we obtain that from the following inequalities:
[TABLE]
∎
Note that the topology induced by may be strictly finer on . Indeed, if , consider , where and fixed. Then while for all . Yet, by combining Propositions 3.1, 3.8 and 3.10, we have the following
Corollary 3.11**.**
Let and be a compact subset. If is -universal, then metrizes the weak- convergence of varifolds on *
In summary, -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 , our goal is now to focus on deformation models for those objects in order to formulate and study the diffeomorphic registration problem on .
4.1. Deformation models
In this section, we discuss different models for how varifolds can be transported by a diffeomorphism of , in other words what are possible group actions of the diffeomorphism group on .
Let us start by considering the case of an oriented rectifiable subset . A diffeomorphism transports 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 , which is defined for all and by:
[TABLE]
in which denotes the determinant of the Jacobian of along (i.e. the change of d-volume induced by along at ) which is given by
[TABLE]
for an orthonormal basis of . One easily verifies that defines a group action which commutes with the action on oriented rectifiable sets, namely
Proposition 4.1**.**
For any oriented rectifiable set and diffeomorphism , .
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 . Indeed if with a measurable density function on and the Lebesgue measure, we can extend to a n-varifold in by taking a constant global orientation in . Then, for any orientation-preserving diffeomorphism , (15) writes in this case: with is the full Jacobian determinant of , leading to the usual action on densities .
However, in contrast with past works on submanifold registration, this is not the only possible group action that could be considered on the space . 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 . 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 be a fixed RKHS of vector fields on and be space of time dependent velocity fields such that for all , belongs to . The flow map is defined for all by and the ODE . If is continuously embedded in , one can show that for all , is a -diffeomorphism of . Moreover, on the subgroup of , one can define the following right-invariant Riemannian metric:
[TABLE]
Let us now consider a source (or template) varifold as well as a target . With the above deformation model and metric, registering to consists in finding a deformation that minimizes with the constraint that is close to in the sense of a kernel metric defined in Section 3.2. This can be reformulated as the following optimal control problem:
[TABLE]
with being the control, the total cost and the state equation is given by 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 and the target varifold . 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 is continuously embedded in , is continuously embedded in and , for some compact subset of , 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 and 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 for instance, it is enough to assume that and are functions such that all derivatives of up to order 2 vanish as .
As an important note, the formulation of (16) extends registration of submanifolds or rectifiable subsets in the sense that if and for two oriented d-rectifiable subsets of then (16) becomes equivalent, thanks to Proposition 4.1, to registering rectifiable subsets, i.e. to the problem
[TABLE]
with , 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 . The state variables are now given by the deformations which we view as elements of the Banach space . Let us denote, for , the mapping . We then introduce the Hamiltonian functional defined by:
[TABLE]
where is the costate variable which is a vector distribution of and denotes the duality bracket in . With the assumptions of Theorem 4.3, it follows from the maximum principle shown in Arguillère \BOthers. (\APACyear2015) that if is a global minimum of the optimal control problem, there exists a path of costates such that the following equations hold:
[TABLE]
with the end time boundary conditions . From the last equation in (18), we can attempt to deduce the form of the optimal . Introducing the Riesz isometry operator and its inverse , we get:
[TABLE]
One additional consequence of (18) is the following conservation of momentum again proved in Arguillère \BOthers. (\APACyear2015): for all and ,
[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 in the cost. This dependency is entirely encompassed by the boundary condition which we may describe a little more precisely based on the following:
Proposition 4.5**.**
The end-time momentum is a vector distribution in of the form
[TABLE]
where , and are continuous fields and for all , and denote the divergence and differential of restricted to .
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 . In that case, one can show, after some derivations, that the above expression of can be rewritten in the form of a vector distribution in with vectors normal to at each . In our more general situation, this is however not possible and is a priori a distribution that involves first order derivatives of the test function .
Now, the conservation law of (20) gives that for all ,
[TABLE]
Using the expression of in Proposition 4.5, and grouping all [math]-th and -st order terms in the resulting expressions, we may write in the general form:
[TABLE]
where and are continuous fields, with and . Furthermore, optimal vector fields satisfy and we have
[TABLE]
Denoting the reproducing kernel of , the reproducing kernel property implies that for all and , . Moreover, the similar property on the kernel first order derivatives Micheli \BBA Glaunès (\APACyear2014) gives that for any ,
[TABLE]
Then, we rewrite the linear maps as for any and where are the component vector fields of . By the above and the linearity of , 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 and 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 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 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 . 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 is a triangulated surface, with being the mesh triangles with specified orientations, one can write and for each approximate by , where is the center of , the oriented plane containing and . This leads to the approximation . As proved in Kaltenmark \BOthers. (\APACyear2017), this approximation provides an acceptable error bound for :
[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 and , 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 with finite total weight, can it be approximated by discrete varifolds and will approximations converge as ? 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 and metrics.
Theorem 5.1**.**
Let
[TABLE]
be the (non-convex) cone of discrete varifolds with at most Diracs. For any oriented varifold with , there exists a sequence such that . Moreover, if has compact support, then we can assume that for all , for some compact set and
[TABLE]
where is a constant that only depends on , and .
Proof.
We first tackle the case of compactly supported . Without loss of generality, we may also assume that is a probability measure. Let and be a closed ball that contains . For brevity, we write . Since we can view as a compact -dimensional submanifold of (using Plücker embedding), M is also regular of dimension (cf. Graf \BBA Luschgy (\APACyear2007)), i.e., and there exist , such that
[TABLE]
Given , by the 5-Times Covering Lemma (cf. Simon (\APACyear1983)), these exists a subset , such that and for all . Therefore,
[TABLE]
We can thus obtain a partition of from the the collection which satisfies and
[TABLE]
Let and and define . For any , with , we have
[TABLE]
Taking the supremum over all with , we obtain . Then for each , we can choose and we obtain such that
[TABLE]
and in particular .
Suppose now that is not compact: we show that for any , there exists a discrete varifold such that . Choose a compact set such that . From the previous case, we can find a discrete varifold such that , and hence . ∎
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 :
Corollary 5.2**.**
With the assumptions from proposition 3.2, one also has . If in addition , an equivalent upper bound as in Theorem 5.1 holds for .
We should point out that the asymptotic convergence rate given by the previous upper bound is rather slow, especially as the dimensions and grow. This is however under very mild assumptions on the varifold . We expect much better convergence properties for certain specific classes of varifolds, for instance assuming Alfors regularity as in 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 and , a natural idea is to look for the optimal discrete varifold in that approximates in terms of the metric . Due to the intricate structure of the set (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 (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 . Then the problem can be recast as the one of finding sparse approximations of in such a dictionary. It remains a NP hard problem but solutions can be approached either through greedy algorithms like orthogonal matching pursuit as proposed in Durrleman \BOthers. (\APACyear2009) or by considering the 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 . Such an increase in dimension makes the approach numerically impractical as soon as and . 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 of Diracs. The approach we develop in this section consists instead in directly tackling the non-convex problem of computing the projection onto for the class of kernel metrics . It shares some connection with the recent work of Chauffert \BOthers. (\APACyear2017) that considers a related problem for standard measures defined on the torus .
Fix a varifold . For any , we seek that is closest to for , namely
[TABLE]
By construction, if then Corollary 5.2 will imply that converges to in the metric . 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 and defining the kernels are non-negative. Then for any and , there exists such that
Proof.
Let be a minimizing sequence, be the dual operator and . Without loss of generality, we may assume that there is a such that remains bounded and tends to as .
Observe that must be bounded. If it’s not bounded, then from the assumptions that , we obtain
[TABLE]
as , which is absurd.
Since , and
[TABLE]
are all bounded sequences of , we may replace them by convergent subsequences, thus we could assume that
[TABLE]
Since , the matrix must has the following form:
[TABLE]
where and are -by- and -by- semi-positive definite matrices. Combining this with the assumption , we obtain
[TABLE]
where , , and . Since is bounded we can assume that . Let , then
[TABLE]
Hence is a minimizer. ∎
However, in general this projection is not unique. We also point out that the existence is a priori not guaranteed if kernels and take negative values. It is so far an open question to determine to what extent one could generalize the result of Proposition 5.3, one particular but important case being the one of current metrics obtained for 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 with respect to the . In particular, we have
[TABLE]
which gives . It shows that for any , .
5.3. -convergence of registration functionals
Ultimately, our purpose is to use the previous approximating discrete varifolds to approximate the diffeomorphic registration problem (16). The natural question that arises is whether replacing the source varifold by its projections 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 -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 with compact support and finite mass and such that , we will need that for some compact set or that . Unfortunately, this does not hold in general since convergence in does not allow to control the support nor the total mass of the sequence .
Yet, provided that , we can actually retrieve the boundedness of the total mass. We assume in what follows that the kernels are such that and .
Lemma 5.4**.**
Let be a sequence of discrete varifolds with finite mass such that there exists a compact with for all . We assume that is bounded. Then is bounded.
Proof.
We prove it by contradiction. Assume that is unbounded. Then, up to extracting a subsequence, we can assume that . Let’s write . Thus . Since and are continuous and , we can find compact subsets and with diameters small enough, so that: , and , where . It follows that, as ,
[TABLE]
which is a contradiction. ∎
Lemma 5.4 suggests that one should enforce the uniform compactness of the supports of the . 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 stays in a compact set containing . We still have to verify the convergence of the resulting sequence:
Proposition 5.5**.**
Let be a varifold with finite mass and be a compact set in which contains . Construct the approximating sequence of by solving the following constrained optimization problem:
[TABLE]
Then converges to in and, if the kernel is -universal, it also converges in .
Proof.
Thanks to Theorem 5.1 and Lemma 5.4, we immediately get that as and . Moreover, if is -universal, then by Proposition 3.10 it implies that . Since and , weak-* convergence implies that converges to in by Proposition 3.1. ∎
We are now able to state the main result of this section. We assume that the source/template varifold is compactly supported and we fix is a compact subset of that contains . Then for any , is defined as in Proposition 5.5 and we introduce the following energy functionals :
[TABLE]
which are the equivalent to the energy of the original problem (16) but replacing the template varifold by its approximations .
Theorem 5.6**.**
With the above notations, we assume that the reproducing kernel of is -universal and satisfies all the conditions of Proposition 5.3. We also assume the continuous embedding . Then, the sequence of functionals -converges to for the weak topology on . Consequently, if is a global minimizer of for each , then is bounded in and every cluster point for the weak topology of is a global minimum of .
Proof.
We first show that whenever converges to weakly in , we have
[TABLE]
Since is lower semicontinuous with respect to the weak topology, we only need to prove the following,
[TABLE]
For all with , we have
[TABLE]
where and , for all . Similar to the computation done in the proof of Theorem 4.3, we see that
[TABLE]
Taking supremum over all with , we obtain the following inequality,
[TABLE]
From Proposition 5.5, converges to in the narrow topology. Hence the right hand side in the equation above converges to [math] as . This proves (29).
Second, we need to show that for each , there exists a sequence converging to weakly such that
[TABLE]
In fact, it suffices here to take to be the constant sequence 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 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 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 may become prohibitively large in practice. We may instead choose to represent an element by an oriented frame of independent vectors for which . Such a representation is of course not unique since elements of are equivalence classes of oriented frames but we leave to the next section the more thorough analysis of the additional invariances that this representation will imply.
We will in fact go one step further by also incorporating the weight of Dirac varifolds in this frame representation itself, which is done as follows. Let be a discrete varifold of the form . For each , we consider a frame such that
[TABLE]
In other words, the oriented space spanned by the frame corresponds to while its -volume matches the weight . Given such a choice of frame for each , we can then identify with the (non-unique) state variable in the vector space . Conversely, such a frame with a matrix of rank for all , corresponds to the (unique) discrete oriented varifold defined by the relations of (31); we will denote it by in what follows.
In this representation, the kernel metrics for discrete varifolds expressed in (23) can be explicitly written as
[TABLE]
where . Note that this expression does not depend on the choice of frames that satisfy the conditions of (31) for (and similarly for ). In the case where is a more general non-discrete varifold in , the computation of involves integrals over 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 metric then formally reduces to an expression equivalent to (32) in which the and 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 . The gradient of can be computed by direct differentiation of expressions like (32) with respect to the and . 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 Diracs composing the varifold . Note that one of the main downside of this projection algorithm, in contrast with the previously mentioned approach of fixing a dictionary and solving a convex sparse decomposition problem, is that we can provide no general guarantees of convergence to a global minimum of (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 be a diffeomorphism of and , the pushforward action 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 . 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 case).
Following once again the Pontryagin maximum principle approach, the Hamiltonian for this discrete representation is given by:
[TABLE]
with denoting respectively the costates for the position and frame vector variables. The PMP then shows that optimal trajectories of the registration problem are governed by the dynamical system:
[TABLE]
while optimal vector fields satisfy
[TABLE]
Plugging the above expression of with respect to in the Hamiltonian (33) gives the reduced Hamiltonian which writes:
[TABLE]
and (38) becomes a coupled system in the variables and called the reduced Hamiltonian equations. Consequently, the set of optimal paths is entirely determined by the initial values and the value of the reduced Hamiltonian 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 , the matrix
[TABLE]
is constant in time.
Proof.
Using the Hamiltonian equations written above, we have for all
[TABLE]
Hence is a constant matrix. ∎∎
Note that, at this point, all those equations are fundamentally modelling the deformation of the frames 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), and thus are independent of the choices of the frame vectors that span the same oriented vector spaces with the same -volumes . This in turn leads to a set of conditions satisfied by the different components of the final costate and, with Lemma 6.1, of the full path . These are summed up by the following result:
Proposition 6.2**.**
Let be optimal trajectory, then for all , the matrices as defined above are constant scalar matrices. In particular, we have for all , and .
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 actually lie in affine subspaces of of lower dimensions , which is precisely the dimension of the ’true’ state space .
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 , which is known, and . One of the standard class of methods in optimal control, known as shooting methods, consist in directly optimizing the cost function over , 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 with respect to the initial costate . The regularization term being equal to thanks to the conservation of the reduced Hamiltonian, its gradient can be obtained by direct differentiation of (6.2). The fidelity term on the other hand depends indirectly on the initial costate 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 with respect to can be computed by flowing backward in time the adjoint Hamiltonian system
[TABLE]
where , the adjoint variables of the system, together with the end-time conditions and . 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 in the directions of and . We then specifically approximate the above by centered finite difference
[TABLE]
for some small , which only requires at each time two evaluations of the same function that appears in the forward reduced Hamiltonian equations.
With the above approach to compute the gradient with respect to , the registration algorithm then consists of essentially the same steps as the aforementioned works:
1:repeat
2: From compute by forward integration of the reduced Hamiltonian system given by (38) and (39).
3: Compute and .
4: Integrate backward the adjoint Hamiltonian equations (42) to obtain .
5: Deduce the gradient of the full cost function with respect to .
6: Update .
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 , where we typically take 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 to the linear subspace . 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 and . In all these experiments, we choose the deformation kernel of to be a diagonal Gaussian kernel . The kernel function is a Gaussian of scale . The choice of these scales is adapted to the sizes of the shapes in each of the experiment. We will not discuss these questions more in detail here, since this is not our main topic and it has been more thoroughly analyzed in previous works such as Bruveris \BOthers. (\APACyear2012); Kaltenmark \BOthers. (\APACyear2017); Hsieh \BBA Charon (\APACyear2019). The function 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 to recover an orientation invariant fidelity metric, or 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 . 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 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 which involve the tangent vectors to the curve that passes through . We then complete it into a 2-varifold by adding a second ”vertical” (i.e. inter-sectional) frame vector , which can be estimated in this case by simply finding the projection of 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 () to the position of each sampled point. A first possible registration approach could be to treat such point clouds as standard measures of (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 is composed of a very high number of Diracs, we compute the compressed varifolds given by the of (24) for increasing values of and evaluate the resulting quantization error in terms of the metric. Then we solve the registration problems to a fixed target from the source varifolds given by the in lieu of , and compare the estimated solutions to the registration of the groundtruth. For comparison, we will evaluate the total energy of the estimated deformation fields 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 of the source curve as a function of (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 , we observe that the error is indeed smaller than with the subsampling approach. We also display a few of the quantized for several values of . 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” 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 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 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 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 with different values of , 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 but that in addition we obtain a very good approximation of with only a small fraction of the initial number of Diracs. Some of the quantized varifolds 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 being a numerical solution for the groundtruth and the solutions based on the quantized source shapes, the third row of Fig. 5 shows the difference of the energies . 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 . 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 , 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 on . 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 . Let us denote by and the RKHS associated to kernels and respectively. Suppose that and are rectifiable sets as above such that and . Without loss of generality, we may assume that . From Lusin’s theorem, there exists a subset of such that is continuous and . Let us denote by , we see that . Since for ,
[TABLE]
(cf Evans (\APACyear2018)), there exists for any .
Let be defined by . Since is continuous on and , there exists such that . Let and , then and . Using the density of in together with the fact that is -universal, there exist and such that in and . Now, since and in , we have
[TABLE]
which is a contradiction. Hence we have
Next, we show that -. Let and assume that . From Lusin’s theorem, there exists subset such that is continuous and . Using the upper density argument as above, we can find such that for all . Since the map restricted to is continuous, there exists a satisfying:
[TABLE]
Define , and . Observe that, from the assumption ,
[TABLE]
From this, we may assume that . Let and be sequences in and such that converges to in and . We obtain
[TABLE]
which is impossible. ∎
Proof of Theorem 4.3 Thanks to the first term in , any minimizing sequence of is bounded in . Let be a subsequence of such minimizing sequence which converges weakly to some in . Using the results of Younes (\APACyear2019) Chapter 7.2, we know that
[TABLE]
Furthermore, for any , we have
[TABLE]
Now, using the embedding
[TABLE]
Taking supremum over all with , we obtain that
[TABLE]
as . Combining this with lower semicontinuity of , we finally obtain that
[TABLE]
and hence is a global minimizer. ∎
Proof of Proposition 4.5 Recall that for all , which we may rewrite as
[TABLE]
Thus, the variation with respect to in the Banach space writes
[TABLE]
where . Moreover
[TABLE]
Taking the variation with respect to along any , we obtain:
[TABLE]
where the last term follows from the differentiation of Gram determinant matrices while the notation 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 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 the integrand involves a linear combination (depending on ) of the partial derivatives of along the subspace i.e. of the elements of the matrix . Thus, without attempting to specify this term explicitly, we can in general write it as where is a continuous map from into giving us
[TABLE]
The result of the theorem then follows by setting , and . ∎
Proof of Proposition 6.2 We can treat the case of each particle separately and thus, without loss of generality, we may directly assume that . We write , for the state and costate variables along an optimal trajectory and
[TABLE]
Consider the group of linear transformations, , i.e., for any ,
[TABLE]
where and are the orthogonal projections of on and , with and . The Lie algebra of is and is the set of all zero trace linear transformations of . Now, consider the action of on defined as:
[TABLE]
for any . We see that for all and therefore .
Now, if we let be a smooth curve in that satisfies and , differentiating the equality shows that for any , we have
[TABLE]
Since , we must have that is a zero trace linear map. For any , we may choose such that and , which leads to . Consequently,
[TABLE]
for some constant . In addition, for any , we can also choose such that and , which gives . It results that .
Finally, since is constant by Lemma 6.1, we obtain that
[TABLE]
for all . ∎
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Allard ( \APA Cyear 1972) \APA Cinsertmetastar allard 1972 first {APA Crefauthors} Allard, W \BPBI K. \APA Cref Year Month Day 1972. \BBOQ \APA Crefatitle On the first variation of a varifold On the first variation of a varifold. \BBCQ \APA Cjournal Vol Num Pages Annals of mathematics 417-491. \Print Back Refs \Current Bib
- 2Allard \B Others . ( \APA Cyear 2012) \APA Cinsertmetastar allard 2012 multi {APA Crefauthors} Allard, W \BPBI K., Chen, G. \BCBL \BBA Maggioni, M. \APA Cref Year Month Day 2012. \BBOQ \APA Crefatitle Multi-scale geometric methods for data sets II: Geometric multi-resolution analysis Multi-scale geometric methods for data sets ii: Geometric multi-resolution analysis. \BBCQ \APA Cjournal Vol Num Pages Applied and Computational Harmonic Analysis 323435–462. \Print Back Refs \Current Bib
- 3Almgren ( \APA Cyear 1966) \APA Cinsertmetastar almgren 1966 plateau {APA Crefauthors} Almgren, F \BPBI J. \APA Cref Year 1966. \APA Crefbtitle Plateau’s problem: an invitation to varifold geometry Plateau’s problem: an invitation to varifold geometry ( \BVOL 13). \APA Caddress Publisher American Mathematical Soc. \Print Back Refs \Current Bib
- 4Ambrosio \B Others . ( \APA Cyear 2000) \APA Cinsertmetastar Ambrosio 2000 {APA Crefauthors} Ambrosio, L., Fusco, N. \BCBL \BBA Pallara, D. \APA Cref Year 2000. \APA Crefbtitle Functions of bounded variation and free discontinuity problems Functions of bounded variation and free discontinuity problems. \APA Caddress Publisher Oxford : Clarendon Press. \Print Back Refs \Current Bib
- 5Ardekani \B Others . ( \APA Cyear 2011) \APA Cinsertmetastar Ardekani 2011 {APA Crefauthors} Ardekani, S., Jain, A., Jain, S., Abraham, T \BPBI P., Abraham, M \BPBI R., Zimmerman, S. \BDBL Younes, L. \APA Cref Year Month Day 2011. \BBOQ \APA Crefatitle Matching Sparse Sets of Cardiac Image Cross-Sections Using Large Deformation Diffeomorphic Metric Mapping Algorithm Matching Sparse Sets of Cardiac Image Cross-Sections Using Large Deformation Diffeomorphic Metric Mapping Algorithm. \BBCQ
- 6Arguillère \B Others . ( \APA Cyear 2016) \APA Cinsertmetastar Arguillere 2016 b {APA Crefauthors} Arguillère, S., Miller, M. \BCBL \BBA Younes, L. \APA Cref Year Month Day 2016. \BBOQ \APA Crefatitle Diffeomorphic Surface Registration with Atrophy Constraints Diffeomorphic Surface Registration with Atrophy Constraints. \BBCQ \APA Cjournal Vol Num Pages SIAM Journal on Imaging Sciences 93975-1003. \Print Back Refs \Current Bib
- 7Arguillère \B Others . ( \APA Cyear 2015) \APA Cinsertmetastar arguillere 14:_shape {APA Crefauthors} Arguillère, S., Trélat, E., Trouvé, A. \BCBL \BBA Younes, L. \APA Cref Year Month Day 2015. \BBOQ \APA Crefatitle Shape deformation analysis from the optimal control viewpoint Shape deformation analysis from the optimal control viewpoint. \BBCQ \APA Cjournal Vol Num Pages Journal de Mathématiques Pures et Appliquées 1041139-178. \Print Back Refs \Current Bib
- 8Aronszajn ( \APA Cyear 1950) \APA Cinsertmetastar Aronszajn 1950 {APA Crefauthors} Aronszajn, N. \APA Cref Year Month Day 1950. \BBOQ \APA Crefatitle Theory of reproducing kernels Theory of reproducing kernels. \BBCQ \APA Cjournal Vol Num Pages Trans. Amer. Math. Soc.68337-404. \Print Back Refs \Current Bib
