Non‐Isometric Shape Matching via Functional Maps on Landmark‐Adapted Bases

Abstract We propose a principled approach for non‐isometric landmark‐preserving non‐rigid shape matching. Our method is based on the functional map framework, but rather than promoting isometries we focus on near‐conformal maps that preserve landmarks exactly. We achieve this, first, by introducing a novel landmark‐adapted basis using an intrinsic Dirichlet‐Steklov eigenproblem. Second, we establish the functional decomposition of conformal maps expressed in this basis. Finally, we formulate a conformally‐invariant energy that promotes high‐quality landmark‐preserving maps, and show how it can be optimized via a variant of the recently proposed ZoomOut method that we extend to our setting. Our method is descriptor‐free, efficient and robust to significant mesh variability. We evaluate our approach on a range of benchmark datasets and demonstrate state‐of‐the‐art performance on non‐isometric benchmarks and near state‐of‐the‐art performance on isometric ones.


Introduction
A common scenario in shape matching is that of very sparse userprovided landmark correspondences that need to be extended to a full map between the considered shapes.The landmarks in question are often of a semantic nature, and thus are very sensitive to exact placement.Consider, for instance the position of the eyes or the nose on a human face (see Fig. 1) that are matched by an artist, e.g., in a texture transfer scenario.In such cases, it is crucial to preserve the landmark correspondences exactly when extending the map.Furthermore, it is desirable for the extension process to be time-efficient and applicable to general, possibly non-isometric shape pairs.submitted to COMPUTER GRAPHICS Forum (6/2022).

arXiv:2205.04800v2 [cs.CV] 22 Jun 2022
Functional map methods [OCB*17] constitute a highly effective shape matching framework, especially when coupled with powerful recent post-processing tools such as ZoomOut and its variants [MRR*19; HRWO20].The existing methods, however, suffer from two major limitations: first, they heavily rely upon the assumption of near-isometry, and second, they typically formulate landmark correspondence via descriptor preservation objectives, combined with other regularizers in the least squares sense.Unfortunately, this implies that the final map is not guaranteed to preserve userprovided landmark correspondences.
In this paper, we propose a novel approach that maintains the efficiency and flexibility of the functional maps pipeline, while overcoming these drawbacks.We organize our proposal in three major stages.First, we introduce a novel functional basis in which to express our map.Crucially, our basis is explicitly adapted to the landmark correspondences, unlike the commonly-used general Laplace-Beltrami eigenbasis.Intuitively speaking, this allows us to enforce landmark preservation by only considering functional maps with a particular (block-diagonal) structure.The design of this landmark-adapted basis is the most technically involved part of our proposal, and relies on solving intrinsic Dirichlet-Steklov and Dirichlet Laplacian eigenproblems.Specifically, we first construct new boundaries at the landmarks, and then formulate and solve the associated boundary value problems.
Second, we remove the assumption of near-isometry by structuring shape matching as a search for bijective near-conformal maps, which are significantly more general than isometries.Following the functional maps pipeline, we express this as a carefully designed energy to be minimized.
Third, we propose an iterative minimization strategy for our energy by following in the footsteps of ZoomOut [MRR*19].In particular, we demonstrate how landmark correspondences can be promoted throughout this iterative refinement.Furthermore, we exploit the landmark-awareness of our basis to provide a simple initial guess of the correspondence.
We test our approach on various benchmark datasets, both isometric and non-isometric.We compare our results to both stateof-the-art functional maps approaches, as well as recent methods that exactly preserve landmark correspondences.We report stateof-the-art accuracy on non-isometric datasets and near state of the art on isometric ones.Meanwhile, the computation time of our approach is significantly lower than that of the competing landmarkpreserving methods.

Contributions. To summarize:
1. We introduce a novel landmark-dependent functional basis by solving an intrinsic Dirichlet-Steklov eigenproblem.2. We formulate a functional maps-based approach to nearconformal shape matching that preserves given landmarks exactly without restrictions on the topology of the shapes.3. We propose an efficient way to both compute the basis and to solve the shape matching problem and report state-of-the-art results on difficult non-isometric benchmarks.

Related Work
Non-rigid shape matching is a well-established research area with a rich history of solutions.Below we review the works that are most closely related to ours, focusing on functional maps and landmarkpreserving methods, and refer the interested readers to recent surveys [vKZHC11;Sah20] for a more in-depth discussion.
Functional Maps Our approach fits within the functional maps framework that was originally introduced in [OBS*12] and extended in many follow-up works, including including [KBB*13; ADK16; RCB*17; EB17; BDK17; NO17; MRR*19; RMOW20] to name a few.An early overview of many functional maps-based techniques is given in [OCB*17].The key idea exploited in all of these techniques is to represent correspondences as linear transformations across functional spaces, which can be compactly encoded as small-sized matrices given a choice of basis.This leads to simple optimization problems that can accommodate a range of geometric objectives such as isometry [OBS*12], accurate descriptor preservation [NO17], bijectivity [ERGB16], orientation preservation [RPWO18] or even partiality [RCB*17] among others.Typically, such objectives are formulated as soft penalties on the functional map and are optimized for in the least squares sense.
Landmarks in functional maps Landmark constraints are commonly used in functional maps-based approaches, especially in an attempt to resolve symmetry ambiguity, present, e.g., when mapping between human shapes.Starting from the segment correspondences advocated in the original approach [OBS*12], and exploited in follow-up works, e.g., [KO19], several techniques also used pointwise landmarks, that were either user-specified [NO17], automatically computed [MMRC18], or even extended to curve constraints [GBKS18].All of these techniques, however, formulate landmark correspondences via functional descriptor preservation, e.g., based on the heat kernel [SOG09; OMMG10] or wave kernel maps [ASC11], which are enforced during optimization only in a least squares sense, alongside other descriptors and regularizers.Therefore, there is no guarantee that the final recovered point-topoint map will satisfy these user-constraints.In contrast, our approach is geared towards preserving the landmark correspondences exactly, while computing a smooth overall map.
Landmark-based matching Landmark-preserving shape correspondence has also been studied in other matching frameworks.
Early methods relied on extrinsic shape alignment, under given constraints, e.g., using thin plate splines [Boo89;CR00] or by extending non-rigid ICP, as done in [SP04] among others.Such approaches, however, rely strongly on the shape embedding and often require a significant number of landmarks to work well in practice.
Another successful class of approaches have aimed to compute correspondences by embedding shapes to a common parametrization domain.This includes powerful approaches based on mapping surfaces to the planar domain, [APL14; WZ14], Euclidean orbifolds [AL15] general flat cone manifolds [APL15] or, more recently, the hyperbolic plane [AL16], which can accomodate an arbitrary number of landmarks.
Finally, recent techniques have also allowed landmarkpreserving shape correspondence by cross-parametrizing the surfaces directly.This includes exploiting direct and inverse averages submitted to COMPUTER GRAPHICS Forum (6/2022).
on surfaces [PBDS13] or finding maps that minimize various notions of distortion, e.g., harmonicity and reversibility (using, first a surrogate high-dimensional embedding) [ESB19] or a related symmetric Dirichlet energy [SCBK20] (via direct optimization on the surface).These recent techniques can lead to accurate results, but are often computationally expensive, and typically place restrictions on the topology of the shape pair, such as having the same genus.In contrast, our method does not suffer from this limitation, as topological stability is one of the features of functional map methods, which is also inherited by our technique.
Basis selection for functional maps Finally, we remark that our construction of landmark-adapted functional bases also fits within the functional map framework, aimed at developing flexibile and effective basis functions.The original article and most follow-up works [OCB*17] have advocated using the eigenfunctions of the Laplace-Beltrami operator (LBO), which are optimal for representing smooth functions with bounded variation [ABK15].However, the Laplace-Beltrami basis has global support and may not be fully adapted to non-isometric shape changes.
The compressed manifold modes [NVT*14; OLCO13; KGB16] have been introduced to offset the global nature of the LBO by promoting sparsity and locality in the basis construction.In a related effort, Choukroun et al. [CSBK18] have proposed to modify the LBO through a potential function, thus defining a Hamiltonian operator, whose eigenfunctions have better localization properties.In [MRCB18], a similar approach was introduced to obtain a basis that is also orthogonal to a given set of functions.The "Coordinate Manifold Harmonics" used in [MMM*20], complement the LBO eigenfunctions with the coordinates of the 3D embedding, allowing to capture both extrinsic and intrinsic information.Finally, a rich family of diffusion and harmonic bases have been proposed in [Pat18], by exploiting the properties of the heat kernel.
While these basis constructions offer more flexibility and have been shown to improve the functional map pipeline in certain cases, e.g., [NVT*14; MMM*20], they nevertheless are typically still geared towards approximate isometries, and only enable approximate constraint satisfaction.In contrast our basis is geared towards landmark-preserving maps during functional map optimization, as well as during refinement.
Dirichlet-Steklov basis Finally, we note that Steklov eigenproblems have been considered within geometry processing [WBPS18] as tool for extrinsic shape analysis.This is achieved by considering the (two-dimensional) surface as the Steklov boundary of its (three-dimensional) interior.In contrast, we consider a fully intrinsic problem by using (one-dimensional) boundaries of small disks centered around the landmarks as the boundary of the remainder of the surface.

Method Overview
In this section, we provide a high-level overview of our approach.Our method takes as input a pair of shapes M, N represented as triangle meshes along with two sets of k landmark vertices We then aim to compute a highquality vertex-to-vertex correspondence ϕ : N → M that preserves the given landmarks exactly.I.e., ϕ(γ N i ) = γ M i for all i.Each of these basis sets is well-suited to describing smooth functions in the vicinity of its corresponding landmark circle.Intuitively these functions complement the Laplace-Beltrami eigenbasis, are harmonic on the interior of the shapes, and are zero at all but one disk boundary: u (i) j | Γ l = 0 for l = i and all j. 4. Compute an optimal functional map by minimizing an energy that promotes near-conformal maps, via an iterative refinement strategy.We split the functional map into k + 1 parts, and separately align the Laplacian eigenbasis and each set of Dirichlet-Steklov ones. 5. Convert the computed functional map to a vertex-to-vertex map between the shapes with the disks cut out.6. Reinsert the landmark correspondences to obtain a landmarkpreserving vertex-to-vertex map between the original meshes.
Our general strategy follows the standard functional map pipeline, especially in its recent variants based on iterative refinement [MRR*19; PRM*21; XLZ21], with several crucial changes.First, our main motivation for introducing disks to represent landmarks in Step 1 is to associate to each landmark a well-defined functional space.In this, we are inspired by techniques that represent landmarks or seed points on a surface via associated harmonic functions [ZRKS05; Pat18] on a mesh.Unlike such harmonic functions, however, our construction is fully justified in the smooth setting.This is because it is impossible to impose boundary conditions on isolated points on a smooth surface.Furthermore, as we submitted to COMPUTER GRAPHICS Forum (6/2022).Secondly, instead of computing a single functional map across Laplace-Beltrami eigenfunctions, we estimate a block-diagonal functional map that aligns each of the k +1 components of the functional space separately.This both improves efficiency and promotes desirable structural map properties.Indeed, we prove that this splitting must hold for conformal maps in the smooth setting, and we observe that it promotes preservation of landmark neighborhoods across a wide range of shape deformations in practice.
Finally, rather than focusing on near-isometries, we build a functional map energy that aims at computing near-conformal maps, and fully avoids the use of descriptor functions.Furthermore, we propose an efficient initialization and an iterative strategy for optimizing this energy, while promoting desirable map properties.This ensures high-quality correspondences even in challenging cases, in which existing functional maps-based methods tend to fail.
In the following sections, we discuss each step of this pipeline in detail.Throughout our discussion related to the basis construction and the structure of pointwise and functional maps, we focus mainly on derivations in the smooth setting, to highlight the theoretically justified nature of our approach.

Functional Basis
Central to our proposal is a careful choice of functional basis for a convenient functional space over the considered shapes.Our basis is built as the union of the solutions to the Dirichlet Laplacian and Dirichlet-Steklov eigenproblems, which we describe in Secs.4.1 and 4.2, respectively.In Sec.4.3 we define the functional space that we use in the rest of the proposal.The constructions described in these sections pertain to manifolds with boundaries and are not yet specialized to our shape matching method, which can be used both for shapes with and without boundaries.The specialization to our case is carried out in Sec.4.4.There, we describe how to create a landmark adapted functional basis by, roughly speaking, treating the landmarks as boundaries.All the constructions discussed in this section are carried out in the smooth setting.Their discretization is discussed in App.B.

Dirichlet Laplacian Eigenproblem
Let M be a smooth, connected, oriented compact Riemannian manifold with metric g and a boundary ∂M.The Dirichlet Lapla- where ∆ denotes the non-negative Laplace-Beltrami operator.Despite the vanishing boundary condition, it can be shown (see [Cha84]) that the eigenfunctions {ψ i } ∞ i=1 can be chosen to form an orthonormal basis for L 2 (M) (square integrable functions of M).Moreover, the eigenvalues and eigenfunctions can be ordered such that 0 < λ 1 ≤ λ 2 ≤ ... → ∞.
In Fig. 3, we illustrate the first few Dirichlet Laplacian eigenfunctions for an annulus in the plane with external radius 1 and internal radius 1/2.We will return to this example of the annulus in our discussion to compare the properties of different bases that we consider.We very briefly discuss the discretization of the Dirichlet Laplacian problem on triangle meshes in App.B.

Dirichlet-Steklov Eigenproblem
Let M be a smooth, connected, oriented compact Riemannian manifold with metric g and a Lipschitz continuous boundary ∂M.Suppose that, up to sets of measure 0, ∂M consists of two disjoint nonempty open sets, denoted D and S. The (mixed) Dirichlet-Steklov eigenproblem is posed as follows: where ∂n denotes the exterior normal derivative.The second and third line of the above are respectively known as the Dirichlet and Steklov boundary conditions, explaining the name Dirichlet-Steklov.
As hinted at in Sec. 3 and explained in detail in Sec.4.4, in our approach, S will be the boundary corresponding to a given landmark, while D will be the union of all other landmark boundaries.
The general theory of the Dirichlet-Steklov and many other similar problems can be found in [Neč12].For a gentle introduction to Steklov eigenproblems, see [Lab17] (in French).
We emphasize that, in contrast to previous uses of the Steklov eigenproblem in [WBPS18], we consider a purely intrinsic problem on the shape surface.I.e., as described in detail in Sec.4.4 our boundaries are one-dimensional, being the boundaries of disks centered at the landmarks.
As it is written above, the Dirichlet-Steklov problem may seem a bit mysterious.However, it becomes much more approachable when written in weak form: In this form, it can be compared to the standard Laplacian eigenproblem: Intuitively, and as we demonstrate in practice, the Dirichlet-Steklov eigenfunctions "focus" on the boundary S and provide detailed information in the vicinity of this boundary.As discussed below, in our method we establish one set of Dirichlet-Steklov eigenfunctions for each landmark, and align those functional spaces across the pair of shapes.
A derivation of the weak form of the Dirichlet-Steklov problem is provided in App. A. The discretization of the resulting problem on triangle meshes is discussed in App.B.
We illustrate some Dirichlet-Steklov eigenfunctions for the annulus in Fig. 4. Notice that the eigenfunctions concentrate on the boundary on which the Steklov boundary condition is imposed.

Functional Space W (M)
In this section, we specify the functional space used for the remainder of our proposal.Recall that our goal is to obtain a variant of the functional maps method suitable for non-isometric shape matching.We propose to search for near-conformal maps.
We thus need to translate the search for near-conformality to the functional maps setting.We do so by building upon the foundations laid in the context of conformal shape differences [ROA*13; CSB*17].Given a pair of surfaces M and N , in [ROA*13], it is observed that to study the deviation from conformality of a map ϕ : N → M, it is useful to consider its pullback F MN as a map between spaces of functions equipped with the Dirichlet form: The Dirichlet form becomes an inner product on the space of smooth functions modulo constant functions.A Hilbert space is then obtained by taking the completion of the space in the induced topology.We denote the space thus obtained by W (M). We remark that this space is different from the standard L 2 space of square integrable functions, due to the additional smoothness conditions.Below we describe both the properties and the utility of this space in the context of our landmark-based shape matching approach.

Landmark-Adapted Basis for W (M)
As highlighted above, a key aspect of our approach is the construction of a novel functional basis that is adapted to the landmarks.
Our main idea is to treat the landmarks as boundaries at which the functional bases satisfy certain boundary conditions.For this, we first slightly modify the shapes under study.Indeed, while advocated in several prior works [ZRKS05;Pat18] in geometry processing, it is not strictly speaking possible to impose boundary conditions at isolated points in the smooth setting.
Let M and N be compact, connected, oriented Riemannian surfaces.For simplicity, we also temporarily assume them to be without boundary.This last assumption is removed later.Let {γ M i } k i=1 ⊂ M and {γ N i } k i=1 ⊂ N be k landmarks in one-to-one correspondence.That is, we will be looking for maps ϕ : N → M such that ϕ(γ N i ) = γ M i for all i = 1...k.Such ϕ are said to be landmark preserving.The functional map representation of ϕ, that is its pullback on functions, will be denoted F MN , as before.
Our first step is to convert the landmarks into proper boundaries.We do so by removing small disks centered at the landmarks and treat the boundaries of these disks as boundaries of the shapes.We make sure that none of the disks intersect.Thus, we end up with a new shape that has k boundary components, one for each landmark.We denote the boundary corresponding to the landmark γ M i by Γ M i .By abuse of notation, we denote the shapes thus modified by M and N , same as their original versions.On triangle meshes, we create boundaries that are fully contained in a one-ring neighborhood of each landmark.This operation is described in detail in App. C. We now use the newly created boundaries to split W (M) into convenient subspaces.These subspaces will be composed of functions satisfying carefully chosen eigenvalue problems and boundary conditions.
We begin by considering the span of Laplace-Beltrami eigenfunctions satisfying Dirichlet boundary conditions on the {Γ i } i : Recall that the eigenfunctions {ψ i } ∞ i=1 form a orthonormal basis for L 2 (M).They remain mutually orthogonal in W (M), but interestingly fail to form a full basis for that space.This counterintuitive behavior is due to the change of topology from L 2 (M) to W (M) and the infinite dimensionality of the functional spaces under consideration.
Let the W (M) closure of the subspace spanned by the {ψ i } ∞ i=1 be denoted by G(M).
Naturally, our next step is to find functions that span the remainder of W (M).This is where the Dirichlet-Steklov eigenfunctions of Sec.4.2 come in.We pose k Dirichlet-Steklov problems, with the j th problem being: This results in k Dirichlet-Steklov eigenbases and spectra denoted {u , respectively.Recall that the {u form an orthonormal basis for L 2 (Γ j ).These functions remain mutually orthogonal in W (M).This follows directly from the weak form of the Dirichlet-Steklov problem (Eq.( 3)).We denote the W (M) closed span of the {u Our key result is that, once put together, the Dirichlet Laplacian eigenfunctions and the k sets of Dirichlet-Steklov eigenfunctions span all of W (M). Lemma 1.The function space W (M) admits the following decomposition: where ⊕ denotes direct sums and ⦹ denotes orthogonal direct sums, and the overline denotes the W (M) closure of the spanned functional space.
Intuitively, the above lemma says that W (M) can be split into a non-harmonic part and k harmonic landmark-associated subspaces, with each landmark getting its own subspace of harmonic functions that are nonvanishing on it.In practice, we always W (M)normalize all of the considered eigenfunctions by dividing each function by its W norm.In all of the following, we use W (M)normalized bases.
The resultant basis is thus normalized.However, it is not quite W (M)-orthogonal, as suggested by the notation used Lemma 1. Specifically, the problem lies in the mutual non-orthogonality of the subspaces H j (M).This is discussed in App.D.
In principle, the energy that we are to minimize (see Sec. 5 below) can be expressed in any basis, even if it is not orthogonal.For our purposes, however, the non-orthogonality of our basis poses a few challenges, which will be detailed later.Fortunately, in practice, our basis can be accurately approximated as orthonormal.A typical matrix of W (M) inner products is shown in Fig. 5 (see App. H for an extended evaluation of this approximation).In Fig. 6 we evaluate the orthonormality in the case of the 2D annulus and observe that it becomes more and more valid as the radius of the inner disk becomes smaller.We will call attention to this approximation when we use it in the implementation of our proposal.
Before proceeding further, we note briefly that on shapes with pre-existing boundaries we impose Neumann boundary conditions (vanishing normal derivatives).The above discussion remains unchanged.Note that imposing Neumann boundary conditions requires no special effort in the discrete setting.We now illustrate our functional basis using the landmark circles and Neumann boundary conditions on both the inner and outer boundary of the annulus in Fig. 7. Notice that as their eigenvalue increases, the Dirichlet-Steklov eigenfunctions rapidly concentrate on the landmark circles.In fact, the eigenfunctions of the closely related Steklov eigenproblem (i.e.without the Dirichlet boundary) are known to decay exponentially with distance from the Steklov boundary, the rate of decay being proportional to the corresponding eigenvalue [PST19].In contrast, the Dirichlet-Laplacian eigenfunctions remain evenly spread in the bulk of the shape.Thus, high eigenvalue Dirichlet-Steklov eigenfunctions are uninformative regarding the bulk of the manifold.Meanwhile, the high eigenvalue Dirichlet Laplacian eigenfunctions remain informative in the bulk even at high eigenvalues.
So far, we assumed that the considered shapes were connected.Our discussion remains unchanged on general shapes, as long as each connected component has at least two landmarks on it, as this is necessary to impose both boundary conditions of the Dirichlet-Steklov eigenproblem.If this is not satisfied for some connected component, at least some of the considered eigenproblems will have eigenfunctions that are piecewise constant per component and correspond to eigenvalue 0. These should not be included in a basis submitted to COMPUTER GRAPHICS Forum (6/2022).for W (M), as they have vanishing W −norm.We avoid this issue by rejecting eigenfunctions with eigenvalues below a certain small threshold.Note that for components with one landmark we only impose the Steklov condition on the corresponding circle, omitting the second line of Eq. (2).

Structure of the Functional Map
Recall that our ultimate goal is to compute a near-conformal diffeomorphism ϕ : N → M that preserves the landmarks.Recall also that we propose to use functional map methods to find it.In this section we translate the structural properties of ϕ into properties of its pullback F MN which helps us to restrict the space of admissible functional maps, which is crucial for our approach.
We begin on a technical note.Since we have replaced landmark points with landmark circles, the notion of landmark preservation has to be slightly adjusted.We no longer can claim something as simple as ϕ(γ N i ) = γ M i for all i, as the landmark points are no longer part of the considered shapes.Instead we impose that ϕ restricts to a diffeomorphism on corresponding landmark circles.That is, Now, suppose that ϕ is indeed a conformal map.Then, F MN satisfies the following lemma.Lemma 2 (Structure of F MN ).Let F MN : W (M) → W (N ) be the pullback of a conformal diffeomorphism that preserves the landmark circles in the sense described above.Then, F MN maps 1. G(M) to G(N ), 2. H j (M) to H j (N ) for all j.
Proof.See App.E.
The above lemma provides necessary, but not sufficient conditions for F MN to be the pullback of a diffeomorphism preserving the landmark circles.Nonetheless, we will use properties (1) and (2) of Lemma 2 to structure our search for F MN .
From now on, we only consider functional maps that satisfy statements (1) and (2) of Lemma 2. This can be seen as k + 1 separate maps, one for each landmark subspace H i and one for the orthogonal complement G, assembled into one block-diagonal functional map.Intuitively, this keeps the overall map tethered to the landmarks.
Landmark preservation At this point, it is worth explaining what we mean when we say that our method preserves the landmark correspondences in the discrete setting.Indeed, the challenge of landmark preservation is to not merely enforce the condition ϕ(γ i , but to also obtain a smooth (or at least continuous) map in the neighborhood of the landmarks (notice that we required ϕ to be a diffeomorphism when discussing the smooth setting).Our method achieves this by using a functional basis whose elements are well suited to describe smooth functions near the landmarks (recall the decay of the Dirichlet-Steklov eigenfunctions away from the Steklov boundary depicted in Fig. 7).By enforcing the functional map structure of Lemma 2 during the entire solution process, we promote vertex-to-vertex maps that smoothly map the neighborhoods of the landmarks of N to the corresponding neighborhoods on M, the smoothness of the map reflecting the smoothness of the functional basis.Furthermore, we reinsert the original pointwise landmarks at the end of the solution process to preserve the initial landmarks exactly.Recall that the landmark vertices are excluded from the meshes the moment the landmark circles are introduced.

Functional Map Energy
The previous section describes our landmark adapted basis construction, and the block-diagonal structure of landmark-preserving conformal maps when expressed in this basis.In this section we specify the optimization problem that we will solve in order to obtain landmark-preserving maps between triangle meshes.
Recall that we propose to look for conformal maps, which can be characterized in terms of the Dirichlet form (W (M) inner product).Theorem 3. Let ϕ : N → M be a diffeomorphism between oriented Riemannian surfaces with pullback F MN : In practice we do not expect to obtain an exact equality of the inner products as described in the previous theorem.Instead, we will search for ϕ and F MN by relaxing the above equality to a minimization problem.Let Φ M and Φ N denote reduced (finite dimensional) functional bases for W (M) and W (N ), respectively.These bases consist of the eigenfunctions of the Dirichlet Laplacian and Dirichlet-Steklov eigenproblems corresponding to small eigenvalues.The precise size of the bases is discussed in App.H.
From now on, we concentrate our attention on the discrete case.
Namely, M and N will now denote oriented manifold triangle meshes.Letting Φ M , Φ M W (M) be the matrix of all inner products of the normalized basis vectors of Φ M , we relax the equality of Theorem 3 to the minimization of the following energy term: We call this the conformal term of the energy.Here, as well as everywhere else in this text, • F denotes the Frobenius norm.
Having covered the conformality of the map, it remains to rephrase the restriction of F MN to pullbacks of landmarkpreserving diffeomorphisms.This assumption cannot be exactly imposed in the discrete case.Still, we would like F MN to exhibit the properties of such a map.In order to do so, we complete our energy by specifying two structural terms.Specifically, the first term promotes F MN being a proper functional map (i.e., the pullback of a vertex-to-vertex map), as recently defined in [RMWO21], and the second promotes the invertibility of F MN [ERGB16].
Let Π N M denote the vertex-to-vertex map from N to M expressed as a matrix (i.e. a binary matrix that contains exactly one 1 per row).Then, F MN should satisfy: where (Φ N ) + denotes the pseudoinverse of Φ N , or in other words, the W (N ) projection onto the reduced basis Φ N .As before, we relax the equality into an energy to be optimized: We call this the properness term of the energy.Notice that we have expressed the energy as a function of both F and Π.We do so as we will have to consider these two objects as independent variables when minimizing the energy.The exact meaning of this is discussed in Sec. 6.
In addition to F MN arising from a point-to-point map, we would also like for it to be invertible.For this, we consider two maps F MN : W (M) → W (N ) and F N M : W (N ) → W (M), the latter arising from a vertex-to-vertex map Π MN : M → N .Thus, in what follows, we will be simultaneously optimizing for maps going in both directions between the shapes.With I being the identity matrix, the invertibility condition is, of course: Once again, we convert the above into minimization form.The invertibility term corresponding to the first line above is The invertibility term E I,N M is defined analogously.
In sum, our search for the correspondence between M and N will involve the joint minimization of the energy and an analogously defined energy E N M .Here, a C , a P and a I are nonnegative tunable weights controlling the relative strength of the conformality, properness and invertibility terms, respectively.Different values of these parameters are explored in App.H.The above energy is conformally invariant in the following sense.Lemma 4 (Energy Invariance).The conformality, properness and invertibility terms of the energy (Eqs.(9), ( 11) and (13)), as well as the energy (their weighted sum, Eq. ( 14)) are invariant under (combinations of) the following transformations: 1. Conformal transformations of the meshes keeping the reduced bases fixed.2. Orthogonal transformations of the reduced bases.
Proof.Let the functional and vertex-to-vertex maps be fixed.Since conformal transformations leave the W inner product invariant, the energy terms are conformally invariant for a fixed choice of functional basis.Statement 1. is now proven.Statement 2. follows from the fact the Frobenius norm is invariant under orthogonal transformations.
Note that, in the lemma above, conformal transformations and changes of basis are treated as independent.In practice, they are not, as reduced bases are usually mesh-dependent.Thus, the change of basis induced by a conformal transformation may fail to be orthogonal and then Lemma 4 will not apply.As long as one works with reduced rather than full bases (i.e.spanning all functions of the mesh), the invariance of the energy under conformal transformations is therefore only approximately guaranteed.
The conformal invariance of the energy can be violated in another way.Suppose that the recipe for constructing the reduced bases produces non-orthogonal bases.Then, the change of basis induced by a conformal transformation may fail to be orthogonal even when full bases are used.
Of course, we raise the previous two issues precisely because our method uses non-orthogonal reduced bases.Thus, Lemma 4 does not offer a full guarantee of conformal invariance for our energy.Still, the bases that we use turn out to be very nearly orthogonal and thus the energy that we employ remains approximately conformally invariant.Obtaining a truly conformally invariant energy (at least up to basis truncation) is a subject for future work.

Solving the Problem
In this section, we propose an efficient approach for the optimization problem posed in Eq. ( 14).Our approach is inspired by a discrete optimization strategy, first suggested in [MRR*19] and recently extended to other general energies [RMWO21].The general idea is to recast the problem in a way that makes every iteration of the optimization into a nearest-neighbor search.The overall process submitted to COMPUTER GRAPHICS Forum (6/2022).
then consists of two qualitatively different parts.First, an initial guess of the correspondence is obtained.Then, the correspondence is refined via the iterative process mentioned above.These steps are explained in Secs.6.1 and 6.2, respectively.

Initial Correspondence
In this section we explain how we obtain an initial guess of the functional maps F MN and F N M .A common way to initialize functional maps with landmarks is via descriptor preservation [OBS*12; RPWO18].However, common descriptors such as HKS or WKS [SOG09; ASC11] strongly rely on the isometry assumption and moreover the initial functional map is not guaranteed to respect landmark correspondences exactly.To overcome this, we propose a simple and lightweight initialization scheme.
Recall that our approach upgrades landmark correspondence to landmark circle correspondence.Moreover, in the smooth setting, we require this correspondence to be a diffeomorphism.We now make the assumption that the correspondence between landmark circles can be seen as a rotation of one circle to match the other.Specifically, we label the vertices of each landmark circle in counter-clockwise order using the outward-facing normal orientation.We can then assign each vertex in a landmark circle coordinates in [0, 1).All that remains is to ensure that the origin of this coordinate system is placed consistently on both shapes.In other words, the matching of two corresponding landmark circles reduces to finding an appropriate shift of one of the parametrizations.
We propose to align the parametrizations of the boundary circles such that the landmark circles are consistently oriented relative to the other landmarks.In order to do so, we construct functions on the landmark circles that have maxima in directions roughly pointing towards other landmarks.Consider the following problem: This results in k harmonic functions, one for each landmark, where each function h i equals 1 on the boundary of landmark circle i, and zeros on the boundaries of other circles.As they stand, these functions are constant on each landmark circle.Their normal derivatives, however, are not.In essence, we use ∂n Γ i h j as an indication of the direction one should take from Γ i to reach Γ j .See Fig. 8 for an illustration.This is similar in spirit to the Geodesics in Heat construction [CWW13], where gradients of solutions to the heat transfer problem are used to construct approximate geodesics.

Denoting the landmark circle coordinates on Γ M
i by θ i , we select the optimal shift α i by solving: In this problem, we consider each landmark i and examine directions to all other landmarks (via normal derivatives).We then align the coordinates of landmarks on M and N so that these directions align in the best possible way.This problem can be solved simply by directly examining all possible shifts and taking the optimum.
In order to gain robustness to changes in triangulation, we first project the normal derivatives (as circular functions ∂n ) onto the reduced basis in order to remove spurious high frequency components.Recall that this makes sense as the Dirichlet-Steklov eigenfunctions belonging to landmark Γ i form a basis for L 2 (Γ i ).
Converting the optimal shifts α i into vertex-to-vertex correspondences on the landmark circles is a matter of a nearest-neighbor search between the circular coordinates of the vertices of Γ M i and the α i -shifted coordinates of the vertices of Γ N i .It remains to convert the resulting vertex-to-vertex map into a functional map.Once again, recall that our reduced basis contains an L 2 basis for each landmark circle.Thus, by using an expression of the form of Eq. (10) we can construct functional maps between H i (M) and H i (N ).Assembling the resulting maps into blockdiagonal matrices gives our initial guesses of F MN and F N M .
In App.G, we compare this approach to two alternative initialization strategies, and demonstrate its relative advantages.

Energy Minimization via Nearest Neighbor Search
The procedure described in Sec.6.1 provides a descriptor-free initial guess for the functional map.In this section we describe a refinement method that significantly improves the map.
Recall that we are looking for a vertex-to-vertex map by minimizing an energy that depends on both the point-to-point and the associated functional map (pullback).In [MRR*19] it is observed that a particular case of such problems can be efficiently solved by considering the two maps as being independent variables.This observation was recently extended to a wide range of energies in [RMWO21].Following this line of work we will move all of the difficult optimization on the side of the vertex-to-vertex map and use Eq.(10) to restore the relationship between the maps.
Our main tool is the following result, standard in functional maps literature [EB17; RMWO21], which allows one to reduce optimization problems of a certain form to nearest neighbor searches.Lemma 5. Let A be a symmetric positive-definite matrix inducing the matrix norm M 2 A = Tr(M T AM).Let Φ be a reduced basis orthogonal with respect to A, that is Φ T AΦ = I.Then, given n pairs of matrices X i and Y i , the following two expressions are equal: Moreover, if A is diagonal, minimizing the above expressions over matrices Π that reflect point-to-point maps (i.e., binary matrices that contain exactly one 1 per row) is equivalent to This problem can be solved via nearest-neighbor search between the rows of the concatenated matrices [X 1 ... Xn] and [ΦY 1 ... ΦYn].
Proof.See [EB17] for a proof of a special case and [RMWO21] (Lemma 4.1) for the general statement.
We first convert the energy of Eq. ( 14) into the form used in the above lemma.For brevity's sake, we will only develop the expression for F MN and Π N M .The expression for the pair F N M and Π MN is analogous.As mentioned in Sec.4.4, we approximate the functional bases Φ M and Φ N to be orthonormal with respect to the Dirichlet form.Then, the energy minimized by the desired F MN and Π N M becomes: Here, we used the approximation of basis orthonormality in two ways.First, we used it to evaluate the inner products in the conformality term (first line of the above equation).Second, we used it to express (Φ N ) + = (Φ N ) T W N , where W N is the so-called cotangent Laplacian on N , which also corresponds to the piecewise linear finite element discretization of the Dirichlet form.We are still a few modifications away from being able to apply Lemma 5 to this problem.
We obtain the desired form for the expression by replacing certain instances of F MN with its expression in terms of the vertexto-vertex map Π N M : (Φ N ) T W N Π N M Φ M .By using the fact that I − F T F 2 F = FF T − I 2 F and making this replacement, we obtain: Now, all of the terms of the above are of the form Φ T AΠX i −Y i 2 F , with W N playing the role of the matrix A. Our energy is thus of the form of line (1) of Lemma 5, up to three terms of the form (I − ΦΦ T A)ΠX i 2 .Notice that (I − ΦΦ T A) is the orthogonal projection onto the orthogonal complement of the reduced (approximately) orthonormal basis Φ.Thus, this term can be seen as a regularizer penalizing solutions lying outside of the considered reduced basis.Indeed, this is how this term is was originally introduced in [EB17].Consequently, by implicitly introducing the appropriate regularizers we can use the first part of Lemma 5 to obtain the following expression for the energy: At this point we are forced to make an approximation.Namely, we assume that the second part of the lemma applies, which would normally require W N to be diagonal.In other words, we convert the problem into a nearest neighbor search without having the guarantee of the equivalence of solutions.Despite this approximation, we have observed that the resulting approach works remarkably well in practice.
This finally brings us to the procedure that we use to minimize the energy.As mentioned above, we will consider the functional and vertex-to-vertex maps as independent variables.Thus, given functional maps F MN and F N M , the point-to-point map Π N M can be found by solving the nearest-neighbor search problem: Here NNS(A, B) denotes a set of nearest neighbor problems: for each row of B among the rows of A. The vertex-to-vertex map Π MN can be obtained analogously.In sum, minimizing the energy with respect to the vertex-to-vertex maps is also a recipe for converting functional maps into vertex-to-vertex maps, while taking into account the original functional map energy.
We are now ready to formulate the optimization algorithm.Following [MRR*19], the overall procedure is based on an iterative spectral upsampling of the functional map.Specifically, we iteratively convert the functional map into a vertex-to-vertex map while increasing the size of the reduced basis.As explained earlier (see Fig. 7), the Dirichlet-Steklov functions are concentrated near the landmark circles.Thus, increasing their number does not provide much additional information about the map in the bulk of the shapes.Thus, we only increase the number of Dirichlet Laplacian eigenfunctions.
Beginning from the initial functional maps F MN and F N M obtained in Sec.6.1, we proceed as follows.
1. Convert F MN and F N M into Π N M and Π MN via Eq.(21).2. Increase the reduced bases Φ M and Φ N by including kstep additional Dirichlet Laplacian eigenfunctions.3. Update the functional maps to the new basis size via F MN = (Φ N ) + Π N M Φ M and F N M = (Φ M ) + Π MN Φ N .4. Iterate steps (1) to (3) until the desired basis size is reached.5. Repeat step (1) using only the original non-landmark vertices.This produces a vertex-to-vertex map between the original meshes, landmarks excluded.6. Insert the landmark correspondence into the vertex-to-vertex map.
A Fast Approximation.We conclude this section by proposing an submitted to COMPUTER GRAPHICS Forum (6/2022).
acceleration strategy to perform the nearest-neighbor search.The method proposed here is unprincipled, but is validated by both the overall quality of our results and explicit tests found in App.G.The method proposed below is the only one used in the main text of this paper.
In the language of Lemma 5, we propose to replace the nearest neighbor search between the concatenated matrices [X 1 ... Xn] and [ΦY 1 ... ΦYn] by a nearest neighbor search between the summed matrices X 1 + ... + Xn and ΦY 1 + ... + ΦYn.This corresponds to solving the following problem: This reformulation helps to decrease the dimensionality of the nearest neighbor searches.Essentially, we assume that the different energy terms will not cancel each other.The payoff for this approximation is that the matrices involved in the nearest-neighbor search become n times smaller.In our case, there are n = 3 energy terms.
The experiments in App.G show that this reduction in matrix size results in a slightly more than threefold speed-up.

Evaluation
We evaluate our method † on standard shape matching datasets, which we describe in Sec.7.1.We first analyze the parameters involved in our computations (Sec.7.2).Second, we conduct an in-depth evaluation to compare our method to state-of-the-art approaches on shape matching benchmarks (Sec.7.3).
For our quantitative evaluation in Fig. 12 (right), Fig. 14, Fig. 16 and Fig. 18, we follow the commonly-used protocol, introduced in [KLF11] by plotting the percentage of correspondences below a certain geodesic distance threshold from the ground truth.

Datasets
We perform all our experiments on the following datasets.
FAUST [BRLB14].This dataset contains models of ten different humans in ten poses each.Despite the variability in the body types of said humans, this dataset is typically considered as nearisometric.We remesh the shapes of the dataset to shapes with approximately 5K vertices and use 300 shape pairs following the procedure of the authors of [RPWO18].Note that the shapes in question are remeshed independently and do not share the same connectivity.
TOSCA [BBK08].This dataset consists of meshes of humans and animals.Following [RPWO18], we split this dataset into isometric and non-isometric shape pairs.We call the resulting datasets TOSCA isometric (284 shape pairs) and TOSCA non-isometric (95 shape pairs) respectively.The shapes of these datasets are remeshed independently to count around 5K vertices per shape.Once again, the remeshed shapes have distinct connectivity.The mean geodesic error is averaged on the 95 shape pairs of the TOSCA non-isometric dataset (remeshed to 5K vertices).Notice how stable our method remains, even for extreme values of r f .SHREC'19 [MMR*19b].This challenging dataset is composed of human shapes with high variability in pose, vertex count (ranging from 5K to 200K vertices) and topology (some shapes are watertight manifold meshes whereas other have holes and other surface noise sources).
FAUST "Wild" [SACO20].This dataset is a variant of FAUST in which challenging differences in connectivity are introduced via remeshing.We use the following types of remeshing of the dataset: a uniform isotropic remeshing (iso), a remeshing where randomly sampled regions are refined (dense), and the remeshing proposed in [GH97] (qes).Finally, we consider correspondences across the 20 template models of the dataset instead of solely considering the initial template shape as the source shape.

SHREC'20 [DLR*20]
. This dataset proposes a collection of 14 animal shapes with a set of landmarks determined by experts.The animal pairs contain parts in correspondence with highly nonisometric deformations.We only consider the correspondences between full shapes for our experiments (test sets 1 to 4).

Parameter Study
We present here the main results concerning the parameters of our method.Other minor experiments on this topic are presented in App.H (influence of the weights in the energy, qualitative illustration of the impact of landmark placement, near-orthogonality assessment for our basis and study of the effect of basis size).

Radius r f
The construction of the landmark boundaries Γ i explained in App.C relies on the user-defined scalar parameter r f ∈ (0, 1).In Fig. 9, we study the influence of r f on the geodesic matching error averaged on the TOSCA non-isometric dataset, with 7 landmark correspondences at their standard locations (see App. I).It demonstrates empirically that this parameter has no significant impact on the matching performance.We therefore set r f = 0.5 in all our other experiments.

Landmark placement
In order to study the influence of landmark placement on our method, we conduct the following experiment on 10 shapes of the TOSCA Isometric dataset (cat category).We consider an increasing number of landmark correspondences, ranging from 3 to 100, submitted to COMPUTER GRAPHICS Forum (6/2022).placed according to four standard surface sampling strategies: (i) random, (ii) euclidean farthest point (iii) geodesic distance farthest point (iv) Poisson disk (as implemented in [Jac*18]).The outcome of these experiments is illustrated in Fig. 10.The farthest point sampling strategies result in the fastest decrease of the error, Poisson disk is slightly slower and random placement is predictably the slowest.This indicates that our method performs best when the extremities of the shapes are prioritized for landmark placement.The landmark placement used in the benchmarks of Sec.7.3 makes use of this observation (see App.I for details).
To complement the above experiment, we show the variance of our method when initializing two sampling strategies with 3 different seeds in Fig. 11 on the full TOSCA non-isomtric dataset.

Remeshing invariance
In order to show that our method remains applicable on shapes with different triangulations, we remesh independently the target pair of each FAUST shape pair and compute the mean geodesic error in Fig. 12 (left).We additionally experiment with the FAUST "Wild" dataset created in [SACO20] to assess invariance to the remeshing proposed by the authors.Fig. 12 (right) and Tab. 1 present the output of this experiment.We observe marginal difference when considering the various remeshing approaches tested, which highlights the insensitivity of the proposed approach to the shape connectivity.Fig. 13 illustrates qualitatively the median transfer obtained on this dataset.Figure 13: Qualitative illustration of the median map quality obtained with our method on three types of remeshing in the FAUST "Wild" dataset (see Sec. 7.1).Despite the great disparity of the underlying meshes, our method provides smooth transfers.

Benchmarks
In this section, we describe the competing state-of-the-art methods that we employ (Sec.7.3.1)and present our main results for shape matching (Sec.7.3.2).

Setup
We compare our method against three competitors that leverage landmark information to compute correspondences between shapes.The detailed setup for each method, including the landmark placement is provided in App.I.The competing methods are:   Weighted Averages (WA) [PBDS13] also defines a parameterization of the input surfaces that preserves landmarks exactly: each point at the surface is expressed as a weighted average of its distance to a set of landmarks.
Functional Maps With ZoomOut Refinement (FMap ZO) [MRR*19] computes correspondences between shapes by leveraging a functional basis defined on the source and target shapes.While the method does not allow to retrieve exact correspondence between user-specified landmarks, it constitutes the current state-of-the-art method for isometric shape matching.

Results
In this section, we present our main results on shape matching.
Isometric shape matching.The evaluation on FAUST and TOSCA Isometric are illustrated in Fig. 14, with averaged errors and runtimes displayed in Tab. 2. On the FAUST data set, our approach remains competitive with a mean geodesic error of 1.40 × 10 −2 and a mean computation time of 8.83 s.On the TOSCA isometric data set, we obtain a slightly better average geodesic error score than competitors.Qualitatively, our method produces smooth texture transfers on both data sets, as highlighted in Fig. 15.
Non-isometric shape matching.We run an evaluation of our method on the TOSCA non-isometric and the SHREC'20 datasets (Fig. 16).The mean error values and timings are showed in Tab. 3. In this challenging setup, our method has the best results in terms of mean geodesic error, while being the second best in terms of computation time.Fig. 17 presents a qualitative evaluation using a texture transfer on a pair of shapes for each data set.
SHREC'19 benchmark.The quantitative evaluation is reported in Fig. 18, with the associated averaged geodesic errors on the right of the figure.Our method obtains the best mean geodesic error score for this difficult benchmark.In addition, a qualitative evaluation via texture transfer is depicted in Fig. 19   Table 3: Quantitative evaluation results on the TOSCA nonisometric (n-i.) and the SHREC'20 lores (without partial shapes) data sets.The average geodesic error (Av.Geo.Err.) and average execution time (Av.Time) on both data sets are displayed for competing approaches and our method.) is displayed for our method and competing approaches.
performance on this dataset is indicative of its stability and applicability across diverse changes in shape topology, such as the introduction of small holes.This is a general feature of the functional maps methods, which our approach inherits.

Source HyperOrb FMap ZO Ours
Figure 19: Qualitative evaluation of our method and competing approaches on a shape pair from the SHREC'19 data set, selected such that the geodesic error of our method is median over the dataset.The best and worst cases are illustrated in App.G.

Conclusion, Limitations and Outlook
We have proposed an efficient functional maps-based shape matching approach that promotes conformal maps and exactly preserves landmark correspondences.This was achieved via the introduction of a novel functional basis and an energy promoting bijective conformal maps.The efficiency of our solution comes from an adaptation of the ZoomOut procedure [MRR*19; RMWO21] using our energy and novel basis.The resulting method exhibits state-ofthe-art performance on non-isometric benchmark datasets and near state of the art performance on isometric ones.
Recall, however, that our usage of the ZoomOut procedure was not fully principled.Indeed, we needed to make some approximations in order to use Lemma 5, which converts certain optimization problems into nearest-neighbor searches.The quality of our results indicates that our approximations were justified, suggesting that Lemma 5 could likely be rigorously extended to suit our needs.In fact, extending Lemma 5 would be of general interest to the functional maps community, as it would enable the efficient minimization of various other energies.
The construction of our landmark-adapted basis required us to upgrade the landmarks to proper boundaries.We did so by cutting submitted to COMPUTER GRAPHICS Forum (6/2022).out small disks centered at the landmarks, resulting in the introduction of landmark circles.The landmark circles offer an intriguing possibility that we have not explored here.Namely, one could augment landmark correspondence to include a user-specified matching of the landmark circles.This could allow for greater semantic or artistic control of the resulting map.Our initialization procedure of Sec.6.1 can be seen as an automated implementation of a similar idea.Furthermore, since our present work has demonstrated the fruitfulness of landmark-adapted bases, it is natural to ask whether better performance can be achieved by improving upon basis construction.In particular, we have noted that the Dirichlet-Steklov eigenfunctions have their amplitude intensely concentrated near the landmark circle equipped with the Steklov boundary condition (see Fig. 7).It seems likely that an analogous basis with less concentrated functions could be better suited to describe the behavior of the functional map near the landmarks.Notice that this dovetails with the idea of user-specified landmark circle correspondence, as the user-provided information would have impact further away from the landmarks.

Acknowledgements
The authors would like to acknowledge the anonymous reviewers for their helpful feedback and suggestions.Parts of this work were supported by the ERC Starting Grants No. 758800 (EXPRO-TEA), the ANR AI Chair AIGRETTE, the Swiss National Science Foundation (SNSF) under project number 188577 and the Association Nationale de la Recherche et de la Technologie (ANRT) via the Convention industrielle de formation par la recherche (CIFRE) grant No. 2019/0433.Finally, the authors wish to thank Prof. Alexandre Girouard for his assistance in navigating the literature on the Steklov eigenproblem, Jing Ren and Simone Melzi for providing tools for shape analysis via functional maps and the TOSCA/FAUST/SHREC'19 datasets, Nicholas Sharp for releasing the "FAUST Wild" dataset and Patrick Schmidt for computing a baseline to his work [SCBK20].
A Word of Warning As a final note on the discretization of the considered eigenproblems, we would like to warn the reader of a small issue one may encounter when numerically solving them.Recall that we want the Dirichlet-Steklov eigenfunctions to be normalized with respect to the boundary mass matrix S M .Solvers for generalized eigenvalue problems, such as Matlab's eigs routine, which we use in our implementation, will typically do so automatically.However, according to our observations, sometimes this automated process will not happen.This seems to be related to the fact that S M is a positive semi-definite matrix rather than a positive definite one.Thus, one needs to explicitly normalize the solutions with respect to S M .In fact, we suggest explicitly normalizing even the Laplacian eigenfunctions, despite the fact that there the mass matrix A M is positive definite on (good quality) triangle meshes.Indeed, A M can fail to be positive-definite on pathological inputs.Consider for instance an otherwise good mesh with an isolated vertex belonging to no triangle.Functions vanishing everywhere except on said vertex have norm 0 with respect to M , despite being nonzero.

Appendix C: Boundary Circles on Triangle Meshes
In Sec.4.4, small disks centered at the landmarks are removed in order to create new boundaries for the shapes under study.Here, we describe in detail how this is achieved on triangle meshes.Crucially, we do not want to unduly disturb the geometry of the shapes.In order to achieve this we construct the new boundaries entirely within the triangles adjacent to the landmarks.
Let's say that we are constructing the boundary circle for the landmark γ i .We begin by selecting the radius r i of the disk to be removed.This is done by finding the length s i of the shortest edge connected to γ i .The minimum is taken over both shapes, which are scaled to be of identical surface area and thus of comparable size.Then, we set r i = r f • s i , where r f ∈ (0, 1) is a user-set parameter.The (surprisingly low) impact of this parameter is studied in Sec.7.2.1.
We are now ready to construct the boundary Γ i .This process is best understood by looking at its illustration in Fig. 20.First, we split each triangle adjacent to the landmark into ns wedges of equal angle, which introduces ns − 1 new vertices at the opposite edge of the original triangle, as well as edges connecting them to the landmark.Then, we introduce ns + 1 new vertices situated on the new edges at a distance r i away from the landmark γ i .We then connect these vertices in a way that creates an approximation of a sector of a disk of radius r i .Doing so produces ns quadrilaterals in the part of the original triangle far from the landmark.We split those quadrilaterals into triangles along their diagonals.This concludes the refinement of the triangles adjacent to the landmark.It remains to refine the triangles adjacent to them across the edges opposite to the landmark.There, the common edges between the triangles contains ns − 1 new vertices.On each triangle, we connect these new vertices to the original vertex not on the common edge.This concludes the refinement process.Note that all of the new triangles are contained within the original ones.An example of a mesh

Input
Step 1 Output with landmark circles constructed in this manner is shown in Fig. 5.
The construction of the boundaries associated to different landmarks is done sequentially over the landmarks.This requires some additional care if the landmarks are placed too close to each other.Indeed, during the construction of Γ i , new faces are created in what was originally the 2−ring neighborhood of the landmark γ i .Thus, if a different landmark γ j is closer than 4 rings away from γ i , there will be overlap between the newly created mesh faces.The resulting mesh will then be dependent upon the order in which the boundary circles Γ i and Γ j are created.In the present paper, we avoid this issue by disallowing such landmark placement.If such landmark placement becomes necessary in a given application, we suggest locally refining the mesh via, say, √ 3−subdivision [Kob00] such that the landmarks are no longer closer than 4 triangle rings from one another.We do not pursue this here.
Appendix D: Proof of Lemma 1 and Discussion on its Meaning Lemma 1.The function space W (M) admits the following decomposition: where ⊕ denotes direct sums and ⦹ denotes orthogonal direct sums.
Proof.Recall that, by construction, W (M) is the completion of smooth functions modulo constants with respect to the Dirichlet form.Thus, we begin our analysis on smooth functions.
Let u be smooth and W (M)-orthogonal to all of the Dirichlet-Laplacian eigenfunctions {ψ i } ∞ i=1 .Then, by Stokes' theorem, submitted to COMPUTER GRAPHICS Forum (6/2022).
See [Jos08] for the relevant theory.
In the discrete setting, we use the same method as in [ESB19] to compute the Dirichlet energy.Namely, the expression becomes where E(M) denotes the edges of the mesh M, w M uv denotes the cotangent weight of the edge (u, v) and D 2 N (•, •) is the matrix of square geodesic distances on N .

Analysis of Alternative Initialization Methods
The iterative optimization procedure detailed in Sec.6.2 requires as an input an initial guess of the functional map.In Sec.6.1, we thus introduce an initialization procedure for this initial guess based on the landmark correspondence and the normal derivatives of certain landmark-dependent harmonic functions.In this section we compare this approach to two alternatives.
For the purposes of this discussion, the approach of Sec.6.1 shall be referred to as the "normal derivatives" method.The two alternatives described below will be termed "trivial" and the "conformal energy", for reasons that should soon become apparent.
The landmark circles can be seen as lists of vertices ordered counter-clockwise as seen from outside the shape.The choice of the first element of this list carries no particular meaning and is left to the whims of the indexing of the faces of the mesh.Thus, the first elements of two corresponding boundary circles need not match.The "trivial" approach consists in assuming that the first elements of the boundary circles do indeed match.This correspondence is then proportionally extended to the rest of the landmark circle.
The "conformal energy" approach stems from the observation that mapping the landmark circles Γ N i → Γ M i induces a restricted functional map H i (M) → H i (N ).The conformal term of the energy (Eq.( 9)) can be easily evaluated on these subspaces.The "conformal energy" approach consists in choosing the shifts {α i } k i=1 (see Sec. 6.1) such that they minimize the conformal energy of the resulting H i (M) → H i (N ) map.Fig. 21 (left) depicts the performance of the three initializations in terms of geodesic error on the SHREC'20 dataset (lores), using 7 landmarks.Tab. 4 provides quantitative evaluations for the same experiment in terms of averaged geodesic error and Dirichlet energy.The "normal derivatives" approach slightly outperforms the other two on all metrics, which is why it is the one used in the main text.

Comparison of the "Principled" and "Fast" Energy Optimization
At the end of Sec.6.2, we introduced an unprincipled way to accelerate the nearest neighbor search used in the solution of our problem.In this section, we quantitatively compare this "fast" method to the "principled" one on the SHREC'20 data set (partial shapes excluded).The output of this evaluation is displayed in Fig. 21 (right) and Tab. 5.While very similar in terms of matching performance, the "fast" method is more than three times faster to compute.We therefore employ it instead of the "principled" approach.Note that the more than threefold speedup is consistent with the fact that the matrices used in the "fast" method are three times smaller.

Complementary benchmark on SHREC'20 lores
As a complement to our main evaluation on SHREC'20 lores, we conducted an evaluation using only 8 pairs from the initial benchmark to compare against the method proposed in [SCBK20] (In-terSurf).InterSurf, WA, HyperOrb FMap ZO and our approach obtain a geodesic error (scaled by a factor ×100) of respectively 11.9, 5.41, 5.99, 8.69 and 5.2.The restricted number of shapes on which we evaluate is due to the fact that InterSurf does not handle shapes with complex topologies well.In particular, the method assumes that the meshes are watertight and share the same genus, in strong contrast to our approach that does not make such assumptions.However, we note that this method was not primarly designed for shape matching.

Additional Qualitative Evaluations
We provide additional qualitative evaluations on isometric and nonisometric shape pairs in order to show best-and worst-case shape matching scenarios for our method.
For isometric shapes, the best pairs are depicted in Fig. 22 and the worst pairs in Fig 23 .For non-isometric shapes, the best pairs are illustrated in Fig. 24 and the worst pairs in Fig. 25. Figure 22: Qualitative evaluation of our method and competitors on isometric shapes from the FAUST dataset (top row) and the TOSCA isometric dataset (bottom row).The shape pair is selected such that the geodesic error of our method is the best over the dataset.

Study of the Weights in the Energy
We define three weights to compute a point-to-point map between two shapes based on the energy (Eq.( 14)): the conformal, the properness and the invertibility weights, denoted respectively a C , a P and a I .Since we normalize the weights, their absolute value is unimportant.
To study how their relative value influences the quality of the output map we conduct a dedicated experiment on the SHREC'20 ergy) when choosing the weight configuration.Roughly speaking, the invertibility and properness terms promote accuracy, while the conformality term promotes smoothness.
Since this trade-off is application-dependent, we leave the finetuning of the energy weights to the end-user and set all weights to 1 in the remaining of our experiments as it provides a satisfactory balance in practice.

Landmark Sampling Qualitative Illustration
We visualize qualitatively the interest of introducing more landmark correspondences in Fig. 29.In this visualisation, since "Hy-perOrb" does not support less than 5 landmark correspondences, no map for 3 and 4 landmark correspondences can be computed for this method.
Note how the regions around the mouth and the eyes are accurately mapped with our approach compared to the two other approaches.

Basis Near-Orthogonality
For each shape M of the SHREC'19 data set [MMR*19a], we compute the matrix with entries m i, j = Φ M i , Φ M j W (M) , where Φ M i designates the i-th basis vector.We use 7 landmarks, 10 Dirichlet-Steklov eigenfunctions, leading to a Dirichlet-Steklov block of size 70 × 70, and 120 Dirichlet Laplacian eigenfunctions.Since we are only interested in the computation of the basis itself in this setup, the landmarks were placed at random locations to maximize the diversity of situations encountered.The average of all matrices is displayed in Fig. 30.Note the clear diagonal behavior, that is in agreement with our observations on a simple sphere shape (Fig. 5).
Figure 29: Qualitative comparison of our method to competitors when increasing the number of landmarks on the same shape pair as for our teaser (Fig. 1).The ground truth landmark locations are denoted by green dots.In the case of FMapZO (no exact landmark preservation), the blue dots indicate the location of the mapped landmarks.We highlight that this computation also sheds light on the robustness of our basis computation to complex triangulation and partiality setups.

Number of Basis Functions
To select the number of basis functions for G(M) and each H j (M) (see Sec. 4.4), we study their respective size N LB and N DS separately, as illustrated in Fig. 31.
Increasing the size of G(M) slightly increases the matching performance up to N LB = 120.In contrast, varying N DS above 10 decreases the quality of the maps.Hence, we fix the following basis sizes throughout the rest of the article: N LB = 120 and N DS = 10.The last landmark is only used on the TOSCA and SHREC'20 data sets.Notice that our landmark placement is reminiscent of farthest point sampling.The landmark placement is common to all considered methods.The other parameters depend on the method used.

Method Configuration
Hyperbolic Orbifold Tutte Embeddings (hyperOrb) and Weighted Averages (WA).These methods do not require any additional parameters.
Functional Maps With ZoomOut Refinement (FMap ZO).A 20 × 20 functional map is computed for each source-target pair in setup 1 and 2, following the setup of [MRR*19].In particular, we use wave kernel signature and wave kernel map functions as descriptors.The descriptor functions are computed at the same ground truth landmark positions used for the other methods.At each landmark location, 12 wave kernel map functions are computed using a basis of 120 LB-eigenfunctions.
The energy employed to compute the functional map leverages the descriptor preservation, descriptor commutativity and LBcommutativity terms.Contrary to [MRR*19], we did not employ the orientation term in the energy.Indeed, with a high number of landmarks as in our setup, the symmetry ambiguities are easily solved by the functional maps pipeline.

Figure 2 :
Figure2: Schematic of the main steps involved in our method to map a source shape (orange) to a target shape (blue) as described in Sec. 3.

Figure 3 :
Figure 3: Several Dirichlet Laplacian eigenfunctions of the annulus in 2D with external radius 1 and internal radius 1/2.Notice that the eigenfunctions concentrate away from the boundary.

Figure 4 :
Figure 4: Some Dirichlet-Steklov eigenfunctions of the same annulus from Fig. 3.The Steklov boundary condition is imposed in turn on the external (top row) and internal (bottom row) boundaries.Notice that the eigenfunctions concentrate on the Steklov boundary.

Figure 5 :Figure 6 :
Figure 5: W -inner products for the first 20 Dirichlet-Steklov eigenfunctions corresponding to six landmark circles on a sphere mesh (left).The first three landmarks are on the top left and the remaining three are on the bottom right.Notice that the different H i subspaces are almost orthogonal.

Figure 7 :
Figure 7: Dirichlet-Steklov (top row) and Dirichlet Laplacian (bottom row) eigenfunctions on an annulus with three landmark circles and Neumann conditions on inner and outer boundaries.The Dirichlet-Steklov eigenfunctions correspond to the landmark on the right.Notice that as the eigenvalues increase, the Dirichlet-Steklov eigenfunctions quickly concentrate around the corresponding landmark, unlike the Dirichlet Laplacian eigenfunctions that remain distributed in the bulk of the annulus.

Figure 8 :
Figure 8: Harmonic function satisfying Eq. (15) corresponding to the top landmark on a disk with three landmark circles.Notice that the gradient of the function roughly points towards the top landmark.Hence, its normal derivatives at the bottom landmark circles can be seen as specifying the direction towards the top landmark.

†Figure 9 :
Figure9: Impact of the r f parameter on the shape matching quality.The mean geodesic error is averaged on the 95 shape pairs of the TOSCA non-isometric dataset (remeshed to 5K vertices).Notice how stable our method remains, even for extreme values of r f .

Figure 10 :
Figure 10: Error summary when increasing the number of landmarks k for different surface sampling strategies.The mean geodesic error on 10 cat shapes of the TOSCA Isometric dataset is reported."FPS" stands for Farthest Point Sampling.
Figure 11: Error summary when increasing the number of landmarks k for two surface sampling strategies.The mean geodesic error on 95 shape pairs of the TOSCA non-isomatric dataset with 3 different seed initializations for each pair is displayed."FPS" and "std dev." respectively stand for Farthest Point Sampling and standard deviation.

Figure 12 :
Figure12: Left: remeshing stability when varying the triangle reduction factor r of the target shape.The geodesic error, averaged over 300 test pairs of the FAUST data set, slightly increases when the target mesh becomes coarse (low value of r).Right: stability of our method when performing resmeshings on the FAUST dataset (Remeshed to 5K vertices and FAUST "Wild" (see Sec. 7.1) ).The geodesic error is measured in mean geodesic distance ×100 after normalizing by the geodesic diameter.The mean values, mean execution times and vertex counts for each remeshing is presented in Tab. 1.

Figure 14 :
Figure 14: Error summary on the FAUST (left) and TOSCA Isometric dataset (right).The geodesic error is measured in mean geodesic distance ×100 after normalizing by the geodesic diameter.

Figure 15 :
Figure15: Qualitative evaluation of our method and competing approaches on isometric shapes.The first row corresponds to shapes from the FAUST data set.The bottom row consists of shapes from the TOSCA isometric data set.The shape pair is selected such that the geodesic error of our method is median over the dataset.The best and worst cases are illustrated in App.G.

Figure 18 :
Figure18: Error summary on 165 shapes of the SHREC'19 data set.The average geodesic error (Av.Geo.Err.) is displayed for our method and competing approaches.

Figure 20 :
Figure 20: Illustration of the creation of a landmark boundary.The landmark position is indicated by a green dot.triangles composing the landmark disk are shown in light red.The boundary circle is highlighted as a red line.Note that a gap of connectivity appears when creating the boundary around the landmark.This gap is closed when the process finishes producing the boundary.

Figure 21 :
Figure21: Left: comparison of initializations for our method, where "Norm.De." and "Conf.En." respectively stand for "Normal Derivative" and "Conformal Energy".Right: comparison of the "fast" and "principled" energy formulations of our method.Both experiments are performed on the SHREC'20 lores dataset (partial shapes excluded).

Finally
, in Fig.26, we show the best and worst pairs for the SHREC'19 benchmark.

Figure 23 :Figure 24 :Figure 25 :Figure 26 :
Figure23: Qualitative evaluation of our method and competitors on isometric shapes from the FAUST dataset (top row) and the TOSCA isometric dataset (bottom row).The shape pair is selected such that the geodesic error of our method is the worst over the dataset.

Figure 27 :
Figure 27: Weight study on the SHREC'20 data set (full shapes remeshed to 1K vertices).The error measure is the mean geodesic error, averaged on the data set.a C , a P and a I are the Conformality, Properness and Invertibility weights.On the left, we fix the conformality weight a C and vary the properness and invertibility weights a P and a I .On the right, we vary one weight a C/P/I and fix the remaining weights either to 0 or to 1.

Figure 28 :
Figure 28: Weight study on the SHREC'20 data set (full shapes remeshed to 1K vertices).The error measure is the Dirichlet energy, averaged on the data set.a C , a P and a I are the Conformality, Properness and Invertibility weights.On the left, we fix the conformality weight a C and vary the properness and invertibility weights a P and a I .On the right, we vary one weight a C/P/I and fix the remaining weights either to 0 or to 1.

Figure 30 :
Figure30: Average of the absolute values of the inner product matrix of each shape in the SHREC'19 data set.Except for the first few Dirichlet-Steklov eigenfunctions, the off-diagonal inner products are negligible.This validates the approximation of orthogonality.We highlight that this computation also sheds light on the robustness of our basis computation to complex triangulation and partiality setups.

Figure 31 :
Figure 31: Effect of varying the size of our basis on the G(M) space (left) and on the H j (M) space (right).Both figures are an average over all pairs of the TOSCA non-isometric dataset.

Table 1 :
Stability of our method when performing resmeshings on the FAUST dataset.The geodesic error (geo.err.) is measured in mean geodesic distance ×100 after normalizing by the geodesic diameter.The corresponding error curves are displayed in Fig.12(right).The execution time (exec.t.) is also reported, along with the mean number of vertices for each remeshing type (nv).submitted to COMPUTER GRAPHICS Forum (6/2022).

Table 2 :
. Our method's strong Quantitative evaluation results on the remeshed FAUST and TOSCA Isometric (TOSCA Iso.) data sets.The average geodesic error (Av.Geo.Err.) and average execution time (Av.Time) on both data sets are displayed for our method and competing approaches.
FAUST TOSCA Source HyperOrb WA FMap ZO Ours

Table 4 :
Quantitative evaluation results on the SHREC'20 lores (without partial shapes) data sets.The average geodesic error (Av. Geo. Err.), the Dirichlet energy (Dir.E.) and average execution time (Av.Time) on both data sets are displayed for the three initialization methods that we tried: Trivial, Conformal Energy ("Conf.En.") and Normal Derivatives ("Norm.De.").Normal Derivatives is the method used in the rest of the paper.

Table 5 :
Average geodesic error (Av.Geo.Err.) and average execution time (Av.Time) associated to the comparison of the "principled" and "fast" computation methods.submitted to COMPUTER GRAPHICS Forum (6/2022).

Table 6 :
-Zero Weight Av.Geo.Err.Dir.E. Quantitative evaluation results on the SHREC'20 data set (full shapes remeshed to 1K vertices) when fixing one weight to 1 (Non-Zero Weight) and setting the remaining weights to 0. The average geodesic error (Av.Geo.Err.) and Dirichlet Energy (Dir. E.) is given for each.