Manish Mandad

Variance-Minimizing Transport Plans for Inter-surface Mapping

Keywords: CCS Concepts:, Computing methodologies → Shape analysis; Inter-surface mapping, optimal transport, Sinkhorn iterations, di usion geometry, area distortion, conformality

We introduce an e cient computational method for generating dense and low distortion maps between two arbitrary surfaces of same genus. Instead of relying on semantic correspondences or surface parameterization, we directly optimize a variance-minimizing transport plan between two input surfaces that de nes an as-conformal-as-possible inter-surface map satisfying a user-prescribed bound on area distortion. The transport plan is computed via two alternating convex optimizations, and is shown to minimize a generalized Dirichlet energy of both the map and its inverse. Computational e ciency is achieved through a coarse-to-ne approach in di usion geometry, with Sinkhorn iterations modi ed to enforce bounded area distortion. The resulting inter-surface mapping algorithm applies to arbitrary shapes robustly, with little to no user interaction.

INTRODUCTION

Finding a map between two surfaces is a recurring task in geometry processing, required by a wide range of applications including shape matching, template tting, attribute and animation transfer, shape completion, blending, morphing, and remeshing. It consists in establishing a meaningful correspondence map between two input surfaces, and has been referred to as inter-surface mapping [START_REF] Schreiner | Inter-surface Mapping[END_REF]], cross-parameterization [START_REF] Kraevoy | Cross-parameterization and Compatible Remeshing of 3D Models[END_REF], consistent parameterization [START_REF] Asirvatham | Consistent Spherical Parameterization[END_REF], and shape alignment in the graphics literature.

If the two surfaces are isometric, e cient algorithms exist to automatically nd a locally isometric map between these shapes. When the two surfaces belong to the same speci c class of shapes, maps based on semantic correspondences (where semantic features in one shape are associated to the corresponding features in the other shape) can also be constructed [START_REF] Oliver Van Kaick | A Survey on Shape Correspondence[END_REF]; recently, approaches that learn what a user considers to be good correspondences have even been proposed [START_REF] Boscaini | Learning Shape Correspondence with Anisotropic Convolutional Neural Networks[END_REF]]. However, in Fig. 1. Inter-surface mapping. Our approach generates a dense, lowdistortion correspondence map between non-isometric surfaces through a geometrically-derived transport map minimizing local variances of associated neighborhoods. Only two user-defined constraints were prescribed in these models, at the tip of diametrically opposed limbs. Colors are used to indicate correspondences between all these models. most practical cases where shapes di er signi cantly, there is no universal consensus on what de nes the best correspondence map. In this more challenging case of arbitrary inputs, one typically invokes a geometric prior instead, such as conformality or area-preservation, to narrow down the type of suitable maps to seek.

In this paper we propose a computational approach for generating inter-surface maps that are bijective and have some degree of continuity, i.e., small neighborhoods should map to bounded neighborhoods. Our method involves the optimization of a transport map between the two surfaces, which de nes an as-conformal-aspossible point-to-point map with bounded area distortion. These resulting geometric properties make our inter-surface maps particularly useful to subsequent geometry processing.

Previous Work

Previous geometric approaches can be classi ed by the type of mapping (point-to-point vs. soft), its properties (bijective, locally injective), its measure of distortion (isometric, conformal) and the optimization strategy employed to nd the mappings. We review the most relevant approaches next.

Embedding spaces. When dealing with near isometric maps, a practical way to simplify matching is to embed the input manifolds into some target space in such a way that corresponding points of di erent isometric shapes are mapped to nearby points in the target space. A popular way of building such embeddings is to use spectral geometry [START_REF] Bronstein | Generalized multidimensional scaling: A framework for isometry-invariant partial surface matching[END_REF][START_REF] Gallot | Embedding Riemannian Manifolds by Their Heat Kernel[END_REF][START_REF] Shtern | Spectral gradient elds embedding for nonrigid shape matching[END_REF][START_REF] Vestner | Product Manifold Filter: Non-rigid shape correspondence via kernel density estimation in the product space[END_REF]. Another approach in this family of methods use a common template mesh as embedding space. The base domain is often chosen to be a coarse mesh obtained [START_REF] Schreiner | Inter-surface Mapping[END_REF]] through mesh decimation , or using overlapping domains [START_REF] Kraevoy | Cross-parameterization and Compatible Remeshing of 3D Models[END_REF] bounded by consistent paths connecting feature points . The mapping is then constructed by optimizing a mutual tessellation of the two surfaces. While this approach can handle non-isometric transformations, decreasing the global non-linear distortion through vertex-by-vertex updates of the map is often very slow. Instead, our approach directly optimizes a transport plan over the two surfaces, bypassing the need for embedding spaces.

Parameterization. Surface mappings can be computed via mesh parameterization as well. Recent approaches o er to construct bounded distortion maps [START_REF] Lipman | Bounded distortion mapping spaces for triangular meshes[END_REF]], angle-preserving maps with as uniform as possible conformal factors [A alo et al. 2013], using cone manifolds [START_REF] Myles | Controlled-distortion constrained global parametrization[END_REF] or orbifolds [Aigerman and Lipman 2015], or even propose to relax the bijective condition of the parameterization to local injectivity [START_REF] Aigerman | Lifted bijections for low distortion surface mappings[END_REF]]. One main limitation of these parameterization-based approaches is that the map distortion is not directly measured between the two input surfaces, but through a common base domain, which leads to suboptimal inter-surface maps. Mapping to a common domain also often requires designing a cut-graph which increases distortion as well (see [Aigerman et al. 2015] for a recent improvement). In contrast, our approach directly optimizes a surface-to-surface map without resorting to an intermediate domain. Moreover, while A alo et al. [2013] sought conformal maps with no (or bounded) area distortion, we look instead for area-preserving (or bounded area distortion) maps that are as conformal as possible.

Soft maps and optimal transport. While point-to-point mapping between surfaces has been favored in the past, soft maps [START_REF] Solomon | Soft Maps Between Surfaces[END_REF]] and functional maps [START_REF] Ovsjanikov | Functional maps: a exible representation of maps between shapes[END_REF]] (where functions or distributions are mapped between surfaces) have been shown relevant to the analysis of correspondences through linear algebra tools. Using a mass transport formulation, Solomon et al. [2012] generate transport maps that are continuous, faithful to geometric descriptors and soft. Continuity of such maps means that nearby regions on one surface should map to nearby mass distributions on the other surface, where "nearby" is to be understood with respect to the geodesic earth mover's distance. In addition, softness of the map is e ciently optimized by L 2 maximization of the spreading of the probabilities in the mapping matrix describing the mass transport plan. The manner in which a soft map transports one distribution to another can also help control map continuity [START_REF] Solomon | Dirichlet Energy for Analysis and Synthesis of Soft Maps[END_REF]]. [START_REF] Corman | Continuous Matching via Vector Field Flow[END_REF] nd inter-surface mappings by integrating vector elds so as to minimize a soft map based energy that favors isometry. Note that softness of transport maps can also be used as a means to accelerate the computation of optimal transport plans: Cuturi devised an entropic regularization of transport which reduces optimal transport computational times by several orders of magnitude, at the price of smoothed-out transport plans [START_REF] Cuturi | Sinkhorn Distances: Lightspeed Computation of Optimal Transport[END_REF][START_REF] Solomon | Convolutional Wasserstein Distances: E cient Optimal Transportation on Geometric Domains[END_REF]. [START_REF] Frogner | Learning with a Wasserstein Loss[END_REF] further extended smoothed transport to unnormalized measures, a concept also referred to as relaxed or unbalanced transport [START_REF] Chizat | Unbalanced Optimal Transport: Geometry and Kantorovich Formulation[END_REF][START_REF] Chizat | Scaling Algorithms for Unbalanced Transport Problems[END_REF]. In sharp contrast, our approach avoids soft maps altogether, and entropic regularization is used sparingly and for numerical stability only; computational scalability is, instead, gained by optimizing the transport plan in a coarse to ne fashion in di usion space (i.e., based on a small fraction of eigenvectors of the Laplace operator).

Gromov-Wasserstein distances. [START_REF] Mémoli | Gromov-Wasserstein Distances and the Metric Approach to Object Matching[END_REF] pioneered the use of the Gromov-Wasserstein distance for deriving a metric approach to the shape matching problem. [START_REF] Solomon | Entropic Metric Alignment for Correspondence Problems[END_REF] computed correspondences through an entropy regularized Gromov-Wasserstein objective function, using the Softassign algorithm [START_REF] Rangarajan | Convergence Properties of the Softassign Quadratic Assignment Algorithm[END_REF]. Their objective function tries to pair points with similar geodesic distances on the two shapes. However, the solver involved in this approach loops through compute-intensive iterations of cubic complexity in the number of points to be matched, and a matrix of pairwise geodesic distances must be stored, rendering the construction of a dense map intractable. In addition, the resulting maps in the Gromov-Wasserstein approach can have one-to-many correspondences (see Fig. 20). While this property may be valuable in speci c applications, it does not provide the type of one-to-one mappings that geometry processing tools typically requires. We deviate from their approach by formulating a very di erent notion of distortion, related to the symmetric Dirichlet energy of transport maps, that naturally induces point-to-point maps. Moreover, our mapping algorithm has nearly linear complexity, thus allowing the construction of inter-surface maps over two orders of magnitude faster in practice than these existing approaches.

Contributions

In this work we introduce an approach for nding inter-surface maps based on a mass transport formulation that minimizes the variance of the images of local neighborhoods. Variance-minimizing transport plans, which we show to be minimizers of the symmetric Dirichlet energy of maps in the continuous limit, inherently favor conformality and induce continuity and sharpness of the resulting map while enforcing bounded area distortion. Although this new formulation becomes non-convex and hence substantially more challenging to minimize, we reformulate it as a biconvex minimization problem. A coarse-to-ne approach based on a hierarchy built in di usion geometry is then devised to e ciently optimize the map through alternating minimizations and reduce the risk of getting trapped in local minima. We demonstrate our robust and scalable algorithm for inter-surface mapping on a series of non-isometric models (see, e.g., Fig. 1).

APPROACH

Our approach stems from the simple idea that given two pointsets X = {x i } and Y = { j } sampling two input surfaces, we should aim at constructing a map between X and Y that associates neighborhoods from one pointset with neighborhoods in the other pointset, and vice-versa. We formulate this inter-surface mapping as an optimal mass transport problem between the two surfaces that minimizes the variance of the image of each neighborhood.

Setup

In order to formulate our map optimization problem, we introduce a few key concepts that we will leverage throughout our exposition. Surfaces as distributions. Our approach considers each pointset as a discrete mass distribution. More speci cally, we discretize the normalized area measures of the two surfaces as weighted sums of Dirac measures centered at input points, i.e., µ = m i δ x i and Variance-Minimizing Transport Plans for Inter-surface Mapping • 39:3 ν = n j δ j , where the masses m i and n j are normalized so that each surface has unit total mass, i.e., i m i = j n j = 1. The mass of each point is assigned based on the estimated local Voronoi area of the surface it covers to account for non-uniform sampling.

Maps as transport plans. Given two distributions µ and ν , a transport plan π between them is by de nition a matrix of size |X|×|Y| whose marginals equal µ and ν . Each entry π i j represents the amount of mass transported from x i ∈ X to j ∈ Y. A transport plan π has two associated transfer operators π X and π Y de ned as

π X (δ x i ) = 1 j π i j j π i j δ j , π Y (δ j ) = 1 i π i j i π i j δ x i .
Note that the denominators are respectively equal to m i and n j . These two linear operators describe how mass distributions on X are mapped to mass distributions on Y with the same total mass, and vice-versa. Therefore, they contain essentially the same information as the transport plan.

Weighting functions. For every point x on a surface, we associate an isotropic weighting function W x centered at x (de ned as either a heat kernel based function in Sec. 3.2 or a Gaussian function in Sec. 4.2) for two purposes: rst, its nite support de nes a notion of local neighborhood; second, weighting functions can be used as local test functions to measure the variance of an inter-surface map, as we explain next.

Variance-minimizing Transport Plans

Our approach consists in computing the transport plan π between µ and ν that minimizes the following variance-based transport cost:

E (π ) = X var π X W x µ mass(W x µ) dµ (x ) + Y var π Y W ν mass(W ν ) dν ( ).
(1)

That is, this cost evaluates the integral of the variance of the image, by the transfer operator and its inverse, of the area measure modulated by local weighting functions (see inset below: (left) original measure, then (middle) modulated by weighting functions, and (right) its corresponding measure through the plan).

• • • • • • X W x µ µ π X Y x x • • • • • • • • • • • • Properties.
The reason for targeting such a variance-minimizing plan is two-fold. First, penalizing the variance of images favors transport plans corresponding to actual point-to-point maps rather than soft, di used maps: any point mass that is mapped to a spread-out distribution incurs a signi cant cost. Second, variance-minimizing maps tend to be bijective, since we also account for the inverse map in the cost evaluation. Bi-continuity is also expected of the minimizers of our cost as we penalize variances of images of small neighborhoods de ned by the support of our weighting functions. Moreover, we will show that our variance-based formulation of the transport cost is amenable to e cient plan optimization.

Discrete symmetric Dirichlet energy. The cost function E (π ) of a plan π can be understood a discrete approximation of the symmetric Dirichlet energy of the map and its inverse: it coincides with the Dirichlet energy for in nitesimal weighting functions with isotropic covariance matrices tending to zero. More speci cally, let f : S → S be a di eomorphism between two surfaces. Choose the weighting function at x ∈ S to be a Gaussian distribution with covariance σ 2 I , where σ is a small positive number. Let be a random vector drawn from this distribution. Then the integrand of E at x is by de nition the variance of f ( ). For small σ , we may further approximate locally the two manifolds by their tangent spaces. Under this approximation f ( ) is a Gaussian with covariance

σ 2 (D f )(D f ) , hence var( f ( )) = 2σ 2 ||D f || 2 F .
As a consequence, the total cost is proportional to the sum of the Dirichlet energy of f and the Dirichlet energy of its inverse. Note that the same conclusion holds for non Gaussian weighting functions, as long as they have isotropic covariances with total variance independent of x, and tending to 0. While generalization of the Dirichlet energy to transport plans have been proposed in the literature (see, e.g., [START_REF] Solomon | Dirichlet Energy for Analysis and Synthesis of Soft Maps[END_REF]]), we are not aware of any existing formulation that also penalizes soft maps, a property that is essential for our purposes as it implies that our optimal plans can be expected to be close to conformal if no other constraints are imposed.

Optimizing Transport Plans

While the formulation of transport cost given above induces good geometric properties for the optimal map, it involves a non-convex functional which may seem di cult to minimize e ciently and robustly. We can, however, introduce auxiliary variables within the variance terms so that nding a variance-minimizing map amounts to solving a biconvex problem, which we achieve through alternating convex minimization.

Reformulation with auxiliary variables. Given a point x and a measure µ = µ i δ x i , we de ne the variance of µ with respect to x as:

var(µ, x ) = i µ i d (x i , x ) 2 , (2)
where distances d (•, •) are computed in a (possibly high-dimensional) Euclidean embedding of the surface X that µ samples, such that the minimum of this variance is found through a simple linear average Fig. 2. Mapping via variance minimization. Le : Two input pointsets forming a line and a curve with a loop respectively (rainbow and grey colormaps depict their geodesic connectivity). Middle le : initial optimal transport map using the 2-Wasserstein distance (correspondences in thin lines). Middle right: our variance-minimizing map a er alternating minimization iterations using ambient, Euclidean distances. Right: result a er minimization using geodesic (i.e., intrinsic) distances instead.

of positions in this embedding. Note that using geodesic distances in the de nition of the variance would be computationally intractable beyond the case of curves, while using R 3 distances would be poorly adapted as it only provides information on the extrinsic structure of the surface; we will thus use a "di usion" embedding (see Sec. 3.1) which is fully intrinsic, but as simple to evaluate as the R 3 distanceand short distances in this embedding will be good approximations of geodesic distances on X, which enforces that our approximation will be accurate for ne enough sampling. As we will see later, these embedding can further be used in a multiscale way, which will lead to better optima and increased e ciency as well. With this notion of variance w.r.t. a point location, we can now introduce a set of point locations η x i and η j , one for each point x i and j on the two surfaces. These auxiliary variables, which we call centers for reasons that will become clear soon, allow us to rewrite our transport cost between X and Y as a function of both the transport plan π and the set of all centers η:

C (π , η) = i m i var π X W x i µ k m k W x i (x k ) , η x i + j n j var π Y W j ν k n k W j ( k ) , η j .
(3)

Note that minimizing this cost functional with respect to η assigns each center η x i (resp., η j ) to the (weighted) center of mass of the image of the normalized weighting functions W x i (resp. W j ) by the transport plan π . For such values of η, the cost functional is precisely equal to the energy de ned in Eq. ( 1).

Minimization through alternating convex problems. This augmented cost functional C (π , η) is now amenable to an alternating minimization similar to Expectation-Maximization (EM): for a xed transport plan π , minimizing C (π , .) is e ciently achieved by relocating the centers to the barycenters of the images of the weighting functions by π (see Sec. 3.2); for xed centers, minimizing C (., η) requires the computation of an optimal transport plan under the usual mass preservation constraints. Alternately minimizing the cost C over π i j and η as described in Alg. 1 is a robust and e cient way to treat the original non-linear problem since both minimization problems are convex. Fig. 2 shows the transport plan after alternating minimization iterations for an example between two curves, where we use a Gaussian weighting function as our local test function to measure variance. Note that the number of iterations needed to properly "unfold" one curve onto the other depends on the width of this function: if the weighting function used as a test function is too narrow, the variance is only measured locally and many iterations are needed as Fig. 3 shows. This observation calls for a coarse to ne treatment where alternating minimizations are performed at multiple scales to improve e ciency, as we detail later.

Additional control over the mapping. The cost we de ned in Eq. ( 3) measures the regularity of a transport plan via its local variance. However, user control over the nal mapping is often desirable: the user may want to prescribe a few correspondences manually. If x i must map to j , we simply constrain the center η x i to be j and η j to be x i . Additionally, one can add to the transport cost a small amount of the standard 2-Wasserstein distance between the surfaces to help the solver disambiguate symmetries, see Sec. 5.

NUMERICAL TOOLS

Our geometric mapping algorithm relies on the frequent evaluation of distances and barycenters on surfaces, as well as the repeated evaluations of optimal transport plans. Before delving into the details of our algorithm, we discuss a few numerical tools we will rely on heavily; in particular, the use of di usion geometry and Sinkhorn iterations will help make our approach both robust and e cient.

Di usion Geometry

The concept of di usion geometry introduced in [START_REF] Coifman | Geometric di usions as a tool for harmonic analysis and structure de nition of data: Di usion maps[END_REF][START_REF] Nadler | Di usion Maps, Spectral Clustering and Eigenfunctions of Fokker-Planck Operators[END_REF] proposes to embed a Riemannian manifold into a "di usion space" whose coordinates are derived from the eigenvectors and eigenvalues of the celebrated Laplace-Beltrami operator. In our context, this intrinsic di usion operator ∆ is approximated on each pointset and we denote its M smallest eigenvalues by {λ k } and its corresponding eigenfunctions by {ϕ k } with ∆ϕ k = λ k ϕ k . Then a point x ∈ X is projected into (truncated) di usion space as x through:

x Φ t -→ x = Φ t (x ) = e -λ 1 t /2 ϕ 1 (x ) e -λ 2 t /2 ϕ 2 (x ) . . . e -λ M t /2 ϕ M (x ) , ( 4
)
where t controls the di usion time scale. For two points x 1 and x 2 , the corresponding heat kernel K t (x 1 , x 2 ) is expressed as:

K t (x 1 , x 2 ) = Φ t (x 1 ) Φ t (x 2 ), (5)
while the di usion distance between these points is computed as:

d di (x 1 , x 2 ) = Φ t (x 1 ) -Φ t (x 2 ) , (6)
or

d 2 di (x 1 , x 2 ) = K t (x 1 , x 1 ) + K t (x 2 , x 2 ) -2K t (x 1 ,
x 2 ) equivalently. For small t, di usion distances can be shown to approximate local geodesic distances very well, while increasing t has a smoothing e ect by diminishing the contribution of high frequency eigenfunctions [START_REF] Coifman | Geometric di usions as a tool for harmonic analysis and structure de nition of data: Di usion maps[END_REF]. In order to keep the computational cost low, we select only the rst M smallest eigenvalues of the Laplace-Beltrami operator (M = 100 in practice, unless indicated otherwise).

Leveraging Di usion Geometry

One of the key points of our approach is the use of di usion geometry to approximate a number of important geometric notions.

Di usion centers. The geodesic barycenter of a weighted neighborhood de ned on a surface is in general not expressible in closed form and may not even be unique if the neighborhood is large. A crucial reason for using di usion distance instead is that we can compute weighted di usion barycenters in closed form as simple weighted averages of points in di usion space as discussed in Sec. 2.3. This property brings signi cant speedup since converting a point to di usion space is done e ciently via Eq. ( 4). Moreover, because of their intrinsic nature, di usion coordinates are much preferable to Euclidean coordinates in R 3 as they are more robust to averaging of distant points. These di usion centers will thus be used in lieu of the true geodesic centers.

Weighting Functions. If distances and variances are computed in di usion space, a consistent choice for the weighting function W x is to rely upon the heat kernel in Eq. ( 5) in order to de ne an approximate geodesic neighborhood. While this choice is nicely isotropic for small di usion times t, its support becomes o -centered for large values of t (see Fig. 4). We use instead the weighting function:

W x (x i ) = 2K t (x, x i ) K t (x, x ) + K t (x i , x i ) , (7)
which remains localized and centered for any di usion time while still satisfying W x (x ) = 1. In addition, we truncate the weighting function below ε = 0.5 (Fig. 4) in order to further localize its associated geodesic neighborhood, which will reduce computational times required to evaluate geodesic variances.

Projection to tangent planes. As part of our algorithm, local projections to tangent planes will be necessary in Section 4.4 to check map consistency. While this is easily achieved for a ne sampling of the surface in Euclidean space, the use of di usion geometry requires the same operation to be performed directly in di usion space. To nd the two tangent vectors forming the di usion tangent plane at a point x on the surface, we proceed as follows. We rst assemble the matrix of size M×3 for which the k-th row is ∇ϕ k (x ). Then, for two approximate tangent vectors u and in Euclidean 3D space (e.g., computed based on two nearest neighbors of x), we multiply this matrix by each of these vectors to get the di usion tangent vectors û and ˆ , that we further orthonormalize. Any point ẑ in di usion space can now be directly projected onto the 2D di usion tangent plane at x through ( ẑ û, ẑ ˆ ). Note that this simple projection to the di usion tangent plane through inner products provides a consistent de nition of the notion of linear approximation of the surface even for coarse sampling, adding signi cant robustness: de ning a tangent space in Euclidean space for a coarse sampling of a surface would be much less reliable.

K t (x,x i ), K t (x,x i )/ K t (x,x )K t (x i ,x i ) and 2K t (x,x i )/(K t (x,x )+K t (x i ,x i ))
for x indicated as the purple vertex. Its neighborhood, i.e., vertices x i for which W x (x i ) > ε, is shown with green spheres proportional to function values.

Sinkhorn Iterations

Since our approach relies on repeated optimizations of transport plans, an e cient solver for optimal transport problems is mandatory. In order to avoid having recourse to Linear Programming, recent work [START_REF] Benamou | Iterative Bregman Projections for Regularized Transportation Problems[END_REF][START_REF] Cuturi | Sinkhorn Distances: Lightspeed Computation of Optimal Transport[END_REF] proposed an entropic regularization of the optimal transport problem which minimizes the Kullback-Leibler (KL) divergence:

min KL(π |K γ ) def. = π ln π K γ -1 dxd , (8)
where γ controls the amount of regularization and K γ denotes the kernel associated to the cost function. This regularized optimal transport map can be e ciently computed using Sinkhorn's matrix row/column normalization algorithm which is several orders of magnitude faster than solving a linear program.

Sinkhorn algorithm. Following [START_REF] Benamou | Iterative Bregman Projections for Regularized Transportation Problems[END_REF][START_REF] Solomon | Convolutional Wasserstein Distances: E cient Optimal Transportation on Geometric Domains[END_REF], given a matrix H discretizing the kernel K γ , the transport plan minimizing Eq. ( 8) between two datasets X and Y is of the form π = D m D v HD w D n , for diagonal matrices D v and D w de ned by vectors v, w that satisfy:

     D v HD w n = 1 |X | , D w H D v m = 1 |Y | , (9)
where H is of size |X| × |Y|, m and n represent the mass (normalized Voronoi area) vectors of X and Y, and 1 k is the vector of ones of size k. Notice that the properties π 1 |Y | = m and π 1 |X| = n are enforced by Eq. ( 9). We nd (v, w) via Sinkhorn iterations [START_REF] Sinkhorn | A Relationship Between Arbitrary Positive Matrices and Doubly Stochastic Matrices[END_REF]] as detailed in Alg. 2, where and ⊗ denote entrywise division and multiplication, respectively. Convergence is assumed when the change in KL divergence is less than 0.1%.

Algorithm 2 Sinkhorn iterations to solve for optimal transport

1: function S 2: v, w ← 1 3: while !converged do 4: v ← 1 |X | H (n ⊗ w) 5: w ← 1 |Y| H (m ⊗ v) 6: return π = D m D v HD w D n
Transport plan Working in log space. The original Sinkhorn approach assembles the kernel matrix H through its elements H i j = exp(-coef(π i j )/γ ), where coef(π i j ) refers to the coe cient of π i j in the cost function (3). The regularization coe cient γ used to compute the optimal transport is recommended to be set such that max{coef π i j }/γ < 200 [START_REF] Cuturi | Sinkhorn Distances: Lightspeed Computation of Optimal Transport[END_REF]] in practice, to prevent numerical blowups due to limited oating-point dynamic range. However, this choice of coefcient usually creates a di use transport plan, con icting with our objective to construct variance-minimizing plans. In our approach, we sidestep this numerical limitation by storing the coe cients of the kernel matrix in log scale and compute the Sinkhorn iterations in log scale directly (see also the concurrent work of [START_REF] Schmitzer | Scaling Algorithms for Unbalanced Transport Problems[END_REF]] for a similar approach). This implies that we now store Working in log scale and subtracting the maximum element drastically reduces the negative e ects of a limited oating-point dynamic range, allowing us to safely set to γ = max{coef π i j }/10 -4 , a 50-fold decrease in regularization. In order to further improve e ciency, we also use the previous values of v and w instead of initializing them to 1 in each successive relaxation iterations.

H i j = -coef(π i j )/γ ,

Mass-Relaxed Sinkhorn Iterations

When dealing with a very non-isometric pair of shapes, enforcing strict mass preservation leads to mappings with high anisotropic di erentials. This anisotropy may become signi cant enough locally to render the optimal map ill-suited for further geometry processing. Relaxing the mass preservation constraint helps remove this issue, while providing a simple way to bound area distortion.

Adapting Sinkhorn iterations.

A way to relax the mass preservation constraint in the framework of entropy-regularized optimal transport is to solve the following optimization problem:

min π ∈Π KL(π |K γ )
where Π denotes the set of matrices satisfying:

π i j ≥ 0, α i m i ≤ j π i j ≤ β i m i , i, j π i j = 1, α j n j ≤ i π i j ≤ β j n j ,
where 0 ≤ α k ≤ 1 ≤ β k and 0 ≤ α k ≤ 1 ≤ β k are the prescribed lower and upper bounds on the allowed mass distortion on each surface.

We write the convex set of allowed solutions as the intersection of two convex sets, namely the set Π 1 of matrices satisfying the rst, second and third conditions above, and the set Π 2 of matrices satisfying the rst, second and fourth conditions. Optimizing over their intersection can be done using Dykstra's algorithm adapted to KL divergence in [START_REF] Heinz | Dykstra's algorithm with Bregman projections: a convergence proof[END_REF][START_REF] Benamou | Iterative Bregman Projections for Regularized Transportation Problems[END_REF]. This algorithm only requires projections of distributions on each convex set (in the sense of KL divergence), which we derive in App. A. It turns out that the KL projections on Π 1 and Π 2 boil down to row and column rescaling respectively; because of these particularly simple forms, the auxiliary variables used in the original Dykstra's algorithm do not play any role, and the overall algorithm reduces to iterated KL projections on Π 1 and Π 2 . Also, since only row and column rescaling is used, we may represent the transport plan π at each iteration as a product of matrices associated to vectors v and w as in the usual Sinkhorn algorithm. Pseudocode for the algorithm is given in Alg. 3, and related functions are in Alg. 4. r ← solution to G _M _ (lb, ub)

6: v ← U _V _G _M (r, lb, ub) 7: lb ← α H (m ⊗ v) 8: ub ← β H (m ⊗ v) 9:
r ← solution to G _M _ (lb, ub) for all i do 3:

if lb[i] ≤ r ≤ ub[i] then x[i] ← r 4: else if ub[i] < r then x[i] ← ub[i] 5: elsex[i] ← lb[i] 6: function r = G _M _ (lb, ub) 7:
x(r ) ←U _V _G _M (r , lb, ub) x(r ) ←U _V _G _M (r, lb, ub)

11: r ← solution to m D v HD x n = 1
Note that each KL projection requires solving an equation to compute parameter r , which must be chosen such that the resulting transport plan has global mass equal to 1. Fortunately, the global mass is piecewise linear and non decreasing as a function of r , with a number of nodes equal to twice the number of points in each dataset. Hence, solving for r can be done by simply sorting these nodes. In the rare case that no solution exists, we set r to a large magnitude, with the appropriate sign.

Three-fold Alternating Minimizations. When bounded area distortion is allowed, our alternating minimization algorithm described in Sec. 2.3 cannot be directly extended to the case where we relax mass preservation: the denominator in our de nition of transfer operators makes the transport plan optimization step more di cult, and seemingly not reducible to a convex problem. To solve this issue, we introduce two additional variables m i and n j , conceptually corresponding to respectively j π i j and i π i j . Rather than enforcing strict equality, we allow slack by setting the following constraints:

     C -1 m i ≤ j π i j ≤ Cm i , C -1 n j ≤ i π i j ≤ Cn j . ( 10
)
for some constant C close to 1 (set to 1.1 in our experiments). We also impose that D -1 m i ≤ m i ≤ Dm i and D -1 n j ≤ n j ≤ Dn j , where D is a user-chosen bound on the allowed mass distortion. We then modify the transfer operators to:

π X (δ x i ) = 1 m i j π i j δ j , π Y (δ j ) = 1 n j i π i j δ x i .
and optimize the functional C (π , η) obtained from C(π , η) in Eq. ( 3) by just replacing π X by π X and π Y by π Y .

We then proceed by alternately optimizing over the three sets of variables π , (m , n ), and η. To solve for π i j , we use the relaxed Sinkhorn iterations described above. Solving for m and n requires the optimization of a positive linear combination of inverses of their coe cients over a hypercube. This convex optimization problem is solved in O(n log n) time using an algorithm similar to the one used above to enforce global mass preservation. Finally, the update of η remains a simple weighted average as before.

OPTIMAL MAP SOLVER

Equipped with the algorithmic and numerical components we described above, we can now formulate an e cient and scalable solver for our inter-surface mapping approach.

Coarse-to-Fine Approach

As demonstrated by the convergence plot in Fig. 3 of the intercurve map shown in Fig. 2, computational times for our varianceminimizing transport plan based approach may depend heavily on the size of neighborhoods we consider. While large neighborhoods make our alternating minimizations converge faster, the number of variables |X|+|Y| for two surfaces is quite large, inducing a high computational cost despite our use of fast transport plan solvers. Instead, we note that nding proper correspondences between arbitrary surfaces lends itself rather naturally to a coarse-to-ne strategy: it is quite simple to nd a rough map between two simpli ed versions of the surfaces as the number of samples involved is low; then building a more precise map from this rough approximation becomes simpler and less prone to being stuck in a local minimum. Therefore, we operate our transport optimization on a hierarchy of approximations of the two surfaces. By starting the map optimization on a coarse sampling of the surfaces (and thus, a small number of variables) and progressively increasing resolution to re ne the map, we signi cantly reduce the total computation time required to nd a variance-minimizing map. As discussed in Sec. 3.2, using di usion for all S 0 , S 1 . . . S do // traverse hierarchy 5:

π ← A M (S i , π ) 6: π ← F _F _A _T (S i , π ) 7:
return π Mapping between X and Y geometry for low sampling of surfaces is also particularly convenient and reliable as it leverages spectral approximations of the ne, original surfaces. This means that any instance of the positions x, , η x and η is replaced by its equivalent in di usion space, i.e., x, ˆ , ηx or η as detailed in Sec. 3.1. When a transport plan has been found at the nest level of our hierarchy, we then switch to a pure Euclidean stage where now the sampling is ne enough to continue re ning the map using more traditional computations of tangent planes and barycenters. Throughout this hierarchical process, we also check that the map is locally well behaved to avoid degeneracies. In this section, we mostly describe the di usion stage, but brie y mention how to adapt this process to the (simpler) Euclidean stage in Sec. 4.5. Pseudocode for our coarse-to-ne strategy is given in Alg. 5.

Constructing Hierarchy in Di usion Space

We begin with the construction of a hierarchical sampling where the cardinality of a neighborhood, de ned via a kernel-derived weighting function W x via N x = {x i ∈ X :W x (x i ) ≥ ε} (as mentioned in Sec. 3.2) remains constant as we go to ner levels in the hierarchy so as to limit computational complexity. We leverage the smoothing properties of the time scale t in di usion geometry (see Sec. 3.1) to construct a hierarchy of dyadic levels, one for each time scale t = 2 -t 0 . Starting from a large di usion time scale t 0 = ln 10/λ 2 [START_REF] Sun | A Concise and Provably Informative Multi-scale Signature Based on Heat Di usion[END_REF]], we build a subsampling S 0 of the input dataset via farthest point insertion; i.e., we add one point at a time from the original pointset such that it is the farthest from the currently selected subset. We stop adding points to this level once the farthest point p just added satis es W x (p) >τ for the nearest neighboring point x currently in S 0 for a value τ satisfying 1 >τ > ε (we take τ = 0.9 in our experiments), since it indicates that we found a dense enough sampling for that level. This approach ensures that the area contributing to the neighborhood of any point x in the level (namely, the set of points x i satisfying ε <W x (x i ) <τ ) is covered by at most a constant number of Voronoi cells, leading Fig. 6. Di usion distance estimates on thin parts. The black dots depict the final sampling during the last, finest di usion stage. Le : with 100 eigenvectors, sampling is insu icient to get satisfactory precision on thin parts (dog's feet, tail, ears). This may translate into twists of the mapping (middle). Right: using 300 eigenvectors significantly improves precision.

to a constant neighborhood cardinality. We then proceed with the construction of the next level using the same furthest point insertion strategy. As the time scale decreases, levels S become more dense until eventually, no additional points can be inserted: since only a limited number of eigenvectors are used to compute di usion distances, our weighting functions will stop being "discerning" past a certain scale (see Fig. 6). Further levels are then computed using Euclidean distances instead, at which point local distances are almost equal to geodesic distances.

Sparsified Sinkhorn Variables

At the coarsest level of the hierarchy, we can a ord to solve the full optimal transportation problem as we have considerably reduced the number of samples. However, this number steadily increases as we proceed to ner levels. In order to maintain a reasonable computational cost at higher resolutions, we limit the number of variables in the Sinkhorn iterations as follows. Consider the transport plan found so far; since our functional favors homeomorphisms, we can expect that the support of the transport plan is sparse since each point is mapped to only a small number of points as a consequence of variance minimization. Now, we consider a slightly enlarged support containing all pairs (x , ) such that there exists a pair (x, ) in the current support with x ∈ N x and ∈ N . For the next round of Sinkhorn iterations, we limit our transport plan to this enlarged support, which limits the number of variables to instantiate (Fig. 7). Note that this process is repeated within each level of the hierarchy several times during our alternating minimization steps, so the support can gradually move to an optimal location as the centers are relocated-thus eliminating the artifacts that a naive reduction of variables in the transport plan could induce. Once the map stops evolving, we move to the next hierarchy level.

Locality. We leverage the locality of the transport map we are solving for to further accelerate the computation involved in Sinkhorn iterations. At a given level of the hierarchy, a naive implementation of Sinkhorn iterations would use a matrix H of size |S | × |S |, and each iteration would have a complexity of |S | 2 . Instead, we exploit our hierarchical framework and the local instantiation of variables used by Sinkhorn iterations by storing H as a sparse matrix that records only the instantiated variables. All other matrix coecients are set to 10 max{coef π i j } to strongly penalize transport. Sinkhorn iterations are then computed e ciently by rst precomputing i m ⊗ v[i] (resp., j n ⊗ w[j]) in Alg. 2, then adjusting accordingly the value of H (m ⊗ v) computed using the sparse matrix H. As we kept the cardinality of neighborhoods constant during the construction of the hierarchy, both storage size and computational complexity involved in the optimal transport solver are linear in the number of samples. Fig. 8 plots the total transport cost throughout iterations while solving for an optimal map. Note that we display the transport cost normalized per level since comparing cost across levels is not meaningful: sample count, di usion time and neighborhoods have changed. Moreover, while each level leads to a seemingly small decrease in energy, it has a signi cant cumulative impact across the hierarchy.

X Y • • • • • • • • • • • • x W x ≥ ε x ∪N π X (N x ) π X

Preventing Map Flips and Twists through Ironing

While our energy was designed to induce well-behaved maps, artifacts can appear during the coarse-to-ne optimization. In particular, a few, localized discontinuities of the map may appear in the form of twists (180 • rotations of the map) and/or ips (mirrored versions of the map), see Fig. 9: they are typically due to thin protrusions and/or low sampling, leading to an inadequate global map during the iterations. While more eigenvectors (Fig. 6) or more iterations with larger weighting functions at a ner scale may resolve a local aw in the map, it would come at a very high computational cost. We thus track potential map aws preemptively at the end of each level of the hierarchy to x the map early if needed. Noting that the local presence of a twist or a ip will engender a sudden change of normal directions in mapped neighborhoods, we check orientation consistency of the local map around each point x ∈ X as follows. We rst project the immediate neighborhood of x on its tangent plane (computed in di usion space, see Sec. 3.2) and build, in this 2D space, a triangulation of the convex hull of this projected neighborhood. This triangulation is then used as reference to verify the consistency of the orientation of its map by now projecting the neighborhood of ηx on Y on its own tangent plane. If any mismatch in orientation is found, we mark this point as inadequate. This process conservatively identi es the regions where the map is inconsistent and divides our data into separate patches delineated with badly oriented regions between them. We then systematically correct the problem areas by what is best described as an "ironing" process: we proceed with an update of the centers where we arti cially increase the weights of neighborhoods where the map is good; that is, we use an asymmetric weighting function through:

W * x (x i ) =      W x (x i ), x i ∈ visited, ωW x (x i ), otherwise,
where ω is set to 0.5 in all our experiments. This e ectively overpowers any twist or ip present, and brings the solution out of the local minimum in which it was. More speci cally, starting from the patch with maximum area on X with a correctly oriented map, we build a strip around a front (i.e., the delineation between already visited and non-visited points) by uniformly and geodesically thickening it. The corresponding strip is built on Y by geodesically thickening the image of the front on X proportional to the mass of the visited and non-visited parts. We then solve for a variance minimizing transport plan satisfying Eq. (1) restricted to these two strips. Once the centers η stabilize, the front is moved outward and the process is repeated until the regions marked as inadequate have been visited, eventually removing ips/twists. This ironing process is very e cient at detecting and removing artifacts of the map and its inverse even if we intentionally ip the map as a stress test, see Fig. 10.

Final Euclidean Levels

When we reach the nal level of the di usion stage, we have a fairly dense correspondence map between the two surfaces. We, however, continue with our hierarchical solver-now using the Euclidean embedding of the surfaces directly. We create ner levels in the hierarchy by reverting to a plain, unnormalized Gaussian function as our weighting function:

W x (x i ) = exp -(x -x i ) 2 /(2σ 2 ) , ( 11
)
where σ is initialized using the average radius r avg of the neighborhoods in the last di usion stage such that ε = exp(-r 2 avg /2σ 2 ).

To avoid any ill-e ects when σ reaches the sampling resolution of the input data, we always have a lower bound on the cardinality of the neighborhood. Furthermore, for each new Euclidean hierarchy level , σ is directly updated to σ /2. We then proceed to the repeated optimizations of the transport plan at each level as before, but now every point x and are considered in Euclidean space directly, and the centers are computed as linear averages of these Euclidean points (which provide, at this scale, a good approximation of geodesic barycenters). It is important to point out that the Euclidean hierarchy levels do improve results quite signi cantly: using di usion geometry with a limited number of eigenvectors basically ignores surface details corresponding to high frequencies.

As Fig. 11 shows on a map between a sphere and a skull, details are properly adjusted once all the Euclidean levels have been treated, resulting in a sharp and variance-minimizing map as expected.

Fig. 11. Hierarchical map optimization. A sphere mesh (65kV) is mapped onto the input skull mesh (35kV) (le ) by relocating the sphere vertices onto the skull using η, across 3 Euclidean steps of the hierarchical solve (di usion steps are not depicted as the corresponding η's live in di usion space).

EXPERIMENTS

We implemented our algorithm in C++, using the Spectra library [START_REF] Qiu | Spectra: A header-only C++ library for large scale eigenvalue problems[END_REF] to compute the eigenvectors of the Laplace operator, the Eigen library for linear algebra operations, CGAL [2016] for mesh data structures, the Intel Threading Building Blocks library for parallel computing and the CUDA library for GPU parallel implementation of the Sinkhorn iterations. Unless otherwise indicated, we used the rst M = 100 eigenvectors of the Laplace operator to create di usion coordinates. For visual evaluation of the maps we use one color per vertex and linear color interpolation per triangle. Each vertex color is set based on a smooth 3D function, sometimes modulated by a 3D axis-aligned checkerboard for enhanced clarity.

In this case we further re ne the reference mesh via mutual tessellation with a strip along the grid edges to improve color sharpness (called strip coloring hereafter).

Initialization. When the input models are properly rigidly aligned in space (up to a translation), minimizing only the Wasserstein-2 transport distance is su cient to initialize η for the rst level of the di usion hierarchy. When the input models are not aligned or very di erent, we initialize our optimization with η = 0 and adding two to four user-de ned point constraints is typically su cient to obtain meaningful maps in these cases.

Sensitivity to parameters. While our exposition mentioned a number of parameters that the implementation depends on, their actual values only a ect e ciency. For instance, our use of di usion times being powers of two may look arbitrary: one could pick bigger increments between di usion times which would lead to a smaller number of di usion levels (as the number of levels is automatically found when one cannot add more samples), or smaller increments which would increase the number of levels in di usion stage. However, the rst case may trigger more alternating minimization iterations at the beginning of a new level as too sudden of a change in di usion may induce a large change of the map; the second case may reduce the number of iterations, but results in more iterations overall as the number of levels is larger. Similarly, di erent values of ε and τ could be used in the di usion stage as well. However, decreasing ε or increasing τ has an exponential impact on the cardinality of the neighborhoods. Our choice of 0.5 and 0.9 respectively has proven robust in all our examples as it roughly enforces the presence of 3-ring neighborhoods for each sample at each level. The parameter M as shown in Fig. 6, a ects the quality of sampling in the di usion stages. Picking only the rst 100 eigenvectors as we do by default may not be appropriate for shapes with very thin features: it could numerically force an excessive amount of map "ironing" otherwise in the Euclidean stages of the hierarchy (as the features appear late), which an increase of M to 300 would prevent quite nicely.

Timings. Computing the eigenvectors using Spectra as a preprocessing step takes up to 2.5 seconds on a mesh with 5k vertices and around 9 seconds on mesh with 100k vertices for 100 eigenvectors. De ning k to be the average neighborhood cardinality and n the total number of samples, the computation cost of each Sinkhorn iteration (including building/storing H) is O(ck 2 n), where c represents the transport plan "spread" (average number of points to which a point is mapped) that depends on γ . Furthermore, we parallelize our computations; consequently, our approach requires 2 to 5 mins for small (20kV) meshes, and up to 15 mins for 100kV. Note that without hierarchy and local variable instantiation, each Sinkhorn iteration would have a time and space complexity of O(n 2 ), with an additional O(ckn 2 ) time for building H. Despite its overall linear complexity, our approach is compute intensive due to three nested loops: the traversal of the hierarchy, the convergence of the centers of mass and the Sinkhorn iterations.

Inter-surface mapping examples. As our approach can nd a map between arbitrary surfaces, we mostly focus on demonstrating our results on non-isometric pairs of shapes. However, the reader is invited to parse through a number of examples of near-isometric maps in Supplemental Material, where we show that our technique compares favorably to previous approaches on the simpler intersurface mapping case.

Fig. 1 shows a baby mapped to a variety of non-isometric shapes. In this example, the user only provided two point-to-point correspondences (for the tip of the left hand/paw and the tip of the right big toe). The rest of the map, displayed through corresponding vertex colors, was found automatically through our coarse-to-ne optimization via alternating minimizations. A mass relaxation of 1.25 was allowed to avoid extreme stretching of the map.

Fig. 12 illustrates the mapping between two human hands. The hands di er in types (left vs. right) and are non-isometric. In addition, one hand extends to the wrist. Only two user-speci ed points are used as constraints: one to guide thumb to thumb and the other on the palm to disambiguate between front and back; the minimization algorithm then automatically nds the expected mapping. We further illustrate the mapping through cross-remeshing, i.e., by relocating the mesh vertices of the rst mesh onto the center of the images (via the optimal mass transport plan) of each vertex neighborhood on the other model (and vice versa).

Fig. 13 depicts a mapping between two obviously non-isometric dinosaur models, with no area relaxation. Even if the shapes di er greatly (in particular, on the tail, limbs, and head), we can automatically nd a smooth map in between them without requiring any user-prescribed correspondences.

Fig. 14 depicts a mapping between two highly non-isometric models, a man and a gorilla, with a maximum area relaxation factor set to 2 (i.e., area distortion is allowed to vary between 0.5 to 2). Four constraints were used (two hands and two feet). The resulting map is natural despite large di erences in limb sizes. Note that if we further increase area relaxation to improve the map, our algorithm starts experiencing slow convergence. Area preservation vs. conformality. Fig. 15 demonstrates the behavior of our approach when the area relaxation factors vary. We compute a map between a canonical torus and a swollen torus using three di erent relaxation factors (from left to right : 1, 2 and 4). As expected the part mapping to the swollen region gets more mass and vice-versa. When using no relaxation, the area is exactly preserved, hence triangles from the thin part of the torus get stretched due to tangential transport. Conversely, triangles are compressed in the opposite direction on the swollen part of the torus. When area distortion is allowed, triangles remain mostly undeformed and conformality is much improved. Note that the amount of relaxation at each location is found automatically by our algorithm based on m i and n j . The user only needs to input the global mass relaxation limit (referred as D in Sec. 3.4). Also note that although the mass distribution is locally non smooth, the resulting map is smooth nonetheless since our transport stencils are a few rings wide.

Robustness. Fig. 16 depicts a mapping between two poses of a kid. While on one pose the arm merges with his chest, the resulting map is only slightly distorted. Through mapping a sphere onto three genus-1 surfaces, Fig. 17 demonstrates the robustness of our approach to topological noise and depicts the mappings when the topology di ers. We modeled (1) a sphere with a small handle, (2) a sphere with a larger handle and (3) a canonical torus to test our approach. As expected the resulting maps are discontinuous but discontinuities are only present locally. No mass relaxation was used in these examples.

Evaluation. We evaluate the performance of our approach on three standard benchmark data sets: models from TOSCA [START_REF] Bronstein | Numerical Geometry of Non-Rigid Shapes[END_REF], SCAPE [START_REF] Anguelov | The Correlated Correspondence Algorithm for Unsupervised Registration of Nonrigid Surfaces[END_REF]] and SHREC'07 [START_REF] Giorgi | Shape retrieval contest 2007: Watertight models track[END_REF]]. Fig. 18 plots the percentage of correct correspondences against the geodesic error, averaged over all models. The units of the plots are identical to the ones used in BIM [START_REF] Kim | Blended intrinsic maps[END_REF]] to facilitate comparison. On near-isometric models from TOSCA and SCAPE our approach recovers the correct isometry in all cases and performs comparably to previous work. On SHREC models, BIM outperforms our approach (although not signi cantly), but our results further improve as we allow for more area relaxation: obviously, an area preserving map necessarily induces high distortion leading to bad semantic correspondences on such data. On these models we only constrained two (usually diagonally opposite) feature points. We, however, note that the curves typically shown for the SHREC models should be taken with a grain of salt: they solely measure mapping quality based on semantic feature points manually selected by a human. This is not telling much about the map: since a map is only useful in further processing task if it is smooth and valid everywhere, measuring performance based on a nite number of semantic points is not very informative, and provided here only for comparison purposes. In fact, we show in Fig. 19 that even for isometric shapes, BIM does not provide a useful map as obvious aws on the ears and folds on shoulders would prevent its use in practical applications. Comparatively our maps (not driven by semantics but by local geometric criteria) are always smooth and non-degenerate. Comparisons. We ran experiments to better compare our approach to BIM. It relies heavily on its initial establishment of pointto-point feature correspondences: if these correspondences fail (which happens increasingly often when shapes become non-isometric and with a wide range of feature sizes), the resulting map is nonsensical. When BIM nds good initial correspondences, their use of blended conformal maps results in visually nice maps. However, closer inspection shows severe limitations: if one relocates the vertices of the source mesh onto the target mesh via their output map, the resulting mapped mesh exhibits several artifacts as shown in Fig. 19. These artifacts often make BIM maps, already limited to genus-0 surfaces, unable to help with further geometry processing operations. Our result on the same example, found without any user-speci ed constraints, de nes a smooth and valid map everywhere. On most of the other examples we show in this paper, BIM (or the more recent non-rigid registration method of [START_REF] Chen | Robust Nonrigid Registration by Convex Optimization[END_REF], which can only handle near-isometric models) fails to provide a usable map. More generally, even if nding sparse correspondences can be achieved by some existing methods, extending these correspondences to a continuous bijection while ensuring low distortion is no small feat (see an e ort in this direction for functional correspondences in [START_REF] Rodolà | Point-wise Map Recovery and Re nement from Functional Correspondence[END_REF]). Instead, our approach strives to compute a dense and smooth map through a coarse to ne variational approach.

The closest approach to ours is arguably the Gromov-Wasserstein (GW) approach recently proposed by [START_REF] Solomon | Entropic Metric Alignment for Correspondence Problems[END_REF]: they also derive a geometry-driven map through optimization as we do. Because they consider both long range and short range geodesic distances in their formulation, one could be led to believe that their method is less myopic than ours which "only" minimizes the Dirichlet energy of the map through variance minimization of increasingly small neighborhoods-albeit, in a coarse to ne fashion. However, it should be pointed out that the GW approach only nds correspondences between a small number of points (typically a few thousands in their experiments): our attempt on a denser mesh (e.g., Fig. 19 with 70k vertices) as needed for a minimally-dense map in geometry processing was already intractable given their O(n 2 ) memory requirements and O(n 3 ) computational complexity. In sharp contrast, both our space and time complexity are near linear in the number of samples; we were thus able to generate examples in this paper using meshes with 20k to 200k vertices. Another drawback of their methodology is that one-to-many maps can be generated as output Fig. 19. Comparison with BIM. The dog model is mapped to another isometric pose with no user-specified constraints. Grey background: results obtained with blended intrinsic maps [START_REF] Kim | Blended intrinsic maps[END_REF]; even the map between two isometric models is, in fact, poor-particularly on the ears and the shoulders. We also depict the meshes relocated via their respective maps. (See Supplemental Material for meshes.) Fig. 20. One-to-many vs. one-to-one maps. Using the example of [START_REF] Solomon | Entropic Metric Alignment for Correspondence Problems[END_REF]] mapping a graph (le ) to a surface, their approach finds correspondences to disconnected regions (center, colored disks to colored regions), while our algorithm (right) enforces locality.

(see Fig. 20), which does not meet the typical objective in geometry processing to generate sharp point-to-point maps.

Limitations. The numerics of our approach works best for uniform isotropic meshes since σ used for weighting function W is set proportionally to the vertex density to ensure that the neighborhood cardinality remains nearly constant. As a consequence, we perform recursive longest edge bisection as a preprocessing stage to obtain a more isotropic vertex distribution so to ensure robustness. Additionally, we are also limited in the mass relaxation we can allow: relaxing the area preservation constraints through unbalanced mass transport increases the number of variables, which may slow down convergence for ratios exceeding a factor 5. In addition, if the mass is relaxed to the extent that a neighborhood maps to a single point, our algorithm will no longer be able to favor conformality locally. An obvious solution is to increase the neighborhood (by adjusting σ ), but this comes again at the cost of increased computation times.

Finally, while our hierarchy with power-of-two di usion time scales have worked reliably on all the examples we tried, a better understanding of the way di usion geometry behaves over the surface may allow for a better choice of the hierarchy, which could improve our runtimes signi cantly.

CONCLUSION

Recent progress in improving the e ciency of optimal transport solvers provides an opportunity to derive inter-surface mapping through direct energy minimization. Our paper proposes a rst variational model of the sort, where transport plans are optimized to be locally sharp so as to yield homeomorphisms. We showed that our approach nds a least-stretched map with exact area preservation or bounded area distortion, and demonstrated its robustness and performance on a range of challenging pairs of models, with di erent and/or non-trivial topology.

One weakness of our approach is that convergence may be slower on models with a wide range of feature sizes. One direction for future work is to use Newton iterations within the nested loops to accelerate convergence of the centroids. We also would like to explore the enduring challenge of partial matching which currently remains a largely unsolved problem. Further relaxing area constraints seems to be a promising tool in this perspective. Finally, exploring other geometric priors or bounds is an interesting research direction.

We wish to thank Vladimir Kim, Justin Solomon and Qifeng Chen for providing us with their respective source codes, and the reviewers for their constructive comments. MD gratefully acknowledges the Inria International Chair program and all the members of the TITANE team for their support.

A KL PROJECTIONS

We derive the KL projection of a matrix ρ on Π 1 , assuming the of ρ are positive. The gradient of the KL divergence to ρ at π is given by log π ρ. By the optimality conditions, the projection π * can be written as: log π * i j /ρ i j = λ + λ i + ξ i j ,

where λ comes from the total mass constraint, λ i from the row mass constraints, and ξ i j from the nonnegativity condition. Parameter λ i is zero when the mass bound is not attained, non negative if the lower mass bound is attained, and non positive if the upper mass bound is attained. Parameter ξ i j is zero when π * i j is positive, and non negative otherwise. Since π * i j is positive by construction, all the ξ i j are null. Let now ϕ i =e λ e λ i and ϕ = e λ . Using these variables, the optimality conditions read:

           ϕ i = ϕ if j ρ i j ϕ i ∈ (α i m i , β i m i ) ϕ i ≤ ϕ if j ρ i j ϕ i = β i m i ϕ i ≥ ϕ if j ρ i j ϕ i = α i m i
So ϕ i is a function of ϕ: ϕ i = f i (ϕ), where f i is the closest projection on the interval [α i m i / j ρ i j , β i m i / j ρ i j ], which is a nondecreasing piecewise linear function with two nodes. In order to determine the correct value for ϕ, it su ces to solve the equation i j ρ i j f i (ϕ) = 1, so that the global mass constraint is satis ed. Taking into account the change of variables due to the area factors m and n yields the algorithm described in Section 3.4.

ηFig. 3 .

 3 Fig.3. Convergence of alternating minimizations. This plot indicates the number of iterations to reach convergence for the curve-to-curve map in Fig.2, against the weighting function's width (indicated as the number of neighbors x i such that W x (x i ) > 0.5) used to measure variance.

Fig. 4 .

 4 Fig.4. Choice of weighting functions. From le to right: K t (x,x i ), K t (x,x i )/ K t (x,x )K t (x i ,x i ) and 2K t (x,x i )/(K t (x,x )+K t (x i ,x i )) for x indicated as the purple vertex. Its neighborhood, i.e., vertices x i for which W x (x i ) > ε, is shown with green spheres proportional to function values.

Fig. 5 .

 5 Fig. 5. Mapping between isometric poses. Result of our algorithm between three di erent dog poses; the reference pose is on the le , and its strip-coloring is displayed on the other two shapes based on resulting maps.

  and perform the inner products a b in Steps 4-6 of Alg. 2 (involving columns or rows of H) by rst computing the maximum element Z = max i {a[i] + b[i]}, and evaluating a b Z + log j exp(a[j]+b[j] -Z ).

Fig. 7 .Fig. 8 .

 78 Fig. 7. Local variable instantiation. For a point x , we instantiate variables only in the neighborhood of the mapping of its neighborhood. Localizing samples it can be mapped to dramatically reduces the search space.

Fig. 9 .

 9 Fig.9. Flips and twists. An ideal map (le ); a mirrored map creating flip (middle); and 180 degree rotation of the map creating a twisted map. In practice, we found that twists rarely occur in our maps and are highly unstable as compared to flips.

Fig. 10 .

 10 Fig. 10. Fixing flips and twists. Starting from a suboptimal map having a discontinuity, ironing proceeds through front propagation to remove artifacts. For each strip, blue indicates visited parts while non-visited parts are in pink. In case of multiple fronts, each front is considered independently.

Fig. 12 .

 12 Fig. 12. Mapping hands. Top le : a source right hand is mapped to a target le hand based on two user-defined point-to-point correspondences. Top right: the map is depicted through a color strip. Bo om le : target mapped onto source mesh. Bo om right: source mapped onto target mesh.

Fig. 13 .

 13 Fig.13. Dinos. While the limbs, tails and heads of these two dinos greatly di er in proportion, our algorithm automatically finds a smooth and intuitive area-preserving map between the two, without any user interaction.

Fig. 14 .

 14 Fig. 14. Mapping man onto gorilla. An area distortion (up to a factor 2) is used, along with 4 constraints (one at the extremity of each limb) to map these two very di erent shapes.

Fig. 15 .

 15 Fig. 15. Area preservation vs. conformality. When mass is relaxed, the transported mass j π i j changes while the original mass m i remains identical. The slackness C helps to gradually adjust the amount of local massrelaxation, globally bounded using the parameter D. Top: colors indicate local relaxed transported mass m i /m i (blue=small, red=large). Top center: colors indicate area change per triangle (blue=shrunk, red=dilated). Bo om center: colors indicate quasi-conformal distortion per triangle (green=conformal, magenta=sheared). Larger relaxation factors allow be er conformality. Bottom: closeups on the triangle meshes to highlight improved conformality.

Fig. 16 .

 16 Fig. 16. Mapping with local contact. While the baby's arm merges with his chest on the second pose, we find a good map with only four userselected constraints (one for each limb extremity).

Fig. 21

 21 Fig. 21 depicts a sample of the various SHREC models we tested and the mappings we obtained, with and without area relaxation.

Fig. 17 .

 17 Fig.17. Mapping for di erent genus. Middle le : a sphere with a protruding handle is mapped to a sphere (le ). The map is discontinuous, but only locally. Middle right: with a bigger handle. Right: a regular torus is mapped to the sphere; the discontinuity is no longer local, but the map is continuous over half of the sphere.

  Fig. Evaluation on benchmark datasets. We plot the percentage of correct correspondences against geodesic error, averaged over all model pairs, for various datasets. Our curves for SHREC'07 show three di erent bounds on allowed mass distortion.

Fig. 21 .

 21 Fig.21. Maps from SHREC'07. We depict mappings between a variety of SHREC models from di erent categories.

ACM Transactions on Graphics, Vol. 36, No. 4, Article 39. Publication date: July 2017.

ACKNOWLEDGMENTS

The hands, skull and dino models are courtesy of the AIM@Shape repository; the dog, man and gorilla models are from TOSCA; the baby models from SHREC'16; Fig. 20 features a mesh from [Solomon et al. 2016], courtesy of the authors; all other models are from the SHREC'07 dataset.

This work has been funded by the European Research Council (ERC grants 257474: "IRON", 339025: "GUDHI" and 340884: "ACROSS"), NSF CCF-1655306 and the German Research Foundation (DFG, Gottfried-Wilhelm-Leibniz Programm).