I Isotopic Approximation within a Tolerance Volume

This thesis is divided into two independent parts.

In the first part, we introduce a method that, given an input tolerance volume, generates a surface triangle mesh guaranteed to be within the tolerance, intersection free and topologically correct. A pliant meshing algorithm is used to capture the topology and discover the anisotropy in the input tolerance volume in order to generate a concise output. We first refine a 3D Delaunay triangulation over the tolerance volume while maintaining a piecewise-linear function on this triangulation, until an isosurface of this function matches the topology sought after. We then embed the isosurface into the 3D triangulation via mutual tessellation, and simplify it while preserving the topology. Our approach extends to surfaces with boundaries and to non-manifold surfaces. We demonstrate the versatility and efficacy of our approach on a variety of data sets and tolerance volumes.

In the second part we introduce a new approach for creating a homeomorphic map between two discrete surfaces. While most previous approaches compose maps over intermediate domains which result in suboptimal inter-surface mapping, we directly optimize a map by computing a variance-minimizing mass transport plan between two surfaces. This non-linear problem, which amounts to minimizing the Dirichlet energy of both the map and its inverse, is solved using two alternating convex optimization problems in a coarse-to-fine fashion. Computational efficiency is further improved through the use of Sinkhorn iterations (modified to handle minimal regularization and unbalanced transport plans) and diffusion distances. The resulting inter-surface mapping algorithm applies to arbitrary shapes robustly and efficiently, with little to no user interaction. v

First of all, I am most thankful to my advisor, Pierre Alliez, for his thorough involvement and patience. His incredible zeal for research and a near supernatural ability of fulfilling multiple deadlines at the same time has been a source of motivation. He pushed me towards ideas which seemed unreasonable at times, but later formed the backbone of the algorithms presented in this thesis.

I am grateful to David Cohen-Steiner for being my co-advisor (though unofficial). Every discussion with him brought me some new insights. His critical thinking and reasoning always amused as well as motivated me to perform better. I owe many thanks to Mathieu Desbrun for helping me out with every minute issue, always excited and being available for discussions, and most importantly, bearing up with my horrible explanation and coding skills. I am also thankful to Leif Kobbelt for his guidance, broad discussions and providing me with a postdoctoral fellowship opportuinity.

I must thank Florent Lafarge, David Bommes, Bishesh Khanal, Yannick Verdie, Sonali, Liuyun Duan, Deepesh Agarwal and Sven Oesau for always being there and helping me out whenever I needed them.

I thank the members of the jury: Marc Alexa, Paolo Cignoni and Marco Cuturi for spending their precious time on my dissertation.

Finally, I attribute some of the shortcomings of this thesis to my family, friends and la belle France who gave me enough reasons to escape work every now and then.

Résumé

Cette thèse comprend deux parties indépendantes.

Dans la première partie nous contribuons une nouvelle méthode qui, étant donnée un volume de tolérance, génère un maillage triangulaire surfacique garanti d'être dans le volume de tolérance, sans auto-intersection et topologiquement correct. Un algorithme flexible est conçu pour capturer la topologie et découvrir l'anisotropie dans le volume de tolérance dans le but de générer un maillage de faible complexité.

Dans la seconde partie nous contribuons une nouvelle approche pour calculer une fonction de correspondance entre deux surfaces. Tandis que la plupart des approches précédentes procède par composition de correspondance avec un domaine simple planaire, nous calculons une fonction de correspondance en optimisant directement une fonction de sorte à minimiser la variance d'un plan de transport entre les surfaces.

Representative Publications

MANDAD, M., Cohen-Steiner, D., and Alliez, P., 2015. Isotopic Approximation within a Tolerance Volume. Proceedings of ACM SIGGRAPH, ACM transactions on Graphics, 34(4):64:1-64:12. ISSN 0730-0301. Cohen-Steiner, D., Desbrun, M., Alliez, P., and Kobbelt L., 2016. Low Distortion Inter-surface Mapping via Variance-Minimizing Mass Transport. (Submitted).

MANDAD, M.,

Chapter 1 Introduction

Faithful approximation of complex shapes with simplicial meshes is a multifaceted problem, involving geometry, topology and their discretization. This problem has received considerable interest due to its wide range of applications and the ever-increasing accessibility of geometric sensors. Increased availability of scanned geometric models, however, does not mean improved quality: while many practitioners have access to high-end acquisition systems, a recent trend is to replace these expensive systems with a combination of consumer-level acquisition devices.

Measurement data generated by, and merged from, these heterogeneous devices are reputedly unfit for direct processing.

Similarly, the growing variety of geometry processing tools often increases the net amount of defects in data: conversion to and from various geometry representations often degrades the input, and rare are the algorithms that have stronger guarantees on their output than they have requirements on their input. As we deal with ever finer discretizations to capture intricate geometric features, this issue of offering strict geometric guarantees to be robust to the occurrence of artifacts is becoming more prevalent.

Geometric guarantees usually refer to upper bounds on the approximation error and to the absence of self-intersections. Topological guarantees refer to homotopy, homeomorphism or isotopy. Surface meshes with such guarantees are required for artifact-free rendering, computational engineering, reverse engineering, manufacturing and 3D printing. While geometric simplification can reduce the number of primitives, topological simplification can repair holes and degeneracies in existing discretizations. Combined, the two may also be used for reconstructing clean shapes from raw geometric data such as point sets or polygon soups with strict guarantees.

Background

Before delving into details, we discuss the existing theory required to understand the topic and our problem statement. As mentioned before, the problem of faithful approximation deals with two sets of theoretical guarantees: Topological and Geometric. We will detail them one by one.

Topology

In the context of this thesis, topological guarantees reflect whether the final (reconstructed) output has same topology as that of the input or not. Intuitively, topology refers to study of the properties that are preserved through deformations, twistings, and stretchings of objects.

Tearing and stitching, however, are not allowed. Based on these properties and how strongly they are preserved, topological guarantees can be classified as homotopy, homeomorphism and isotopy.

Homotopy. Mathematically, two given spaces 𝑋 and 𝑌 , are homotopy equivalent or of same homotopy type if there exist continuous maps 𝑓 ∶ 𝑋 → 𝑌 and 𝑔 ∶ 𝑌 → 𝑋 such that 𝑓 ∘ 𝑔 is the identity map on 𝑌 , and, 𝑔 ∘ 𝑓 is the identity map on 𝑋. Having the same homotopy type enforces the existence of isomorphisms between respective homology and homotopy groups.

In other words, as explained by [START_REF] Chazal | Smooth manifold reconstruction from noisy and nonuniform approximation with guarantees[END_REF], having the same homotopy type enforces a one-to-one correspondence between cycles, holes, connected components, tunnels, cavities and/or higher dimensional topological features of the two sets 𝑋 and 𝑌 such that one set can be continuously deformed into the other set. These correspondences along with the relation between these features is defined by the maps 𝑓 and 𝑔. For example, a solid disk or solid ball is homotopy equivalent to a point. A torus and a trefoil knot are homeomorphic. Note that there is no need for the existence of a continuous deformation for the two spaces to be homeomorphic. Furthermore, they are not isotopic as the knots will have to pass through each other to form an unknot and it will not be possible to maintain a homeomorphism throughout the deformation.

Homeomorphism. Being a stronger condition, homeomorphism is generally considered as the notion of equality in topology. Two spaces with a homeomorphism between them are called homeomorphic. For example, a classic example in topology suggests that an unknot (torus) and a knot (trefoil knot) are homeomorphic (Figure 1.1). This is because one of the geometric objects can be stretched and bent continuously to form the other object. Mathematically, given two topological spaces 𝑋 and 𝑌 , a function 𝑓 ∶ 𝑋 → 𝑌 between these spaces is a homeomorphism if 1. 𝑓 is a bijection and is continuous, and 2. the inverse of 𝑓 exists and is continuous.

Isotopy. Isotopy is a finer relation than homeomorphism and the strongest condition for topological similarity between two objects. In our context, isotopy means that there exists a smooth deformation that maps one shape to another while maintaining a homeomorphism between the two. Note that maintaining a homeomorphism during deformation is not a necessary criteria for two objects to be homeomorphic. As a result, a coffee mug and donut are isotopic (Figure 1.2), but they are not isotopic to a trefoil knot (Figure 1.1).

Figure 1.2: Isotopy. A coffee mug can be continuously deformed into a donut (torus). Intuitively, two surfaces are isotopic if they can be continuously deformed into the other without introducing self-intersections. Note that all the objects including the intermediate deformations are homeomorphic.

Geometry

Geometric guarantees mainly deal with upper bounds on the approximation error of the reconstructed output and with the absence of self-intersections in it. In this thesis we look at these guarantees only in the context of triangle meshes.

Approximation

It is usually desirable to significantly reduce the complexity (mainly redundancies) of the discrete surface representation while preserving, as much as possible, the geometry (features) of the object. This results in a tradeoff between the mesh complexity and mesh quality as illustrated in Figure 1.3. While the mesh complexity can be easily defined in terms of the vertex count, the mesh quality has many aspects and is, in general, computed in terms of Hausdorff distance.

Along with Hausdorff distance, other fairness criteria for example based on normal deviation, triangle shape, dihedral angle, valence etc. can also be used.

A good approximation can either be achieved by reconstructing an anisotropic mesh directly or via mesh simplification. Such an approximation can drastically reduce the size of an oversampled 3D data (for example, scan data or an isotropic mesh) and the output mesh can be efficiently used for geometry processing, multi-resolution hierarchies, level-of-detail (LOD) rendering and adaptation to hardware capabilities. Hausdorff Distance. The Hausdorff distance is a measure of (dis-)similarity between two sets and is widely used in various applications. It is measured as the maximum distance of a set to the nearest point in the other set. More formally, given two sets 𝐴 and 𝐵, the (one-sided) Hausdorff distance 𝑑 𝐻 (𝐴 → 𝐵) from a set 𝐴 to set 𝐵 is a maximin function (Figure 1.4), defined as:

𝑑 𝐻 (𝐴 → 𝐵) = sup ∀𝑎∈𝐴 inf ∀𝑏∈𝐵 ‖𝑎 -𝑏‖.
(1.1)

B A d H (B → A) d H (A → B)
Figure 1.4: Hausdorff Distance. The one-sided Hausdorff distance is defined as the maximum distance of a set to the nearest point in the other set. Note that this distance is not a metric because it is not always symmetric, i.e. 𝑑 𝐻 (𝐴 → 𝐵) ≠ 𝑑 𝐻 (𝐵 → 𝐴) and 𝑑 𝐻 (𝐴 → 𝐵) = 0 does not imply that both sets are identical.

Such distance between a given pair of surfaces can be easily computed via densely sampling them and finding the nearest neighbours [START_REF] Cignoni | Metro: Measuring error on simplified surfaces[END_REF].

On a general note, asymmetry is a property of maximin functions, while minimin functions are symmetric. As a result, the two-sided Hausdorff distance is symmetric (assuming that the two sets do not have any repetitive elements), which can be calculated as the maximum of the two one-sided distances: max{𝑑 𝐻 (𝐴 → 𝐵), 𝑑 𝐻 (𝐵 → 𝐴)}.

(1.2) Mesh Simplification. Given a mesh 𝑀 = (𝑉 , 𝐹 ), where 𝑉 and 𝐹 denote the vertices and faces respectively, a mesh approximation problem can be defined as finding a new mesh 𝑀 ′ = (𝑉 ′ , 𝐹 ′ ) such that either 1. |𝑉 ′ | = 𝑛 < |𝑉 | and ‖𝑀 -𝑀 ′ ‖ is minimal, or 2. ‖𝑀 -𝑀 ′ ‖ < 𝛿 and |𝑉 ′ | is minimal.

B A A

Figure 1.5: Edge collapse. Given an edge 𝐴𝐵 (middle), a half-edge collapse (left) locates the target vertex at one of the vertices of the edge. A general edge collapse (right), relocates the target vertex based on the user defined criteria. Note that an edge collapse always has a unique triangulation (combinatorially, depending on its one-ring neighbours), although the final vertex location may change.

A Figure 1.6: Vertex removal. From left to right: Given a vertex 𝐴, it is eliminated leading to a hole in the triangulation which is then filled via re-triangulation. Note that unlike an edge collapse operator (which has vertex location as the degree of freedom), a vertex removal operator has degrees of freedom in terms of triangulation.

In general, incremental decimation algorithm is used to implement mesh decimation. Using a greedy optimization approach, a decimation operator is applied on a mesh region until no further reduction is possible. Half-edge collapse is commonly used as a decimation operator due to its simplicity, but, other methods like general edge collapse (Figure 1.5) and vertex removal with re-triangulation also exist (Figure 1.6). Note that both of them are Euler operations as they do not change the genus of the mesh.

Self-Intersection

A mesh with self-intersections does not necessarily mean that its topology (homeomorphism) has been altered. In fact, most of the time, purely based on connectivity, the mesh does not lose its homeomorphic property at all. However, the mesh seen as a function loses the bijective property and behaves as a one-to-many relation. Such a mesh cannot be used for practical applications, for example, 3D printing, as at self-intersections, the Euclidean distance is zero for points with non-zero geodesic distance. Note that an isotopic mesh is always free from (extra, if any) self-intersections. This, in fact, bridges this geometric property back to topological properties mentioned above. Flips. An edge collapse may not always result in a valid triangulation. For a triangulation to be valid, the target location of vertex must lie in the visibility kernel (green shaded region) of the one-ring. Left: a valid triangulation in 2D. Right: A self-intersection can also occur in the form of flips due to invalid triangulation. The same can also be extended to higher dimensions.

Related Work

A vast array of methodologies have been proposed for surface reconstruction and shape approximation over the years, ranging from Delaunay triangulation to implicit surfaces and decimation to optimization through clustering and refinement. Fewer, however, provide error bounds. In addition, they only apply to specific types of input geometry, and often fail to satisfy geometric and/or topological guarantees as we now review.

Reconstruction with Topological Guarantees

A great variety of methods have been proposed for surface reconstruction. Among them, while only a limited come with the theoretical guarantee of homeomorphism, only a few recent articles tackle the more difficult problem of ensuring isotopy.

The credit for the first reconstruction algorithm to provide theoretical guarantees goes to Amenta for Crust [Amenta et al., 1998;Amenta and Bern, 1998]. This algorithm exploits the structures of the Voronoi diagram of the input point set to reconstruct the surface, and converges to the original surface as the sampling density increases. The output mesh of this combinatorial algorithm interpolates the input points rather than approximating them to form a surface and hence, can handle variable sampling density. A 2-dimensional curve reconstruction example is illustrated in Figure 1.9. The resulting mesh is homeomorphic and geometrically close to the original (smooth) surface when the prescribed sampling condition is satisfied. Furthermore, they proposed the definition of poles, formed by a subset of Voronoi vertices of the input point set which can represent an approximated medial axis when input data set is sufficiently dense.

Homeomorphism. [START_REF] Amenta | A simple algorithm for homeomorphic surface reconstruction[END_REF] introduced the Cocone algorithm which improved the previous algorithm [Amenta and Bern, 1998] both in theory and in practice.

Following Crust, [START_REF] Boissonnat | Smooth surface reconstruction via natural neighbour interpolation of distance functions[END_REF] designed another algorithm based on natural neighbors and reconstruction using signed distance functions. The algorithm works in any dimension. Another Delaunay based approach [START_REF] Cohen-Steiner | A greedy delaunay-based surface reconstruction algorithm[END_REF] provided a greedy algorithm for surface reconstruction form unorganized point sets. At each iteration, the algorithm ensured that the reconstructed surface is an orientable manifold, possibly with a boundary.

Note that the fundamental property used for proofs is that a closed surface is homeomorphic (and in fact isotopic) to the restricted Delaunay triangulation of a sufficiently dense point sample.

However, the same does not hold for surfaces with boundaries. Although all these algorithms work nicely on dense data sets, they have difficulty when the data is undersampled resulting into a reconstructed surface that may not even be a manifold. In fact, in most such cases, these algorithms produce holes or other artifacts in the regions of undersampling. The reason for these artifacts can be attributed to relying totally on Delaunay triangulation or Voronoi diagram which, in practice, fails to form a topologically correct surface due to noise, undersampling or sharp surface features.

Watertight Reconstruction. To deal with undersampling, Amenta et al. [2001a,b] provided the Power Crust algorithm based on approximating the medial axis transform of the object and then using an inverse transform to produce the surface. The algorithm uses a weighted Voronoi diagram of the poles, referred to as a power diagram, instead of the Delaunay triangulation (as in Amenta and Bern [1998]) of the point set. The algorithm is guaranteed to produce a watertight surface (along with producing an approximate medial axis) and does not depend in any way on the quality (sampling density) of the input points. However, it introduces many extra points in the output and does not produce a triangulated surface. [START_REF] Edelsbrunner | Surface Reconstruction by Wrapping Finite Sets in Space[END_REF] proposed a Morse theoretic approach to determine the surfaces which also produces watertight surfaces. It also provides an algorithm to reconstruct the surface based on Delaunay complexes and extracting a subcomplex through repeated collapsing. All ideas and the results generalize to arbitrary dimensions. However, the algorithm may not recover small features of the surface even if it is densely sampled. [START_REF] Adamy | Surface reconstruction using umbrella filters[END_REF] introduced an umbrella filter algorithm coupled with a linearprogramming based topological post-processing step designed for topologically correct watertight triangle mesh reconstruction. [START_REF] Dey | Tight cocone: A water-tight surface reconstructor[END_REF] extended [START_REF] Amenta | A simple algorithm for homeomorphic surface reconstruction[END_REF] to provide the Tight Cocone algorithm. The algorithm first computes a preliminary surface using Cocone which almost completes the reconstruction except in the vicinity of undersampled areas. A subsequent marking and peeling phase completes the reconstruction by filling holes and computes the output surface as the boundary of the union of Delaunay tetrahedra. The algorithm guarantees that the surface cannot have any holes. However, it relies on the locality of undersampling and may produce a non-manifold (though watertight) or empty output if this property is violated. Furthermore, it cannot reconstruct surfaces with boundaries.

Another class of approaches based on implicit surfaces have the advantage that the output surface always remains watertight and they work well on noisy dataset. [START_REF] Hoppe | Surface reconstruction from unorganized points[END_REF] pioneered an approach to explicitly construct the geometry of an implicit manifold. The algorithm computes a signed distance field to the underlying surface using local tangent plane construction, and extracts an isosurface from the distance field using the marching cubes algorithm. [START_REF] Mullen | Signing the unsigned: Robust surface reconstruction from raw pointsets[END_REF] proposed a method that improves the correctness of the local surface orientation in order to fill gaps in regions where data is missing. [START_REF] Schnabel | Completion and reconstruction with primitive shapes[END_REF] (Image taken from [START_REF] Alliez | Voronoi-based variational reconstruction of unoriented point sets[END_REF])

Isotopy. Given a nonsingular compact 2-manifold surface 𝑆 without boundary, [START_REF] Sakkalis | Ambient isotopic approximations for surface reconstruction and interval solids[END_REF] introduced algorithm for creating a family of approximating surfaces to 𝑆 while ensuring that each approximant (𝑆 ′ ) is ambient isotopic to the original. The algorithm proceeds by constraining the approximant to lie within a bounded offset of 𝑆. While most existing reconstruction methods provide only for piecewise-linear approximations, this method is also suitable for higher order approximations. The main result of [START_REF] Sakkalis | Ambient isotopic approximations for surface reconstruction and interval solids[END_REF] is that 𝑆 and 𝑆 ′ are isotropic if projection on 𝑆 defines a homeomorphism from 𝑆 ′ to 𝑆 (Figure 1.11).

S S 2δ

x

n x x -δn x x + δn x

S(δ)

Figure 1.11: Conditions for ambient isotopy. Assuming 𝑆 to be a compact nonsingular 2-manifold without boundary and a normal tubular neighborhood 𝑆(𝛿) for a given 𝛿 (refer to [START_REF] Sakkalis | Ambient isotopic approximations for surface reconstruction and interval solids[END_REF] for details). If (a) 𝑆 ′ ⊂ 𝑆(𝛿), and (b) ∀𝑥 ∈ 𝑆, the line segment (𝑥 -𝛿𝑛 𝑥 , 𝑥 + 𝛿𝑛 𝑥 ) intersects 𝑆 ′ precisely at one point, then, 𝑆 and 𝑆 ′ are ambient isotopic. [START_REF] Amenta | Computational topology: Ambient isotopic approximation of 2-manifolds[END_REF] showed that a specific piecewise linear approximation of a closed surface is isotopic, using the same condition as in [START_REF] Sakkalis | Ambient isotopic approximations for surface reconstruction and interval solids[END_REF] indirectly. Both these approaches involve not only the topology of the surfaces, but also the geometry in providing topological guarantees. Furthermore, since projection is involved, the conditions cannot be met when 𝑆 is not smooth. [START_REF] Chazal | A condition for isotopic approximation[END_REF] proposed a purely topological condition which is sufficient to ensure isotopy but a practical algorithm is missing for this theoretical result.

Another work [START_REF] Chazal | Smooth manifold reconstruction from noisy and nonuniform approximation with guarantees[END_REF], provided reconstruction with topological guarantees using unions of balls centered at the data set. [START_REF] Attali | Vietoris-rips complexes also provide topologically correct reconstructions of sampled shapes[END_REF] showed that Rips complexes can also be used to provide topologically correct reconstruction of shapes which might be of some computational interest in higher dimensions.

Isotopic surface reconstruction from noisy point sets, provided the sampling is dense enough, has also been investigated well [START_REF] Dey | Curve and Surface Reconstruction: Algorithms with Mathematical Analysis[END_REF]. Another approach which can handle noisy data set by estimating normals from big Delaunay balls [START_REF] Dey | An adaptive mls surface for reconstruction with guarantees[END_REF] (Figure 1.12), extended the classical Moving Least Squares (MLS) [START_REF] Shepard | A two-dimensional interpolation function for irregularly-spaced data[END_REF] (out-of-core extension [START_REF] Fiorin | Out-of-core mls reconstruction[END_REF]) which is a popular method for functional approximation. A recent approach, [START_REF] Dey | Isotopic reconstruction of surfaces with boundaries[END_REF] showed that a restricted Delaunay triangulation of a sufficiently dense point sample of a surface is isotopic to a collar extension. Although this approach is able to deal with boundaries and even allows reconstruction of non-orientable surfaces, it cannot handle noisy data.

Isotopic approximation in the context of implicit curves and surfaces has also been studied.

Lately, [START_REF] Burr | Complete subdivision algorithms, ii: Isotopic meshing of singular algebraic curves[END_REF]; [START_REF] Lin | Adaptive isotopic approximation of nonsingular curves: The parametrizability and nonlocal isotopy approach[END_REF] introduced new algorithms by combining the parametrizability (as in [START_REF] Snyder | Interval analysis for computer graphics[END_REF]) and the nonlocal isotopy (as in [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF]). (Image taken from [START_REF] Dey | An adaptive mls surface for reconstruction with guarantees[END_REF])

Discussion. Most of the algorithms discussed above cannot handle noisy data and, in fact, none can handle outliers. Although some of them provide isotopic guarantee, they fail to provide an approximation with a tradeoff between complexity and quality. All these methods generate isotropic meshes, overly complex, which would require another algorithm for simplification with geometric guarantees.

Other Reconstruction/Mesh-Repair Algorithms

Apart from the algorithms discussed above, there are many reconstruction algorithms that do not provide theoretical guarantees over the topology but we highlight them to complete the state of the art in the context of this thesis.

Heterogeneous Data. [START_REF] Hornung | Robust reconstruction of watertight 3d models from non-uniformly sampled point clouds without normal information[END_REF] used unsigned distance function to avoid the topological noise artifacts caused by misalignment of 3D scans, which are common to most volumetric reconstruction techniques. The proposed volumetric algorithm estimates local surface confidence values within a dilated crust around the input samples and the surface which maximizes the global confidence is then extracted by computing the minimum cut of a weighted spatial graph structure. the method generates a watertight mesh even for noisy and highly irregular data containing large holes, without loosing fine details in densely sampled regions (Figure 1.13).

Lately, [START_REF] Argudo | Biharmonic fields and mesh completion[END_REF] proposed a volumetric method based on bi-harmonic fields. The algorithm operates locally, within an expanded bounding box of each hole, and therefore scales well with the number of holes in a single, complex model. (Image taken from [START_REF] Hornung | Robust reconstruction of watertight 3d models from non-uniformly sampled point clouds without normal information[END_REF]) [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] casted the surface reconstruction from oriented points as a spatial Poisson problem. The approach considers all the points at once, without resorting to heuristic spatial partitioning or blending, and hence, is highly resilient to data noise. Recently, Giraudot et al.

V crust V int V ext
[2013] extended the robust distance function [START_REF] Chazal | Geometric Inference for Measures based on Distance Functions[END_REF] and proposed a noise-adaptive robust distance function and a surface reconstruction algorithm. The resulting algorithm is resilient to variable noise, outliers and missing data as illustrated in Figure 1.14. The algorithm only assumes that the inferred shape is a smooth closed submanifold of known dimension.

Reconstruction and simplification algorithms based on optimal transportation [START_REF] Goes | An optimal transport approach to robust reconstruction and simplification of 2d shapes[END_REF][START_REF] Digne | Featurepreserving surface reconstruction and simplification from defect-laden point sets[END_REF] has also been studied as they are robust to both noise and outliers.

However, they lack topological guarantees and are too compute intensive to be practical.

Methods involving repairing using a volumetric approach [START_REF] Ju | Robust repair of polygonal models[END_REF] have also been targeted. Other methods involve identifying consistent sub-meshes in the original triangle data which are then merged together by snapping corresponding boundary segments or by stitching small patches into the remaining gaps and holes [START_REF] Bohn | Automatic CAD-model Repair[END_REF][START_REF] Guéziec | Converting sets of polygons to manifold surfaces by cutting and stitching[END_REF][START_REF] Liepa | Filling holes in meshes[END_REF][START_REF] Attene | A lightweight approach to repairing digitized polygon meshes[END_REF]. Other local inconsistencies, for example small handles or tunnels can also be repaired [START_REF] Guskov | Topological noise removal[END_REF]]. We refer the readers to [START_REF] Attene | Polygon mesh repairing: An application perspective[END_REF] for a comprehensive survey on polygon mesh repairing. [START_REF] Bischoff | Automatic restoration of polygon models[END_REF] used a combination of a volumetric geometry representation and the original triangle data in order to exploit the advantages of both (Figure 1.15). [START_REF] Shen | Interpolating and approximating implicit surfaces from polygon soup[END_REF] introduced a method for building interpolating or approximating implicit surfaces from polygonal data which can preserve sharp features while smoothing noise, generate envelopes around input data and produce mesh with low polygon count. Anisotropic Reconstruction. Instead of a two step process -reconstruction followed by simplification -direct anisotropic reconstruction has been targeted via particle-based meshing [START_REF] Bossen | A pliant method for anisotropic mesh generation[END_REF][START_REF] Zhong | Particlebased anisotropic surface meshing[END_REF], Riemannian Voronoi diagrams [START_REF] Leibon | Delaunay triangulations and voronoi diagrams for riemannian manifolds[END_REF], Centroidal Voronoi Tessellation (CVT) [START_REF] Du | Anisotropic centroidal voronoi tessellations and their applications[END_REF][START_REF] Sun | Obtuse triangle suppression in anisotropic meshes[END_REF][START_REF] Richter | Mahalanobis centroidal voronoi tessellations[END_REF], Delaunay Triangulation [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications. part i algorithms[END_REF] and parameterization [START_REF] Zhong | Anisotropic surface meshing with conformal embedding[END_REF]. Meshing in higher dimensional space has also been explored by embedding in a 6D space and efficiently computing CVT [START_REF] Lévy | Variational Anisotropic Surface Meshing with Voronoi Parallel Linear Enumeration[END_REF]. Based on Delaunay refinement [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF][START_REF] Boissonnat | Anisotropic Delaunay Mesh Generation[END_REF] made the star around each vertex to be consisting of the triangles that are exactly Delaunay for the given metric. These stars are later stitched together gradually using a refinement algorithm. As a totally different approach, [START_REF] Azernikov | Anisotropic meshing of implicit surfaces[END_REF] proposed a grid-based approach to mesh implicit surfaces according to their curvature tensor. The method uses a deformed grid and produces quad-meshes. Recently, [START_REF] Fu | Anisotropic simplicial meshing using local convex functions[END_REF] transformed the anisotropic meshing problem into a functional (convex) approximation problem which generalizes optimal Delaunay triangulation [START_REF] Long | Optimal delaunay triangulations[END_REF][START_REF] Chen | Optimal anisotropic meshes for minimizing interpolation errors in 𝑙 𝑝 -norm[END_REF].

Discussion. When dealing with imperfect or heterogeneous data, various methods involving repairing, conversion, or reconstruction are designed to generate clean meshes but they do not yield low-polygon-count approximations with bounded error. Moreover, the topological guarantees are missing, omitting a few which provide watertight meshes as output.

Anisotropic reconstruction algorithms are too slow and are not effective enough to provide low-polygon count as compared to other state of the art approaches, for example, mesh simplification algorithms. Furthermore, handling noisy or undersampled data is often out of scope.

Guaranteed Error

Hausdorff. In order to build approximation with bounded error, a large collection of simplification algorithms exist. Approximation with bounded error has been targeted through clustering [START_REF] Kalvin | Superfaces: Polygonal mesh simplification with bounded error[END_REF] and mesh decimation [START_REF] Bajaj | Error-bounded reduction of triangle meshes with multivariate data[END_REF][START_REF] Cohen | Simplification envelopes[END_REF][START_REF] Klein | Mesh reduction with error control[END_REF][START_REF] Guéziec | Surface simplification inside a tolerance volume[END_REF][START_REF] Ciampalini | Multiresolution decimation based on global error[END_REF][START_REF] Cignoni | Simplification of tetrahedral meshes with accurate error evaluation[END_REF][START_REF] Cohen | Successive mappings: An approach to polygonal mesh simplification with guaranteed error bounds[END_REF][START_REF] Botsch | GPU-based tolerance volumes for mesh processing[END_REF][START_REF] Ovreiu | Mesh simplification using a two-sided error minimization[END_REF]. A combination of both, clustering and mesh decimation has also been studied [START_REF] Zelinka | Permission grids: Practical, error-bounded simplification[END_REF], where the tolerance volume was discretized to accelerate the validity tests for the atomic simplification operations. Typically, the methods only check for topology preservation and the normal-flip test [START_REF] Ronfard | Full-range approximation of triangulated polyhedra[END_REF] which only prevents local self-intersection. The simplification algorithms use a distance metric to estimate the error, [START_REF] Garland | Surface simplification using quadric error metrics[END_REF]; [START_REF] Hoppe | New quadric metric for simplifying meshes with appearance attributes[END_REF] being the popular ones. We further refer the readers to [START_REF] Cignoni | A comparison of mesh simplification algorithms[END_REF] and [START_REF] Luebke | A developer's survey of polygonal simplification algorithms[END_REF] for comprehensive surveys on simplification algorithms.

Other Fairness Schemes. In general, the error metric considered is the one-sided Hausdorff distance to the input mesh, but the normal deviation has also been considered [START_REF] Borouchaki | Simplification of surface mesh using hausdorff envelope[END_REF]. Although the distance metric is very efficient in measuring geometric errors, it has difficulties in distinguishing important shape features such as high-curvature regions. Measuring geometric error for such features has been targeted via curvature tensor [START_REF] Coll | Adaptive simplification of huge sets of terrain grid data for geosciences applications[END_REF] and discrete curvature norm [START_REF] Kim | Surface simplification using a discrete curvature norm[END_REF]. [START_REF] Kobbelt | A general framework for mesh decimation[END_REF] generalized the error metric by combining the local surface properties -the approximation error (function value), the local distortion (1st order derivatives) and the local curvature (2nd order derivatives) -into one linear functional. Recently, [START_REF] Morigi | Multilevel mesh simplification[END_REF] presented new approach based on the evolution of surfaces under p-Laplacian operators which provides a natural geometric clustering.

Discussion. All these approaches, in general, need a clean mesh as an input. Furthermore, most of them are not generic enough to handle heterogeneous input data, and they are not designed to guarantee a valid, intersection-free output.

Self-intersection Free

By Prevention. Most of the algorithms provide guarantees for intersection-free output during mesh decimation by building an envelope volume around the surface. [START_REF] Cohen | Simplification envelopes[END_REF] provided a generic genus-preserving simplification algorithm which guarantees a two-sided Hausdorff error bound along with self-intersection prevention and sharp feature preservation. The envelope is created by simply offsetting each vertex of the original surface in the direction of its normal vector to transform the fundamental triangles into those of the envelope. Offset vertices are not allowed to go beyond the Voronoi regions of its adjacent fundamental triangles to avoid selfintersection in within the envelope (Figure 1.16). The original mesh is then decimated within the envelope.

The main disadvantage of the tolerance volume based approaches is that the space needs to be split. As a result, we cannot cross the tolerance even if the operation corresponds to an optimal target location with no self-intersection. [START_REF] Gumhold | Intersection free simplification[END_REF] proposed a method based on preventing and avoiding self-intersections during the mesh decimation. However, searching for the locus of points in space that avoids intersections when applying a decimation operator, and tunneling out of situations where every operator is forbidden is often too labor-intensive to be considered a practical solution.

By Construction. Another class of approaches based on Delaunay filtering and refinement, instead, provide intersection-free approximations by construction. [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF] n 1 n 2 n 3

Voronoi edge

Figure 1.16: Creating envelope using offset vertices. Left: The vertices are offsetted in the normal direction by the desired Hausdorff error to create the Simplification envelope [START_REF] Cohen | Simplification envelopes[END_REF]. Right: To avoid self-intersection within the envelope, the offset vertices are bounded by the Voronoi regions. This however, might reduce the level of approximation in certain regions. Note that the boundaries of this envelope will be isotopic to original surface.

presented an algorithm that constructs provably good surface samples and meshes. A notable feature of the algorithm is that the surface needs only to be known through an oracle that, given a line segment, detects whether the segment intersects the surface and, in the affirmative, returns the intersection points making the algorithm useful in a wide variety of contexts and for a large class of surfaces. Note that all the methods discussed above in Section 1.2.1 which provide topological guarantees obviously fall in this catogery.

Discussion. The algorithms which provide self-intersection free guarantee by prevention need a clean mesh to start with. Moreover, reconstructing a clean surface first and then using these algorithms will accumulate the Hausdorff error and the output will not be within the desired bounds. Other algorithms which guarantee intersection-free by construction, unfortunately, generate only isotropic meshes (as we saw before) and do not target very coarse approximations and, as such, cannot be used for shape simplification.

Positioning

In ℝ 3 , Chazal and Cohen-Steiner [2004] showed that when seeking a homeomorphic approximation 𝑆 ′ of a connected surface 𝑆, a simple topological condition is sufficient to guarantee that the two surfaces are isotopic (Figure 1.17). 𝑆 and 𝑆 ′ are isotopic iff:

1. 𝑆 ′ is contained in a topological thickening of 𝑆 and separates the boundary components of this thickening, and, 2. 𝑆 and 𝑆 ′ are homeomorphic.

We contribute next a constructive approach for this theoretical result in the form of an algorithm that matches these conditions in order to ensure that the output surface mesh is an isotopic approximation.

Positioning

S Ω S Ω Figure 1.17: Isotopic conditions. Two connected and closed surfaces 𝑆 and 𝑆 ′ are isotopic iff (𝑎) 𝑆 ′ is included in a topological thickening of 𝑆 and separates its boundary components, and (𝑏) the genus of 𝑆 ′ does not exceed the genus of 𝑆. Note that if 𝑆 ′ is closed and separates boundary components of topological thickening of 𝑆, then its genus cannot be lower than the genus of 𝑆.

Problem Statement

We state the problem as follows. The input is a tolerance volume 𝛺 (Figure 1.17) that is a topological thickening of a surface 𝑆 which we want to approximate. By topological thickening of 𝑆 we mean a compact subset of ℝ 3 homeomorphic to 𝑆 × [0, 1]. Our goal is to generate as output, a surface triangle mesh:

1. located within 𝛺, 2. isotopic to the boundary components of 𝛺, and 3. with a low triangle count.

Note that our problem statement contains all the features which were missing from the current state of the art. The first criteria ensures that the approximation error is bounded by the Hausdorff distance of the topological thickening. It also ensures that the output separates the boundary components of 𝛺. This will later help us in providing theoretical guarantees over its topology (isotopy). Although we have explicitly mentioned only Hausdorff distance in the problem statement, other fairness schemes, for example, normal deviation (as discussed in Section 3.1) can always be added. Secondly, by definition, the boundary components of a topological thickening are always isotopic to the original surface. Furthermore, assuming that the input is provided in the form of tolerance volume ensures that original surface 𝑆 need not be provided in a specific form, and hence, our algorithm can handle heterogeneous data. The second criteria also ensures that our output is free from self-intersections. Finally, the third criteria ensures that we solve the mesh simplification problem (Section 1.1.2) and the output size is as low as possible.

This approximation problem was originally stated by Klee for polytopes in arbitrary dimensions. In 2D, the problem is commonly referred to as the minimum nested polygon problem, and has been investigated well [START_REF] Aggarwal | Finding minimal convex nested polygons[END_REF]. The 3D instance of this problem, referred to as minimum nested polyhedron problem has been shown to be NP-hard [START_REF] Agarwal | Surface approximation and geometric partitions[END_REF]].

Discussion

Despite being a long standing problem, there is still no robust and practical solution to this enduring scientific challenge. Yet, it is both relevant to, and timely for, the increasing variety of industrial applications that involve raw geometric data. In this part of the thesis, we develop an algorithm for the above problem that yields approximations with very low triangle count, while enjoying topological guarantees under relatively mild assumptions on the tolerance volume.

Note that while the assumption that 𝛺 is a proper thickening makes the analysis easier, it is not always necessary and our approach may also work when boundary components of 𝛺 have, for instance, additional spurious handles. We also extend our algorithm to non-closed and non-manifold surfaces. If 𝛺 is not provided as input, we may generate it from a possibly defectladen approximation of 𝑆 (𝛴, e.g., a point cloud or a polygon soup) using either simple offsets in the noise-free case, or sublevel sets of a robust distance function (e.g. [START_REF] Chazal | Geometric Inference for Measures based on Distance Functions[END_REF]).

Hence, under relatively mild conditions, our algorithm is able to solve the problem of robust reconstruction, repair and simplification concurrently.

Chapter 2 Algorithm and Guarantees

Algorithm Overview

Initialization

For initialization, we generate a 𝜎-dense set 𝒮 sampled on the tolerance boundary 𝜕𝛺, 𝜎 being typically set to a fixed fraction of the minimum separation 𝛿 between the 𝜕𝛺 𝑖 . That is, the balls of radius 𝜎 centered on 𝒮 cover 𝜕𝛺.

For the base algorithm we assume that 𝜕𝛺 has only two components 𝜕𝛺 1 and 𝜕𝛺 2 . We assign to each sample 𝑠 of 𝒮 a function value:

ℱ(𝑠) = +1 if 𝑠 ∈ 𝜕𝛺 1 , and ℱ(𝑠) = -1 if 𝑠 ∈ 𝜕𝛺 2 .
We then construct an initial 3D Delaunay triangulation (𝒯) with the eight corners of a loose bounding box of 𝒮. We assign to these eight vertices the same function value as that of the samples of the outer boundary of 𝛺. We maintain a piecewise-linear function 𝑓 interpolated on 𝒯, and its zero-set, denoted by 𝒵. During refinement of 𝒯 with a subset of 𝒮 (described next), the classification of 𝒮 provides us with a means to detect when 𝒵 lies within 𝛺, with a safety margin (defined in Section 2.6).

Figure 2.2 illustrates in 2D a refinement on 𝒯, the corresponding zero-set 𝒵 and the classification of 𝒮.

Refinement

We now refine the triangulation 𝒯 through inserting Steiner points selected from 𝒮, until the correct topology is met. More specifically, we insert one sample point at a time into 𝒯 and update the Delaunay property, until 𝒵 classifies all samples of 𝒮 as good, or equivalently, until 𝒵 separates the boundaries 𝜕𝛺 𝑖 of 𝛺.

Greedily inserting the sample 𝑠 with maximum error at each step is a natural idea for achieving the above goal with few samples. For each tetrahedron we maintain a list of sample points (⊂ 𝒮) contained in it, and a global modifiable priority queue during refinement with the maximum error points of these tetrahedra. Figure 2.3 illustrates several steps of a refinement sequence in 2D, until complete classification of 𝒮.

Unfortunately, the above basic refinement algorithm is not sufficient for at least two reasons. The first reason relates to the fact that we are dealing with a finite sample of 𝜕𝛺. Even if all sample points end up being well classified, this still leaves the possibility that 𝒵 crosses 𝜕𝛺 in-between the samples. To prevent this from happening, we enforce that all samples are well classified with an 𝛼 margin, as well as an upper bound on the Lipschitz constant of the piecewise-linear function.

The second reason relates to the quality of normals. In certain configurations (e.g., Figure 2.4), the normal directions are grossly wrong even in locally smooth areas. To alleviate this issue we detect so-called misoriented tetrahedra by checking that the piecewise linear function they define is locally well adapted to the geometry of 𝛺 (condition 3 below). We note that this condition is not required for the topological correctness of the algorithm. Our modified refinement algorithm iteratively refines the triangulation until all the following criteria are met in order:

1. For some given 0 < 𝛼 < 1 ∶ ∀𝑠 ∈ 𝒮, 𝜖(𝑠) ≤ 1 -𝛼 (𝛼 is set to 0.2 in all experiments).

2. The height of every tetrahedron contributing to 𝒵 is at least 2𝜎/𝛼. The height is defined as the distance between the supporting lines or planes of the maximal faces with different labels (Figure 2.5).

3. The piecewise-linear function defined by each tetrahedron 𝑡 classifies well the samples of 𝒮 on both 𝜕𝛺 𝑖 that are nearest to the vertices of a shrunk copy of 𝑡. The size of this shrunk copy is set to 70% of the size of 𝑡 in all experiments.

The term "in order" herein means that at each iteration, we look at the first condition that is violated and attempt to satisfy it by inserting a Steiner point as described below. If the condition is not satisfied after exhausting all candidate Steiner points, we move to the next condition.

Since we are dealing with a 𝜎-dense sample, for an 𝛼 margin (condition 1), an upper bound of 𝛼/𝜎 on the Lipschitz constant of the piecewise-linear function suffices to ensure that the zero-set does not cross 𝜕𝛺. Noticing that the Lipschitz constant is nothing but twice the inverse height of a tetrahedron, we get an easy-to-check criterion (condition 2). The first criterion is met by adding the sample point with maximum error while the two other criteria are met by adding the sample point nearest to the circumcenter of a bad tetrahedron. The full refinement algorithm incorporates one more condition (condition 4): while the output 𝒵 of the above algorithm does not have the expected genus, we refine the heterogeneous tetrahedron with the largest circumradius by adding the sample point closest to its circumcenter.

This additional layer is needed to get topological guarantees on the result. However, in practice, we did not encounter a single case where the genus was not correct after the first iteration. Also note that the circumradius criterion for refining tetrahedra is blind in the sense that it does not necessarily refines the mesh where topological defects are present. While it is possible to improve the criterion from this point of view, we did not pursue this goal since it has no practical relevance.

Upon termination of the refinement step, the union of all tetrahedra of 𝒯 which contribute to 𝒵, bound a simplicial tolerance volume (𝛤 ), seen as an approximation of 𝛺. The boundary facets of 𝛤 are denoted by 𝜕𝛤 .

Simplification

The zero-set 𝒵 is now topologically correct. In the simplification step we reduce its complexity via decimation of 𝒯 combined with a mutual tessellation with 𝒵. Note that we stop enforcing that 𝒯 is a Delaunay triangulation, which allows for increasingly anisotropic triangulations.

Simplification is achieved through performing a series of edge-collapse operators on 𝒯. These The validity of 𝒯 requires checking for two conditions. The combinatorial topology of 𝒯 is preserved via the link condition [START_REF] Dey | Topology preserving edge contraction[END_REF]] (Figure 2.6).

The valid embedding of 𝒯 is preserved by computing the visibility kernel (𝒦 𝒯 )(𝑃 𝑄) of the polyhedron formed by the one-ring of the edge 𝑃 𝑄 (Figure 2.7). If the visibility kernel is non empty then locating the target vertex into this kernel preserves a valid embedding. Preserving the classification of 𝒮 requires further restricting the visibility kernel of edge. As this problem is non-convex, we resort to a point sampling of the kernel during the simulation of each edge-collapse operator. To obtain faithful normals locally in a smooth area, we use the same method as before (Figure 2.4) and check in advance whether the final solution is locally well adapted to the geometry of 𝛺 or not.

P Q P Q K T (P Q)
In order to improve efficiency we always perform simpler halfedge collapse before general edge collapse operators. A halfedge collapse operator locates the target vertex at one of the vertices of the edge. In addition, we adopt a multi-staged decimation approach with the following steps:

1. Collapse edges of 𝜕𝛤 .

2. Mutual Tessellation of 𝒵 into 𝒯.

Collapse edges of 𝒵.

4. Collapse edges between 𝛤 and 𝒵, which may induce further edge collapses of 𝒵 (previous step).

Intuitively, we perform the steps in increasing order of computational complexity: first the operations with low number and discrete degrees of freedom, then with higher or continuous degrees of freedom. As for other decimation algorithms we need to define an error to sort the operators and to optimize the target vertex placement when performing a general edge collapse operator.

In order to preserve fidelity to the initial zero-set, we use as error the sum of square distances between the target vertex and the set of supporting planes of the zero-set facets located in the 2-ring of the collapsed edge. The edge collapse operators are sorted via a priority queue sorted by increasing error.

Simplicial Tolerance

In this step we collapse only a subset of the edges of the simplicial tolerance boundary 𝜕𝛤 .

Denote by 𝑃 𝑄 such an edge (Figure 2.8). We select as target vertex a sample point (1) from 𝒮

(2) located within the visibility kernel 𝒦 𝒯 (𝑃 𝑄) of 𝑃 𝑄, (3) inducing a zero-set that preserves the classification of 𝒮 along with normals and (4) that minimizes the aforementioned error. Denote by 𝑆 |𝒦 𝒯 (𝑃 𝑄) the initial set of candidate sample points from 𝒮 located within the visibility kernel 𝒦 𝒯 (𝑃 𝑄). To avoid exhaustive search, we discard the sample points leading to errors in the classification of 𝒮, as located in invalid regions, denoted 𝛹 . Figure 2.9 illustrates 𝛹 = 𝑎 ∩ 𝑏 ∩ 𝑚 where the point with maximum error is chosen only over 𝜕𝛺 1 . A similar invalid

P Q K T S |K T (P Q) Ω ∂Ω 1 ∂Ω 2
A B T m n b a E X Y A B E T n m b a X Y = = Ψ Ψ ∂Ω 1 ∂Ω 2 ∂Ω 2 ∂Ω 1 = = Figure 2
.9: Invalid region. Assume an edge of 𝜕𝛤 is collapsed into the target point 𝑇 . Line segment 𝑋𝑌 denotes the zero-set of △𝐴𝐵𝑇 after collapse. Line 𝑛 represents the extreme zero-set of △𝐴𝐵𝑇 which preserves the classification of the point with maximum error 𝐸. Line 𝑚 delineates the corresponding locus for 𝑇 . The intersection of the two half-spaces delineated by 𝑎 and 𝑏 represents the locus of 𝑇 which keeps 𝐸 within △𝐴𝐵𝑇 . If 𝑇 is located in the invalid region 𝛹 (gray) then the classification of 𝐸 is not preserved. Left: case where 𝐴 and 𝐵 belong to the same 𝜕𝛺 𝑖 . Line 𝑛 is parallel to 𝐴𝐵 and passes through 𝐸. Right: case where 𝐴 and 𝐵 belong to two different 𝜕𝛺 𝑖 . Notice that 𝑌 is fixed and 𝑛 is the supporting line of 𝐸𝑌 .

region is computed by considering the point of maximum error on 𝜕𝛺 2 within △𝐴𝐵𝑇 . We then collapse iteratively all edges of 𝜕𝛤 to the point which exhibits the minimum error, as discussed above. When no more edges of 𝜕𝛤 are collapsible, we perform a mutual tessellation between 𝒵 and 𝒯 by inserting all vertices and faces of 𝒵 into 𝒯. The newly inserted vertices are assigned the function value ℱ = 0. We then label all tetrahedra of 𝒯 in accordance to their associated tolerance boundary component 𝜕𝛺 𝑖 . This provides us with a means to preserve the classification in the next simplification steps. A sample 𝑠 ∈ 𝜕𝛺 𝑖 is constrained to lie within a tetrahedron with label 𝑖. Intuitively, this step implements a transition from a function embedded in a volume mesh of the tolerance, to a surface mesh embedded within the 3D triangulation. Figure 2.10 illustrates such mutual tessellation in 2D.

Mutual Tessellation

Zero-set

After mutual tessellation we collapse the edges of 𝒵. Figure 2.11 illustrates the visibility kernel of an edge that preserves a valid embedding upon a collapse operator. To accelerate the computations, we compute invalid regions as described above. Figure 2.12 illustrates the invalid region by considering the point of maximum error on both 𝜕𝛺 𝑖 for a

P Q K T (P Q)
△𝐴𝐵𝑇 when it contains one zero-set vertex. The invalid regions 𝛹 are constructed similarly when △𝐴𝐵𝑇 contains several zero-set vertices. To further reduce the computational time when simulating general edge collapse operators, we use an octree for hierarchical sampling of 𝒦 𝒯 to find the best target location and ignore further sampling of 𝒦 𝒯 for the octree cells lying inside 𝛹 . Figure 2.13 provides a summary of the conditions required for a valid edge collapse.

P Q ∂Ω2 ∂Ω1 Z Q P KT Q P Q P Figure 2
.13: Conditions for valid edge collapse. From left to right: To collapse an edge 𝑃 𝑄, we check for its link condition and compute the visibility kernel (middle left). Note that, all the points within this visibility kernel may not preserve the classification of sample points. The region (subset of visibility kernel) that preserves the classification of sample points is depicted in green (middle right). Since, finding such a kernel is a non-convex problem, we resort to hierarchical sampling of the visibility kernel to find the target point which preserves classification and minimizes the sum of square distances between the target vertex and the set of supporting planes of the zero-set facets located in the 2-ring of the collapsed edge as depicted in green (right). Euclidean distance function to the input triangle soup. Note that until mutual tessellation the zero-set is made up of triangles and quadrangles before being converted into a pure triangle mesh after mutual tessellation.

All Edges

Extensions

The algorithm above can be extended in order to deal with boundaries and non-manifold surfaces.

While the guarantees can be extended to non-closed surfaces, for non-manifold cases, they are more difficult to formulate and are probably beyond the reach of existing tools. Still, the output is guaranteed to have the correct homotopy type in these cases also.

Non-closed Surfaces

Our algorithm deals with surfaces with boundaries when 𝛴 is provided as input along with 𝛺.

We first detect sample points corresponding to boundaries, referred to as ℬ, using 𝛾𝛿 ′ -radius balls centered at 𝛴 (Figure 2.17), where 𝛾 denotes a user defined parameter derived from the reach of input data. Parameter 𝛿 ′ is defined as the minimum distance between the current center of ball and 𝒮.

Non-manifold Surfaces

To handle non-manifold surfaces when computing the error of a sample, we evaluate 𝑓 with respect to each component of 𝜕𝛺. More specifically, to evaluate the error 𝜖 at a sample 𝑠 ∈ 𝜕𝛺 𝑖 we define

∀𝑝 ∈ 𝒮, ℱ(𝑝) = ⎧ { ⎨ { ⎩ +1, 𝑝 ∈ 𝜕𝛺 𝑖 -1, 𝑝 ∉ 𝜕𝛺 𝑖 .
(2.

2)

The zero-set is ignored when its end points lie outside 𝛺; this configuration occurs when the input geometry is made up of several components (Figure 2.19). where several surfaces of 𝜕𝛺 meet. These artifacts are located inside tetrahedra containing vertices from three or more 𝜕𝛺 𝑖 . In addition, such tetrahedron may contain several zero-sets corresponding to the total number of possible permutations when assigning function values to its vertices. We remove these artifacts by joining all zero-set edges to the centroid of the zero-set vertices located on the edges of this tetrahedron. Figure 2.20 illustrates our algorithm at work on a non-manifold geometry.

Guarantees

We first derive geometric conditions under which the first three conditions of the refinement algorithm are met upon termination. Denote by 𝜀 the radius of the largest ball that can fit within 𝛺, and by 𝛿 the minimum separation between the two boundary components of 𝛺.

Condition 1 (2.3-1) is necessarily satisfied at the end of the refinement process since any sample not meeting the condition will be added to the triangulation.

Assume now that Condition 2 (2.3-2) is not satisfied upon termination. This means that there exists a heterogeneous tetrahedron 𝑡 with height lower than 2𝜎/𝛼. Denote by 𝐵 its circumscribed ball and 𝑟 its radius. Ball 𝐵 cannot contain any sample point from 𝒮, else that sample would have been added by the algorithm. Since 𝒮 is a 𝜎-sample of 𝜕𝛺, we get that the shrunk ball 𝐵 -𝜎 does not intersect 𝜕𝛺. Hence it is either empty, within 𝛺, or outside 𝛺. In the first case, 𝑟 ≤ 𝜎.

In the second case, we have that 𝑟 -𝜎 ≤ 𝜀. In the third case, because 𝑡 is heterogeneous, 𝐵 meets both boundary components of 𝛺, hence 𝜎 ≥ 𝛿. As a partial conclusion, upon termination, and assuming 𝜎 < 𝛿, the circumradius 𝑟 of a tetrahedron violating condition 2 cannot exceed

𝜀 + 𝜎. ρ ∂Ω 2 ∂Ω 1 ≥ h x
We now formulate our condition. Given two subsets 𝐴 and 𝐵 of ℝ 3 , define the margin of (𝐴, 𝐵) to be the maximum thickness of a slab separating 𝐴 and 𝐵. If no such slab exists then the margin is set to zero. We say that a tolerance volume 𝛺 is (𝜌, ℎ)separated if for all 𝑥 ∈ ℝ 3 , the margin of (𝜕𝛺 1 ∩ 𝐵(𝑥, 𝜌), 𝜕𝛺 2 ∩ 𝐵(𝑥, 𝜌)) is at least ℎ.

From the above discussion, if 𝜎 < 𝛿, and if 𝛺 is (𝜀 + 𝜎, 2𝜎/𝛼)-separated, condition 2 will be satisfied at the end of the algorithm. Similarly, condition 3 (2.3-3) will ultimately hold assuming a stronger local separation assumption on 𝛺. However, since this condition is not essential for the topological correctness of our algorithm, we do not elaborate further on expliciting the required separation constants.

Finally, concerning condition 4, we note for future reference that if the correct genus is not met upon termination, we can bound the circumradius 𝑟 of any heterogeneous element as above. That is, assuming 𝜎 < 𝛿, we have that 𝑟 ≤ 𝜀 + 𝜎.

Theorem 2.6.1. Let 𝜅 be such that for any two points 𝑥 and 𝑦 in 𝜕𝛺 𝑖 at distance at most 2(𝜀 + 𝜎), the geodesic distance between 𝑥 and 𝑦 is at most 𝜅 times their Euclidean distance.

Assuming 𝛺 is a ((5 + 𝜅)(𝜀 + 𝜎)/2, 2𝜎/𝛼)-separated topological thickening of a surface 𝑆, the output of the algorithm is isotopic to 𝑆.

Proof. Our proof for isotopy utilizes the same criterion as [START_REF] Chazal | A condition for isotopic approximation[END_REF].

This result states that two connected 3-dimensional surfaces 𝑆 ′ and 𝑆 are isotopic if:

1. 𝑆 ′ is included in a topological thickening of 𝑆 and separates its boundary components.

2. 𝑆 ′ is connected and its genus does not exceed the genus of 𝑆.

In our setting, we take 𝛺 to be a topological thickening of 𝑆. Being the zero-set of a piecewiselinear function, 𝒵 is a 2-manifold. Because condition 2 of the base algorithm is satisfied, all points in 𝜕𝛺 are well classified, meaning that 𝒵 is contained inside 𝛺 and separates its boundary components.

It remains to show that the second condition holds when the refinement algorithm terminates.

First, because 𝑆 ′ is a zero-set of a piecewise linear function, any component of 𝑆 ′ must enclose a vertex of the Delaunay triangulation. As there are no vertices in the interior of 𝛺, these components must induce non-zero homology classes in 𝛺. These components are included in the simplicial tolerance, which is fibered by line segments where the piecewise linear function is monotone. This implies that 𝑆 ′ is connected.

S

Now assume that the genus of 𝑆 ′ exceeds the one of 𝑆. This means that 𝑆 ′ contains a spurious handle. Because sub and superlevel sets of piecewise linear functions are homotopy equivalent to subcomplexes of the background triangulation, we get that the two components of 𝛺\𝑆 ′ contain linked homogeneous polygonal cycles in the Delaunay triangulation. For each edge in these polygonal cycles, we may form an elementary cycle by stitching the edge with a geodesic shortest path drawn on the appropriate boundary component of 𝛺 and joining the two endpoints of the edge. Because the two polygonal cycles are linked, there must be two linked elementary cycles with different labels.

Assume the correct genus is not met upon termination of the algorithm. Then we may assume that the Delaunay edges in the linked polygonal cycles are edges of heterogeneous tetrahedra.

Hence their length is at most 2(𝜀 + 𝜎). Hence the length of the elementary cycles are at most 2(1 + 𝜅)(𝜀 + 𝜎). Because two of them are linked, the tolerance volume 𝛺 cannot be ((5 + 𝜅)(𝜀 + 𝜎)/2, 0)-separated (Appendix A.0.1).

The above conditions are clearly met when the tolerance volume is a sufficiently small offset of a smooth surface, for example. Also, the algorithm can also be shown to work in situations not covered by the above theorem, as for instance tolerances bounded by convex surfaces.

Ω ∂Ω 1

∂Ω 2

The numerical constants in the theorem may be further optimized, and other types of conditions, e.g. based on the separation 𝛿 can also be proved to be sufficient. In particular, if 2(𝜀+𝜎) √ 𝜅 2 -1 < 𝛿, the algorithm is correct. The inset figure depicts an example of tolerance volume where the algorithm would fail. It is apparent from the proof above that such configurations are essentially the only way the algorithm can fail. Note that even in such situations, the output of the algorithm will be a manifold surface, albeit possibly with a too large genus.

Finally, we note that if the only pursued goal was to provide topologically guaranteed output for tolerances that are topological thickenings, then a trivial solution would be to output one of the tolerance boundary surface. However, such methods would be limited to topological thickenings, while our approach may work in less favorable situations. In addition, algorithms inspired by the trivial idea above do not seem to allow successful subsequent simplifications, even in favorable cases.

All atomic operations performed over the tetrahedra or sample points are easily parallelized as they are independent. All experiments are performed on an Intel 2.4GHz 16-core machine with 128 Gb RAM. The tolerance errors are specified as a percentage of the longest edge of the bounding box of the input data. Margin 𝛼 is set to 0.2 in all experiments.

Though for simplicity of exposition, we assigned ℱ at bounding box vertices as that of the outer boundary. In practice, we found that by multiplying this assignment by the distance to shows outputs of our algorithm on a smooth surface (fertility) with 𝛿 ranging from 0.15% to 8%.

𝛺
The one-sided Hausdorff is measured from the output to the input and is bounded in all cases. Note that when 𝛿 is large (bottom right), the tolerance volume is not a topological thickening anymore as the topology of the inner boundary of the tolerance changes. We also run our algorithm on Armadillo (Figure 3.3) and Vase Lion (Figure 3.4) -more general meshes made of smooth and flat parts.

Experiments

To evaluate our algorithm, we perform several experiments and compare it with other state of the art algorithms as follows:

Vertex count over time. Comparisons. A strict qualitative comparison with previous work is not possible as our problem statement differs. The one-sided Hausdorff distance preserved in general mesh decimation algorithms is measured from the input to the output mesh, while we guarantee the other side of the Hausdorff distance. Nevertheless, we plot our one-sided Hausdorff distance against the number of vertices of the final output mesh, for five other mesh approximation algorithms: simplification envelopes [START_REF] Cohen | Simplification envelopes[END_REF], a decimation algorithm from Lindstrom-Turk [START_REF] Lindstrom | Evaluation of memoryless simplification[END_REF] (without error bounds), the MMGS remeshing algorithm [START_REF] Borouchaki | Simplification of surface mesh using hausdorff envelope[END_REF] (with Hausdorff error bound), MMGS with the mesh anisotropy option and a decimation algorithm with error bound implemented in OpenFlipper [START_REF] Möbius | Openflipper: An open source geometry processing and rendering framework[END_REF]. 0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 9,000

0 1 2 3 Vertex count 𝐻 |𝑂→𝐼 (%)

Our algorithm Simplification Envelopes

Lindstrom-Turk MMGS MMGS -A OpenFlipper 0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 9,000

0 1 2 Vertex count 𝐻 |𝐼→𝑂 (%)

Our algorithm Simplification Envelopes

Lindstrom-Turk MMGS MMGS -A OpenFlipper Albeit none of the other approaches except [START_REF] Cohen | Simplification envelopes[END_REF] target an intersection-free output, we indicate with a dot a self-intersection of the output mesh (Figure 3.6). For completeness we also plot the other-sided Hausdorff distance over the same input and output datasets (Figure 3.7). In both cases we achieve a lower vertex count for a given tolerance error, at the price of higher computational times. We also compare our algorithm with [START_REF] Cohen | Simplification envelopes[END_REF] on Armadillo (Figure 3.8). For a very large tolerance, the vertex count for our algorithm and simplification envelopes is comparable. However, on most part of the curve, our algorithm generates on average 10% fewer vertices, for a given tolerance error. Also note that the simplification envelopes require a manifold mesh as input. In addition, they cannot simplify the geometry of highly undulating surfaces beyond a certain limit, due to the specific type of tolerance volume used (Figure 3.9).

Normals.

We do not provide quantitative guarantees of faithful approximation of normals in all cases. Instead, our proof (Section 2.6) yields good normals in smooth areas, when the sampling 𝜎 is dense enough. Figure 3.10 plots the distribution of normal deviation with respect to different shrinkage factors used for condition 3 for the Fertility model. Dropping this condition may cause in practice the normal deviation to exceed 90 ∘ . Figure 3.11 provides on the Lucy model a visual comparison of normals of our output to the input surface triangle mesh (left).

Dealing with sharp creases subtending small angles may require an extremely dense 𝜎, which also translates into dense refinement. When using a too large value for 𝜎, activating the preser- vation of normals during refinement may further translate into dense refinement as the normals are ill-defined locally. We cannot always deduce whether overly dense refinement comes from the preservation of inferred normals or from the recovery of the topology. As the decimation preserves the inferred normals we may end up with overly complex meshes. However, relaxing constraint 3 during halfedge collapses can alleviate this issue as illustrated in Fig 3 .12. We plot the distribution of normal deviation of our output for the fertility model (𝛿 = 0.2%) when condition 3 is not used, and when a 30%, 70% and 100% shrinkage factor is used. The final vertex count of the output is 3, 498, 3, 538, 3, 624 and 5, 715 respectively. Robustness. A primary virtue of our algorithm is its resilience to the type and defects of the input dataset. More specifically, and as our algorithm takes a tolerance volume as input, the robustness of our algorithm is delegated to the construction of a robust tolerance volumethen the output is guaranteed to be homotopy-equivalent to the given tolerance volume. Said differently, our algorithm is oblivious to the dataset inside the tolerance as long as the tolerance is well-behaved. Figure 3.13 illustrates the robustness of our algorithm on point sets and defectladen triangle soups sampled on the elephant model, with two levels of noise. We use as tolerance volume a sub-level of the robust distance function based on distances between measures [START_REF] Chazal | Geometric Inference for Measures based on Distance Functions[END_REF].

Limitations. Despite its guarantees and qualitative performances, our algorithm is computeintensive, especially when setting a small tolerance. In Figure 2.16 the tolerance is set to 𝛿 = 0.6% and our algorithm runs for approximately 3h and consumes 2.1𝐺𝑏 of peak RAM memory. The time complexity is dominated by the simplification step. Table 3.1 lists the time taken by each step of the algorithm against the vertex count of 𝒵.

Chapter 4 Conclusion

We introduced a novel approach to the problem of isotopic approximation within a tolerance volume. We depart from common approaches by leveraging a dense point set sampled on the boundary of the input tolerance volume. We designed a pliant meshing algorithm that first proceeds by Delaunay refinement in order to recover the correct topology through classifying the samples, and then by topology-preserving simplification in order to discover the anisotropy within the tolerance volume. A distinctive feature of our approach is its robustness to input data sets.

As our approach is oblivious to the type and defects of the datasets within the tolerance volume, it can reconstruct, repair and simplify concurrently. Compared to error-driven simplification algorithms, with or without error bounds, our approach makes full use of the tolerance volume and achieves lower vertex counts for a given tolerance error, in addition to intersection-free outputs. Such lower vertex counts, however, come at a price: our current implementation is compute-intensive, and the computational times escalate when the tolerance decreases.

Future Work. The current version of the algorithm is highly parallelizable: in practice we observed that the running times are inversely proportional to the number of processors (Fig 3.14).

As future work we wish to extend our approach so as to make it out-of-core. This will help us tackle the problem of both space and time complexity which are currently governed by the sampling density and simplification step. A natural direction is to cut the tolerance volume into sub-parts before stitching, but it requires another line of work to preserve the guarantees during refinement and simplification. We also observe that as the tolerance decreases, the vertex count at the end of refinement stage increases exponentially (Fig 3.5). While this is partly due to the increase in complexity (and number of features) of the shape, it is mainly due to the fact that the Delaunay triangulation is intrinsically isotropic. This not only increases the running times of the refinement step, but also has a direct impact on that of the simplification step: Over-refinement requires more edge-collapse operations for simplification. An interesting topic for future work is to better analyze this behavior and restrict over-refinement possibly via relaxing the Delaunay property of the triangulation or via a constrained Delaunay triangulation.

Note that our algorithm assumes that a valid tolerance -which is a topological thickeningis provided as input. However, we do not provide an algorithm which is guaranteed to generate such a valid tolerance. Some previous work such as the medial axis transform [START_REF] Moon | Mathematical theory of medial axis transform[END_REF]Amenta et al., 2001b] and the concept of big Delaunay balls [START_REF] Dey | An adaptive mls surface for reconstruction with guarantees[END_REF] may provide a solution to this issue. They may also help us providing a guarantee over the two-sided Hausdorff distance which is currently not provided.

Finally, another stimulating direction is the concept of a progressive approximation algorithm, in which we could guarantee that every additional CPU cycle spent by the algorithm is making progress toward the optimal solution that matches the global minimum vertex count. Finding a mapping between two discrete surfaces is a recurring problem in geometry processing, motivated by a wide range of applications including template fitting, attribute and animation transfer, shape completion, blending, morphing, and remeshing. It consists in establishing a meaningful correspondence map between two input surfaces [START_REF] Van Kaick | A survey on shape correspondence[END_REF],

Chapter 5 Introduction

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], or shape alignment in the graphics literature.

Given two arbitrary shapes, an infinite number of mappings exist between them. If the two shapes are isometric, efficient algorithms exist to establish a good map between these shapes.

However, in most practical cases where shapes may differ significantly, there is no universal consensus on what defines the best correspondence map, and one typically invokes a geometric prior such as conformality or area-preservation to narrow down the type of suitable maps to look for. For inter-surface maps to be useful to subsequent geometry processing, they should be bijective (one-to-one) and have some degree of continuity, i.e., small neighborhoods should map to bounded neighborhoods. Another common mapping requirement is semantic meaningfulness:

semantic features in one shape should map to the corresponding features in the other shape.

Background

Before detailing the exact problem formulation and our algorithm, we mention the numerical tools that will help us to solve the problem robustly and efficiently.

Diffusion geometry.

For a Riemannian manifold 𝑀 , its Laplace-Beltrami operator admits eigenvalues 𝜆 𝑖 and corresponding eigenfunctions 𝜙 𝑖 that satisfy △ 𝑀 𝜙 𝑖 = 𝜆 𝑖 𝜙 𝑖 . This operator can be well approximated for both triangulated [START_REF] Meyer | Discrete Differential-Geometry Operators for Triangulated 2-Manifolds[END_REF][START_REF] Desbrun | Discrete differential forms for computational modeling[END_REF] and polygon surfaces [START_REF] Wardetzky | Discrete laplacians on general polygonal meshes[END_REF], as well as pointsets [START_REF] Belkin | Constructing laplace operator from point clouds in rd[END_REF][START_REF] Liu | Point-based manifold harmonics[END_REF] and has been extensively used in geometry processing [START_REF] Sorkine | Differential representations for mesh processing[END_REF]. [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps[END_REF]; [START_REF] Nadler | Diffusion maps, spectral clustering and eigenfunctions of Fokker-Planck operators[END_REF] proposed the notion of diffusion distance which has also been lately introduced in geometry processing [START_REF] Rustamov | Laplace-Beltrami eigenfunctions for deformation invariant shape representation[END_REF][START_REF] Curran Associates | A hierarchical segmentation of articulated bodies[END_REF]. If we denote the heat kernel 𝐾 𝑀,𝑡 as 𝐾 𝑀,𝑡 (𝑥, 𝑦) = 𝛷(𝑥) ⊤ 𝛷(𝑦), (5.1)

where 𝛷(𝑥) = (𝑒 -𝜆 1 𝑡/2 𝜙 1 (𝑥), 𝑒 -𝜆 2 𝑡/2 𝜙 2 (𝑥), … ) for any point 𝑥 ∈ 𝑀 and 𝑡 controls the diffusion time scale, then the diffusion distance between two points 𝑥 and 𝑦 on manifold 𝑀 can be computed as 𝑑 diff (𝑥, 𝑦) = ‖𝛷(𝑥) -𝛷(𝑦)‖.

(5.2)

Note that the heat kernel and diffusion distances are related via 𝑑 diff (𝑥, 𝑦) 2 = 𝐾 𝑀,𝑡 (𝑥, 𝑥) + 𝐾 𝑀,𝑡 (𝑦, 𝑦) -2𝐾 𝑀,𝑡 (𝑥, 𝑦).

For small 𝑡, the diffusion distance approximates the local geodesic distances well, while increasing 𝑡 generates smooth long range geodesic distances [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps[END_REF]. Moreover, by selecting only the first 𝑛 smallest eigenvalues of the Laplace-Beltrami operator (which can be solved efficiently [START_REF] Vallet | Spectral geometry processing with manifold harmonics[END_REF]), the diffusion distances can be approximated at low computational cost. Note that the diffusion distance Eq. ( 5.2) isometrically embeds in (Euclidean) diffusion space, i.e. the geodesic distance between two samples 𝑥 and 𝑦 in ℝ 3 is approximated by the Euclidean distance between 𝛷(𝑥) and 𝛷(𝑦) in ℝ 𝑛 .

Optimal transportation

The problem of optimal transportation was originally proposed by Monge in 1781 [START_REF] Monge | Mémoire sur la théorie des déblais et des remblais[END_REF]:

given two mass distributions m and n ≥ 0 on ℝ 𝑑 , with ∫ m = ∫ n = 1, the aim is to find a map where, 𝐴 denotes any Borel set in ℝ 𝑑 . Kantorovich also studied the problem independently and proposed a suitable framework [START_REF] Kantorovich | On the translocation of masses[END_REF] by transforming it into a linear problem on a convex set, which later led to the proof of existence of a solution. The main idea was to look for optimal transport plans in the product space ℝ 𝑑 × ℝ 𝑑 instead of optimal transport maps. Since then, the problem has been widely generalized to general cost functions (beyond the Euclidean distance) and over general measures, and advances have been made to solve it efficiently [START_REF] Smith | Note on the optimal transportation of distributions[END_REF][START_REF] Cuesta-Albertos | A characterization for the solution of the monge-kantorovich mass transference problem[END_REF][START_REF] Villani | Topics in Optimal Transportation[END_REF][START_REF] Carlier | From knothe's transport to brenier's map and a continuation method for optimal transport[END_REF][START_REF] Igbida | A monge-kantorovich mass transport problem for a discrete distance[END_REF]. One of such extension, the Wasserstein metric, also commonly referred to as the earth mover's distance, is the measure of distance between two given probability distributions on a given metric space.

In the context of this thesis, we look at the optimal transportation problem of the form inf 𝜋∈𝛱(𝜇 0 ,𝜇 1 )

∬ 𝑐(𝑥, 𝑦)𝑑𝜋(𝑥, 𝑦), (5.5) between two probability distributions 𝜇 0 and 𝜇 1 where, 𝑐(𝑥, 𝑦) denotes the transport cost between 𝑥 and 𝑦 and 𝛱(𝜇 0 , 𝜇 1 ) denotes the set of all couplings of 𝜇 0 and 𝜇 1 , i.e., the collection of all measures with source and target marginals as 𝜇 0 and 𝜇 1 respectively (Fig. 5.2).

x c(x, y) µ 0 µ 1 y Figure 5.2: Optimal transportation problem. We search for a transport plan such that the total cost Eq. (5.5) is minimized. Note that, being a probability distribution, the total mass of both 𝜇 0 and 𝜇 1 is normalized to 1.

Entropic Regularization. While computing an optimal transport plan for a given linear cost functional can be achieved via linear programming, dealing with large variable counts negatively affects computation times. Recently, [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] proposed an entropic regularization of the optimal transport problem to smooth the conventional optimal transportation problem as

inf 𝜋∈𝛱 [∬ 𝑐(𝑥, 𝑦)𝜋(𝑥, 𝑦)𝑑𝑥𝑑𝑦 -𝛾𝐻(𝜋)] , (5.6)
where 𝐻(𝜋) refers to the entropy of the coupling, namely,

𝐻(𝜋)

def.

= -∬ 𝜋(𝑥, 𝑦)ln 𝜋(𝑥, 𝑦)𝑑𝑥𝑑𝑦, (5.7)

and 𝛾 controls the amount of regularization or smoothing. The cost function 𝑐 can be associated to a kernel 𝒦 𝛾 (𝑥, 𝑦) = exp(-𝑐(𝑥, 𝑦)/𝛾), so that we can compute the entropy-regularized transport plan from the smallest Kullback-Leibler (KL) divergence as

𝛾 [1 + min 𝜋∈𝛱 KL(𝜋|𝒦 𝛾 )] ,
(5.8)

where,

KL(𝜋|𝒦 𝛾 )

def.

= ∬ 𝜋 [ln 𝜋 𝒦 𝛾

-1] 𝑑𝑥𝑑𝑦.

(5.9)

Sinkhorn Iterations. This regularized optimal transport plan can be efficiently computed using Sinkhorn's matrix normalization algorithm [START_REF] Sinkhorn | A relationship between arbitrary positive matrices and doubly stochastic matrices[END_REF][START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF].

Following [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF]; [START_REF] Solomon | Convolutional wasserstein distances: Efficient optimal transportation on geometric domains[END_REF], given a matrix H discretizing the kernel 𝒦 𝛾 , the transport plan minimizing Eq. ( 5 as its diagonal, and 1 𝑘 is the vector of ones of size 𝑘. We find (v, w) via Sinkhorn iterations as detailed in Alg. 1, where ⊘ and ⊗ denote the entry-wise division and multiplication, respectively.

The Sinkhorn algorithm is known to converge at a linear rate [START_REF] Franklin | On the scaling of multidimensional matrices[END_REF][START_REF] Knight | The sinkhorn-knopp algorithm: Convergence and applications[END_REF]] and hence, is several orders of magnitude faster than solving a linear program. Note that the result, ⎧ { ⎨ { ⎩ 𝜋1 |n| = m, and (5.11) leading to the two distributions can be easily deduced using Eq. (5.10).

𝜋 ⊤ 1 |m| = n,
Algorithm 1 Sinkhorn iterations to solve for optimal transport

1: function Sinkhorn 2: v, w ← 1 3: while !converged do 4: v ← 1 |m| ⊘ H(n ⊗ w) 5: w ← 1 |n| ⊘ H ⊤ (m ⊗ v) 6: return 𝜋 = D m D v HD w D n ▷ Transport plan
Dykstra's Algorithm. In general, when convex sets of constraints are not affine subspaces, iterative Bregman projections [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF] (or Sinkhorn iterations, being a specific case) on these sets do not converge to the KL projection on the intersection. In such cases, the more general Dykstra's algorithm [START_REF] Dykstra | An algorithm for restricted least squares regression[END_REF] can be used as long as the constraints are convex sets. This algorithm is proven to converge [START_REF] Bauschke | Dykstra's algorithm with bregman projections: a convergence proof[END_REF]]. [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF] extended the Dykstra's algorithm to KL settings which can be used to solve such problems, for example, a relaxed transport plan. Starting with an initial solution 𝜋 (-1) = H and 𝑙 temporary variables 𝑞 (𝑘) = 1 -one corresponding to each of the 𝑙 convex constraint set 𝐶 𝑘 -we iteratively solve for

⎧ { ⎨ { ⎩
𝜋 (𝑛) = 𝑃 𝐾𝐿 𝐶 (𝑛 mod 𝑙) (𝜋 (𝑛-1) ⊗ 𝑞 (𝑛 mod 𝑙) ), 𝑞 (𝑛 mod 𝑙) = 𝑞 (𝑛 mod 𝑙) ⊗ 𝜋 (𝑛-1) ⊘ 𝜋 (𝑛) , (5.12) where, 𝑃 𝐾𝐿 𝐶 𝑘 denotes the projection according to the Kullback-Leibler divergence Eq. (5.9) over the constraint set 𝐶 𝑘 . The algorithm converges as 𝑛 → ∞.

Positioning

In our work, we adopt a mass transport formulation to seek inter-surface maps. However, instead of favoring softness and continuity as in most previous work, we promote sharpness of the map by penalizing the variance of the image of local neighborhoods. We show that our variance-based formulation favors continuity, conformality and sharpness of the map all at once. Although the problem becomes non-convex and hence substantially more challenging to minimize, we provide a formulation that turns this non-linear problem into a biconvex problem where geodesic variance is minimized in diffusion distance space, leading to a robust and efficient algorithm for inter-surface mapping. Additionally, we directly optimize a surface-to-surface map without resorting to an intermediate domain like many previous works. Our approach notably contrasts with Aflalo et al.

[2013] which sought conformal maps with no (or bounded) area distortion: we look instead for area-preserving maps (or maps with bounded area distortion) that are as conformal as possible.

Contributions. We propose an approach to find a dense map between two arbitrary surfaces.

We formulate this inter-surface mapping problem as a minimization of a Dirichlet energy for transport maps that we define in terms of the variances of the images of local neighborhoods. This energy penalizes stretching and favors conformal, point-to-point homeomorphisms. Although we cannot guarantee the existence of homeomorphisms, the idea behind our formulation, however, is to map neighborhoods to small geodesic neighborhoods, where small means low variancehence aspiring to form a homeomorphism.

We further provide an algorithm to compute such a variance-minimizing inter-surface map via alternating convex optimizations. A coarse to fine approach based on a hierarchy built in diffusion space [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps[END_REF] is developed to both improve efficiency and reduce the risk of getting trapped in local minima. In addition, we modify and extend the Sinkhorn iteration algorithm [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF][START_REF] Solomon | Convolutional wasserstein distances: Efficient optimal transportation on geometric domains[END_REF] to handle low regularization and unbalanced optimal transport [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF], where mass preservation can also be relaxed to enable further minimization of the map distortion.

Related Work

Over the last decade, there has been significant amount of work on conformal mapping. However, apart from some recent advancements to generate area-preserving maps, as of now, the problem of area-relaxed maps which are as conformal as possible remaines algorithmically infeasible. We discuss next these algorithms, along with other approaches in order to broaden the discussion to the general topic.

Area-preserving Mapping. [START_REF] Angenent | On Area Preserving Mappings of Minimal Distortion[END_REF] proposed an algorithm for area-preserving map with minimal distortion in the context of 3𝐷 medical visualization, for visualizing functional magnetic resonance imaging data of neural activity in the brain. They define a variant of the Dirichlet energy and minimize it over the space of area-preserving diffeomorphism via steepest descent. Monge-Kantorovich theory of optimal mass transportation has also been used to build area-preserving mapping due to its mass preservation property [START_REF] Zhu | Area-Preserving Mappings for the Visualization of Medical Structures[END_REF][START_REF] Haker | Optimal mass transport for registration and warping[END_REF]. Recently, [START_REF] Su | Area preserving brain mapping[END_REF] proposed an alternate method based on Monge-Brenier theory which leads to a convex energy minimization and power Voronoi diagram construction.

The algorithm generates a unique solution which is invariant under isometric transformations and with an improved complexity, i.e. with a linear instead of quadratic number of variables.

Mapping based on area-distortion minimization has also been studied [START_REF] Desbrun | Intrinsic parameterizations of surface meshes[END_REF] leading to a quartic polynomial which can be efficiently solved using a simple root finder.

Apart from them, other previous approaches mainly target conformal maps and can be classified 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 find the mappings. We review the most relevant approaches next, and refer the reader to van Kaick et al. [2011] for a comprehensive survey on shape correspondence.

Embedding. When dealing with near isometric maps, embedding of the input manifolds into some target (Euclidean) space [START_REF] Gallot | Embedding Riemannian manifolds by their heat kernel[END_REF][START_REF] Bronstein | Generalized multidimensional scaling: A framework for isometry-invariant partial surface matching[END_REF][START_REF] Shtern | Spectral gradient fields embedding for nonrigid shape matching[END_REF]] is a practical way to simplify matching since corresponding points of different isometric shapes are mapped to nearby points in the target space. An extension of this family of approaches use as embedding space a common template mesh. The base domain is often chosen to be a coarse mesh obtained through mesh decimation [START_REF] Schreiner | Inter-surface mapping[END_REF], or by defining overlapping domains bounded by consistent paths connecting feature points [START_REF] Kraevoy | Cross-parameterization and compatible remeshing of 3D models[END_REF]]. 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. Mesh Parameterization. Mesh parameterization was originally introduced in computer graphics by [START_REF] Bennis | Piecewise surface flattening for non-distorted texture mapping[END_REF] for texture mapping and has been extensively used for geometry mesh processing. Since it is impossible to avoid both angle and area distortion for non-zero Gaussian curvature surfaces, [START_REF] Degener | An adaptable surface parameterization method[END_REF] proposed parameterization based on global area and angle distortion. Other approaches based on minimizing the distortion of different intrinsic measures [START_REF] Desbrun | Intrinsic parameterizations of surface meshes[END_REF], conformal mapping via least squares approximation [START_REF] Lévy | Least squares conformal maps for automatic texture atlas generation[END_REF][START_REF] Liu | A local/global approach to mesh parameterization[END_REF] and globaly conformal parameterization based on gradient fields [START_REF] Gu | Global conformal surface parameterization[END_REF] have also been studied.

Recent approaches construct maps with bounded distortion [START_REF] Lipman | Bounded distortion mapping spaces for triangular meshes[END_REF], use cone manifolds [START_REF] Myles | Controlled-distortion constrained global parametrization[END_REF] or orbifolds [Aigerman and Lipman, 2015], and relax the bijective condition of the parameterization to local injectivity [START_REF] Aigerman | Lifted bijections for low distortion surface mappings[END_REF]. A recent approach contributed by [START_REF] Aflalo | Conformal mapping with as uniform as possible conformal factor[END_REF] seeks angle-preserving maps with as uniform as possible conformal factors, i.e., as area-preserving as possible conformal maps. 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 intersurface maps. Mapping to a common domain also often requires designing a cut-graph which can increase distortion as well (see Aigerman et al. [2015] for a recent improvement).

Soft Maps. While point-to-point mapping between surfaces has been favored for the past few decades, recent notions of soft maps [START_REF] Solomon | Soft maps between surfaces[END_REF] and functional maps [START_REF] Ovsjanikov | Functional maps: a flexible 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, [START_REF] Solomon | Soft maps between surfaces[END_REF] generated 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 distributions on the other surface, where "nearby" is to be understood with respect to the geodesic earth mover distance. In addition, softness of the map is achieved by 𝐿 2 penalization, spreading the probabilities in the mapping matrix describing the mass transport plan: this quantity is convex and amenable to efficient optimization. The manner in which a soft map transports one distribution to another was also pointed as relevant to better understand and control map continuity [START_REF] Solomon | Dirichlet energy for analysis and synthesis of soft maps[END_REF]. A recent work, [START_REF] Solomon | Entropic metric alignment for correspondence problems[END_REF], further demonstrates the power of entropic regularization to render challenging problems computationally tractable, but acknowledges that limiting the induced oversmoothing is still an open problem. While soft maps are an elegant and powerful paradigm, sharper maps, from which point-to-point correspondences can be inferred, remain a necessity for most geometry processing algorithms.

Finally, a recent trend is to adopt a machine learning approach in order to alleviate the need to define a prior on the type of deformation [START_REF] Litman | Learning spectral descriptors for deformable shape correspondence[END_REF][START_REF] Rodolà | Dense non-rigid shape correspondence using random forests[END_REF].

Chapter 6 Algorithm

Optimal Transport Approach

Given two input pointsets 𝑋 = {𝑥 𝑖 } and 𝑌 = {𝑦 𝑗 } sampling two surfaces, we aim to construct a map between 𝑋 and 𝑌 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 input pointsets for which the transport plan minimizes the variance of the image of each small neighborhood, where the variance of the plan is efficiently measured in diffusion space.

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 specifically, we discretize the normalized area measures of the two surfaces as weighted sums of Dirac measures centered at input points, i.e., 𝜇 = ∑ 𝑚 𝑖 𝛿 𝑥 𝑖 and 𝜈 = ∑ 𝑛 𝑗 𝛿 𝑦 𝑗 , where the masses 𝑚 𝑖 and 𝑛 𝑗 are normalized so that each surface has total mass ∑ 𝑖 𝑚 𝑖 = ∑ 𝑗 𝑛 𝑗 = 1. The mass of each point can be assigned based on its estimated local Voronoi area of the surface it samples to account for irregular sampling, or simply to a uniform constant if such an estimate is unavailable.

Maps as transport plans. Given two distributions 𝜇 and 𝜈, a transport plan 𝜋 between them is by definition a matrix of size |𝑋| × |𝑌 | whose marginals equal 𝜇 and 𝜈. Each entry 𝜋 𝑖𝑗 represents the amount of mass transported from 𝑥 𝑖 ∈ 𝑋 to 𝑦 𝑗 ∈ 𝑌 . A transport plan 𝜋 has two associated transfer operators 𝜋 𝑋 and 𝜋 𝑌 defined as

𝜋 𝑋 (𝛿 𝑥 𝑖 ) = 1 ∑ 𝑗 𝜋 𝑖𝑗 ∑ 𝑗 𝜋 𝑖𝑗 𝛿 𝑦 𝑗 , 𝜋 𝑌 (𝛿 𝑦 𝑗 ) = 1 ∑ 𝑖 𝜋 𝑖𝑗 ∑ 𝑖 𝜋 𝑖𝑗 𝛿 𝑥 𝑖 .
Note that the denominators are respectively equal to 𝑚 𝑖 and 𝑛 𝑗 . These two linear operators describe how mass distributions on 𝑋 are mapped to mass distributions on 𝑌 with the same total mass, and vice-versa. Therefore, they contain essentially the same information as the transport plan.

Neighborhoods as weighting functions. For a given point 𝑥 on a surface, we denote by 𝑊 𝑥 a weighting function centered at 𝑥. A natural choice for 𝑊 𝑥 is a truncated Gaussian function, but any near-isotropic bump function that decreases as a function of the geodesic distance from 𝑥 may do. Such weighting functions will be used over the two input surfaces as test functions to measure the variance of an inter-surface map at different locations since their finite supports define geodesic neighborhoods.

Regularity measure for transport plans

Our approach consists in minimizing the following variance-based functional over all possible plans 𝜋 between 𝜇 and 𝜈:

𝐸 (𝜋) = ∫ 𝑋 var 𝜋 𝑋 ( 𝑊 𝑥 𝜇 mass(𝑊 𝑥 𝜇)
) 𝑑𝜇(𝑥) + ∫ 𝑌 var 𝜋 𝑌 ( 𝑊 𝑦 𝜈 mass(𝑊 𝑦 𝜈)

) 𝑑𝜈(𝑦), (6.1)

• • • • • • X W x µ µ π X Y x x • • • • • • • • • • • •
that is, we integrate the variances of images by the transfer operators of area measures modulated by local weighting functions (see inset). Note that since we consider our surfaces as discrete measures, these integrals reduce to weighted sums over data points.

Soft vs. sharp maps. The reason for choosing this cost functional is two-fold. First, penalizing the variance of the image of the map favors transport plans corresponding to actual pointto-point maps rather than soft (diffused) maps: any point mass that is mapped to a spreadout distribution incurs a significant cost. Second, these variance-minimizing maps tend to be bijective, since we also account for the inverse map in the cost evaluation. Note that bi-continuity may also be expected for the minimizers of our cost as we penalize variances of images of small neighborhoods defined by our weighting functions.

Link to Dirichlet energy. We show in App. A that for infinitesimal weighting functions with isotropic covariance matrices tending to zero, our cost converges to the sum of the Dirichlet energies of the map and its inverse. Hence, our cost may be viewed as a symmetric Dirichlet energy for transport plans. While generalization of the Dirichlet energy to transport plans have already 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 penalizing soft maps, a property that is essential for our purpose. Moreover, we show next that our particular variance-based formulation is amenable to the use of efficient optimization methods.

Minimization via alternating convex optimizations

While the formulation above formalizes the notion of a good map between two surfaces, it involves a non-linear and in general non-convex functional. As such, it may seem difficult to minimize efficiently and robustly. We can, however, introduce auxiliary variables within the variance terms so that finding a variance-minimizing map amounts to solving a biconvex problem, which we achieve through alternating convex minimization. Improving robustness. The regularization coefficient 𝛾 used to compute the optimal transport is recommended to be set such that 𝛾 > max{coef (𝜋 𝑖𝑗 )}/200 [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] in practice, to prevent numerical blowups due to limited floating-point dynamic range. However, this choice of coefficient usually creates quite a diffuse transport plan, conflicting with our objective to construct varianceminimizing plans. We bypass this numerical limitation by storing the coefficients of the kernel matrix in log scale and compute the Sinkhorn iterations in log scale directly. This implies that we now store 𝐻 𝑖𝑗 = -coef(𝜋 𝑖𝑗 )/𝛾, and perform the inner products 𝑎 ⊤ 𝑏 in Steps 4-6 of Alg. 1 (involving columns or rows of 𝐻) by first subtracting the maximum element 𝑀 = max 𝑖 {𝑎[𝑖] + 𝑏[𝑖]} from each element, and evaluating

𝑎 ⊤ 𝑏 ≔ 𝑀 + log ∑ 𝑗 exp(𝑎[𝑗] + 𝑏[𝑗] -𝑀 ).
Working in log scale along with the subtraction of the maximum element drastically reduces the negative effects of a limited floating-point dynamic range, allowing us to safely set to 𝛾 = max{coef (𝜋 𝑖𝑗 )}/10 4 , a 50-fold decrease in regularization.

As 𝛾 → 0, Sinkhorn's algorithm can require many iterations to converge. In practice we thus perform only 100 iterations between 𝜂 updates. Once the 𝜂 updates are small enough, we perform annealing by first increasing 𝛾 by 50%, followed by 5 extra 𝜂 updates with 1,000 Sinkhorn iterations each. Parameter 𝛾 is then reverted back to its original value. This scheduling experimentally improves the convergence of Sinkhorn iterations, and we observe that mass preserving constraints are nearly satisfied upon completion (Fig. 6.3). In order to further improve efficiency, we also use the previous values of 𝑣 and 𝑤 instead of initializing them to 1 in each successive relaxation iterations.

Barycenters in diffusion space

Computing the geodesic barycenters 𝜂 described in Sec. 6.1.3 is not tractable (Alg. 2 Step 3).

Observe however, that we do not need to compute them on the surfaces: we only require an estimate of the (geodesic) variance. We thus resort to a surrogate version through diffusion distances using the first few eigenvectors of the Laplace-Beltrami operator precomputed on the input datasets. Since diffusion distances behave like Euclidean distances in diffusion space, we can explicitly solve for 𝜂 in diffuson space by minimizing Eq. ( 6

Coarse to Fine Solver

Despite our use of a fast method to compute optimal transport plans, the number of variables |𝑋|+ |𝑌 | can be quite large, slowing down convergence of our alternating minimization.

A common approach to both accelerate convergence and prevent the alternating minimization from getting stuck in local minima of the energy is to use a coarse to fine strategy. By starting the map optimization on a coarse sampling of the surfaces (and thus, a small number of variables) and progressively increasing resolution to refine the map, we significantly reduce the total computation time required to find a variance-minimizing map. Our approach consists in first establishing a hierarchical sampling where the cardinality of a neighbourhood, defined as 𝒩 𝑥 = {𝑥 𝑖 ∈ 𝑋 ∶ 𝑊 𝑥 (𝑥 𝑖 ) ≥ 𝜀} as mentioned in Sec. 6.2.1, remains constant as we go to finer levels in the hierarchy, so that the number of variables needed at each level remains small.

Hierarchical Sampling. We leverage the smoothing properties of the time scale 𝑡 in diffusion distance to construct a hierarchy of dyadic levels, one for each time scale 𝑡 ℓ = 2 -ℓ 𝑡 0 . Since our weighting function is defined according to Eq. ( 6.4), reducing the time scale shrinks the support of 𝑊 𝑥 to a smaller neighborhood. Starting from a large diffusion time scale 𝑡 0 , we build a sampling 𝒮 0 of the input dataset via farthest point insertion until the farthest point 𝑝 satisfies 𝑊 𝑥 (𝑝) > 𝜏 ∀𝑥 ∈ 𝒮 0 for some 𝜏 > 𝜀 (we take 𝜏 = 0.9 in our experiments), indicating a dense enough sampling for that level. In practice, this approach ensures that the area contributing to the neighborhood of any point 𝑥 in the level (namely, the set of points 𝑧 satisfying 𝜀 < 𝑊 𝑥 (𝑧) < 𝜏 )

is covered by at most a constant number of Voronoi cells, leading to a constant neighborhood cardinality. As the time scale decreases, the levels 𝒮 ℓ become more dense until, eventually, no additional points can be inserted: since only a limited number of eigenvectors are used to compute diffusion distances, our weighting functions will stop being "discerning" past a certain scale. Further levels can be computed using Euclidean distances instead: locally, geodesic and Euclidean distances behave the same at such scales.

In these Euclidean stages of the hierarchy, we thus revert to a plain (unnormalized) Gaussian function for the weighting function: we use

𝑊 𝑥 (𝑥 𝑖 ) = exp (- (𝑥 -𝑥 𝑖 ) 2 2𝜎 2 ) , (6.6)
where 𝜎 is initialized using the average radius 𝑟 avg of the neighborhoods in the last diffusion stage such that 𝜀 = exp(-𝑟 2 avg /2𝜎 2 ). To avoid any ill-effects when 𝜎 tends to input data discretization limits, we always have a lower bound on the cardinality of the neighbourhood. Furthermore, for each new hierarchy level ℓ in Euclidean space, 𝜎 is directly updated to 𝜎/2. Note that the Euclidean hierarchy stages further improve results: using only the first few eigenvectors of the Laplace-Beltrami operator for diffusion distances would fail to capture local surface details corresponding to high frequencies.

Local Variable Instantiation. At the coarsest level of the hierarchy, we can afford to solve the full optimal transportation problem. 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 (𝑥 ′ , 𝑦 ′ ) such that there exists a pair (𝑥, 𝑦) in the current support with 𝑥 ′ ∈ 𝒩 𝑥 and 𝑦 ′ ∈ 𝒩 𝑦 . For the next Sinkhorn iterations, we limit our transport plan to this enlarged support, which limits the number of variables to instantiate (Fig. 6.4). 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 during center relocation-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.

X Y We plot the relative decrease of the total cost (6.3) (normalized per level) at every alternating minimization iteration while optimizing the map between two dogs (Fig. 7.2). Blue depicts the diffusion distance phase (10 levels) and red depicts the Euclidean phase (6 levels). The dashed lines represent a change of hierarchy level.

• • • • • • • • • • • • x W x ≥ ε N x ∪N π X (N x ) π X
Sinkhorn iterations are then computed efficiently by first precomputing ∑ 𝑖 m ⊗ v[𝑖] (resp.

∑ 𝑗 n ⊗ w[𝑗]) in Alg. 1 and then accordingly adjusting the value of H(m ⊗ v) computed using the sparse matrix H. As we keep the cardinality of neighborhoods constant during the construction of the hierarchy, both storage and computation involved in the optimal transport solver become in practice linear in the number of samples. Fig. 6.5 shows the result obtained via hierarchical sampling, while Fig. 6.6 plots the total cost (normalized per level) over iterations while solving for a map. The total cost was normalized per level as it is hard to compare the cost between two levels as sample count, diffusion time and neighborhood evolve. Note that while each level corresponds to a decrease of only 4% of the total energy, each reduction has exponential impact across the hierarchy.

Flips and twists

We achieve dramatic speedup by approximating geodesic distances via diffusion distances, reducing the neighborhood over the hierarchy, and limiting the instantiation of variables to a small neighborhood. However, these algorithmic and numerical shortcuts lead to a practical drawback on surfaces with thin protrusions: a few, localized discontinuities of the map may appear in the form of twists (180 degree rotations of the map) and/or flips (mirrored versions of the map), see Fig. 6.7. While increasing the number of eigenvectors used for diffusion distances and increasing mesh density prevent these issues, it comes at the cost of increased running times. (left, cactus (top) to itself (bottom) built by intentionally placing bad constraints), ironing proceeds through front propagation to remove all artifacts. For each strip, blue indicates the visited parts while the non visited parts are in pink. In case of multiple fronts, each front is considered independently. Note that the main idea behind whole approach is to slightly distort the neighborhood and favor already visited parts. Using high value for 𝜔 (close to 1) cannot unfold the mapping, while low values (close to 0) seriously pull 𝜂's towards already-visited parts.

Instead, we found more efficient to run the algorithm as is, and simply fix these small and local map artifacts at the end through what can be best described as "ironing": we correct problem areas by increasing the weights of correct neighborhoods during the 𝜂 update, effectively overpowering these twists and flips and bringing the solution out of the local minimum at which it was. More specifically, starting from a user defined point or a patch on 𝑋 with a correct map in its neighborhood, 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 𝑌 by geodesically thickening the image of the front on 𝑋 proportional to the mass of the visited and non-visited parts. We now consider these two strips as two new (small) datasets, and solve for a variance minimizing transport plan satisfying Eq. (6.1) while keeping 𝜂 for the visited part constant. The update 𝜂 in non-visited part is achieved using an asymmetric weighting function:

𝑊 * 𝑥 (𝑦) = ⎧ { ⎨ { ⎩ 𝑊 𝑥 (𝑦), 𝑦 ∈ visited, 𝜔𝑊 𝑥 (𝑦), otherwise, (6.7)
where, 𝜔 set to 0.5 in all our experiments. Once these 𝜂 stabilize, the front is moved outward and the process is repeated until covering the whole surface, eventually removing flips/twists as shown in Fig. 6.8 and the video in the additional material. This final ironing process is very efficient, and gracefully removes all remaining artifacts of the map and its inverse, even if we intentionally flip the map significantly as in Fig. 6.8.

Relaxing Area Preservation

When dealing with non isometric shapes, enforcing strict mass preservation leads to mappings with anisotropic differentials. Depending on the shapes, this anisotropy may become significant.

Relaxing the mass preservation constraint helps remove this issue.

Sinkhorn Iterations for Relaxed Optimal Transport

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(𝜋|𝒦 𝛾 )
where 𝛱 denotes the set of matrices satisfying:

⎧ { { { ⎨ { { { ⎩ 𝜋 𝑖𝑗 ≥ 0, ∑ 𝑖,𝑗 𝜋 𝑖𝑗 = 1, 𝛼 𝑖 𝑚 𝑖 ≤ ∑ 𝑗 𝜋 𝑖𝑗 ≤ 𝛽 𝑖 𝑚 𝑖 , 𝛼 ′ 𝑗 𝑛 𝑗 ≤ ∑ 𝑖 𝜋 𝑖𝑗 ≤ 𝛽 ′ 𝑗 𝑛 𝑗 , (6.8)
where, 0 ≤ 𝛼

( ′ ) 𝑘 ≤ 1 ≤ 𝛽 ( ′ )
𝑘 are the prescribed lower and upper bounds on the allowed mass distortion. In practice we use 𝛼 𝑘 = 𝛽 𝑘 = 1 and 𝛼 ′ 𝑘 = 1/𝛽 ′ 𝑘 , which is a relaxation of the set of diffeomorphisms whose Jacobian and inverse Jacobian are bounded by a constant. We may write the convex set of allowed solutions as the intersection of two convex sets, namely the set 𝛱 1 of matrices satisfying the first, second and third condition above, and the set 𝛱 2 of matrices satisfying the first, second and fourth condition. Optimizing over their intersection can be done using Dykstra's algorithm as adapted to KL divergence in [START_REF] Bauschke | 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. B.

It turns out that the KL projections on 𝛱 1 and 𝛱 2 respectively boil down to row and column rescaling; because of these particular 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. A pseudocode for the algorithm is given below (Algo. 3). 𝑟 ← solution to Global_Mass_v(lb, ub)

Algorithm 3 Relaxed Sinkhorn iterations

6: v ← Update_Vector_Global_Mass(𝑟, lb, ub) 7: lb ← 𝛼 ′ ⊘ H ⊤ (m ⊗ v) 8: ub ← 𝛽 ′ ⊘ H ⊤ (m ⊗ v) 9:
𝑟 ← solution to Global_Mass_w(lb, ub) which must be chosen so that the resulting transport plan indeed has global mass equal to 1.

Fortunately, the global mass is piecewise linear and non decreasing as a function of 𝑟, with a number of nodes equal to twice the number of points in each dataset. Hence, solving these equations can be done by sorting these nodes, in time linear in the dataset size up to a logarithmic factor. Also, it may happen that the equations have no solution, in which case we set 𝑟 to a large magnitude, with the appropriate sign. x(𝑟) ←Update_Vector_Global_Mass(𝑟, lb, ub)

14:

𝑟 ← solution to m ⊤ D v HD x n = 1

Three-fold Alternating Minimizations

Our alternating minimization algorithm cannot be directly extended to the case where we relax the mass preservation constraint. Indeed, the denominator in our definition of transfer operators makes the transport plan optimization step more difficult, and seemingly not reducible to a convex problem. To solve this issue, we introduce additional variables 𝑚 ′ 𝑖 and 𝑛 ′ 𝑗 , which ideally would correspond respectively to ∑ 𝑗 𝜋 𝑖𝑗 and ∑ 𝑖 𝜋 𝑖𝑗 . Rather than enforcing strict equality, we allow some slack and set the constraints

⎧ { ⎨ { ⎩ 𝐶 -1 𝑚 ′ 𝑖 ≤ ∑ 𝑗 𝜋 𝑖𝑗 ≤ 𝐶𝑚 ′ 𝑖 , 𝐶 -1 𝑛 ′ 𝑗 ≤ ∑ 𝑖 𝜋 𝑖𝑗 ≤ 𝐶𝑛 ′ 𝑗 .
(6.9)

for some constant 𝐶 close to 1 (set to 2 and gradually decreased to 1.1 in our experiments).

Additionally, we impose that 𝐷 -1 𝑚 𝑖 ≤ 𝑚 ′ 𝑖 ≤ 𝐷𝑚 𝑖 and 𝐷 -1 𝑛 𝑗 ≤ 𝑛 ′ 𝑗 ≤ 𝐷𝑛 𝑗 , where 𝐷 is a userchosen bound on the allowed mass distortion. We can then form the following modified transfer 𝑌 . We then proceed by alternately optimizing over the three variables 𝜋, (𝑚 ′ , 𝑛 ′ ), and 𝜂. To solve for 𝜋 𝑖𝑗 , we use the relaxed Sinkhorn iterations described above. Solving for 𝑚 ′ and 𝑛 ′ requires the optimization of a positive linear combination of inverses of their coefficients over a hypercube. This convex optimization problem is solved in near linear time using an algorithm similar to the one used in the previous paragraph for enforcing global mass preservation. Finally, the update of 𝜂 remains a simple weighted average as before. Algorithm 5 provides a pseudocode for the overall algorithm. A result of this approach is depicted in Fig. 6.9. Chapter 7

Algorithm 5 Overall Algorithm

Results

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 heat diffusion operator, the Eigen library for linear algebra operations, CGAL Initialization. When the input models are properly aligned in space (up to a translation), minimizing only the Wasserstein-2 transport distance is sufficient to initialize 𝜂 for the first level of the hierarchy. Unless otherwise specified a very small additional Wasserstein-2 distance term (coefficient 1𝑒 -5) is then added to our cost function in order to help tunneling out of local minima. When the input models are not aligned, we initialize with 𝜂 = 0 and use no Wasserstein-2 distance term. A few (two to four) user-defined point constraints are often sufficient to obtain meaningful maps in these cases.

Experiments

Comparison against most of the methods which target conformality will show that they are more conformal while our method is able to produce area-relaxed maps which are as conformal as possible. Although such a direct comparison might not be meaningful in most of the cases, nevertheless, we evaluate our algorithm with the state of the art algorithms as follows: Isometric Case. We begin our experiments on isometric models from the TOSCA dataset [START_REF] Bronstein | Numerical Geometry of Non-Rigid Shapes[END_REF]. Fig. 7.1 shows a centaur in three isometric poses. All maps are obtained without any user interaction in these cases. Non-isometric Case. We also ran our algorithm on models that are increasingly non isometric. We challenge our algorithm in Fig. 7.5 by mapping between two highly non-isometric models of genus 1, a torus and a kitten, with no area relaxation. As expected, the map gets highly non-conformal in this example, but the resulting map is nice nonetheless. Fig. 7.6 depicts a mapping between two non-isometric models, a man and a gorilla, with a maximum area relaxation factor set to 10 (i.e., area distortion allowed from 0.1 to 10). Top left: a sphere with a small handle is mapped to a sphere. The mapping is discontinuous, but the discontinuity is local. Top right: same with a mediumsize handle. Bottom: a regular torus is mapped to a sphere. The discontinuity is no longer local, but the map is continuous over half of the sphere.

Four constraints have been used (two hands and two feet). The resulting map is nice despite large differences in the proportion of the limbs. Fig. 7.7 depicts a mapping between a man onto Homer, with no area relaxation. Four constraints are specified: two on the hands and two on the feet. As expected, the distortion is very large due to vast differences in shape, especially on the feet and head. However, semantic correspondences are well captured almost everywhere.

We further challenge our algorithm by mapping a sphere onto three genus-1 surfaces (Fig. 7.8), in order to evaluate the robustness to topological noise and analyze the mappings when the topology differs. We thus modeled (1) a sphere with a small handle, (2) a sphere with Chapter 8

Conclusion

Recent advances in improving the efficiency of optimal transport solvers is providing an opportunity to derive inter-surface mapping through direct energy minimization. Our paper proposes a first variational model of the sort, where transport plans are optimized to be locally "sharp" so as to yield homeomorphisms. We showed that our approach finds a least-stretched map with exact area preservation or bounded area distortion, and demonstrated its performance on a range of challenging models, with possibly different or non-trivial topology.

Limitations.

Although our algorithm provides a novel formulation, there is no direct practical application where it can be applied. It cannot be used directly for shape matching as our cost function does not target semantically correct map. For other applications involving finding a balance between conformality and area preservation, several methods which outperform ours (in terms of running times) already exist.

Furthermore, the algorithm is not robust to thin features and depends on ironing (Section 6.2.5) for the continuity of the map. However, it does not guarantee continuity in all the cases and may fail in practice. Note that the strip in ironing is generated by isotropically extending the front in the both directions proportional to the area. Hence, it might not work when area-relaxation is involved. Moreover, the ironing step mainly works only in cases where the topology of the front or the strip on both meshes is identical: appearance of any significant enough extra feature (thin protrusions) in one mesh will break it.

The numerics of our approach works best for uniform isotropic meshes since 𝜎 used for weighting function 𝑊 is set proportionally to the vertex density, thus ensuring that the neighborhood cardinality remains nearly constant. As a consequence, we locally refine meshes with anisotropic vertex distribution as a preprocessing stage to ensure robustness. As relaxing the area preservation constraints through unbalanced mass transport increases the number of variables, we cannot allow too large of a mass relaxation. In addition, if the mass is relaxed to the extent that a neighborhood maps to a single point, our algorithm is no longer 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.

Future Work. Along with Sinkhorn iterations, other compute-intensive parts can be transferred to GPU. This will significantly improve the running time of our algorithm. An obvious direction will be to incorporate shape signatures in our algorithm. We tried Heat Kernel Signature [START_REF] Sun | A concise and provably informative multi-scale signature based on heat diffusion[END_REF], but it does not significantly improve the results because it is also based on the surface metric which is already captured by our objective function. One weakness of our approach is that convergence may be slow, especially on areas with continuous symmetries such as cylinders, and on models with a wide range of feature size. This is also probably relating to the convergence properties of Sinkhorn iterations for different sparse matrices. It will be interesting to better analyze and understand this behavior. One direction left for future work is to use a Newton iteration within the nested loops to accelerate convergence of the centroids. Recall that the centroids are just an approximation in Euclidean stages and do not lie on the surface. Another stimulating direction is implementing the concept of variable neighborhood size. This will allow us to directly proceed with our algorithm on anisotropic meshes (to some extent) without the need for isotropic refinement.

As future work we wish to explore the enduring challenge of partial matching which is currently not achievable through area relaxation. Exploring other geometric priors or bounds is also a stimulating avenue that we plan to tackle.

Figure 1 . 1 :

 11 Figure 1.1: Unknot and knot.A torus and a trefoil knot are homeomorphic. Note that there is no need for the existence of a continuous deformation for the two spaces to be homeomorphic. Furthermore, they are not isotopic as the knots will have to pass through each other to form an unknot and it will not be possible to maintain a homeomorphism throughout the deformation.

Figure 1 . 3 :

 13 Figure 1.3: Size vs Quality tradeoff. As we gradually decrease the number of vertices, the quality of mesh (features) degrades leading to coarser approximation of original mesh.

Figure 1 . 7 :

 17 Figure1.7: Self-Intersection. Intersection with the nearby mesh can be caused due to vertex removal or edge collapse operator in regions where the feature size is relatively smaller as compared to the desired Hausdorff distance error in the approximation.

Figure

  Figure1.8: Flips. An edge collapse may not always result in a valid triangulation. For a triangulation to be valid, the target location of vertex must lie in the visibility kernel (green shaded region) of the one-ring. Left: a valid triangulation in 2D. Right: A self-intersection can also occur in the form of flips due to invalid triangulation. The same can also be extended to higher dimensions.

Figure 1 . 9 :

 19 Figure 1.9: Reconstruction using CRUST algorithm. From left to right: input point set; its Voronoi diagram along with Voronoi vertices and; the Delaunay triangulation of the union of input points and Voronoi vertices. The crust edges connecting the input points are depicted in black. In general, all Delaunay based algorithms attempt to identify or extract a correct subset from Delaunay triangulation for the desired curve or surface through various filtering methods or geometric heuristics.

  used implicit surfaces and combined this with graph-cuts to get a watertight model. Similarly, Shalom et al. [2010] used generalized cones with a data point at their apex to indicate exterior regions. Alliez et al. [2007] created an implicit surface using the anisotropy of the Voronoi cells to estimate the normals of the surface (Figure 1.10). However, such implicit surface methods have problems in reconstructing sharp features and are in general, not suited to deal with massive data sets.

Figure 1 .

 1 Figure 1.10: Voronoi-based variational reconstruction of unoriented point sets. From left to right: input point set; its Voronoi diagram; covariance matrices of the cells; Steiner points added through Delaunay refinement (isotropic tensors are assigned to Steiner points); piecewise linear function f (solution of a generalized eigenvalue problem).(Image taken from[START_REF] Alliez | Voronoi-based variational reconstruction of unoriented point sets[END_REF])

Figure 1 .

 1 Figure 1.12: Adaptive moving least squares surface. From left to right: To handle noisy datasets, outward normals are estimated from big Delaunay balls at a subset of sample points and then the points are projected using these normals.(Image taken from[START_REF] Dey | An adaptive mls surface for reconstruction with guarantees[END_REF])

Figure 1 .

 1 Figure 1.13: Hole filling. From left to right: Given an input point cloud, an unsigned distance function by volumetric diffusion is computed. The output surface is supposed to lie in the voxel crust 𝑉 𝑐𝑟𝑢𝑠𝑡 between the outer and the inner boundary. A spatial graph structure is then embedded within the voxel grid, with small edge weights for high confidence voxels and vice versa. The boundaries are connected to a sink and a source node, respectively. The output surface is finally computed through the min-cut of this graph. Note that the algorithm does not require any information about the local surface orientation.(Image taken from[START_REF] Hornung | Robust reconstruction of watertight 3d models from non-uniformly sampled point clouds without normal information[END_REF])

Figure 1 .

 1 Figure 1.14: Robustness to noise and outliers. Input point set and the reconstructed curve are depicted in black and red respectively.(Image taken from[START_REF] Giraudot | Noise-adaptive shape reconstruction from raw point sets[END_REF])

Figure 1 .

 1 Figure 1.15: Hole-fixing by morphological operations. From left to right: The boundary cells are dilated (magenta) into the empty cells (green) and the outside component (blue) is determined.The outside component is then dilated back into the already dilated cells (magenta). This results in a clean separation of outside (blue) and inside cells (white) from which a surface topology of the output mesh can be deduced.(Image taken from[START_REF] Bischoff | Automatic restoration of polygon models[END_REF])

Figure 2 . 1 :

 21 Figure 2.1: Overview of our algorithm. Top: input tolerance 𝛺, sampling of 𝜕𝛺, mesh refinement by inserting a subset of the sample points, and topology condition met. Samples that are well classified are depicted in green, and in red otherwise. The boundary of the simplicial tolerance volume 𝜕𝛤 is depicted with blue edges. Bottom: simplification of 𝜕𝛤 , mutual tessellation of zero-set, simplification of zero-set, and final output.

Figure 2 .

 2 Figure 2.1 depicts the three main steps of our approach: First, the initialization step generates a dense point sample 𝒮 on the boundary of the tolerance volume 𝜕𝛺. Second, we proceed coarseto-fine through refinement of a 3D Delaunay triangulation by inserting one sample of 𝒮 at a time, and while maintaining a piecewise-linear function interpolated on the triangulation. The function value at the triangulation vertices is set in accordance to the index of each boundary component 𝜕𝛺 𝑖 (+1 or -1). The term zero-set refers to the isosurface where the interpolated function evaluates to zero. Refinement is performed until the zero-set is entirely contained into 𝛺 and

Figure 2 . 2 :

 22 Figure 2.2: Classification of 𝒮. The black solid edges depict the zero-set 𝒵 of 𝑓. A sample classified as good is depicted in green, and in red otherwise.

Figure 2 . 3 :

 23 Figure 2.3: Refinement of 𝒯. Top: initial triangulation and one Steiner point inserted. Middle: more Steiner points inserted. Bottom: more Steiner points inserted, and complete classification of samples.The zero-set is depicted with black solid edges. Samples classified as good are depicted in green, and in red otherwise. A Steiner point to be inserted at the next iteration is depicted in red. Upon termination the edges of 𝜕𝛤 are depicted in blue.

Figure 2 . 4 :

 24 Figure 2.4: Misoriented element. Left: The edges of 𝒵 are depicted with solid black lines. The zeroset of △𝐴𝐵𝐶 (red) has an incorrect normal. Right: The piecewise-linear function defined on △𝐴𝐵𝐶 should classify well the samples of 𝒮 (on both 𝜕𝛺 𝑖 forming △𝐴𝐵𝐶) which are nearest (orange) to the vertices of a shrunk triangle (green).

Figure 2

 2 Figure 2.5: Height of a tetrahedron contributing to 𝒵. The height is defined as the distance between the supporting primitive of the maximum dimension simplices formed by the tetrahedron vertices with common labels. Left: distance between point 𝐴 and supporting plane of △𝐵𝐶𝐷. Right: distance between supporting lines of edges 𝐴𝐵 and 𝐶𝐷.

Figure 2 . 6 :

 26 Figure 2.6: Link condition in 2D. Left: the edge 𝐴𝐵 is collapsible as 𝐿𝑘(𝐴) ∩ 𝐿𝑘(𝐵) = 𝐿𝑘(𝐴𝐵). Right: the edge 𝐴𝐵 is not collapsible as 𝐿𝑘(𝐴) ∩ 𝐿𝑘(𝐵) ≠ 𝐿𝑘(𝐴𝐵).

Figure 2 . 7 :

 27 Figure 2.7: Visibility kernel condition in 2D.Left: the edge 𝑃 𝑄 is collapsible and a valid embedding is preserved when the target vertex is located within the kernel 𝒦 𝒯 (𝑃 𝑄) (orange) of the polygon formed by the one-ring of the edge. Right: the kernel is empty and hence the edge 𝑃 𝑄 is not collapsible.

Figure 2 . 8 :

 28 Figure 2.8: Edge 𝑃 𝑄 of 𝜕𝛤 . 𝑃 𝑄 is candidate to be collapsed. The edges of 𝒯 are depicted with dashed black lines. Visibility Kernel 𝒦 𝒯 (𝑃 𝑄) is depicted (partially) in orange. The edges of 𝒵 (black solid lines) are not part of 𝒯. The green dots (𝑆 |𝒦 𝒯 (𝑃 𝑄) ) depict the subset of samples from 𝒮 located within 𝒦 𝒯 (𝑃 𝑄).

Figure 2 .

 2 Figure 2.10: Mutual tessellation. Left: before mutual tessellation. Middle: after mutual tessellation. Right: classification of tetrahedra in accordance to 𝜕𝛺 𝑖 . The edges of 𝒵 and 𝜕𝛤 are depicted with solid black and blue lines, respectively.

Figure 2 .

 2 Figure 2.11: Kernel of an edge 𝑃 𝑄 of 𝒵. 𝑃 𝑄 is candidate to be collapsed. The edges of 𝒵 are depicted with solid black lines. The visibility kernel 𝒦 𝒯 (𝑃 𝑄) of 𝑃 𝑄 is depicted in orange.

Figure 2

 2 Figure 2.14: Inaccessible regions. The grayed area of 𝛺 is inaccessible due to the restrictions imposed by 𝒦 𝒯 (𝑃 𝑄) to the simplicial tolerance.

Figure 2 .Figure 2

 22 Figure 2.16 illustrates all steps of our algorithm on a mechanical part. The input is a raw triangle soup (20k triangles). The tolerance volume is computed as the sub-level [0-0.6] of the

Figure 2 .

 2 Figure 2.16: Blade. From left to right: Input tolerance (𝛿 = 0.6%); 𝒵 after refinement (20.4𝑘 vertices); simplification of 𝜕𝛤 (5.3𝑘𝑣); mutual tessellation and simplification of 𝒵 (1.01𝑘𝑣); and the final output (752𝑣).

Figure 2 . 17 :

 217 Figure 2.17: Boundary sample points. The intersection between 𝒮 and ball located on the input surface 𝛴 is composed of a single connected component near a boundary.

Figure 2 .

 2 Figure 2.18 depicts a range scan of the Kitten point cloud with boundaries due to missing data. The holes are preserved by the non-closed variant of our algorithm (middle). Note that by not ignoring the parts of the zero-set outside 𝛺, we can also fill the large hole on the nearly flat area. Nevertheless, more work is needed to reliably deal with more complicated holes.

Figure 2 .

 2 Figure 2.18: Preserving or repairing holes. Left: Range scan of a kitten. Middle: Our output as non-closed surface (1.3𝑘 vertices). Right: Our output with holes filled (1.2𝑘 vertices).

Figure 2 .

 2 Figure 2.19: Several Components. From left to right: sampling of 𝜕𝛺, refinement until matching the topology and final output. The zero-set outside tolerance volume is depicted with dashed black lines and will be ignored.

Figure 2 . 20 :

 220 Figure 2.20: Dealing with a non-manifold geometry. From left to right: refinement until matching the topology, mutual tessellation and final output.

  significantly reduces over-refinement. However this choice might lead to large interpolated values at samples of the outer tolerance boundary. Since this does not hinder classification, we do not further refine in such cases. One way to implement this idea, is to replace the error 𝜖(𝑠) for classification by 𝜖(𝑠) = 1 -𝑓(𝑠)/ℱ(𝑠).

Figure 3 . 2 :

 32 Figure 3.2: Fertility. From top left to bottom right: 𝛿 is set to 0.15, 0.25, 0.4, 0.8, 2.0 and 8.0. The final vertex count is 4, 642, 2, 768, 1, 484, 767, 417 and 120, respectively. The input model is a surface triangle mesh of the fertility model with 14𝑘 vertices.

Figure 3 .

 3 Figure 3.1 illustrates our algorithm at work on a mechanical part (blade), for several separation distances between the boundaries of the tolerance volume. The overall time consumed by the algorithm ranges from 34 minutes for 𝛿 = 1.5% to around 7 hours for 𝛿 = 0.15%.Figure 3.2

Figure 3 . 3 :

 33 Figure 3.3: Armadillo. 𝛿 is set to 0.1% and 0.9%. The final vertex count is 26, 189 and 1, 518. The input model is a surface triangle mesh of the armadillo model with 173𝑘 vertices.

Figure 3 . 4 :

 34 Figure 3.4: Vase Lion. From top left to bottom right: 𝛿 is set to 0.2, 0.3, 0.5, 2.5, and 6.0. The final vertex count is 17, 723, 11, 042, 6, 127, 637 and 109, respectively. The input model is a surface triangle mesh of the vase lion model with 200𝑘 vertices.

Figure 3 . 5 :

 35 Figure 3.5: Evolution of mesh complexity over time for different 𝛿(%). Each mark depicts the completion of: refinement; simplification via halfedge collapses of 𝜕𝛤 ; edge collapses of 𝜕𝛤 ; mutual tessellation and halfedge collapses of 𝒵; edge collapses of 𝒵; and simplification of all edges. The input is a raw surface mesh of the fertility model (14kv).

Figure 3 . 6 :

 36 Figure 3.6: Comparisons. We plot the one-sided Hausdorff distance (output to input) (𝐻 |𝑂→𝐼 ) against the final number of vertices. A dot indicates a self-intersection of the output mesh. The input mesh is a clean surface triangle mesh of the fertility model (14kv).

Figure 3 . 7 :

 37 Figure 3.7: Comparisons. We plot the other-sided Hausdorff distance (input to output) (𝐻 |𝐼→𝑂 ) against the final number of vertices over the same data.

Figure 3 . 8 :

 38 Figure 3.8: Comparison with the Simplification envelopes [Cohen et al. 1996]. We plot the final number of vertices against the one-sided Hausdorff distance (output to input) (𝐻 |𝑂→𝐼 ). The input is a clean surface triangle mesh of Armadillo (173𝑘 vertices).

Figure 3 . 9 :

 39 Figure3.9: Aggressive simplification. As the tolerance in the Simplification Envelopes is generated by offsetting the vertices along the normals, this approach cannot simplify the geometry of highly undulating surfaces beyond a certain Hausdorff limit. Left: input mesh. Middle: output from the Simplification Envelopes (𝐻 |𝑂→𝐼 = 60%). Right: output from our algorithm (𝐻 |𝑂→𝐼 = 10%).

Figure 3

 3 Figure3.10: Distribution (in log scale) of normal deviation. We plot the distribution of normal deviation of our output for the fertility model (𝛿 = 0.2%) when condition 3 is not used, and when a 30%, 70% and 100% shrinkage factor is used. The final vertex count of the output is 3, 498, 3, 538, 3, 624 and 5, 715 respectively.

Figure 3 .

 3 Figure 3.11: Lucy. Left: Input model (50𝑘𝑣). Right: Output from our algorithm (16.7𝑘 vertices, 𝛿 = 0.13%).

Figure 3 .

 3 Figure 3.12: Screw. Left: when using a too large 𝜎 the decimation step terminates prematurely (2.49𝑘𝑣) due to normal constraints which are not satisfied during refinement. Right: Relaxing the normal constraint during halfedge collapse enables further decimation (1.35𝑘𝑣).

Figure 3 .

 3 Figure 3.13: Robustness to noise and type of input datasets. From left to right: point set, triangle soup with low noise, noisy point set, and triangle soup with high noise. The corresponding 𝛿 and output vertex count are 0.9, 1.2, 2, 5% and 2, 191, 1, 897, 1, 082, 502 respectively.

Figure 5 . 1 :

 51 Figure 5.1: Inter-surface mapping. We generate a dense, low-distortion map between two surfaces of arbitrary genus through the computation of a transport map minimizing the local variance of geodesic neighborhoods.

  𝜋 ∶ ℝ 𝑑 → ℝ 𝑑 between m and n such that it minimizes ∫ |𝜋(𝑥) -𝑥|n(𝑥)𝑑𝑥,

  .8) between two discrete probability distributions m and n is of the form 𝜋 = D m D v HD w D n , with unique diagonal matrices D v and D w defined by vectors v, w satisfying: w n = 1 |m| , D w H ⊤ D v m = 1 |n| , (5.10) where, both 𝜋 and H are of size |m| × |n|, D x represents the diagonal matrix with the vector x

  Figure 6.1: Mapping via variance minimization in 2D. An 8-shaped pointset 𝑋 and a squareshaped pointset 𝑌 (left); initial and intermediate solution (middle); and convergence after 2 iterations of a map between the two input shapes (right).

Figure

  Figure 6.2: Weighting functions. From left to right are depicted 𝐾(𝑥, 𝑦), 𝐾(𝑥, 𝑦)/√𝐾(𝑥, 𝑥)𝐾(𝑦, 𝑦) and 2𝐾(𝑥, 𝑦)/(𝐾(𝑥, 𝑥) + 𝐾(𝑦, 𝑦)). The neighborhood of the pink vertex corresponding to {𝑦 ∶ 𝑊 𝑥 (𝑦) > 𝜀} is shown with green vertices.

Figure 6

 6 Figure 6.3: Annealing. Percentages of mass densities outside the interval [0.999, 1.001], minimum (blue) and maximum (red) densities. From left to right: After 𝜂 stabilizes, 95%, [0.415, 2.41]; First step of annealing: 90%, [0.54, 2.1]; Second step: 90%, [0.97, 1.2]; Final map: 0.01%, [0.9902, 1.009].

Figure 6 . 4 :

 64 Figure6.4: Local variable instantiation. For a point 𝑥, we instantiate variables only in the neighborhood of the mapping of its neighborhood. By localizing the samples it may be mapped to, we dramatically reduce the search space, without adverse effects.

Figure 6 . 5 :

 65 Figure6.5: Hierarchical solving. An isotropic sphere mesh (65𝑘𝑣) is mapped to the input skull mesh (35𝑘𝑣) (left) by relocating the vertices of the sphere mesh onto the skull using 𝜂, across 3 Euclidean steps of the hierarchical solve (diffusion steps are not depicted as the corresponding 𝜂's live in diffusion space).

Figure 6 . 7 :

 67 Figure 6.7: Diffusion distance estimates on thin parts. Colored spheres on the dog depict the final sampling during the last, finest diffusion stage. Top: with 100 eigenvectors, sampling is insufficient to get satisfactory precision on thin parts of the dog(feet, tail, ears). This may translate into twists of the mapping (right). Bottom: using 300 eigenvectors significantly improves precision.

Figure 6 . 8 :

 68 Figure 6.8: Fixing flips and twists. Starting from a suboptimal map having discontinuity(left, cactus (top) to itself (bottom) built by intentionally placing bad constraints), ironing proceeds through front propagation to remove all artifacts. For each strip, blue indicates the visited parts while the non visited parts are in pink. In case of multiple fronts, each front is considered independently. Note that the main idea behind whole approach is to slightly distort the neighborhood and favor already visited parts. Using high value for 𝜔 (close to 1) cannot unfold the mapping, while low values (close to 0) seriously pull 𝜂's towards already-visited parts.

  D m D v HD w D n ▷ Relaxed transport plan Note that each KL projection requires solving an equation for determining parameter 𝑟,

  ] ← lb[𝑖] 9: function 𝑟 = Global_Mass_v(lb, ub) 10: x(𝑟) ←Update_Vector_Global_Mass(𝑟, lb, ub) 11: 𝑟 ← solution to m ⊤ D x HD w n = 1 12: function 𝑟 = Global_Mass_w(lb, ub) 13:

  𝒮 0 , 𝒮 1 … 𝒮 ℓ do // traverse hierarchy 3: while !converged do 4: 𝜋 ← min 𝒞(., 𝜂, (𝑚 ′ , 𝑛 ′ )) // solve relaxed optimal transport problem 5: (𝑚 ′ , 𝑛 ′ ) ← min 𝒞(𝜋, 𝜂, .) // update local mass relaxation 6: 𝜂 ← min 𝒞(𝜋, ., (𝑚 ′ , 𝑛 ′ )) // relocate centers to geodesic barycenters 7: return [𝜋, 𝜂, (𝑚 ′ , 𝑛 ′ )] ▷ Area-relaxed Mapping

Figure 6 . 9 :

 69 Figure6.9: Relaxing area preservation. We compute a map between a canonical torus and a swollen torus using three different relaxation factors (from left to right : 1, 2 and 4). Top: colors indicate the relaxed mass transported to this location (blue=small, red=large). As expected the part mapping to the swollen region gets more mass and vice-versa. Middle: we depict the quasi-conformal distortion for each triangle, ranging from 0.0 (degenerated) to 1.0 (perfect). Larger relaxation factors get better conformality. 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. Bottom: closeups on the triangle meshes to highlight improved conformality. Note that although the mass distribution is locally non smooth, the resulting map is smooth nonetheless since our transport stencils are several rings wide.

[

  The CGAL Project, 2016] for geometric computations. Parallelization was acheived using Intel Threading Building Blocks library. Unless otherwise indicated, we used the first 100 eigenvectors of the heat diffusion operator to approximate geodesic distances. For visual evaluation of the maps we use one color per vertex on the triangle meshes. Each vertex color is set based on a smooth 3D function, sometimes modulated by a 3D axis-aligned checkerboard (called grid coloring hereafter) to improve legibility.

Figure 7 . 1 :

 71 Figure 7.1: Centaur in three isometric poses. A dense map is successfully established between three isometric poses of the centaur model. No user interaction was required.

Figure 7 . 2 :

 72 Figure 7.2: Isometric deformation of a dog. Our approach successfully establishes a proper dense map between three isometric poses of a dog; the left pose is chosen as reference for coloring.

  Fig. 7.2 shows a dog in three isometric poses. Except for the initialization stage we solve for mappings that minimize only the variance of the transport plans. The mapping to the middle dog is obtained without any user interaction, the Wasserstein-2 term used during initialization being sufficient to disambiguate the intrinsic left/right symmetry.The mapping to the right dog however requires 4 user-specified constraints, one per leg, to avoid left/right flipping. Grid coloring of the source dog (left) is mapped to the two target dogs (middle, right) using the variance-minimizing transport plan computed by our algorithm. For both the dog and centaur models the isometric maps are found without any area relaxation.

Figure 7 . 3 :

 73 Figure 7.3: Mapping with user interaction. Top left: a source right hand is mapped to a target left hand using two user-defined constraints. Top right: the map is depicted through a color grid. Bottom left: source mapped onto the target mesh. Bottom right: target mapped onto the source mesh.

Fig. 7 .

 7 Fig. 7.3 illustrates the mapping between two human hands. The hands differ in types (left vs.

Figure 7 . 4 :

 74 Figure 7.4: Non-isometric mapping. The dog model is mapped to the horse model without any user-specified constraints.

Figure 7 . 5 :

 75 Figure 7.5: Highly non-isometric mapping. A kitten is mapped to a torus. Left: grid on kitten mapped to the torus. Middle: kitten mesh mapped to the torus. Right: torus mesh mapped to the kitten. Closeups highlight the angle distortion of the mapping.

Fig. 7 .

 7 Fig. 7.4 illustrates a fully automatic mapping between a dog and a horse. No point constraints were used for this example. Ears map to ears, legs to legs, and so on, solely due to our computation of a variance-minimizing map without any mass relaxation. We used 300 eigenvectors to get a reliable estimate of diffusion distances on the thin features -legs and tail -of the dog model: a smaller number of eigenvectors would lead to unreliable geodesic distance estimates, slowly translating into twists of the map along thin protrusions in the Euclidean stages (see Fig. 6.7) which would then rely totally on ironing.

Figure 7 . 6 :

 76 Figure 7.6: Non-isometric mapping. An area relaxation (factor of 10) is used, along with 4 constraints (one at the extremity of each limb) to map the two very different shapes of a gorilla and a man. Requiring more relaxation to better improve the map further would face very slow convergence.

Figure 7 . 7 :

 77 Figure 7.7: Non-isometric mapping. For this challenging example, no mass relaxation was used. Four constraints, one at the extremity of each limb were used to guide the map.

Figure 7 . 8 :

 78 Figure7.8: Mapping between different genus. Top left: a sphere with a small handle is mapped to a sphere. The mapping is discontinuous, but the discontinuity is local. Top right: same with a mediumsize handle. Bottom: a regular torus is mapped to a sphere. The discontinuity is no longer local, but the map is continuous over half of the sphere.

a

  larger handle and (3) a canonical torus. As expected the maps are discontinuous but the discontinuities are only local. No mass relaxation was used in these examples. Fig. 7.9 depicts a mapping between two poses of a kid. While on one pose the arm merges with his chest, the map is only little distorted.

Figure 7 . 9 :

 79 Figure 7.9: Mapping with local contact. While the arm of the baby touches (i.e., merges with) his chest on the second pose, we find a good map with 4 constraints (one at the extremity of each limb).

Table 3 . 1 :

 31 Timing (in seconds) for each step of our algorithm for the blade model depicted in Figure 2.16.

	Stage	# Vertices (𝒵)	Time (s)	Time per iteration
	Refinement	20, 447	655	0.0319
	halfedge -𝜕𝛤	10, 217	326	0.0318
	general -𝜕𝛤	5, 346	4, 658	0.956
	halfedge -𝒵	2, 292	153	0.050
	general -𝒵	1, 015	1, 478	1.157
	All edges	752	4, 537	17.185
	Total		11, 807	0.2941
	Most of the time spent by the algorithm is in the exhaustive search to find the best point
	location for an edge collapse operator, and this time escalates as the decimation proceeds. Un-
	like other mesh decimation algorithms, the running time of our algorithm is decreasing with
	parameter 𝛿 of the specified tolerance. A small tolerance requires dense sampling and hence |𝒮|
	increases together with the time consumed to classify the samples. Another dominating factor is
	the sampling density used to probe a kernel when searching for the best point of a general edge
	collapse operator. The halfedge collapse operators are on average two orders of magnitude faster,
	but are not sufficient to generate coarse meshes. On the positive side, each operation performed
	over all sample points and tetrahedra of the triangulation is parallelizable. Figure 3.14 plots the
	speed up in the run-time of our algorithm versus the number of CPU cores.

Figure

3

.14: Speed up. We plot the speed up in the run-time of our algorithm versus the number of CPU cores (input: Fertility, 𝛿 is set to 0.4%).

  Reformulation with auxiliary variables. Given a point 𝑥 and a measure 𝜇 = ∑ 𝜇 𝑖 𝛿 𝑥 𝑖 , we denote by var(𝜇, 𝑥) the variance of 𝜇 with respect to 𝑥: to be understood as geodesic distances on the surface that 𝜇 samples and that 𝑥 lies on. Using this notion of variance w.r.t. a position, we can now introduce a set of positions 𝜂 𝑥 𝑖 and 𝜂 𝑦 𝑗 , one for each point 𝑥 𝑖 and 𝑦 𝑗 . These auxiliary variables, which we call centers for reasons that will become clear soon, allow us to write our transport cost between 𝑋 and 𝑌 as a function of both the transport plan 𝜋 and the set of all centers 𝜂: Note that this cost functional equals the energy defined in Eq. (6.1) if each center 𝜂 𝑥 𝑖 (resp., 𝜂 𝑦 𝑗 ) corresponds to the (weighted) center of mass of the image of the normalized weighting functions 𝑊 𝑥 𝑖 (resp. 𝑊 𝑦 𝑗 ) by the transport plan 𝜋. 𝜋, minimizing 𝒞 (𝜋, .) is efficiently achieved by relocating the centers to the geodesic barycenters of the images of the normalized weighting functions by 𝜋 (see Sec. 6.2.3); for fixed centers, minimizing 𝒞 (., 𝜂) requires the computation of an optimal transport plan under the usual mass preservation constraints. Alternately minimizing the cost 𝒞 over 𝜋 𝑖𝑗 and 𝜂 as described in Alg. 2 is a robust and efficient way to treat the original non-linear problem since both minimizations are convex problems.

	Minimization through alternating convex problems. This augmented cost functional 𝒞 (𝜋, 𝜂)
	is now amenable to an alternating minimization similar to Expectation-Maximization (EM): for
	a fixed transport plan
		var(𝜇, 𝑥) = ∑ 𝜇 𝑖 𝑑(𝑥 𝑖 , 𝑥) 2 ,	(6.2)
	where distances are 𝒞 (𝜋, 𝜂) = ∑ 𝑖 ∑ 𝑗	𝑚 𝑖 var [𝜋 𝑋 ( 𝑛 𝑗 var [𝜋 𝑌 ( ∑ 𝑘 𝑛 𝑘 𝑊 𝑦 𝑗 (𝑦 𝑘 ) 𝑊 𝑥 𝑖 𝜇 ∑ 𝑘 𝑚 𝑘 𝑊 𝑥 𝑖 (𝑥 𝑘 ) ) , 𝜂 𝑦 𝑗 ] . ) , 𝜂 𝑥 𝑖 ] + 𝑊 𝑦 𝑗 𝜈	(6.3)

  .3), i.e., 𝑥 𝑖 ) denotes the value of the weighting function 𝑊 𝑥 at 𝑥 𝑖 . We thus use the following closed form formula as an update rule for 𝜂 once a map 𝜋 is known:

	arg min 𝜂 𝑥	⎣ ⎡ ⎢	𝑥 𝑖 ∈𝑋 ∑	𝑦 𝑗 ∈𝑌 𝑊 𝑥 (𝑥 𝑖 ) ∑	⎦ [𝜋 𝑖𝑗 ‖𝛷(𝑦 𝑗 ) -𝜂 𝑥 ‖ 2 ] ⎤ ⎥	,
	where 𝑊 𝑥 (𝜂 𝑥 =	∑ 𝑖,𝑗	𝑊 𝑥 (𝑥 𝑖 ) 𝜋 𝑖𝑗 𝛷(𝑦 𝑗 ) ∑ 𝑖,𝑗 𝑊 𝑥 (𝑥 𝑖 ) 𝜋 𝑖𝑗	.	(6.5)

  and optimize the functional 𝒞 ′ (𝜋, 𝜂) obtained from 𝒞(𝜋, 𝜂) in (6.3) by replacing 𝜋 𝑋 by 𝜋 ′ 𝑋 and 𝜋 𝑌 by 𝜋 ′

	operators	𝜋 ′ 𝑋 (𝛿 𝑥 𝑖 ) =	1 𝑚 ′ 𝑖	∑ 𝑗	𝜋 𝑖𝑗 𝛿 𝑦 𝑗 , and, 𝜋 ′ 𝑌 (𝛿 𝑦 𝑗 ) =	1 𝑛 ′ 𝑗	𝑖 ∑	𝜋 𝑖𝑗 𝛿 𝑥 𝑖 ,

Acknowledgments

Our algorithm is implemented in C++ using the CGAL library for the triangulation data structures and the Intel Threading Building Blocks library for parallelization. 3D tolerance volumes are rendered via 3D texture mapping using pixel shaders from the NVIDIA Cg Toolkit.

Interlocked loops

Theorem A.0.1. Assume 2 interlocked loops, each formed by joining a segment of length 𝑙 𝑒 with a continuous curve of length 𝑙 𝑐 . Then the continuous curves of the two loops cannot be ( 𝑙 𝑒 +𝑙 𝑐 4 + 𝑙 𝑒 , 0)-separated.

Proof. Each loop is contained within a ball of radius 𝑙 𝑒 +𝑙 𝑐 4 . Let 𝐵 be such a ball for the first loop 𝐶 1 . Because the two loops are linked, the part of the second loop 𝐶 2 lying in 𝐵 cannot be linearly B separated from the first loop. If this part only consists of the curve part, the conclusion follows. Else, let 𝐵 ′ be the ball obtained by enlarging 𝐵 by 𝑙 𝑒 . Now 𝐵 ′ ∩ 𝐶 2 must contain the two endpoints of the segment part of 𝐶 2 . Also, 𝐵 ′ ∩ 𝐶 2 and 𝐶 1 are not linearly separable. Since the segment lies in the convex hull of the curve part of 𝐵 ′ ∩ 𝐶 2 , the conclusion follows.

Part II

Low Distortion Inter-surface Mapping via

Variance-Minimizing Mass Transport

Additional control over the mapping. The cost we defined in Eq. ( 6.3) measures the regularity of a transport plan via its local variance. However, user control over the final mapping is often desirable: the user may want to prescribe a few correspondences manually. If we wish 𝑥 𝑖 to map to 𝑦 𝑗 , we constrain the center 𝜂 𝑥 𝑖 to be 𝑦 𝑗 and 𝜂 𝑦 𝑗 to be 𝑥 𝑖 . Additionally, the user may want to favor mappings that match points with similar attributes such as colors or shape signatures. This can be done by adding to our cost functional an additional transport cost, with a cost matrix given by distances between feature vectors associated with each input point. Finally, adding a small amount of the usual 2-Wasserstein distance between the surfaces can help finding good initial solutions or help the solver disambiguate symmetries in the models, and thus finding better maps automatically, see Sec. 7.

Algorithms and Numerics for Scalability

We designed several additional algorithmic and numerical components to offer an efficient and scalable solver for our approach.

Weighting function

A natural (and consistent) choice as weighting function 𝑊 is to rely upon the heat kernel from Eq. (5.1) in order to define an approximate geodesic neighborhood. While this choice is nicely isotropic for small diffusion times 𝑡, we noticed that its support becomes off-centered for large values of 𝑡 (see Fig. 6.2). We use instead the following weighting function:

that remains localized and centered for all diffusion times while enforcing 𝑊 𝑥 (𝑥) = 1. In addition, we truncate the weighting function below 𝜀 = 0.5 (Fig. 6.2) in order to further localize its associated geodesic neighborhood, which will reduce computational times required to evaluate the geodesic variance.

Fast and robust optimal transport

Instead of solving Step 4 of Alg. 2 via linear programming, we use Sinkhorn iterations (Alg. 1) to compute the optimal transport plan efficiently. While the original approach requires assembling the kernel matrix H through its elements 𝐻 𝑖𝑗 = exp(-coef(𝜋 𝑖𝑗 )/𝛾), where coef(𝜋 𝑖𝑗 ) refers to the coefficient of 𝜋 𝑖𝑗 in the cost function (6.3), we proceed in log space instead to improve robustness.

Comparisons. We first compare our approach on an isometric model from the TOSCA dataset, with the blended intrinsic maps (BIM) [START_REF] Kim | Blended intrinsic maps[END_REF] and a recent continuous matching approach (referred to as CM) contributed by [START_REF] Corman | Continuous matching via vector field flow[END_REF]. In Fig. 7.10 we plot the quality of correspondences for our approach, with or without ironing. Our method outperforms BIM and underperforms CM.

-2 ⋅ 10 -2 0 2 ⋅ 10 -2 4 ⋅ 10 -2 6 ⋅ 10 -2 8 ⋅ 10 -2 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 In Fig. 7.11 we evaluate our approach to our main target problem, i.e., non-isometric mappings, on six models from the body model dataset [START_REF] Yang | Semantic parametric reshaping of human body models[END_REF], with area relaxation (factor 5). One man model is chosen as reference. We depict for each model the xyz mapping function as well the relaxed area. All maps are point-to-point and smooth, and the relaxed area reflects the differences in proportions of the limbs and improves the quality of the maps under the area relaxation bound. We also evaluate the correspondence quality and compare it to BIM, while taking as ground truth the semantic-based correspondences available from the dataset. Our method performs on average similarly to, or slightly better than, BIM. Note however that our approach is neither driven by semantic nor guided by feature points, hence we are not really competing on semantic-based correspondences. Our approach is instead devised to generate maps that are as conformal as possible within area preserving or area bounding constraints.

Timings.

Computing the eigenvectors using spectra as a preprocessing step takes up to 2.5 seconds for 100 eigenvectors on a mesh with 5k vertices and 236.4 seconds for 200 eigenvectors on mesh with 150k vertices. Assuming 𝑘 be the average neighborhood cardinality and 𝑛 the total number of samples, the computation cost of each Sinkhorn iteration (including building/storing H) is 𝒪(𝑐𝑘 2 𝑛), where 𝑐 represents the transport plan spreading (average number of points to which a point is mapped) that depends on 𝛾. Furthermore, we parallelize computations, requiring 2 to 5 mins for small (20kV) meshes, and up to 40 mins for 200kV. Note that without hierarchy and local variable instantiation, each Sinkhorn iteration would have a time and space complexity of 𝒪(𝑛 2 ), with an additional 𝒪(𝑐𝑘𝑛 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 geodesic centers of mass and the Sinkhorn iterations.

Convergence to Dirichlet Energy

Let 𝑓 ∶ 𝑆 ↦ 𝑆 ′ be a diffeomorphism between two surfaces. Take for the weighting function at 𝑥 ∈ 𝑆 a Gaussian with covariance 𝜎 2 𝐼, where 𝜎 is a small positive number. Let 𝑣 be a random vector drawn accordingly. Then in the cost we define, the integrand at 𝑥 is by definition the variance of 𝑓(𝑣). For small 𝜎, we may approximate locally the two manifolds by their tangent spaces. Under this approximation 𝑓(𝑣) is a Gaussian with covariance 𝜎 2 𝐷𝑓𝐷𝑓 * , hence var(𝑓(𝑣)) = 2𝜎 2 ||𝐷𝑓|| 2 𝐹 . As a consequence, the total energy is proportional to the sum of the Dirichlet energy of 𝑓 and the Dirichlet energy of its inverse.

The same conclusion holds for non Gaussian weighting functions, as long as they have isotropic covariances with total variance independent of 𝑥, and tending to 0.

Appendix B KL projections

Let us derive the KL projection of a matrix 𝜌 on 𝛱 1 , assuming the entries 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: In particular, 𝜙 𝑖 is determined by 𝜙. More specifically, 𝜙 𝑖 = 𝑓 𝑖 (𝜙), where 𝑓 𝑖 is the closest point projection on interval [𝛼 𝑖 𝑚 𝑖 / ∑ 𝑗 𝜌 𝑖𝑗 , 𝛽 𝑖 𝑚 𝑖 / ∑ 𝑗 𝜌 𝑖𝑗 ], which is a non decreasing piecewise linear function with two nodes. In order to determine the correct value for 𝜙, it suffices to solve the equation ∑ 𝑖𝑗 𝜌 𝑖𝑗 𝑓 𝑖 (𝜙) = 1, so that the global mass constraint is satisfied. Taking into account the change of variables due to the normalized masses m and n yields the algorithm described in Section 6.3.