Elastic Shape Analysis of Surfaces with Second-Order Sobolev Metrics: A Comprehensive Numerical Framework

This paper introduces a set of numerical methods for Riemannian shape analysis of 3D surfaces within the setting of invariant (elastic) second-order Sobolev metrics. More specifically, we address the computation of geodesics and geodesic distances between parametrized or unparametrized immersed surfaces represented as 3D meshes. Building on this, we develop tools for the statistical shape analysis of sets of surfaces, including methods for estimating Karcher means and performing tangent PCA on shape populations, and for computing parallel transport along paths of surfaces. Our proposed approach fundamentally relies on a relaxed variational formulation for the geodesic matching problem via the use of varifold fidelity terms, which enable us to enforce reparametrization independence when computing geodesics between unparametrized surfaces, while also yielding versatile algorithms that allow us to compare surfaces with varying sampling or mesh structures. Importantly, we demonstrate how our relaxed variational framework can be extended to tackle partially observed data. The different benefits of our numerical pipeline are illustrated over various examples, synthetic and real. Supplementary Information The online version contains supplementary material available at 10.1007/s11263-022-01743-0.


Motivation
Over the past decades, advances in imaging techniques and devices have led to significant growth in the quantity and quality of "shape data" in several fields, such as biomedical imaging, neuroscience, and medicine. By "shape data", we mean objects whose predominantly interesting features are of geometric and topological nature; examples of which include functions, curves, surfaces or probability densities. Naturally, this prompted the emergence of new mathematical and algorithmic approaches for the analysis of such objects, which led to the development of the growing fields of geometric shape analysis and topological data analysis, see e.g Younes (2010), Srivastava and Klassen (2016), Kendall et al. (1999), Edelsbrunner and Harer (2022), Carlsson (2014), Bronstein et al. (2008Bronstein et al. ( , 2021. In this paper, we will focus on 3D surface data, which is becoming increasingly prominent in several areas due to the emergence of high accuracy 3D scanning devices. The domain of geometric shape analysis has produced several mathematical frameworks and numerical algorithms for the comparison and statistical analysis of 3D surfaces that have proven to be useful in numerous applications Srivastava & Klassen, 2016;Younes, 2010). In the context of shape analysis of surfaces, we distinguish between two fundamentally different scenarios: the analysis of surfaces with known point correspondences (parametrized surfaces), and that of surfaces where the point correspondences are unknown (unparametrized surfaces). In the discrete case, i.e. for simplicial meshes, working with parametrized surfaces thus involves having known one-to-one correspondences between the surfaces' vertices, which implies in particular that their mesh structures are required to be consistent. There exists a plethora of different numerical frameworks for shape analysis of parametrized surfaces, see e.g. Jermyn et al. (2017), Kilian et al. (2007), Rumpf and Wirth (2015a), Iglesias et al. (2018) and Pierson et al. (2022) and the references therein. Nevertheless, one rarely ever encounters 3D surface data with consistent mesh structures in practical applications, and thus methods designed for shape analysis of parametrized surfaces are severely limited when used for applications with real data. This motivates the need for registering unparametrized surfaces, i.e., finding the unknown point correspondences between them, as well as the development of tools to compare unparametrized surfaces, i.e. to quantify similarity between them, and for the statistical analysis of unparametrized surfaces.
While algorithms that are designed for the comparison and statistical analysis of unparametrized surfaces usually lead to a counterpart for dealing with parametrized surfaces, the converse is unfortunately far from being true. In the past, this difficulty has been approached by first registering the data in a pre-processing step under a certain metric or objective function, before subsequently comparing and perform statistical analysis of the data independently of this registration metric (Audette et al., 2000). This practice is, however, being increasingly questioned as it is easy to construct examples where it leads to a severe loss in the data structure, see e.g. Srivastava and Klassen (2016) and the references therein. This motivates the need for a more comprehensive solution, where the registration, comparison and statistical analysis are performed jointly. In the field of geometric shape analysis, this is achieved by viewing the space of parametrized surfaces as an infinite-dimensional manifold and equipping it with a Riemannian metric that is invariant under the action of the reparametrization (registration) group. This invariance implies that the Riemannian metric descends to a metric on the quotient space of unparametrized surfaces, which consequently allows us to perform the registration, surface comparison and ensuing statistical analysis in a unified framework.

Related Work in Riemannian Shape Analysis
The Riemannian approach to shape analysis has several benefits. First, a Riemannian metric models a very natural notion of similarity: a Riemannian metric measures the cost of deformations, and can thus be used to define the distance (similarity) between two surfaces as the cost of the cheapest deformation that transforms one surface (the source) onto the other (the target). Furthermore, a Riemannian framework not only leads to a notion of similarity between pairs of surfaces, but also allows one to compute optimal point correspondences and optimal deformations (called geodesics) between the (aligned) surfaces, cf. Figs. 1 and 2 for various examples that are computed with the framework of the present paper. Finally, the Riemannian approach directly allows one to apply the methods of geometric statistics (Pennec, 2006;Pennec et al., 2019) to develop a comprehensive statistical framework for shape analysis, cf. the algorithms developed in Sect. 5.
Riemannian metrics on spaces of surfaces come in two flavors: intrinsic metrics, which are defined directly on the surface, and extrinsic metrics, which are inherited from right invariant metrics on the diffeomorphism group of R 3 . Intuitively the first approach corresponds to deforming the surface only, while the second approach applies a deformation of the whole ambient space in which the surface is embedded. The latter, which is inherited from the principles of Grenander's pattern theory (Grenander, 1996), has notably led to the celebrated LDDMM framework, for which powerful numerical toolboxes have been developed (Beg et al., 2005;Charlier et al., 2021).
The present paper operates in the intrinsic setup. More specifically, we deal with the class of reparametrization invariant Sobolev metrics on spaces of surfaces. This class of metrics was first introduced in the context of spaces of curves by Michor and Mumford (2007), and Mennucci et al. (2008). While these two articles focused on the theoretical properties of these Riemannian metrics, several ensuing numerical frameworks have been developed, see e.g. Srivastava et al. (2011), Bauer et al. (2019a), Nardi et al. (2016) and the references therein. Subsequently, this framework has been generalized to the space of surfaces by the last author and collaborators (Bauer et al., 2011. While several of the theoretical results for these metrics on the space of curves have been generalized for the space of surfaces (e.g., local well-posedness of the geodesic equation, non-vanishing geodesic distance), a comprehensive numerical framework is largely missing.
The most popular numerical approach for shape analysis of surfaces is based on the Square Root Normal Field (SRNF) framework proposed in Jermyn et al. (2012). This framework defines the SRNF transformation φ, which is a mapping from the space of surfaces I that takes values in Fig. 1 Examples of optimal deformations (geodesics) between different types of data with unknown point correspondences: genus zero surfaces (line 1 and 2), higher genus surfaces with boundaries and inconsistent topologies (line 3 and 4), shape complexes (line 5 and 6), partial matching (line 7). Animations of the obtained surface deformations can be found in the supplementary material and on the github repository Fig. 2 Point correspondences obtained after matching two unparametrized surfaces: the coloring of the two surfaces encode the obtained point correspondences. In addition, we highlight the obtained matching for selected points by displaying connecting lines L 2 (S 2 , R 3 ), where S 2 is the unit sphere. This mapping can then be used to define a (pseudo) distance function on I via the pullback of the L 2 distance. This framework is related to intrinsic Riemannian metrics on surfaces as the resulting (pseudo) distance function is a first-order approximation of the geodesic distance of a particular (degenerate) Sobolev metric of order one (Jermyn et al., 2012). The simplicity of the computation of this pseudo-distance has led to several implementations Laga et al., 2017), which have been shown to be effective in applications, see e.g. Kurtek et al. (2014), Joshi et al. (2016), Matuk et al. (2020) and Laga et al. (2022).
However, the SRNF framework has several theoretical shortcomings: first, the non-injectivity of φ implies that the pullback of the L 2 metric by φ is degenerate. Consequently, there arises the phenomenon that distinct shapes are indistinguishable by the SRNF shape distance. This behavior was originally studied in (Klassen & Michor, 2020) and was further discussed in (Bauer et al., 2022b), where it was shown that for each closed surface there exists a convex surface which is indistinguishable by the SRNF distance. Moreover, the image of I via φ is not convex, which implies that the SRNF distance is indeed only a first-order approximation of a geodesic distance function rather than a true geodesic distance on I, i.e., the SRNF distance does not come from geodesics (optimal deformations) in I. Furthermore, the problem of inverting the SRNF transformation φ to recover an optimal deformation in I from a geodesic in L 2 (S 2 , R 3 ) is highly ill-posed.
Consequently, to overcome the theoretical challenges of the SRNF pseudo-distance, it is natural, instead, to consider the reparametrization invariant Sobolev metrics mentioned previously. In Su et al. (2020), a first step towards obtaining a more general numerical framework was achieved: the authors proposed a numerical framework for a family of firstorder Sobolev metrics. The main drawback of this framework is the requirement for the data to be given by a spherical parametrization, which severely limits its applicability in practical contexts (see the comments below). In addition, numerical experiments suggested that it would be beneficial to consider metrics that involve higher-order terms to prevent the occurrence of certain numerical instabilities. Prior to the present paper, and to the best of the authors' knowledge, no implementation of more general (higher) order metrics was available.
The major difficulty in the implementation of these Riemannian frameworks is the discretization of the reparametrization group. In Laga et al. (2017) and Su et al. (2020), the authors used a discretization via spherical harmonics. These methods provide a relatively fast and stable approach for solving the registration of two spherical surfaces but requires that the surfaces are of genus-zero and are given by their spherical parametrizations. However, in most applications, data is typically given as triangular meshes that are not a priori homeomorphic to S 2 . As the reparametrization problem is highly non-trivial and computationally expensive, a better approach for working with real data consists of developing methods that deal directly with triangular meshes.
Inspired by the use of tools from geometric measure theory and in particular by varifold norms with the LDDMM model, Bauer et al. (2021) proposed a varifold matching framework to register surfaces with respect to the SRNF pseudo-distance. This approach provides several benefits: notably, the reparametrization group does not need to be discretized and its action on I does not need to be implemented. This allows one to work with simplicial meshes without having to first produce spherical parametrizations. Moreover, it extends to the analysis of surfaces with more general topologies with or without boundaries. Yet, this framework still suffers from the theoretical disadvantages of the SRNF pseudo-distance discussed above and it has been observed that the degeneracy of the distance can also lead to important numerical artefacts (c.f. Fig. 7 below). Consequently, it seems natural to combine this framework with more general Riemannian metrics on I, which is one of the main contributions of the present article, as we explain in the following section.

Contributions
The central contribution of the present paper is the development of an open-source numerical framework for the statistical shape analysis of surfaces (triangular meshes) under second-order reparametrization invariant Sobolev metrics. In addition, our framework allows one to deal with topologically inconsistent and/or partially observed data. The code is available on github: https://github.com/emmanuel-hartman/H2_SurfaceMatch Towards this end, we extended the relaxed varifold matching framework of Bauer et al. (2021) to compute the geodesic distance with respect to a reparametrization invariant second-order Sobolev metric on I and introduce a natural discretization of this metric for triangular meshes. This framework is the first implementation of higher-order Sobolev metrics on parametrized and unparametrized surfaces. In contrast with Bauer et al. (2021), our framework directly produces geodesics (i.e. the optimal deformation path), and the addition of higher-order terms prevents the formation of numerical artefacts as mentioned above. By splitting the metric into separate terms, we are also able to control the geometric changes penalized by the metric. This allows us to model different deformations as well as control the regularizing effects of the higher-order terms, thus making our framework versatile for a variety of applications.
In addition to providing a framework for surface matching, we develop a comprehensive statistical pipeline for the computation of Karcher means, tangent principal component analysis, and parallel transport. As an application of the latter, we demonstrate how it can be used for motion transfer between surfaces. Thus, our framework is well adapted to the statistical analysis of populations of shapes such as the ones appearing in biomedical applications. To further improve the robustness of our proposed methods, we also implement a weighted varifold matching framework by extending the idea proposed in the context of curves and shape graphs by Sukurdeep et al. (2022). The joint estimation of weights on the source surface enables this augmented model to deal more naturally with partial matching constraints or missing parts in the target shape, or differences in topology between the two shapes.

Outline
Our paper is structured as follows: In Sect. 2, we introduce the family of H 2 -metrics (second-order Sobolev metrics) on the space of parametrized and unparametrized surfaces. In Sect. 3, we formulate a varifold-based relaxed matching problem that allows us to estimate geodesics and distances induced by H 2 -metrics on the shape space of unparametrized surfaces. We then describe a set of numerical approaches for the computation of these geodesic and distance estimates in Sect. 4, before leveraging these algorithms for the development of more general tools for the statistical shape analysis of sets of surfaces in Sect. 5. Finally, we extend our secondorder elastic surface analysis framework to the setting of surfaces which may have incompatible topological properties or exhibit partially missing data in Sect. 6.

Sobolev Metrics on Surfaces
In this section, we introduce the theoretical background on second-order elastic Sobolev Riemannian metrics for spaces of parametrized and unparametrized surfaces, which will provide the key ingredient of our statistical framework for shape analysis of surfaces.

Metrics on Spaces of Parametrized Surfaces
We begin by introducing the main definitions and known theoretical results on second-order Sobolev metrics for spaces of parametrized surfaces in R 3 , which we shall rely on for the remainder of this paper. Let M denote a 2-dimensional compact manifold, possibly with boundary, whose local coordinates are denoted by (u, v) ∈ R 2 . A parametrized immersed surface in R 3 is an oriented smooth mapping q ∈ C ∞ (M, R 3 ), which in addition, we assume to be regular, i.e., we require its differential dq to be injective at every point of M. The set of all such parametrized surfaces, which we denote by I, is itself an infinite-dimensional manifold, where the tangent space at any q ∈ I, denoted T q I, is given by C ∞ (M, R 3 ). Any such tangent vector h ∈ T q I can be thought of as a vector field along the surface q.
Next we introduce the reparametrization group D, which is the group of orientation-preserving diffeomorphisms of M, i.e., the space of all ϕ ∈ C ∞ (M) such that det (Dϕ(u, v)) > 0 for all (u, v) and ϕ −1 ∈ C ∞ (M), where Dϕ denotes the differential (or Jacobian) of the diffeomorphism ϕ. For any immersed surface q ∈ I and ϕ ∈ D, we say that q • ϕ ∈ I is a reparametrization of q by ϕ.
Our goal is to equip the manifold I with a Riemannian metric that will subsequently enable us to develop a framework for the comparison and statistical shape analysis of surfaces. Recall that any Riemannian metric G on I induces a (pseudo) distance on this space, which is given for any two parametrized surfaces q 0 , q 1 ∈ I by (1) with the infimum being taken over the space of all paths of immersed surfaces connecting q 0 and q 1 , which we write as: with ∂ t q(t) denoting the derivative with respect to t of this path. In finite dimensions this distance, which is called the geodesic distance, is always non-degenerate, i.e., a true distance. In our infinite-dimensional setting it can, however, be degenerate (Michor & Mumford, 2005). As our main goal will be the analysis of unparametrized surfaces, we will require our Riemannian metric to be invariant under the action of the aforementioned reparametrization group D, i.e., we require G to satisfy for all q ∈ I, h, k ∈ T q I and ϕ ∈ D, which will imply that the induced geodesic distance as defined in (1) satisfies for all q 0 , q 1 ∈ I and ϕ ∈ D. This will later allow us to consider the induced Riemannian metric (and distance function) on the quotient space of unparametrized surfaces, cf. Sect. 2.2. The simplest and potentially most natural such metric is the reparametrization invariant L 2 -metric, which is given by where vol q is the surface area measure of the immersion q, which in local coordinates (u, v) is given by where the subscripts denote partial derivatives, × denotes the cross product on R 3 , and | · | denotes the norm on R 3 . This Riemannian metric is, however, not useful for any application in shape analysis, as it results in vanishing geodesic distance on both the spaces of parametrized and unparametrized surfaces (Bauer et al., 2012;Michor & Mumford, 2005). Vanishing geodesic distance refers to the phenomenon where the geodesic distance induced by the L 2 -metric between any two surfaces is zero. Consequently, we are interested in stronger Riemannian metrics that induce meaningful distances. A natural approach to strengthen the metric consists of incorporating derivatives of the tangent vector, leading to the class of first-order Sobolev metrics. Therefore we let g q = q * ·, · be the pullback metric of the Euclidean metric on R 3 , see Fig. 3 for an explanation of this construction.
A first-order Sobolev metric is then given by To interpret the first-order term g −1 q (dh, dh), we view the differential dh as a vector valued one form, i.e., as a map from T M to R 3 . Then the inverse of the pullback metric g −1 q can be used to pair such mappings. To understand this pairing better, we can fix a set of coordinates and view all the involved objects as matrix fields. Then we have where dh T denotes the point wise transpose of the matrix field dh. By the results of Bauer et al. (2011), we know that this metric indeed overcomes the degeneracy of the L 2metric, i.e., the corresponding geodesic distance function is non-degenerate. Next we further decompose the first-order term into four different terms which each have a geometric interpretation. Therefore, we write where A straight-forward calculation shows that these terms are orthogonal with respect to the inner product see (Su et al., 2020). Consequently we have: The geometric meaning of the first three terms becomes clear in the following result: Remark 1 (Su et al. 2020) Let q ∈ I and h ∈ T q I. The term measures the change of the pull-back metric g q while keeping the volume form constant (shearing). The second term M g −1 q (dh + , dh + ) vol q measures the change of the volume density vol q (scaling), while the third term measures the change in the normal vector n q (bending).
The interpretation of the last summand is less clear: it can be thought of as measuring the deformation of the local parametrization by a rotation in the parameter space M.

Remark 2
The class of first-order Sobolev (pseudo-) metrics, i.e. metrics obtained as weighted combinations of the four first-order terms discussed above, have often been referred to as elastic metrics in the shape analysis literature (Jermyn et al., 2012. Beyond just the mere high level analogy of these metrics measuring some form of bending or stretching energies, it turns out that these can in fact be more precisely connected to classical linear elasticity, specifically as the thin shell limit of the elastic energy of a layered isotropic material. For the purpose of concision, we will not elaborate on this particular point in this paper, but this connection will be highlighted in more details in an upcoming preprint. The above considerations suggest that metrics of this form provide a meaningful class of metrics for shape analysis of surfaces: they overcome the degeneracy of the L 2 -metric and admit a physical interpretation of the different terms involved. There is, however, numerical evidence that these first order metrics are still too weak for our targeted applications, see the experiments in Fig. 7. Thus we will augment the metric with a further higher-order term involving the Laplacian q induced by the immersion q, which using Einstein summation is given in local coordinates (u, v) by where |g q | denotes the determinant of the pullback metric in the local coordinate frame. This allows us to define a second-order term via M q h, q h vol q .
By adding up all the zero, first and second-order terms, we arrive at the main object of the present article: the family of H 2 -metrics (second-order Sobolev Riemannian metrics) for surfaces, which is given by: Here a 0 , a 1 , b 1 , c 1 , d 1 , a 2 are non-negative weighting coefficients for the different terms in the metric. Note that this family incorporates the Riemannian metric corresponding to the SRNF (pseudo) distance  and the families of elastic Riemannian metrics as proposed by Jermyn et al. (2012) and Su et al. (2020). For a general treatment of properties of Sobolev metrics we refer to the article (Bauer et al., 2011), and for a detailed explanation of the influence of these coefficients on our numerical experiments, see the discussion in Sect. 4.7.
The following result, which summarizes the invariances of our family of metrics, ensures that the metric descends to quotient spaces with respect to the corresponding group actions: Lemma 1 The family of H 2 -metrics G is invariant under the action of the group of reparametrizations D, the group of rotations Rot(R 3 ) and the group of translations R 3 , i.e., for any q ∈ I, h, k ∈ T q I, R ∈ Rot(R 3 ) and τ ∈ R 3 we have It follows that geodesic distance is also preserved by these transformations.
Proof The invariance to the finite dimensional groups of rotations and translations follows by the fact that all the terms of the metric are invariant under this action. The invariance under the action of the infinite-dimensional group of reparametrization D follows from an application of the substitution formula for integration.
Having defined the class of H 2 -metrics, we can now formulate the two main building blocks for our framework for the comparison and statistical shape analysis of surfaces: the geodesic boundary problem and the geodesic initial value problem. Geodesic Boundary Value Problem Given two parametrized surfaces q 0 and q 1 , the geodesic boundary value problem consists in finding paths of shortest length that connect the given surfaces q 0 and q 1 , i.e., calculating the geodesic distance between q 0 and q 1 . Here the Riemannian length of a path q : [0, 1] → I is defined as Paths of minimal length are called minimizing geodesics. By a standard result in Riemannian geometry (Lang, 2012), finding minimizing geodesics is equivalent to minimizing the Riemannian energy: In Sect. 4 we will explain how to discretize this functional for discrete meshes, which will in turn allow us to solve the minimization problem using standard finite dimensional optimization methods. Note that the solution of the geodesic boundary value problem gives rise to both optimal (i.e., energy-minimizing) deformations as well as a notion of a distance between the given shapes. Thus, this operation will be the main building block of all our algorithms. Geodesic Initial Value Problem While the geodesic boundary value problem searches for the shortest path between two given surfaces, the geodesic initial value problem searches for the optimal deformation path of a given surface in a given initial deformation direction. Solving the geodesic initial value problem amounts to solving the geodesic equation, which is the first-order optimality condition of the energy functional defined above. In our situation, the geodesic equation will be a non-linear partial differential equation that is of second-order in time and fourth-order in the two-dimensional space coordinates. As this equation is rather lengthy and not particularly insightful, we refrain from formulating it here and instead refer the interested reader to the article , where the geodesic equation is derived for a general class of Riemanniann metrics on I that are induced by abstract pseudo-differential operators and thus include in particular the class of metrics studied in the present work. In addition,  established local well-posedness (existence and uniqueness) of the corresponding (geodesic) initial value problem. To circumvent dealing directly with the intricacies and difficulty of solving highly non-linear and higher-order partial differential equations, we instead calculate the solution to the initial value problem using the methods of discrete geodesic calculus as developed in Rumpf and Wirth (2015b); see Sect. 4 for a detailed description.
In the context of our statistical shape analysis framework for surfaces, the geodesic initial value problem will be of importance for calculating shape averages, for principal component analysis and in our motion transfer applications.

Metrics on Unparametrized Surfaces
In the previous section we introduced a class of Riemannian metrics on the space of parametrized surfaces. Our main goal is, however, to compare surfaces regardless of how they are parametrized. To this end, we introduce the space of unparametrized surfaces, which is defined as the quotient space of parametrized immersed surfaces modulo the reparametrization group, i.e., S = I/D, and refer to it as shape space. This space consists of equivalence classes Since the family of H 2 -metrics on I introduced in (10) is reparametrization invariant, cf. Lemma 1, it induces a corresponding family of Riemannian metrics on the quotient space; such a construction is referred to as a Riemannian submersion, see Bauer et al. (2011) for a detailed explanation in the context of Sobolev metrics on surfaces. Consequently the geodesic distance dist G given in (1) corresponding to an H 2 -metric G descends to a distance function on the quotient shape space which is given as follows: By expanding the expression above, one notes that for given surfaces [q 0 ] and [q 1 ], computing the geodesic distance can be written as the following constrained minimization problem: where the space of paths of immersed surfaces P q 1 •ϕ q 0 is defined in (2). In practice, computing this distance thus requires solving a matching problem that consists of finding the optimal reparametrization and optimal path of immersions between the surfaces. We refer to the constrained minimization problem in (15) as the geodesic boundary value problem on shape space. Compared to the matching problem for parametrized surfaces, the main difficulty in terms of numerically solving this problem consists of discretizing the action of the reparametrization group D. We will circumvent this issue by introducing a relaxed version of (15), which will make use of methods from geometric measure theory, cf. Sect. 3.
While the geodesic boundary value problem on shape space is significantly more challenging than its counterpart on parametrized surfaces, it turns out that the geodesic initial value problem for these two spaces is essentially equivalent: solving the geodesic initial value problem on the space of parametrized surfaces for an initial condition that is in the tangent space of shape space, called a horizontal initial condition, gives rise to a solution in the space of unparametrized surfaces. This observation follows from powerful results in Riemannian geometry and in particular from the conservation law for the horizontal initial momentum, which stems from the reparametrization invariance of the Riemannian metric (Bauer et al., 2011). Consequently, this will allow us to use the same methods for solving the initial value problem on the space of parametrized and unparametrized surfaces, which will be described in Sect. 4.

Remark 3
We can also consider the space of unparametrized surfaces modulo rotations and translation. Since the class of H 2 -metrics is also invariant with respect to these finite dimensional group actions, cf. Lemma 1, it descends to a class of Riemannian metrics on this quotient space. Computing the induced geodesic distance on this quotient space involves an additional minimization over the rotation group Rot(R 3 ) and over the translation group R 3 in addition to minimizing over the reparametrization group and over the space of paths of immersed surfaces.

Relaxed Matching Problem
We now focus our attention on the actual computation of the geodesic distance on the shape space of surfaces. As outlined in the previous section via Eq. (15), this computation involves a joint optimization over paths of immersed surfaces and reparametrizations. In practice, the space of parametrized surfaces I, and hence the path of immersions, can be discretized by considering a triangular mesh as the domain M of the function space I, and considering piece-wise linear functions defined on M, which gives rise to triangulated surfaces, as outlined in Bauer et al. (2021). More general discretizations schemes, such as spline discretizations, could be used as well. Discretizing the paths of immersed surfaces implies that the minimization over those paths can be framed quite naturally as a standard finite dimensional optimization problem, as we will outline in Sect. 4. However, dealing with the minimization over the infinite-dimensional reparametrization group is typically more difficult, and discretizing such a group and its action on surfaces is not straightforward. Recently, an alternative approach was proposed in  where this minimization over reparametrizations of q 1 is dealt with indirectly by instead introducing a relaxation of the end time constraint using a parametrization blind data attachment term. Broadly speaking, this approach consists in considering the relaxed matching problem: where the minimization occurs over paths of immersed surfaces q(·) ∈ C ∞ ([0, 1], I) that satisfy the initial constraint q(0) = q 0 only, and where (q(1), q 1 ) is a term that measures the discrepancy between the endpoint of the path q(1) and the true target surface q 1 , with λ > 0 being a balancing parameter. If we choose a discrepancy term that is independent of the parametrization of either of the two surfaces, then solving the relaxed problem above would yield (q(1), q 1 ) ≈ 0, which yields q(1) ≈ q 1 • ϕ. Thus, this approach allows us to approximate the end time constraint in (15) without the need to explicitly model the reparametrization itself. Furthermore, this relaxed matching framework allows for inexact matching when computing the distance, which will turn out to be crucial when extending this framework to surfaces that can exhibit different topologies, such as surfaces with different genuses, and to surfaces with partial correspondences, as we shall outline in Sect. 6.

Varifold Representation and Distance
We now describe how to construct the key ingredient in the relaxed model outlined above, namely, an effective and simple to compute data attachment term which gives a notion of discrepancy between unparametrized surfaces. Among different possible approaches, we will rely specifically on methods derived from geometric measure theory which have been used for that particular purpose in several past works on surface registration Charon & Trouvé, 2013;Feydy et al., 2017;Roussillon & Glaunès, 2019;Vaillant & Glaunès, 2005), see also the recent survey . In this paper, we adopt the framework of oriented varifolds introduced in Kaltenmark et al. (2017), following an approach similar to the authors' previous works (Bauer et al., 2019aSukurdeep et al., 2022). We point out, however, that the majority of the present work could be adapted without much difficulty to some of the other types of data attachment terms developed in the aforementioned papers. Given any parametrized surface q ∈ I, the varifold μ q associated to q is a positive Radon measure on the product space R 3 × S 2 , where S 2 is the unit sphere. More specifically, μ q is the image measure (q, n q ) * vol q where n q is the unit oriented normal field of q, and vol q is the area form on M induced by q. In other words, for any Borel set A fundamental property is that this varifold representation does not depend on the parametrization of q. Namely, for any ϕ ∈ D, one has μ q•ϕ = μ q , and thus it induces a well-defined mapping on the quotient space S which can be further shown to be an embedding of S into the space of varifolds.
Then, any norm · on the space of varifolds should induce a distance on S, given for any [q 0 ], [q 1 ] ∈ S by μ q 0 − μ q 1 , where we again emphasize that this expression does not depend on the choice of parametrizations for q 0 and q 1 in the respective equivalence classes [q 0 ] and [q 1 ]. While there are many possible metrics that one can introduce on spaces of measures, norms defined from positive definite kernels on R 3 × S 2 have been shown to lead to particularly advantageous expressions for numerical computations. Specifically, following the setting of Kaltenmark et al. (2017), we consider the class of norms · V * , where V is a reproducing kernel Hilbert space of functions on R 3 × S 2 , whose kernel is of the form in which and are two functions defining a radial kernel on R 3 and a zonal kernel on S 2 , respectively. We discuss specific choices for and when presenting our numerical approach for solving the relaxed matching problem in Sect. 4.4. Following from the particular form of μ q 0 and μ q 1 as well as the reproducing kernel property in V , the inner product of the two varifolds in V * can be explicitly derived as: Consequently, the squared varifold kernel distance between μ q 0 and μ q 1 is obtained as follows: It can be shown that, under the right regularity and density assumptions on the kernel k V (c.f. Proposition 4 in Kaltenmark et al., 2017), the above leads to a true distance when restricting to embedded unparametrized surfaces. However, note that there is no notion of geodesics in S corresponding to the varifold distance, as the straight path (1 − t)μ q 0 + tμ q 1 in V * is not associated to a corresponding path in the space of surfaces, due to the non-surjectivity of the mapping q → μ q . Yet, the squared varifold distance μ q 0 − μ q 1 2 V * still provides a valid discrepancy term for the relaxed matching problem in (16), which has the additional important advantage of being simple to discretize and evaluate numerically, as we shall detail in Sect. 4.4.
Lastly, we point out that the above varifold discrepancy metrics are also equivariant to the action of rigid motions. Specifically, for any q 0 , q 1 ∈ I and any R ∈ Rot(R 3 ), τ ∈ R 3 , we have which follows directly from the form of the kernel (17).

Relaxed Surface Matching
The squared varifold distance is ideally suited for use as the discrepancy term in the relaxed matching problem outlined in (16) due to its reparametrization invariance, which finally allows us to formulate the varifold-based relaxed matching problem for surfaces: Given q 0 , q 1 ∈ I, we consider the variational problem: where the minimization occurs over paths of immersed surfaces q(·) ∈ C ∞ ([0, 1], I) that satisfy the initial constraint q(0) = q 0 , and where λ > 0 is a balancing parameter.
Note that the (relaxed) endpoint constraint q(1) ≈ q 1 • ϕ for some ϕ ∈ D is encoded in the varifold attachment term. The interpretation of the two terms in the relaxed energy (20) is as follows: the first term (the energy of the path of immersed surfaces) measures the cost of the optimal deformation, whereas the second term is merely a data attachment term that enforces the endpoint constraint. In this relaxed surface matching framework, we refer to q 0 as the source, q 1 as the target and q(1) as the deformed source.
We note, however, that the model formulated in (20) is asymmetric in the sense that interchanging q 0 and q 1 will affect the obtained minimizer. Although this is a common phenomenon for relaxed optimization problems, we next propose a symmetric formulation of the varifold-based relaxed geodesic boundary value problem. To do so, we will lift the constraint of q(0) being q 0 and instead add a second varifoldbased data attachment term which measures the similarity of q(0) and q 0 : Given q 0 , q 1 ∈ I, we consider the variational problem: where the minimization is performed over paths of immersed surfaces q(·) ∈ C ∞ ([0, 1], I), and where λ 0 , λ 1 > 0 are balancing parameters.
This symmetric formulation of the relaxed matching problem has several advantages which we will leverage in the implementation and simulations presented in the next sections. First, for λ 0 = λ 1 , the variational problem (21) is indeed symmetric in the sense that for any path t → q(t), the time reversed path t → q(1 − t) has the same energy value for the problem of matching q 1 onto q 0 and thus the value of the infimum is the same for both matching problems. More importantly, the introduction of q(0) allows us to decouple the topological or mesh properties of the immersions in the path q(·) with those of the source shape q 0 . As we shall explain more in details in Sect. 4, this allows us to select the vertex sampling and mesh structure of the surfaces in the geodesic path independently of that of the source q 0 , which can be used to adapt the efficient multiresolution scheme of Bauer et al. (2021) for numerically solving the matching problem.

Numerical Optimization Approach
In this section, we will present a set of numerical approaches for solving the geodesic boundary value problem for parametrized surfaces introduced earlier in Sect. 2.1, the varifoldbased relaxed matching problem for unparametrized surfaces introduced in Sect. 3.2, as well as the geodesic initial value problem introduced in Sect. 2. Our source code is openly available at: https://github.com/emmanuel-hartman/H2_SurfaceMatch First, we describe how to discretize parametrized surfaces. We will do so by considering oriented triangulated surfaces, which are also called oriented triangular meshes, that are represented by a set of vertices, edges, and faces. We view the vertices V of a mesh as an ordered set of points in R 3 , i.e., where n is the number of vertices in the mesh. Occasionally we may want to view V equivalently as a single point in R 3n .
The edges E of a triangular mesh are subset of N 2 where (i, j) ∈ E if and only if there is an oriented edge from v i to v j . Similarly, we view the faces F of a triangular mesh as a subset of points in the direction of the oriented normal vector. Canonically, we choose to use only the representative (i, j, k) of a face such that i < j, k.
In the context of the geodesic boundary value problem for parametrized surfaces, the relaxed matching problem for unparametrized surfaces, and the initial value problem, we are required to solve optimization problems over paths of immersed surfaces. In the discrete setting, we will solve these minimization problems by searching over paths of meshes that each lie in a solution space, M, defined as the set of meshes with a fixed combinatorial structure, i.e., the set of meshes with a fixed number of vertices and a fixed set of edges and faces. Thus, each q ∈ M is determined precisely by the locations of the vertices and it is natural to consider M ∼ = R 3n . However, we can equivalently view q ∈ M as a piecewise linear surface determined exactly by the vertices. Therefore we view q as the map q : where for each f ∈ F, σ 2 f is the simplex given by This interpretation of a mesh will prove useful for defining the geometric quantities used in the definition of the H 2 -metric.

The H 2 -Metric on the Space of Triangular Meshes
To establish a numerical framework based on the class of H 2metrics defined in (10), we must first establish a discretization of each of the components that appear in its definition. The field of discrete differential geometry establishes discrete counterparts to smooth geometric quantities such as volume forms, derivatives and the Laplacian. A review of the derivations of these discrete quantities can be found e.g. in Crane (2018). We will either discretize these quantities per face or per vertex of a given triangular mesh depending on the context in which they will be used in our computation.
Recall that q ∈ M is entirely determined by the vertices V ∈ R 3n . Thus, it is natural to discretize tangent vectors on the vertices of the mesh, i.e., a tangent vector h is viewed as a set of vectors in R 3 assigned to each vertex v ∈ V of the mesh. Therefore, Next, we will explain how we discretize the terms that appear in the definition of the H 2 -metric, i.e., the volume form, the pullback metric, the normal vector, and the surface Laplacian. For a graphic explanation of our discretization, we refer to Fig. 4.
Recall that when we view a mesh as a map q, as in (22), it is affine on the simplex corresponding to each face. Therefore, it is natural to discretize the first-order terms on each face. Given a face f ∈ F, where we will assume f = (0, 1, 2) for simplicity of notation, with vertices Given a tangent vector h, we can compute its differential on the face f as Consequently a discrete version of the pullback metric g q , the volume density vol q and the normal vector n q are given by: We denote these discrete versions by g f , vol f and n f to emphasize that they are defined on the faces. Given that the faces are affine, it is somewhat "unnatural" to discretize the Laplacian, a second-order quantity, on the faces of a mesh. Rather, the natural place to discretize the surface Laplacian is on the dual cells of a mesh. Each such Fig. 4 Defining H 2 -metrics using discrete differential geometry. The cell dual to the vertex v is shown in blue dual cell corresponds to a vertex of the mesh; see Fig. 4 for an illustration. Given a tangent vector h = {h v ∈ R 3 |v ∈ V } to a mesh q = (V , E, F) ∈ M, the Laplacian q applied to h at a vertex v ∈ V is given by where α vw and β vw are angles as shown in Fig. 4. This discretization can be derived using either finite element methods as in Crane (2018) or discrete exterior calculus as in Crane et al. (2013).
The zeroth and second-order terms of the metric contain also the volume form of our mesh, which was previously defined for each face. In order to assign this volume form to a vertex v (instead of to the faces), we sum up one third of the volume of each face containing v. Thus, the volume form at a vertex v is given by Thus we have derived discrete versions of all terms that appear in the definition of the H 2 -metric. Therefore, given a mesh q ∈ M and a pair of tangent vectors h, k ∈ T q M we arrive at the following expression for the discrete version of the family of H 2 -metrics:

Discretizing the H 2 Path Energy
Having discussed how to compute the Riemannian metric at a triangular mesh, we now explain how to discretize the Riemannian energy of a path of meshes. Indeed, given a path of triangular meshes in the solution space, denoted by V : whereV (t) denotes the derivative with respect to time of the path. We re-emphasize that each mesh in the path has the same, fixed combinatorial structure, implying that the path is entirely determined by the locations of the vertices of the meshes, hence the notation V above. Furthermore, we note that a further discrete approximation is required to compute the energy of the path, namely we have to discretize the time interval [0, 1]. To that end, we consider piecewise-linear (PL) approximations for paths of meshes. Given a PL path with N + 1 evenly spaced breakpoints 0 = t 0 < t 1 < · · · < t N = 1, we can compute the tangent vector for the first N points via finite differences. Thus for i ∈ {0, 1, ..., N − 1}, we havė As a result, the energy of a PL path in M reduces to

Solving the Geodesic Boundary Value Problem for Parametrized Surfaces
We are now able to formulate our numerical approach for solving the geodesic boundary value problem (BVP) between parametrized surfaces. Given source and target surfaces q 0 , q 1 ∈ I respectively, whose discretized versions are determined by their vertices V 0 and V 1 respectively, our goal will be to approximate solutions to the geodesic boundary value problem in M by minimizing the energy in (23) over all PL paths with fixed endpoints, those being V 0 and V 1 respectively. In doing so, we have reduced the boundary value problem to a finite dimensional, unconstrained minimization problem on R 3n(N −1) ; the free variables being the vertices of the interpolating meshes between V 0 and V 1 . We implement the discrete energy functional (23) using pytorch, which allows us to take advantage of the automatic differentiation functionality to calculate the gradient of this energy with respect to the vertices of the interpolating meshes. We then use the L-BFGS algorithm, as introduced in Liu and Nocedal (1989), to minimize the energy. We describe this process in Algorithm 1 below.

Algorithm 1 Geodesic BVP for Parametrized Surfaces
procedure Parametrized_Geodesic_BVP(V 0 , V 1 , V ) V 0 , V 1 : vertices of the given source and target surfaces V : initial guess for vertices of the interpolating meshes of the PL path To speed up computations (convergence), we implemented a multi-resolution method in time, i.e., we iteratively refine the temporal discretization of the path and repeat Algorithm 1, where we initialize at each iteration with an up-sampled version of the previous solution. An example of a solution to the boundary value problem for parametrized surfaces can be seen in Fig. 5.

Discretizing the Varifold Norm
In order to tackle the varifold-based relaxed matching problem for unparametrized surfaces introduced in (20), we must discuss the discretization of the varifold data attachment term μ q(1) − μ q 1 2 V * introduced in Sect. 3.1. We specifically need to compute the squared kernel distance between the two varifolds μ q(1) and μ q 1 associated to the piecewise linear surfaces given by the two triangular meshes (V (1), E 0 , F 0 ) and (V 1 , E 1 , F 1 ) respectively. The power of the varifold framework is that it applies equally well to this case and allows us to compare discrete shapes with significantly different mesh structures, including those with different topologies.
Indeed, we note that an efficient discretization of the kernel inner product of (18) consists in approximating the integral of the kernel over each pair of faces from F 0 and F 1 respectively by using its value at those faces' centers. In other words, we consider the following approximation: where c f denotes the barycenter of the face f given by c f = 1 3 v|v∈ f v. We emphasize that the quantities c f 0 , n f 0 and vol f 0 are here calculated based on the vertices V (1) of the endpoint of the path of meshes, with the edges and faces E 0 and F 0 being the same as for the initial mesh in the path. The full discrepancy term μ q(1) − μ q 1 2 V * is then once again calculated as in (19), i.e., through a squared expansion of the norm induced by the kernel inner product (18), where each of the inner products is approximated as above. We emphasize that if the two meshes are exactly aligned, then the discrepancy term μ q(1) − μ q 1 2 V * will be minimized, while its value will be larger if the two meshes are highly misaligned.
Although several choices of kernels are available (cf. Charon et al., 2020;Kaltenmark et al., 2017 for more detailed presentations), in all the numerical simulations of this paper, ), a Gaussian kernel of width σ > 0, for the radial kernel on R 3 . The value of this scale parameter σ is typically adapted to the size of the surfaces to be matched. As for the zonal kernel on S 2 , we take (n f 0 · n f 1 ) = (n f 0 · n f 1 ) 2 , which is known as the Cauchy-Binet kernel on the sphere.
Since the calculation of the varifold metric involves a number of kernel evaluations that is quadratic in the number of faces, it typically represents the bulk of the numerical cost of the proposed matching algorithm. For this reason, in our implementation, we rely on the pykeops library (Charlier et al., 2021), which provides efficient GPU routines to compute such large sums of kernel functions and enables the automatic differentiation of those expressions.

Solving the Geodesic Boundary Value Problem for Unparametrized Surfaces
Using the discretization of the H 2 -path energy described in Sect. 4.2 and the discretization of the varifold norm described in Sect. 4.4, we can reduce both the non-symmetric (20) and the symmetric (21) relaxed surface matching problem to a finite dimensional, unconstrained minimization problem. Note, that free variables for the non-symmetric problem are the vertices at time t i for i ≥ 1, while the free variables in the symmetric version include the vertices at time t 0 . The main difference between these two algorithms is, however, that the mesh structure in the non-symmetric version is prescribed by the given data, i.e., the mesh structure (topology) in the solution space is given by the mesh structure (topology) of the source q 0 . In the symmetric version the mesh structure of the solution is a user input and can be different from the mesh structure of both the source and the target. We describe this process below in Algorithm 2. To speed up convergence, we implemented a multi-resolution method in both time and space, i.e., we iteratively refine the temporal discretization of the path and the mesh discretization of the surfaces in the path and repeat Algorithm 2, where we initialize at each iteration with an up-sampled version of the previous solution.

Solving the Initial Value Problem
We now turn our attention to a numerical approach for solving the geodesic initial value problem (IVP) on the space of parametrized surfaces. We re-emphasize, as noted in Algorithm 2 Relaxed Matching for Unparametrized Surfaces procedure Relaxed_Matching(V 0 , V 1 , V ) V 0 , V 1 : triangular meshes for the source and target. V : initial guess for a PL path in M.
Sect. 2.2, that the geodesic initial value problem on the spaces of parametrized and unparametrized surfaces are essentially equivalent. As a result, the procedure we describe in this section gives rise to a solution in the space of unparametrized surfaces as well.
To solve the geodesic initial value problem, we follow the variational discrete geodesic calculus approach developed in Rumpf and Wirth (2015b). Given a surface q ∈ M and a tangent vector h ∈ T q M (which is assumed to be horizontal for unparametrized surfaces) our method involves approximating the geodesic in the direction of h with a PL path V having N + 1 evenly spaced breakpoints. To simplify notation, we will denote surfaces in the PL path at time t i = i N for i = 0, . . . , N by V (t i ) := V i . At the first step, we set V 0 = q and V 1 = q + 1 N h, and find V 2 such that V 1 is the geodesic midpoint of V 0 and V 2 , i.e., we solve for V 2 such that Differentiating with respect toṼ and evaluating the resulting expression at V 1 , we obtain the system of equations where B i is the i-th basis vector of R 3n . We denote the system of equations in (24) by F(V 2 ; V 1 , V 0 ) = 0, where we stress again that V 0 and V 1 are here fixed. We solve this system of equations for V 2 using a nonlinear least squares approach, i.e., by computing via the L-BFGS algorithm, where we again take advantage of the automatic differentiation capabilities of pytorch in our implementation. We then iterate this process step by step to compute V 3 , V 4 , . . . , V N . We summarize our approach via the pseudocode in Algorithm 3. An example of a solution to the initial value problem can be seen in Fig. 5. These results show excellent consistency between the solutions of the corresponding boundary and initial value problems.

Algorithm 3 Geodesic Initial Value Problem
procedure Geodesic_IVP(q, h, N ) q : a surface in M h : a tangent vector in T q M N : number of time steps

Influence of the Metric Coefficients
In this section we present examples detailing the influence of the choice of constants in the H 2 -metric on the geodesics obtained after matching parametrized and unparametrized surfaces via Algorithms 1 and 2 respectively. We also report the influence of the constants on the corresponding computation times, see Table 1.
A synthetic example of a geodesic boundary value problem for a variety of choices of constants can be seen in Fig. 6. We note that the zeroth-order term weighted by a 0 corresponds to the invariant L 2 -metric and penalizes how far the vertices move weighted by their corresponding volume forms. In Fig. 6 on the second row, we see an example where a 0 dominates the other coefficients and as a result, the further a vertex moves, the more shrinking we observe for faces incident to these vertices. The second-order term weighted by a 2 penalizes paths through meshes with high local curvature. In the third line of Fig. 6, we present a path where a 2 is chosen Fig. 6 Influence of constants. An example of the same boundary value problem with different choices for the H 2 -metric coefficients (a 0 , a 1 , b 1 , c 1 , d 1 , a 2 ). First row: (1, 1, 1, 1, 1, 1), second row: (10, 1, 1, 1, 1, 1), third row: (1, 1, 1, 1, 1, 0.1), fourth row: (1, 10, 10, 1, 1, 0.1), fifth row: (1, 1, 10, 0, 1, 10), sixth row: (1, 100, 1, 1, 1, 1) to be small relative to the other coefficients and as a result the midpoints of the geodesic with respect to this choice of coefficients have points with higher local curvature. As noted in Remark 1, the terms corresponding to the weighting coefficients a 1 , b 1 , and c 1 measure the shearing of faces, stretching of faces, and the change in the normal vector, respectively. In the fourth row of Fig. 6, we choose a 1 and b 1 to be large and a 2 to be small. As a result, the geodesic with respect to this choice of coefficients passes through meshes where portions of the pipe are flattened, which produces vertices with higher local curvature without shearing or stretching the faces of the mesh.
In Fig. 7, we highlight the importance of the second order term for complex matching problems. In this figure we consider a matching problem between two surfaces undergoing strong deformations, which in addition, have inconsistent topologies. In previous work of two of the authors , the same example has been considered for the SRNF pseudo-distance: in this framework the obtained result admitted significant singularities in the form of thin spikes that were appearing in areas of high deformations, cf. the (1, 1, 1, 1, 1, 2) (turquoise), an H 2 -metric with coefficients allowing for partial matching: (1, 1, 1, 1, 1, 2) (violet). The target is displayed on the right. One can observe the regularizing effect of the second-order terms (turquoise and violet) and, in addition, how topological inconsistencies (such as the thin arc near the left ear) are correctly removed in the partial matching framework (violet) instead of getting shrunk to almost zero volume (turquoise) yellow skull in Fig. 7. We repeated this experiment using the metrics implemented in this article; in the second figure (the skull in green) one can see the resulting endpoint of an H 1metric. While the resulting match is slightly superior to the one of the SRNF framework, it still exhibits some of the spike singularities. A theoretical explanation for the appearance of these singularities can be found in the observation that the H 1 -metric is not strong enough to control the L ∞ -norm-by the Sobolev embedding theorem the H 1 -metric is exactly at the critical threshold. Consequently small areas can move far with a limited cost, which can potentially lead to these spike type singularities. This observation suggests that this behavior should not occur for higher-order metrics and, indeed, this is also reflected in our experiment: the turquoise skull, which was obtained using an H 2 -metric, does not exhibit any spike singularities and leads to an overall superior matching. Note, that the thin arc in the right ear region of the skull is not a spike, but stems from the inconsistent topology of the shapes, cf. the arc at the right ear of the animal skull. We tackle these topological inconsistencies in Sect. 6, where we introduce a partial matching framework which would automatically erase such regions.

Remark 4
With our relaxed matching framework, one can obtain an adequate matching even when the meshes under consideration are of low quality, i.e., even if they include a certain level of degradation caused by topological or geometric noise, such as the presence of holes or degenerate triangles for instance. Indeed, our approach avoids the need for an exact matching of the source and target meshes (thanks to the varifold relaxation term in our variational formulation of the matching problem), which helps us avoid instances where enforcing an exact matching would lead to e.g. overfitting parts of the source to noisy parts in the target, thus leading to inaccurate results; see the experiments in Sukurdeep et al. (2019) for an illustration in the context of planar curves. Moreover, in our framework, the mesh structure for solutions to the geodesic boundary value problem is userdefined, i.e., meshes in the geodesic path can be prescribed to have any desired topology or resolution (independently of the mesh structure and resolution of the boundary shapes (i.e., the source and/or target)). In the case where the mesh quality of the boundary shapes is low, issues might arise as the varifold term is most sensitive to large discrepancies in the match. Yet, with proper initialization, this issue can be mitigated and we still obtain desirable results. This is depicted in Fig. 7, where we are able to match two skulls despite the presence of topological noise in the data.

Comparison with Other Shape Matching Frameworks
While the previous section outlines how our approach for surface matching compares with other intrinsic Riemannian frameworks based on the SRNF and H 1 -metrics, in this section, we present a (theoretical) comparison of our approach with a wider class of numerical frameworks for surface matching. For an overview on a variety of surface matching frameworks, we refer the interested reader to the survey article (Biasotti et al., 2016). This survey includes, amongst others, an overview of methods based on the metric (measure) space approach for shape matching, where one considers geometric objects (e.g. point clouds or meshes) as metric spaces (possibly with a probability distribution defined on them), and in which objects are compared via Gromov-Haussdorf distances, or via its extensions like the Gromov-Wasserstein distance (Mémoli, 2011). Here, we will focus mainly on extrinsic Riemannian models for shape analysis (Beg et al., 2005;Younes, 2010), the functional map framework (Ovsjanikov et al., 2012;Ren et al., 2018) and the related evolutionary non-isometric geometric matching (ENIGMA) approach (Edelstein et al., 2019) for obtaining shape correspondences, as well as shape interpolation methods based on as-isometric-as-possible or as-rigid-as-possible deformations (Kilian et al., 2007), Hamiltonian dynamics (Eisenberger & Cremers, 2020) or thin shell models (Iglesias et al., 2018), and finally deep learning based approaches for shape registration (Cosmo et al., 2020;Huang et al., 2021;Trappolini et al., 2021). First, we discuss the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework of Beg et al. (2005) and Younes (2010), which stands as the main alternative Riemannian framework for shape analysis. In contrast to the intrinsic metrics considered in this paper, the LDDMM model consists in building shape metrics extrinsically via rightinvariant metrics on a given subgroup of Diff(R 3 ), the group of diffeomorphisms of the ambient space. Then the distance and geodesic between two surfaces is essentially computed by looking for a diffeomorphism of the whole space that warps the source onto the target surface while minimizing the kinetic energy as defined by the metric on Diff(R 3 ). There are several fundamental differences between the intrinsic and extrinsic frameworks, which have already been emphasized in previous publications (Bauer et al., 2019a, b). To give a brief summary of those differences, a first important distinction is that the LDDMM approach imposes more constraints on the regularity of the surface transformation as it must be induced by a smooth deformation of the ambient space itself. A direct consequence is that this approach guarantees diffeomorphic evolution of the source shape and thus prevents the formation of singularities or self-intersection along geodesics, which is in general not the case with the intrinsic H 2 -metric framework (see Bauer et al., 2019a for examples of such phenomenon in the space of curves). On the other hand, geodesics in the diffeomorphic model are only welldefined between surfaces that belong to the same orbit for the action of the specific subgroup of diffeomorphisms and thus relaxing the problem using e.g. varifold distances is a necessity in practice. From a numerical point of view, LDDMM registration is also an optimal control problem, and it is typically solved based on a geodesic shooting scheme (Vialard et al., 2012). The Hamiltonian dynamical equations generally require evaluating kernel functions between all pairs of vertices in the source surface. Thus the complexity for the integration of these systems is quadratic in the number of vertices, which is an important difference with the linear complexity one gets with intrinsic metrics. This implies that for surfaces with a large number of vertices, the complexity of each iteration of the matching optimization scheme is dominated by the computation of the varifold term and its gradient in the intrinsic framework of this paper while it becomes dominated by the integration of the Hamiltonian system in the case of LDDMM. This is illustrated in Table  1 that shows the time per iteration of the optimization algorithm for the different models.
Another popular method for the analysis of unregistered surfaces is the functional map framework (Ovsjanikov et al., 2012;Ren et al., 2018), which allows one to find optimal maps (optimal point-to-point correspondences) between pairs of surfaces by finding optimal pairings between realvalued functions defined on the surfaces. This is done by using a least squares approach to solve a linear system based on the Laplace-Beltrami operator and Wave (or Heat) Kernel Signature descriptors. These quantities describe the local geometry of surfaces and the method is largely successful at matching regions with similar local geometries. However the global matching of the framework benefits significantly from a good prior selection of landmarks, and extensions of the method, such as the evolutionary non-isometric geometric matching (ENIGMA) approach (Edelstein et al., 2019), have been proposed to obtain point-to-point correspondences in a fully automatic way, even in the context of surfaces with different topologies. Moreover, such methods struggle with topological or geometric noise, such as the presence of holes, degenerate triangles or thin spikes in the meshes, which is a difficulty that our intrinsic H 2 -metric framework handles well, as demonstrated in Fig. 7. Furthermore, methods for shape matching using optimal transport techniques have been proposed. In particular, the Gromov-Wasserstein distance for object matching (Mémoli, 2011) treats surfaces as metric-measure spaces and solves for a probabilistic coupling that preserves pairwise distances of points in the metricmeasure spaces. From such couplings, one can produce approximate point-to-point correspondences of metric measure spaces. The computation of pairwise distance matrices limits the effectiveness of these methods for high resolution meshes as these computations are quadratic with respect to the number of vertices. Furthermore, approaches like functional maps, Gromov-Wasserstein, and ENIGMA allow us to obtain approximate point-to-point correspondences, but they do not provide an optimal deformation between the registered shapes. Thus the optimal deformations have to be calculated in a second (independent) post processing step, using e.g. the as-isometric-as-possible, as-rigid-as-possible framework of Kilian et al. (2007), the Hamiltonian dynamics method given by Eisenberger and Cremers (2020) or a thin shell model as presented in Iglesias et al. (2018). Consequently, in this setup, the registration, deformation, and statistical shape analysis are performed separately, which has been shown to be less desirable as it can introduce a significant bias in the resulting statistical analysis (Srivastava & Klassen, 2016). As discussed in Sect. 1, one major advantage of Riemannian frameworks, including intrinsic frameworks like the one presented in this paper or extrinsic frameworks like LDDMM, is that they do not suffer from this shortcoming as the registration, geodesic interpolation and statistical analysis are all performed under the same metric setting.
More recently, several deep learning methods for the registration and analysis of surfaces have emerged (Cosmo et al., 2020;Huang et al., 2021;Trappolini et al., 2021). Such methods are exciting as they can provide significant computational gains by allowing one to register, or even interpolate between surfaces, via a simple forward pass through a pre-trained network, which can be highly desirable in practice especially when working with very large and high dimensional datasets of surfaces. Nevertheless, the quality of point-topoint correspondences or optimal deformations obtained via these deep learning methods relies on having access to a very large database of ground truths to train the network, which in practice is difficult and costly to obtain. As a result of this lack of good training data, these deep learning methods are thus susceptible to having poor generalization capabilities, resulting in situations where the method performs poorly on data that is significantly different or of significantly worse quality than the training data. One potential future application of our intrinsic H 2 -metric framework is that it could be used to generate high quality training data (in the form of geodesic distances, point-to-point correspondences or optimal deformations) for deep learning methods for the analysis of surfaces. Such ideas have recently been introduced in the case of functional data (Chen & Srivastava, 2021;Nunez et al., 2021) and in the setting of planar curves (Hartman et al., 2021;Nunez et al., 2020), where early results have been encouraging.

Statistical Shape Analysis of Surfaces
Beyond the comparison of two surfaces, the mathematical and numerical framework developed in the previous sections can be used as building blocks for the development of more general tools for the statistical shape analysis of sets of surfaces. In this section, we discuss in particular how to extend our approach to calculate sample averages, perform principal component analysis, and approximate par-allel transport between parametrized and unparametrized surfaces. Our implementation of these statistical shape analysis methods are available in our open source code 1 .

Karcher Mean
A central tool in any statistical shape analysis toolbox is the notion of a Karcher mean. Let (M, G) be a (possibly) infinite-dimensional Riemannian manifold with corresponding geodesic distance function dist G . Given data x 1 , . . . , x K ∈ M, the Karcher meanx is the minimizer of the sum of squared distances to the given data points, i.e., Note that the existence and uniqueness of the Karcher mean is a priori not guaranteed, but requires that the data is sufficiently concentrated, i.e., belongs to a ball in the geodesic distance whose radius depends on the curvature of the manifold M. The Karcher mean can be computed by a gradient descent method, as proposed e.g. in Pennec (2006). This method requires the computation of K geodesic boundary value problems at each gradient step, where usually a relatively large number of iterations (gradient steps) is necessary. For computational efficiency, we instead implemented an alternative algorithm to approximate the Karcher mean based on the iterative geodesic centroid procedure as proposed e.g. in Ho et al. (2013): Given data points x 1 , . . . x K ∈ M and an initial guess x 0 , we generate a sequence of estimates for the Karcher mean, namelyx i for i = 0, . . . , N iter where N iter = O(K ), by settingx 0 = x 0 , and iteratively defininĝ x i = x(1/(i + 1)), with x(t) being the geodesic connectinĝ x i−1 to a data point x k which has been uniformly chosen at random (with replacement) from the dataset. Thus one only has to calculate N iter = O(K ) geodesics in total, which is linear in the number of data points.
A pseudo-code of this method is presented in Algorithm 4. For parametrized surfaces, we can initialize this algorithm with the Euclidean mean (average) of the vertices of the surfaces in our sample, assuming of course that they have been centered. We then iteratively solve the geodesic boundary value problem for parametrized surfaces using Algorithm 1, whose inputs at the ith iteration are the current Karcher mean estimateV as the source, the randomly chosen surface V k as the target, and a linearly interpolated path between the source and target as initial guess for the PL path V . For unparametrized surfaces, there is the additional difficulty that the data might have inconsistent mesh structures. In order to extend the computation of the Karcher mean to this situation, one needs to initialize the Karcher mean estimate (V ,Ē,F) to some user-defined mesh, which will determine the connectivity and topology of the Karcher mean, and then iteratively solve the relaxed matching problem for unparametrized surfaces using Algorithm 2. Note that, as inputs for the relaxed matching problem at the i th iteration, we can use the current Karcher mean estimate (V ,Ē,F) as the source, the randomly chosen surface (V k , E k , F k ) as the target, and a constant path of the Karcher mean estimate for the initial PL path V . An example of a population of unparametrized shapes can be seen in Fig. 10, together with their Karcher mean which has been computed via Algorithm 4.

Algorithm 4 Karcher Mean
. . , V K : triangular meshes from the samplē V : initial guess for Karcher mean

Dimensionality Reduction
Dimensionality reduction is a key tool in modern statistics and machine learning. We illustrate how to construct two popular dimensionality reduction tools for the statistical shape analysis of surfaces using our framework, namely data visualization through multidimensional scaling, and principal component analysis.

Visualizing the Distance Matrix Using Multidimensional Scaling
Multidimentional scaling (MDS) is a well-known procedure for mapping K points in a high (or infinite) dimensional space into a lower dimensional space, while maintaining information about the pairwise distances between these K points.
More specifically, given a (possibly) infinite-dimensional Riemannian manifold (M, G) with corresponding geodesic distance function dist G , with data x 1 , . . . , x K ∈ M, the goal of MDS is to find pointsx 1 , . . . ,x K ∈ R d for some d > 0 Fig. 8 Visualizing the distance matrix between ten human body shapes using multidimensional scaling. The geodesic distance naturally clusters the population into male and female shapes such that: In the context of statistical shape analysis, one can use MDS to project a dataset of surfaces as points in Euclidean space for data visualization purposes, or as an intermediary step in clustering applications with sets of surfaces, see Fig. 8.

Tangent PCA
Principal component analysis (PCA) is an important dimensionality reduction technique in statistics for analyzing the variability of data in Euclidean space. More precisely, given data points x 1 , . . . , x K ∈ R d having zero mean, the goal of PCA is to produce a sequence of linear subspaces {W } d =1 that maximizes the variance of the data when it is projected onto those subspaces (Fletcher et al., 2004). This sequence of linear subspaces W = span({w 1 , . . . , w }) for = 1, . . . , d are constructed by finding an orthonormal basis {w 1 , . . . , w d } of R d which can be computed as the set of ordered eigenvectors of the sample covariance matrix of the data. Thus, PCA amounts to finding the vectors {w } d =1 , which are called the principal components of the data.
Extending PCA to manifolds, even in a finite dimensional setting, is not straightforward nor canonical due to the dif-ference in tangent space at each point of the manifold. As a result, several different models and heuristics have been proposed for manifold PCA. Among those, tangent PCA (Fletcher et al., 2004) is probably the simplest as it relies on directly linearizing the problem around a single point (the Karcher mean). More specifically, let (M, G) be a (possibly) infinite-dimensional manifold. Consider data points x 1 , . . . , x K ∈ M, and a reference pointx ∈ M, for which a natural choice is e.g. the Karcher mean of the data points. The goal of tangent PCA is to find a set of principal component geodesics for the data. By principal component geodesics, we mean a set of geodesics all starting atx whose initial velocities are given by tangent vectors {w } d =1 ∈ Tx M that are computed as the principal components of the data in the linear space Tx M. Thus, tangent PCA amounts to performing standard PCA in Tx M, which can be interpreted as finding the "principal tangent vectors" for the data, i.e., the initial velocities which uniquely determine the geodesics starting at the reference pointx along which one has to move on M in order to maximize the "variability" of the data.
We implemented an algorithm for performing tangent PCA when given a set of K surfaces and a reference point, with details given in Algorithm 5. Our method consists of solving K geodesic boundary value problems via Algorithm 1 (for parametrized surfaces) or Algorithm 2 (for unparametrized surfaces, resp.), using the reference point as the source and each surface in our dataset as respective targets. This produces K geodesics, which we use to estimate K tangent vectors {h k } via finite differences, i.e., by taking the difference between the vertices in the geodesic paths at the first two time points. We then perform PCA on these tangent vectors with respect to the metric GV at the reference point. This is specifically done by computing We then recover the principal component vectors w = K k=1 v ,k (h k −h) and the principal component geodesics by solving initial value problems starting atV in the direction of λ w using Algorithm 3. Note that we solve these IVPs in the positive and negative principal directions ±λ w respectively. While we only write the pseudocode for tangent PCA on parametrized surfaces in Algorithm 5, the method works verbatim for the case of unparametrized surfaces, except that the relaxed matching algorithm (Algorithm 2) is used to solve the K geodesic BVPs.
To illustrate the effectiveness of tangent space PCA, cf. Figs. 9 and 10, we first display the principal component geodesics for an unparametrized dataset of surfaces in Fig. 10. As a second, more large scale experiment, we analyze the faces of the CoMA dataset (Ranjan et al., 2018). As this data comes with known point correspondences, we are able to interpret the data as parametrized surfaces. To evaluate our method we separate the data into a testing set Fig. 9 Tangent PCA for a set of parametrized surfaces. On the left we display the first three principal component geodesics of a training set. On the right, we display a reconstruction of two elements from a separate testing set, where each vertex is colored based on the Euclidean error of the reconstruction of ∼ 2000 meshes and a training set of ∼ 700 meshes. In Fig. 9, we illustrate the principal component geodesics of the training set computed using our method. To reconstruct a target mesh, we then perform an unparametrized geodesic matching from a template to the target with respect to the first 40 tangent PCA basis vectors. In particular, we optimize the relaxed matching energy over all paths where the tangent vectors of the path can be written as a linear combination of the tangent PCA bases. In Fig. 9, we also display such a reconstruction of two surfaces from the testing set. When we reconstruct the entire testing set in this way we achieve 75% of all vertices within a Euclidean error of 1mm. For comparison, the percentage of vertices within 1mm accuracy is 47% when using traditional PCA and 72% when using the Mesh Autoencoder methods of Ranjan et al. (2018).

Parallel Transport
Parallel transport is a method of transporting geometric data (tangent vectors) between different points in a manifold. In our situation this concept has a natural application to motion transfer, as shown in Fig. 11. Given a geodesic (i.e., a motion) between a source and target surface (e.g. the two cats in Fig. 11), we can transfer the motion to a new source shape (e.g. the lioness in Fig. 11) by parallel transporting the initial velocity from the geodesic motion to the new source shape, and then solving an initial value problem starting at the new source shape with initial velocity given by the parallel transported tangent vector. This procedure requires that we approximate parallel transport of tangent vectors on M. We use an implementation of Schild's ladder to produce a firstorder approximation of parallel transport (Guigui & Pennec, 2021;Kheyfets et al., 2000). Given a Riemannian manifold M, with x 0 , x 1 ∈ M, h ∈ T x 0 M, and letting V be a geodesic such that V (0) = x 0 and V (1) = x 1 , the calculation of parallel transport using Schild's ladder requires one to iteratively compute several small geodesic parallelograms with one side corresponding to a small step along V and the other side being a small step in the direction of h. The transport of h for this small step along V is defined to be the log map of the side opposite of h. One then repeats the computation of these rungs until reaching x 1 . An algorithmic explanation of this method is given in Algorithm 6 below.
Algorithm 6 Parallel Transport procedure Parametrized_Parallel_Transport(V , h, N ) V : geodesic to transport the tangent vector along h : tangent vector to be transported N : number of iterations for Schild's ladder

Partial Matching
In this final section, we further extend the second-order elastic surface analysis framework introduced in the previous sections by augmenting the surface matching model with the estimation of spatially-varying weights on the source shape. As we will show, this approach will enable us to compare and perform statistics on sets of surfaces which may have incompatible topological properties or exhibit partially missing data.

Limitations of the Previous Framework
We start by motivating the need for this extended approach. Indeed, the relaxed matching framework presented so far in (21) (as well as its non-symmetric version (20)) is pri- Fig. 11 Example of parallel transport using Schild's ladder. We compute the initial tangent vector in the direction of the top geodesic, use Schild's ladder to transport the tangent vector along the geodesic between the leftmost surfaces, and finally compute the geodesic on the the bottom as an IVP. Animations of the obtained motion transfer can be found in the supplementary material and on the github repository marily designed for the comparison of complete surfaces with consistent topology, as illustrated by the examples in Fig. 12.
Although the matching obtained from (21) is inexact and may in practice be able to handle small inconsistencies including topological noise, it remains ill-suited for datasets involving surfaces with significant missing parts or important topological differences (either artifactual or not). Attempting to compare two such surfaces based on model (21) is likely to lead to highly singular behaviour in the estimated geodesics and distances. This was already emphasized in the case of planar shapes (such as curves and shape graphs) in the authors' previous publication (Sukurdeep et al., 2022), and can be further observed in the case of 3D surfaces, as seen e.g. in Fig. 7 with the formation of geometric artifacts such as the thin arc around the ear of the skull, and in Fig. 12 with phalanges that shrink to almost zero volume.
To address this shortcoming in our model, we propose to incorporate partial matching capabilities in our framework. Extending the idea introduced in Sukurdeep et al. (2022), we do so indirectly by considering surfaces augmented with a weight function defined on their support, leveraging the flexibility of the varifold representation for that purpose. This will lead to a new matching formulation between pairs of weighted surfaces, where, in combination to the geometric matching process, one can vary the weights assigned to different components or parts of the source surface. In particular, this allows us, by setting weights to 0 in specific areas, to remove parts of the source when they have no corresponding parts in the target surface, as shown in Figs. 12 and 13.

The Varifold Norm on the Space of Weighted Surfaces
We first define a parametrized weighted surface as a couple (q, ρ), where q ∈ I is a parametrized surface as previously defined and ρ : M → [0, 1] is a function on the parameter space M. For each (u, v) ∈ M, one can interpret ρ (u, v) as the weight assigned to the point q(u, v) on the surface. The primary reason to assume the values of ρ in the interval [0, 1] is that we are focusing on the issue of partial matching. In such a scenario, it is indeed natural to impose this constraint, with the interpretation being that the weight function to be estimated in the matching problem introduced below should vanish on parts of the transformed surface that need to get erased to adequately match the target, while remaining roughly equal to 1 on the other parts. Note that in other situations such as shapes with multiplicities, one could consider more general R + -valued weight functions. Any such weighted surface (q, ρ) can still be represented as the varifold that we write μ q,ρ := ρ · μ q , which is defined as the image measure (q, n q ) * (ρ vol q ), where n q is the unit oriented normal field of q, as defined earlier in Sect. 2.1, and ρ vol q is the area form on M induced by q rescaled by the weight function ρ. With this definition, the kernel metrics on varifolds outlined earlier in Sect. 3.1 immediately induce a fidelity metric between weighted surfaces. Specifically, the kernel inner product in V * between two weighted varifolds μ q 0 ,ρ 0 and μ q 1 ,ρ 1 is given explicitly by: where we have dropped the coordinates (u 0 , v 0 ) and (u 1 , v 1 ) in the above expression for concision. This simply amounts to a weighted version of (18). The squared weighted varifold kernel distance μ q 0 ,ρ 0 − μ q 1 ,ρ 1 2 V * can again be obtained via a quadratic expansion, exactly as in (19).

Relaxed Surface Matching with Weights
We are now able to formulate the extension of the matching problem of Sect. 3.2 to weighted surfaces: Given a pair of weighted surfaces (q 0 , ρ 0 ) and (q 1 , ρ 1 ), we consider the variational problem: where the infimum is taken over paths of immersed surfaces q(·) ∈ C ∞ ([0, 1], I), and also over all weight functions ρ : M → [0, 1], with λ 0 , λ 1 > 0 being balancing parameters.
In this framework, we refer to (q 0 , ρ 0 ) as the source, (q(1), ρ) as the transformed source, and (q 1 , ρ 1 ) as the target. Note that in addition to relaxing the end time constraint, we have also relaxed the initial constraint of q(0) being q 0 via a second varifold fidelity term in the model above. Similarly to what was explained in the context of (21), this allows us to choose the topological and mesh properties of the path q(·) independently of those of the source q 0 , once again paving the way for the use of the efficient multiresolution scheme from Bauer et al. (2021) to numerically solve this matching problem between weighted surfaces. We note that one could also formulate an asymmetric version of problem (27) by instead enforcing the initial constraint q(0) = q 0 as in (20). Furthermore, aside from the constraint of ρ taking its values in [0, 1], the variational problem (27) does not involve any cost penalty on the weight function. Yet it would be possible to add regularizers for the weight function to the functional, including for instance the total variation norm of ρ − ρ 0 as done in Sukurdeep et al. (2022) for planar shapes so as to promote piecewise constant weight functions. In the context of partial matching, it may also be relevant to enforce ρ to take values close to 0 or 1, which can be achieved e.g. by adding a double well pointwise penalty of the form M (ρ(u, v)(ρ(u, v) − 1)) 2 vol q (u, v). A clear downside to including extra regularizers is the added layer of complexity to the matching model due to the presence of extra terms and balancing parameters. For that reason, we decided to focus this work on the above unpenalized formulation.

Remark 5
We emphasize that in (27), we only allow for weight variations on the transformed source, which lets us model in particular the erasure of parts of q(1) so as to match the target. This is useful e.g. in the context of partial matching problems with missing data in the target shape, see Fig. 12. One can easily adapt the model to allow for weight estimation on the target by minimizing over a weight functionρ defined on q 1 , with the weights on the transformed source shape being kept fixed. More generally, one could technically model weight variations on both shapes, by jointly optimizing over two weight functions ρ andρ. However, the latter case requires careful regularization on those functions in order to prevent the trivial solution of setting all weights to 0. We will thus leave the study of this case to future work.

Numerical Optimization with Weights
We now discuss our approach for numerically solving the matching problem between weighted surfaces, whose discretization can be performed in similar fashion as previously. A discrete weighted surface (q, ρ) is once again represented as a triangular mesh (V , E, F) as in Sect. 4, while the weight function ρ shall be modelled by its discrete set of values at the center c f of each face f ∈ F of the mesh, i.e., by the vector in [0, 1] |F| with entries ρ f := ρ(c f ).
Then, letting (V , E, F, ρ) and (Ṽ ,Ẽ,F,ρ) denote the discretizations of two weighted surfaces (q, ρ) and (q,ρ), we can first approximate the varifold inner product: where n f , nf and vol f , volf are the unit normals and volume forms that have been discretized over the faces f ∈ F and f ∈F of the meshes, as outlined in Sect. 4.1. The full varifold fidelity term μ q,ρ − μq ,ρ 2 V * is then obtained as in (19), via the quadratic expansion of the squared norm.
Equipped with the discretizations of the H 2 -path energy described in Sect. 4.2, of the varifold norm described in Sect. 4.4, and of the weighted varifold norm described above, we are led to numerically solve (27) as finite dimensional optimization problem, where the minimization occurs jointly over the vertices of the discretized piece-wise linear path of meshes V : [0, 1] → M and over the discretized weight function ρ ∈ [0, 1] |F| . In order to deal with the box constraints on the values of ρ, we minimize the discretized matching functional using the bound constrained limited memory BFGS (L-BFGS-B) algorithm (Byrd et al., 1995), whose implementation is available through scipy. We summarize the weighted surface matching approach in Algorithm 7 below.

Algorithm 7 Relaxed Matching for Weighted Surfaces
procedure Weighted_Matching((V 0 , ρ 0 ), (V 1 , ρ 1 ), V , ρ) V 0 : triangular mesh for the source ρ 0 : weights on the source V 1 : triangular mesh for the target ρ 1 : weights on the target. V : initial guess for a PL path in M. ρ : initial guess for weights on the transformed source

Partial Matching Experiments
To illustrate the capabilities of the weighted surface framework for partial matching, we performed several numerical experiments. In all the figures, we compute the linear interpolation (1 − t)ρ 0 + tρ between the initial and estimated weight function and show this interpolated weight function along the geodesic through a transparency map in order to highlight in a more visual way the effect of weight variations.
First, we demonstrate the benefits of weight estimation when comparing surfaces with missing parts, see the last row of Fig. 1 where we perform a matching with a partially-observed femur bone, and Fig. 12 where we use an incomplete set of phalanges. Partially observed or incomplete data is a common occurrence in practice and can be due to several factors, including segmentation issues, inconsistent field of views or occlusions during the data acquisition process. Typically, matching surfaces with missing parts using standard elastic surface matching techniques will result in the transformed source getting bent, stretched or compressed in an attempt to fill in some of those missing parts. This can result in unnatural deformations (see the fairly extreme Fig. 12 Matching with missing data. We use a complete set of phalanges (i.e., hand bones) as the source, and a different set of phalanges as the target, where some bones on the index finger and thumb were artificially removed. Top row: We matched the surfaces without weight estimation using Algorithm 2. The parts of the transformed source that are getting matched to the removed bones from the target get shrunk to almost zero volume. The estimated geodesic distance is 117.006. Bot-tom row: We augment the surfaces with weights and use Algorithm 7 to match them. Our model correctly "erases" (i.e., estimates vanishing weights) the appropriate parts of the transformed source to account for the corresponding missing bones on the target. This produces a natural looking geodesic between the source and target, without the production of singularities, with a lower estimated geodesic distance of 114.564

Fig. 13
Splitting into multiple components. We match a single sphere with two disconnected spheres using Algorithm 7. The transformed source q(1) contains a "bridge" between the two spheres in the target where the algorithm estimates zero weights shrinking of phalanges in Fig. 12) and in turn in an overestimation of the geodesic distance between these surfaces, making any subsequent statistical shape analysis spurious for datasets of partially observed surfaces. As evidenced by the second row of Fig. 12, the proposed approach overcomes this difficulty through the automatic estimation of vanishing weights at the location on the parts of the source shape corresponding to the missing ones of the target.
Second, we computed the matching of surfaces with completely different topologies, as shown with the example of the skulls in Fig. 7 and the synthetic examples from Figs. 13 and 14. The estimation of vanishing weights indirectly allows us to recover several useful types of transformations that are otherwise not achievable in the original model. For instance, it enables the model to erase the thin arc near the left ear of the turquoise skull in Fig. 7 as opposed to geometrically shrinking it. It further allows for the splitting of a surface into several connected components, as shown in Fig. 13, as well as the creation of holes when matching surfaces with different genuses, as illustrated by Fig. 14. It should be noted, however, that this approach does not directly model topological changes in the mesh of the transformed source (which remains the same along the geodesic), but rather allows us to compare objects with different topologies by erasing parts of the transformed source via the weight function.

Fig. 14
Matching with highly inconsistent topological structures. We match a sphere (genus zero surface) and a torus (genus one surface) via Algorithm 7. Our model artificially accounts for the creation of a hole, i.e., the change in topology, via the estimation of vanishing weights Lastly, we considered a case of Karcher mean estimation under partial observations. As a proof of concept, we computed the Karcher mean of a set of five distinct hands, each missing a different finger which was artificially removed, see Fig. 15. Following the same principle as the algorithm for Karcher mean estimation presented in Sect. 5.1 (Algorithm 4), we applied Algorithm 7 to iteratively solve weighted matching problems from the current Karcher mean estimate to a randomly chosen surface from our dataset. As the initial guess for the Karcher mean, we used a complete hand (i.e. a closed mesh with five fingers) from a different subject. While other choices for the initial guess, e.g. an ellipsoid, are possible, poorly chosen initializations will result in slower convergence to the Karcher mean and potentially to a lower mesh quality of the estimated Karcher mean. The joint estimation of weights at each successive matching prevents the geometric shrinking of one of the fingers and ultimately results in the realistic looking Karcher mean displayed in Fig. 15. In Fig. 16, we also show the computed geodesics from the Karcher mean to each subject. We also report the Riemannian energy of the geodesic path for each of these geodesics in Table 2. As a point of comparison, we also ran the non-weighted Algorithm 2 between the Karcher mean estimate and each of the corresponding complete hands (i.e., without the artificially removed fingers), and report the Riemannian energy of the resulting geodesics in the last column of Table 2. We observe that the geodesic distance estimates  Table 2 are comparable and quite consistent in both scenarios. This point highlights the reliability of the distance estimates obtained with weight estimation and hints at the potential viability of this approach for statistical shape analysis of datasets of partially observed surfaces.

Data Sources
The mesh data for our numerical simulations in this paper was obtained from several sources, including the meshes made available by Robert Sumner and Jovan Popovic from  (Bronstein et al., 2008), the dynamic FAUST dataset (Bogo et al., 2017), the CoMA dataset (Ranjan et al., 2018), and the dataset of faces from Vlasic et al. (2004).

Conclusion and Future Work
In this paper, we introduced a mathematical framework and several numerical algorithms for the estimation of geodesics and distances induced by second-order elastic Sobolev metrics on the space of parametrized and unparametrized surfaces. We leveraged our surface matching algorithms to develop a comprehensive collection of routines for the statistical shape analysis of sets of 3D surfaces, which includes algorithms to compute Karcher means, perform dimensionality reduction via multidimensional scaling and tangent PCA, and estimate parallel transport across surfaces. We also proposed to resolve the issue of partial matching constraints in the situation of missing data and inconsistent topologies through the additional estimation of a weight function defined on the source shape. A more in-depth quantitative evaluation and comparison of our method against other approaches for partial shape matching, such as Cosmo et al. (2016) and Rodolà et al. (2017) is left for future work. Finally, we want to mention several limitations of the method presented in this article. First, parameter selection may be an important issue in this framework if no reasonable priors are available for the choice of the H 2 -metric coefficients or the kernel scale used to compute the varifold relaxation term. As we illustrated in the numerical experiments, those can all have significant influence on the quality of registration and on the behavior of geodesics. While we use different practical strategies such as multiscale schemes to mitigate this issue, a subject of active current investigation is precisely to develop effective approaches to obtain such parameter estimates in a data-driven way, see (Bauer et al., 2022a) for a recent preprint on this topic in the context of elastic metrics on curves.
A second potential limitation is the choice of initialization for the geodesic path, as the variational problems we tackle are non-convex. In this work, we typically initialize algorithms using a time constant path with either the source or target mesh. However, one could expect extra robustness if more adapted initializations are chosen, which could be computed for instance as the output of some other fast surface matching procedure.
Lastly, while the computational cost of our approach is favorable when compared to other Riemannian frameworks for shape analysis, it still involves running an optimization procedure with quadratic complexity at each iteration. As a result, the numerical pipelines of this paper might become somewhat impractical when working with high resolution meshes (e.g. O(10 6 ) vertices or more) or for populations with a large number of subjects. One way to overcome this issue would be to leverage deep learning architectures, similar to related work in the case of functional data (Chen & Srivastava, 2021;Nunez et al., 2021) or planar curves (Hartman et al., 2021;Nunez et al., 2020). This would reduce the computation of quantities such as distances and geodesics to a simple forward pass through a neural network trained from supervised data obtained from our intrinsic H 2 -metric framework. Conversely, many existing neural network methods applied to surface mesh processing, such as the autoencoder models of Cosmo et al. (2020), Ranjan et al. (2018) and Huang et al. (2021), attempt to learn a latent space representation of the surface dataset using various metric priors as regularization. We also plan to investigate in the future the effectiveness of second-order Sobolev metrics to regularize latent space representation learning of mesh autoencoders.
Funding Open access funding provided by Austrian Science Fund (FWF).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.