Massive amounts of data are being generated, collected and processed all the time. A considerable portion of them are sampled from objects with geometric structures. Such objects can be tangible and ubiquitous in our daily life. Inferring the geometric information from the data, however, is not always an obvious task. Moreover, it's not a rare case that the underlying objects are abstract and of high dimension, where the data inference is more challenging.

This thesis studies two problems on geometric data analysis. The first one concerns metric reconstruction for filamentary structures. We in general consider a filamentary structure as a metric space being close to an underlying metric graph, which is not necessarily embedded in some Euclidean spaces. Particularly, by combining the Reeb graph and the mapper algorithm, we propose a variant of the Reeb graph, which not only faithfully recovers the metric of the filamentary structure but also allows for efficient implementation and convenient visualization of the result.

Then we focus on the problem of shape comparison. In this part, we study the stability properties of some recent and promising approaches for shape comparison, which are based on the notion of functional maps. Our results show that these approaches are stable in theory and potential for being used in more general setting such as comparing high-dimensional Riemannian manifolds.

Lastly, we propose a pipeline for implementing the functional-maps-based frameworks under our stability analysis directly on point cloud data. Though our pipeline is experimental, it undoubtedly extends the range of applications of these frameworks.

Chapter 1. Résumé

Bien que cette généralisation ne bénéficie pas des garanties théoriques, elle permet d'étendre le champ d'application des méthodes fondées sur les applications fonctionnelles. En effet, du point de vue plus général de l'analyse de données, les résultats expérimentaux que nous obtenons en étudiant les formes 3D laissent à penser que l'approche fonctionnelle appliquée à l'analyse géométrique de données plus abstraites et/ou de grande dimension a du potentiel.

Globalement, dans cette thèse, nous présentons une analyse théorique dans deux contextes différents, cette analyse consiste à inférer des informations à partir de données géométriques. Nous avons proposé des garanties théoriques variées, mais avons également apporté des contributions dans le domaine des applications : un nouvel algorithme de reconstruction de structure filamentaire et quelques extensions des approches fonctionnelles existantes sont proposées dans les première et seconde parties.

CHAPTER 1

Résumé

Des quantités importantes de données sont générées, collectées et analysées en permanence. Dans de nombreux cas, ces données sont échantillonnées sur des objets possédant une structure géométrique spécifique. Certains de ces objets sont omniprésents dans notre vie quotidienne. L'inférence de leur structure géométrique à partir de ces données est une tâche souvent ardue. Cette tâche est d'ailleurs rendue plus difficile si les objets sous-jacents sont abstraits ou de grande dimension.

Dans cette thèse, nous nous intéressons à deux problèmes liés à l'analyse géométrique de données. Dans un premier temps, nous nous penchons sur l'inférence de métriques pour des structures filamentaires. Les structures filamentaires apparaissent naturellement dans les données du quotidien telles que les collections de traces GPS collectées par des véhicules sillonnant le réseau routier, les distributions d'épicentres de tremblements de terre se concentrant le long de failles géologiques ou encore les réseaux de vaisseaux sanguins. Cette liste n'étant pas exhaustive. Dans de nombreuses applications pratiques, les données apparaissent sous la forme d'une matrice de distances entre points, ce qui permet de mettre sur les données une métrique discrète. Par conséquent, il est intéressant de construire une inférence des données à partir de distances.

En général, nous représentons des structures filamentaires par des espaces métriques que nous supposons être proches d'un graphe métrique. La méthode que nous proposons est une variante des graphes de Reeb, elle combine les graphes de Reeb et l'algorithme Mapper. Cette variante permet non seulement d'approcher de façon judicieuse la métrique associée aux structures filamentaires, mais elle permet également d'implémenter et de visualiser le résultat aisément. En particulier, nous fournissons des garanties théoriques fortes pour cet algorithme.

Nous nous concentrons ensuite sur le problème de la comparaison de formes. Nous étudions un ensemble de méthodes récentes et prometteuses pour la comparaison de formes qui utilisent la notion d'applications fonctionnelles. Depuis les travaux précurseurs sur les applications fonctionnelles, de nombreux travaux ont suivi traitant de l'exploration de formes, de la segmentation d'images, de traitement de champs de vecteurs etc. Cependant, malgré le succès de telles méthodes, beaucoup moins d'attention a été portée à l'étude des propriétés théoriques fondamentales de ces approches fonctionnelles. Nos résultats théoriques montrent que ces approches sont stables et peuvent être utilisées dans un contexte plus général que la comparaison de formes comme par exemple la comparaison de variétés Riemanniennes de grande dimension.

Enfin, en nous reposant sur notre analyse théorique, nous proposons une généralisation des applications fonctionnelles aux données de type nuages de points. Dans le domaine de l'analyse de forme, la représentation des formes 3D par des nuages de points est sans doute la plus fréquente et primitive. Cependant, la plupart des approches fonctionnelles en traitement de données géométriques utilisent comme données des formes discrétisées par des maillages triangulaires. Nos travaux produisent des résultats raisonnables par rapport à ceux obtenus avec les maillages, malgré le fait que les nuages de points contiennent une information géométrique intrinsèque beaucoup plus pauvre que celle fournie par des maillages.

CHAPTER 2

Introduction

Data play a role that is becoming increasingly important nowadays. Massive amounts of data are being generated, gathered and processed in various areas such as scientific research, economic activities, manufacturing productions all the time.

Perhaps data come with geometric structures are the most familiar to us. As a matter of fact, our perception system takes in and analyzes geometric data everyday. For example, looking at the scattering blue points in Figure 2.1(a), probably people would outline the skeleton of the point cloud and draw the red curves in mind. If the perception system is fed multiple objects like in Figure 2.1(b), then a more complicated task is accomplished-people would easily discriminate the two shapes and highlight the red bump of the second one as a prominent difference in mind. a sphere and a deformed one, the deformation in-between can be visually captured.

However, these geometric data inferences are limited in several aspects. First of all, our perception system outputs only qualitative results. When differentiating shapes, we can tell where the difference is, but can't describe quantitatively how large it is. Besides, though we are able to comfortably deal with low-dimensional objects such as those depicted in Figure 2.1, we have trouble even visualizing objects of dimension more than 3, not to mention outlining or comparing. Moreover, our perception system can not process abstract inputs directly. If the plots in Figure 2.1 are replaced with tables listing the coordinates of the points, we can hardly see anything meaningful there.

Fortunately, with the aid of computers, we are able to quantitatively process complex and abstract data. However, computers are not gifted with a perception system. It is then critical to develop methods that guide the computers to efficiently infer information from data. Approaches have been taken for data inference from various perspectives. For example, dimension reduction algorithms [START_REF] Belkin | Laplacian Eigenmaps for Dimensionality Reduction and Data Representation[END_REF], Lafon 2004, Tenenbaum 2000b] attempt to project high-dimensional data into low-dimension spaces for visualization; clustering algorithms [Planck 2006, Lloyd 1982[START_REF] Ester | A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise[END_REF], Comaniciu 2002] separate data into groups of a relatively small number; and more recent frameworks of topological data analysis [START_REF] Edelsbrunner | [END_REF][START_REF] Zomorodian | Computing persistent homology[END_REF], Carlsson 2009, Chazal 2012] enable extracting data that are not embedded in an Euclidean space.

In this thesis, we study two problems on geometric data inference respectively related to the perception tasks illustrated in Figure 2.1.

The first problem is metric reconstruction for filamentary structures. An intuitive example of data sampled from a filamentary structure is the blue point cloud in Figure 2.1(a), which looks like some graph, say, the red curves. We assume that the given filamentary structure X is a metric space (not necessarily embedded in some Euclidean space as the example) that is close with respect to the so-called Gromov-Hausdorff distance to an unknown metric graph 1 . And our goal is to construct a metric graph that is close to X in the Gromov-Hausdorff distance. We propose a new method to address this question in a rigorous mathematical framework. This method relies on the Reeb graph and so-called Mapper algorithm, which are soon to be introduced.

As a result of our investigation to the first problem, we develop a method to construct a metric graph G that is guaranteed to be close to X. In general, The metric of X is impossible to be fully reconstructed by G. How to evaluate and visualize the parts of X which are less well approximated by the resulting graph is then a natural follow-up problem. Though there is a natural map between X and G constructed by our method, it is still a challenging problem. For instance, naive point-to-point comparison is unstable due to noisy data in practice. Motivated by this, the second problem is about comparing two geometrical objects associated by a map. This is indeed a vast and difficult problem. Our contribution to it is different than that to the first one: instead of proposing a new model/method, we consider some recent and promising methods on comparing 3D shapes that are potentially applicable in the general setting, say, comparing Riemannian manifolds of higher dimensions. With the understandings obtained from analyzing the chosen methods, we take some experimental approach based on them to the general problem. Especially, our main effort is devoted to analyzing the stability properties of the chosen methods. More precisely, the goal is to understand how the outputs of the method change with respect to perturbations on the inputs and/or the parameters. Investigating the stability of a method not only verifies its robustness in theory, but also provides valuable insights for implementations.

In the following we give a brief overview of the methods involved in this thesis. 1 A metric graph is a graph endowed with a metric.

Reeb Graphs

Let f be a real-valued continuous function defined on a topological space X. We then define a relation on X with respect to f : for x, y ∈ X, we let x ∼ f y if and only if f (x) = f (y) and x, y are in the same connected component of the pre-image f -1 (f (x)). The Reeb Graph (originally proposed in [Reeb 1946]) of f is then defined as the quotient space X\ ∼ f . For the point of view of level-set of f , for any a ∈ R, we decompose the subset f -1 (a) = {x ∈ X : f (x) = a} into one or more connected components, and then collapse each component into a point in the quotient space (see Figure 2.2 for an illustration).

f f

Reeb graph Mapper X a Figure 2.2: Illustrations of the Reeb graph (on the left of X) and the output of the mapper algorithm (on the right of X) with respect to the same function f . The cover used in the mapper algorithm is the 4 intervals colored differently to the right-most.

Under some regularity conditions, the Reeb graph is a 1-dimensional structure, making interpreting results convenient. Besides, the Reeb graph is defined upon a function on X instead of upon X itself. This property allows the users to view X through the lens of different functions. In the recent decades, its capability of encoding both the geometrical and the topological features of shapes is appreciated in the area of computer graphics, say, shape comparison [Hilaga 2001[START_REF] Escolano | [END_REF], skeleton extraction [Ge 2011] (see [Biasotti 2008] for a survey).

Mapper Algorithm

Many works studying variations on the Reeb graph have been proposed, such as the extend Reeb graphs [START_REF] Biasotti | [END_REF][START_REF] Escolano | [END_REF]], the contour trees [Pascucci 2003, Van Kreveld 2006, Carr 2003] and so on. A generalization of the Reeb graph-the mapper algorithm [Singh 2007] also falls into this category. The mapper algorithm takes in a continuous function f on X as an input. Instead of directly collapsing the connected components of the level-set of each single point a ∈ f (X), the algorithm considers a cover of f (X), the image of f . By a cover we mean a collection of open sets

{U i } i∈Γ such that f (X) ⊂ i∈Γ U i . For any U i , let f -1 (U i ) = α∈Λ i V α
i , where V α i 's are the connected components of f -1 (U i ). In the output of the mapper algorithm, each set V α i is represented as a node. An edge connects two nodes if the intersection of the corresponding sets is non-empty. A triangle is added among three nodes if the three corresponding sets are intersecting, likewise a tetrahedron is created in the case of four sets intersecting. In general, a k-simplex is added to the output if there exist k -1 sets have a common non-empty subset. Thus the output of the mapper algorithm is not necessarily a graph. In general it is a so-called abstract simplicial complex (see the definition in Section 3.3.2).

The mapper algorithm is more of a qualitative analysis tool for high-dimensional data set. The main aim of this algorithm is to robustly capture and easily visualize the topological structure of the input data.

Functional Maps

In [START_REF] Ovsjanikov | [END_REF]], the authors propose a functional representation of a map between two shapes. Given two shapes M and N and a map T : M → N , a functional map, T F , is a pull-back induced by T . Namely, given a real-valued function w ∈ F(N ), we define T F (w) = w • T ∈ F(M ). Therefore T F is a map from F(N ) to F(M ). The first benefit of this representation is that T F is a linear operator between the two function spaces. In fact, by definition of T F , T F (αf + βg) = (αf + βg) • T = αf • T + βg • T = αT F (f ) + βT F (g), where α, β ∈ R and f, g ∈ F(N ).

This linearity property indicates that T F admits a (potentially infinite) matrix representation. Suppose that {ϕ N i }, {ψ M j } form a basis of F(N ) and of F(M ) respectively. Let w = i=1 a i ϕ N i and accordingly

T F (w) = j=1
b j ψ M j , then there exists a unique matrix C T such that C T a = b holds for any w, where

a = (a 1 , a 2 , • • • , a n , • • • ) and b = (b 1 , b 2 , • • • , b n , • • • ).
We then call C T the matrix representation of T F with respect to basis {ϕ N i } and {ψ M j }. Especially, the authors take the first k M and k N eigenfunctions of the Laplace-Beltrami operators on M and N respectively as approximations to the full bases of the function spaces, and obtain a compact matrix representation of C T which is of dimension k M × k N .

Geometrically, truncating the eigenbasis in this way amounts to putting a low-frequency filter on the function space. In discrete setting, M and N are usually represented by m and n discrete points. Thus the function spaces are respectively vector spaces of dimension m and n, and T F is a m × n matrix. Using above scheme, T F is then approximated by a k M × k N matrix C T written in the truncated bases. It has been validated empirically that for a pair of shapes with m, n ∼ 10 4 , k M , k N ∼ 102 are sufficient to obtain a reasonable approximation.

Shape Difference Operators

Based on the functional map, in [Rustamov 2013] the authors propose a framework to encode the differences between shape with so-called shape difference operators. Let ⟨•, •⟩ F (M ) and ⟨•, •⟩ F (N ) be two inner products on the function spaces F(M ) and F(N ) respectively. In general, the equality ⟨f, g⟩ F (N ) = ⟨T F (f ), T F (g)⟩ F (M ) does not hold for all f, g ∈ F(N ). The key observation of [Rustamov 2013] is that by the Riesz representation theorem and the linearity of T F , there always exist a self-adjoint linear operator on the function space F(N ), S such that ⟨f, S(g)⟩ F (N ) = ⟨T F (f ), T F (g)⟩ F (M ) .

This operator compensating the differences between the two inner products in some sense also captures the differences between the shapes induced by T F (and equivalently T ), and is called the shape difference operator. Particularly, in [Rustamov 2013], for each of the shapes M and N , two typical function spaces and the associated inner products are discussed, resulting in two types of shape difference operators that capture differences between M and N from separated aspects.

On the other hand, it is also interesting to study the gap between ⟨f, g⟩ F (N ) and ⟨T F (f ), T F (g)⟩ F (M ) . The framework of [START_REF] Ovsjanikov | [END_REF]] constructs a functional to measure the gap at any function f :

F (f ) = ⟨T F (f ),T F (f )⟩ F (M ) ⟨f,f ⟩ F (N )
for a specified inner product, which is one of the two considered in [Rustamov 2013]. By maximizing the functional, the authors identify a function at which the gap is the largest. Moreover, a multiscale framework is presented there by putting an extra constraint on the maximization. This constraint is related to the truncated function basis mentioned before. Particularly, the authors force f to be spanned by the first k eigenfunctions of the Laplace-Beltrami operator on N and generate a collection of maximizers corresponding to different k's.

Contributions

This thesis consists of two major parts: the first part (chapter 2) investigates the problem of metric reconstruction for filamentary structures and the second part (chapter 3 and chapter 4) studies the stability properties of the shape difference operators and presents an empirical method to perform the functional-map-based frameworks directly on point cloud data.

Chapter 2: Particularly, given a geodesic space 2 (X, d X ), we study the Reeb graph of the distance function to a point. A distance function to x ∈ X is defined as d(y) = d X (x, y), which is the distance in X from x to y. The same Reeb graph has been studied in [Ge 2011] as a tool for capturing a 1-dimension skeleton of some general data. We further define a metric on the Reeb graph of the distance function so that it's as well a metric graph. We prove that this metric graph is well-approximating X in the sense of the Gromov-Hausdorff distance.

From the practical point of view, data usually come as discrete points sampled from the underlying metric space X. We incorporate the idea of the mapper algorithm to construct a variant of the Reeb graph called α-Reeb graph, with which it is more convenient to deal with the discrete sampling points in practice. And theoretically, the α-Reeb graph not only enjoys similar guarantees on the metric approximation as the Reeb graph, but also comes with certain topological guarantees.

Overall, we conclude that the α-Reeb graph is a reliable approximation of X in both geometrical and topological sense. The above are results of a collaboration with F. Chazal and J. Sun, which have been published in Discrete and Computational Geometry [START_REF] Chazal | [END_REF]].

Our framework assumes that the input X is already a (discrete) metric space. In many cases, the raw data are discrete isolated points sampled from some underlying metric space X. It's then appealing to study how accurately we can recover the metric of X with the sampling points. In this chapter we study in a special case where X is a metric graph embedded in an Euclidean space. We assure that under some technique conditions on X and the sampling density, the metric we recover from the sampling points is well-approximating the ground truth.

Chapter 3: As mentioned above, the second problem considered in this thesis is about comparing geometrical objects associated by a map. We notice that similar approaches have been made in the area of computer graphics. Particularly, the map-based shape analysis frameworks based on the functional map seem promising approaches to this problem. Though the original formulations are proposed for shapes, i.e., 2-dimensional Riemannian manifolds, they can naturally be defined for Riemannian manifolds of arbitrary dimension.

In this chapter, we concentrate on stability analysis of two frameworks-the shape difference operators and the map analysis and visualization from [START_REF] Ovsjanikov | [END_REF]]. More generally, we assume that the inputs are two n-dimensional Riemannian manifolds M, N and a map T : M → N . We verify two type of stability properties: one is with respect to perturbations on the input manifolds and the other is with respect to the changing scale, which is peculiar to the latter framework.

The results of this chapter are obtained in collaboration with F. Chazal and M. Ovsjanikov. The manuscript is soon to be submitted.

Chapter 4: The functional-map-based frameworks are usually proposed to analyze 3D shapes, i.e., 2dimensional Riemannian manifold embedded in R 3 . The 3D shapes can be nicely approximated and represented by polygon meshes. Implementing these frameworks on meshes are obvious and efficient. In a more general setting, where the input are Riemannian manifolds of dimension higher than 2 or embedded in a Euclidean of dimension higher than 3, the implementation is not straightforward any more.

In this chapter, we propose a pipeline for dealing with data in the most primitive form-point cloud data (PCD). The idea is to construct counterparts in the PCD setting of the ingredients necessary in the mesh setting. To test our method, we perform the functional-map-based frameworks with our pipeline on point clouds sampled from 3D shapes. The results are compared to the ones obtained from the mesh setting. Empirically, the pipeline works well with the test data and shows robustness with respect to noisy data.

At the end, we emphasize that the work in this chapter is exploratory and experimental. The experimental results suggest that this direction is worth further exploration and requires corresponding theoretical analyses.

The work presented in this chapter is obtained in collaboration with F. Chazal and M. Ovsjanikov. The manuscript is in preparation.

CHAPTER 3

Approximation for Filamentary Structures

Using Reeb-Type Graphs

Introduction

With the advance of sensor technology, computing power and the Internet, massive amounts of geometric data are being generated and collected in various areas of science, engineering and business. As they are becoming widely available, there is a real need to analyze and visualize these large scale geometric data to extract useful information out of them. In many cases these data are not embedded in Euclidean spaces and come as (finite) sets of points with pairwise distance information, i.e. (discrete) metric spaces. A large amount of research has been done on dimensionality reduction, manifold learning and geometric inference for data embedded in (possibly high dimensional) Euclidean spaces and assumed to be concentrated around low dimensional manifolds [START_REF] Belkin | Laplacian Eigenmaps for Dimensionality Reduction and Data Representation[END_REF], Lafon 2004, Tenenbaum 2000b]. However, the assumption of data lying on a manifold may fail in many applications. In addition, the strategy of representing data by points in an Euclidean space may introduce large metric distortions as the data may lie in highly curved spaces, instead of in flat Euclidean space raising many difficulties in the analysis of metric data. In the past decade, with the development of topological methods in data analysis, new theories such as topological persistence (see, for example, [START_REF] Edelsbrunner | [END_REF][START_REF] Zomorodian | Computing persistent homology[END_REF], Carlsson 2009, Chazal 2012]) and new tools such as the Mapper algorithm [Singh 2007] have given rise to new algorithms to extract and visualize geometric and topological information from metric data without the need of an embedding into an Euclidean space. In this chapter we focus on a simple but important setting where the underlying geometric structure approximating the data can be seen as a branching filamentary structure i.e., more precisely, as a metric graph which is a topological graph endowed with a length assigned to each edge. Such structures appear naturally in various real-world data such as collections of GPS traces collected by vehicles on a road network, earthquakes distributions that concentrate around geological faults, distributions of galaxies in the universe, networks of blood vessels in anatomy or hydrographic networks in geography just to name a few. It is thus appealing to try to capture such filamentary structures and to approximate the data by metric graphs that will summarize the metric and allow convenient visualization.

Overview

In this chapter we address the metric reconstruction problem for filamentary structures. The input of our method and algorithm is a metric space (X, d X ) that is assumed to be close with respect to the so-called Gromov-Hausdorff distance d GH to a much simpler, but unknown, metric graph (G ′ , d G ′ ). Our algorithm outputs a metric graph (G, d G ) that is proven to be close to (G ′ , d G ′ ) in both geometry and topology. Our approach relies on the notion of Reeb graph (and some variants of it introduced in Section 3.4) and our main theoretical results are stated in the following two theorems.

Theorem 3.5 [Recovery of Geometry]. Let (X, d X ) be a compact connected geodesic space, let r ∈ X be a fixed base point such that the metric Reeb graph (G, d G ) of the function d = d X (r, .) : X → R is a finite graph. If for a given ε > 0 there exists a finite metric graph

(G ′ , d G ′ ) such that d GH (X, G ′ ) < ε then we have d GH (X, G) < 2(β 1 (G) + 1)(17 + 8N E,G ′ (8ε))ε
where N E,G ′ (8ε) is the number of edges of G ′ of length at most 8ε and β 1 (G) is the first Betti number of G, i.e. the number of edges to remove from G to get a spanning tree. In particular if X is at distance less than ε from a metric graph with shortest edge larger than 8ε then d GH (X, G) < 34(β 1 (G) + 1)ε.

Note that β 1 (G) ≤ β 1 (X) and thus d GH (X, G) is upper bounded by the quantities depending only on the input X.

Theorem 3.7 [Recovery of Topology]. Let (X, d X ) be a compact connected path metric space and

(G ′ , d G ′ ) is a metric graph so that d GH (X, G ′ ) < ε. Let r ∈ X, α > 60ε and I = {[0, 2α), (iα, (i + 2)α)|1 ≤ i ≤ m} covers the segment [0, Diam(X)] such that the 2α-Reeb graph G associated to I and the function d = d X (r, .) : X → R is a finite graph.
If no edges of G ′ are shorter than L and no loops of G ′ are shorter than 2L with L ≥ 32α + 9ε, then we have G and G ′ are homotopy equivalent.

Approximating the Reeb graph (G, d G ) from a neighborhood graph is usually not obvious. If we compute the Reeb graph of the distance function to a given point defined on the neighborhood graph we obtain the neighborhood graph itself and do not achieve our goal of representing the input data by a simple graph. See Table 3.1. It is then appealing to build a two dimensional complex having the neighborhood graph as 1-dimensional skeleton and use the algorithm of [START_REF] Harvey | [END_REF], Parsa 2013] to compute the Reeb graph of the distance to the root point. Unfortunately adding triangles to the neighborhood graph may widely change the metric between the data points on the resulting complex and significantly increase the complexity of the algorithm. We overcome this issue by introducing a variant of the Reeb graph, the α-Reeb graph, inspired from [Singh 2007] and related to the recently introduced notion of graph induced complex [Dey 2013a], that is easier to compute than the Reeb graph but also comes with approximation guarantees (see Theorem 3.6). As a consequence our algorithm runs in almost linear time (see Section 3.8).

Raw data usually do not come as geodesic spaces. They are given as discrete sets of points (and thus not connected metric spaces) sampled from the underlying space (X, d X ). Moreover in many cases only distances between nearby points are known. A geodesic space (see Section 3.3.1 for a definition of geodesic space) can then be obtained from these raw data as a neighborhood graph where nearby points are connected by edges whose length is equal to their pairwise distance. The shortest path distance in this graph is then used as the metric. In our experiments we use this new metric as the input of our algorithm. Particularly, in Section 3.7 we study a special case that the underlying metric space (X, d X ) is a metric graph embedded in a d-dimensional Euclidean space R d and data are sampled from (X, d X ). Our result shows that the metric of a certain neighboring graph built on the sampling points well approximates d X under proper conditions on X and the sampling density.

Related Works

Approximation of data by 1-dimensional geometric structures has been considered by different communities. In statistics, several approaches have been proposed to address the problem of detection and extraction of filamentary structures in point cloud data. For example Arial-Castro et al [Arias-Castro 2006] use multiscale anisotropic strips to detect linear structure while [START_REF] Genovese | [END_REF], Genovese 2012] and more recently [Genovese 2014] base their approach upon density gradient descents or medial axis techniques. These methods apply to data corrupted by outliers embedded in Euclidean spaces and focus on the inference of individual filaments without focus on the global geometric structure of the filaments network.

In computational geometry, the curve reconstruction problem from points sampled on a curve in an Euclidean space has been extensively studied and several efficient algorithms have been proposed [START_REF] Amenta | [END_REF], Dey 2000, Dey 2001]. Unfortunately, these methods restrict to the case of simple embedded curves (without singularities or self-intersections) and hardly extend to the case of topological graphs. In a more intrinsic setting where data come as finite abstract metric spaces, [Aanjaneya 2012] propose an algorithm that outputs a topologically correct (up to a homeomorphism) reconstruction of the approximated graph. However this algorithm requires some tedious parameters tuning and relies on quite restrictive sampling assumptions. When these conditions are not satisfied, the algorithm may fail and not even outputs a graph. Compared to the algorithm of [Aanjaneya 2012], our algorithm not only comes with metric guarantees but also whatever the input data is, it always outputs a metric graph and does not require the user to choose any parameters. Closely related to our approach is the data skeletonization algorithm proposed in [Ge 2011] that computes the Reeb graph of an approximation of the distance function to a root point on a 2-dimensional complex built on top of the data whose size might be significantly larger than a neighboring graph. The algorithm of [Ge 2011] also always outputs a graph but it does not come with metric guaranties. Recently, Bauer, Ge and Wang [Bauer 2014] define a metric based on the function for Reeb graph and show it is stable under Gromov-Hausdorff distance. The implementation of our algorithm relies on the Mapper algorithm [Singh 2007], that provides a way to visualize data sets endowed with a real valued function as a graph, where the considered function is the distance to the chosen root point. However, unlike the general Mapper algorithm, our methods provides an upper bound on the Gromov-Hausdorff distance between the reconstructed graph and the underlying space from which the data points have been sampled.

In theoretical computer science, there is much of work on approximating metric spaces using trees [START_REF] Bǎdoiu | [END_REF], Abraham 2007[START_REF] Chepoi | [END_REF] or distribution of trees [Dhamdhere 2006[START_REF] Fakcharoenphol | [END_REF] where the trees are often constructed as spanning trees possibly with Steiner points. Our approach is different as our reconstructed graph or tree is a quotient space of the original metric space where the metric only gets contracted (see Proposition 3.2). Finally we remark that the recovery of filament structure is also studied in various applied settings, including road networks [Chen 2010[START_REF] Tupin | Detection of linear features in SAR images: application to road network extraction[END_REF]], galaxies distributions [Choi 2010].

Preliminaries

Metric Spaces

A metric space is a pair (X, d X ) where X is a set and d X : X × X → R is a non-negative map such that for any x, y, z ∈ X, d X (x, y) = 0 if and only if

x = y, d X (x, y) = d X (y, x) and d X (x, z) ≤ d X (x, y) + d X (y, z).
A continuous path (or curve) σ, is a continuous map from I, a real interval, to X. A path is called simple if it is not self-intersecting, or equivalently σ is an injective map. The length of a path can be induced by the metric space (X, d X )

Definition 3.1 Let I = [a, b] be an interval. A partition of I is a finite sequence of point in I such that a = t 0 ≤ t 1 ≤ t 2 ≤ • • • ≤ t k = b.
Let Γ be the set of all partitions of I. Then the length of σ, L(σ) is defined as

L(σ) = sup {t 0 ,t 1 ,...,t k }⊂Γ k i=1 d X (σ(t i-1 ), σ(t i ))
We call (X, d X ) a path metric space if the distance between any pair of points is equal to the infimum of the length of the continuous curves joining them. Hereinafter we only consider compact path metric spaces, which are geodesic according to the following theorem (see also [Gromov 2002]).

Theorem 3.1 [Hopf-Rinow Theorem] If (X, d X
) is a complete, locally compact path metric space, then:

• Closed balls are compact, or ,equivalently, each bounded, closed domain is compact.

• Each pair of points can be joined by a minimizing geodesic.

A non-geodesic metric space is illustrated in Figure 3.1.

(0, 0) (1.0) Figure 3.1: A metric space that is not geodesic: let X = R 2 \( 1 2 , 0)
, and d X be the euclidean metric. The distance between point (0, 0) and (1, 0) is 1, however there does not exist a continuous curve joining them of lenth 1.

Two compact metric spaces (X, d X ) and (Y, d Y ) are isometric if there exits a bijection φ : X → Y that preserves the distances, i.e., for any x, x ′ ∈ X, d Y (φ(x), φ(x ′ )) = d X (x, x ′ ). The set of isometry classes of compact metric spaces can be endowed with the Gromov-Hausdorff distance that can be defined using the following notion of correspondence (see [START_REF] Burago | [END_REF]).

Definition 3.2 Let (X, d X ) and (Y, d Y ) be two compact metric spaces. Given ε > 0, an ε-correspondence between (X, d X ) and (Y, d Y ) is a subset C ⊂ X × Y such that: i) for any x ∈ X there exists y ∈ Y such that (x, y) ∈ C; ii) for any y ∈ Y there exists x ∈ X such that (x, y) ∈ C; iii) for any (x, y), (x ′ , y ′ ) ∈ C, |d X (x, x ′ ) -d Y (y, y ′ )| ≤ ε. Definition 3.3 The Gromov-Hausdorff distance between two compact metric spaces (X, d X ) and (Y, d Y ) is defined by d GH (X, Y ) = 1 2 inf{ε ≥ 0 :
there exists an ε-correspondence between X and Y } It follows obviously from the definitions that two metric spaces are isometric if and only if the Gromov-Hausdorff distance between them is zero.

Neighboring Graphs on Point Clouds in R d In many applications, many metric spaces are isometrically embedded in R d . Sampling from an embedded metric space results in a collection of points

X = {x 1 , x 2 , • • • , x n } ⊂ R d .
No information about the underlying metric space can be inferred from solely the isolated points, thus a first step towards processing X is to build connections between points. Usually, we build a neighboring graph on top of X by connecting points being close to each other in a certain criterion. Now that X is a subset of R d , the Euclidean distance is available for measuring closeness between points. Hereinafter we denote the Euclidean distance between x and y by ∥x -y∥. Based on the Euclidean distances, two typical neighboring graph constructions are defined as the following Definition 3.4 Given an integer k, we denote the set of the nearest k neighborhoods of x i among X\x i by N (x i , k). The k-nearest graph on top of X is then constructed by connecting x i to x j if and only if

x i ∈ N (x j , k) or x j ∈ N (x i , k).
Additionally, the mutual k-nearest graph is constructed by connecting x i to x j if and only if x i ∈ N (x j , k) and x j ∈ N (x i , k). Definition 3.5 Given a positive constant δ, the δ-neighborhood graph of X is constructed by connecting x i to x j if and only if ∥x ix j ∥ ≤ ε.

There certainly exists many other graph constructions for a point cloud. Starting with a graph R = ( X, E) embedded in R d , a metric on X is induced by the graph and we denote it by d R , whose definition is actually a discrete version of the definition of curve length 3.1.

Definition 3.6 Given a pair of points x, y ∈ X ⊂ R d , let P = (x 1 , x 2 , • • • , x i ) be a path in R = ( X, E) such that x 1 = x, x i = y and (x j , x j+1 ) ∈ E, ∀1 ≤ j ≤ i -1. A metric d R : X × X → [0, +∞) is then defined as: d R (x, y) = inf P ∈Γ i-1 j=1 ∥x j -x j+1 ∥
where Γ is the set of all possible paths connecting x to y in R.

It follows from the above definition and the triangle inequality that d R (x, y) ≥ ∥x -y∥ for any x, y ∈ X.

Metric Graphs A graph G = (V, E) is obtained by taking a finite set of vertices, V , and joining some of them by edges, which are elements of the edge set E. A metric segment of length a (a > 0) is a metric space isometric to an interval [0, a] ⊂ R. A metric graph is then a graph whose edges are metric segments. Equivalently, a metric graph is also seen as a graph with a length assigned to each of its edges. Such a length assignment naturally induces a metric on G, we then denote the metric by d G : G × G → R. For any two points g 1 , g 2 ∈ G (they are not necessarily elements of V ), d G (g 1 , g 2 ) is the length of the shortest path between them along edges of G.

As we will mention soon, the goal of this part of our work is to approximate filamentary metric spaces with certain graphs. Especially, we consider filamentary metric spaces that are close to metric graphs in Gromov-Hausdorff distance (see Definition 3.3).

Topology

After introducing the geometrical notions, we give a brief review of topological notions that will be involved in Section 3.6.

It is well-known that a metric on a space induces topological properties such as open sets. In fact, let B(x, r) be the open ball centered at x ∈ X with radius r, i.e., B(x, r) = {y : d X (x, y) < r}. The set of all the open balls in X forms a basis of a topology on space X. Denote the topology by τ X , then (X, τ X ) is a topological space induced by the metric d X .

Homotopy Equivalence We first introduce the homotopic relationship between functions from a topological space (X, τ X ) to (Y, τ Y ).

Definition 3.7 Two continuous functions f 1 , f 2 from X to Y are homotopic if and only if there exists a continuous function

H : X × [0, 1] → Y , such that H(x, 0) = f 1 (x) and H(x, 1) = f 2 (x)
Then we define homotopy equivalence between topological spaces. Definition 3.8 Two topological spaces X and Y are homotopy equivalent if and only if there exists two continuous functions f :

X → Y and g : Y → X such that f • g (resp., g • f ) is homotopic to id Y (resp., id X ), which is the identity map.
Especially, if a topological space X is homotopy equivalent to a one-point space Y = {y}, then we call X a contractible space. From the point view of homotopy, cycles are the only non-trivial topological features in a graph. In fact, one can classify graphs with the number of circles. A simple illustration of the following theorem is in Figure 3 The number of cycles of a graph, |E| -|V | + 1, is also known as the first betti number β 1 (G) of G ,i.e., the rank of the first homology group of G. Roughly speaking, one need to remove at least β 1 (G) edges from G to obtain an acyclic graph, i.e., a tree.

Abstract Simplicial Complex An abstract simplicial complex ∆ is a family of non-empty sets that satisfies if A ∈ ∆ and B ⊂ A, then B ∈ ∆.

Covering and Nerve Let V = {V j } j∈J be a finite open covering of space X, meaning that X = j∈J V j and J is a finite indexing set.

The nerve of the covering V is a simplicial complex N (V) whose vertex set is the indexing set J, and

vertices {j 1 , j 2 , • • • , j k+1 } span a k-simplex in N (V) if and only if V j 1 ∩ V j 2 ∩ • • • ∩ V j k+1 ̸ = ∅.
The following Nerve lemma asserts that the homotopy type of X is possibly recovered by a good cover.

Lemma 3.1 [Nerve Lemma] Let V = {V j } j∈J be a finite open covering of X.
If the intersection of any subset of V is either empty or contractible, then the nerve of V, N (V), is homotopy equivalent to the union of X.

Reeb-type Graph

In this section, we describe a construction to build a Reeb-type graph for approximating a metric space (X, d X ). Let (X, d X ) be a compact geodesic space and let r ∈ X be a fixed base point. Let d : X → R be the distance function to r, i.e., d(x) = d X (r, x). The set of points in X with distance to r equal to d(x) is then the pre-image of d(x), i.e., d -1 (d(x)). G is defined by: for any g, g ′ ∈ G, g ≤ G g ′ if and only if d * (g) ≤ d * (g ′ ) and there exist a continuous path γ in G connecting g to g ′ such that d • γ is non decreasing, makes G a partially ordered set.

The

The α-Reeb graphs. Computing or approximating the Reeb graph of (X, d) from a finite set of point sampled on X is usually a difficult task. To overcome this issue we also consider a variant of the Reeb graph that shares very similar properties to the Reeb graph. This variant is closely related to the graph constructed by the Mapper algorithm introduced in [Singh 2007] making its computation much easier than the Reeb graph (see Section 3.8).

Let α > 0 and let I = {I i } be a covering of the range of d by open intervals of length at most α. The transitive closure of the relation x ∼ α y if and only if d(x) = d(y) and x, y are in the same path connected component of d -1 (I i ) for some interval I i ∈ I is an equivalence relation that is also denoted by ∼ α . The quotient space G α = X/ ∼ α is called the α-Reeb graph 1 of d.

We denote by π : X → G α the quotient map. Notice that π is continuous and as

X is path connected, G α is path connected. The function d induces a function d * : G α → R + that satisfies d = d * • π. The relation on G α is defined by: for any g, g ′ ∈ G α , g ≤ Gα g ′ if and only if d * (g) ≤ d * (g ′ ) and there exist a continuous path γ in G α connecting g to g ′ such that d • γ is non decreasing, makes G α a partially ordered set.
Notice that without making assumptions on X and d, in general G and G α are not finite graphs. However when the number of path connected components of the level sets of d is finite and changes only a finite number of times then the Reeb graph turns out to be a finite directed acyclic graph. Similarly, when the covering of X by the connected components of d -1 (I i ), I i ∈ I is finite, the α-Reeb graph also turns out to be a finite directed acyclic graph. This happens in most applications and for example when (X, d X ) is a finite simplicial complex or a compact semialgebraic (or more generally a compact subanalytic space) with d being semi-algebraic (or subanalytic).

All the results and proofs presented in Section 3.5 are exactly the same for the Reeb and the α-Reeb graphs. In the following paragraph and in Section 3.5, G denotes indifferently the Reeb graph or an α-Reeb graph for some α > 0. We also always assume that X and d (and α and I) are such that G is a finite graph.

Metrics on Reeb and α-Reeb Graphs.

Let us define the set of vertices V of G as the union of the set of points of degree not equal to 2 with the set of local maxima of d * over G, and the base point π(r). The set of edges E of G is then the set of the connected components of the complement of V . Notice that π(r) is the only local (and global) minimum of d * : since X is path connected, for any x ∈ X there exists a geodesic γ joining r to x along which d is increasing; d * is thus also increasing along the continuous curve π(γ), so π(x) cannot be a local minimum of d * . As a consequence d * is monotonic along the edges of G. We can thus assign an orientation to each edge: if e = [p, q] ∈ G is such that d * (p) < d * (q) then the positive orientation of e is the one pointing from p to q. Finally, we assign a metric to G. Each edge e ∈ E is homeomorphic to an interval to which we assign a length equal to the absolute difference of the function d * at two endpoints. The distance between two points p, p ′ of e is then |d * (p)d * (p ′ )|. This makes G a metric graph (G, d G ) isometric to the quotient space of the union of the intervals isometric to the edges by identifying the endpoints if they correspond to the same vertex in G. Note that d * is continuous in (G, d G ) and for any p ∈ G, d * (p) = d G (π(r), p). Indeed this is a consequence of the following lemma. Proof: As d * • δ is strictly increasing, when δ enters an edge e by one of its end points, either it exits at the other end point or it stops at p ′ if p ′ ∈ e. Moreover δ cannot go through a given edge more than one time. As a consequence δ can be decomposed in a finite sequence of pieces e 0 = [p,

p 1 ], e 1 = [p 1 , p 2 ], • • • , e n-1 = [p n-1 , p n ], e n = [p n , p ′
] where e 0 and e n are the segments joining p and p ′ to one of the endpoint of the edges that contain them and e 1 , • • • , e n-1 are edges. So, the length of δ is equal to

(d * (p 1 ) -d * (p)) + (d * (p 2 ) - d * (p 1 )) + • • • + (d * (p ′ ) -d * (p n )) = d * (p ′ ) -d * (p) and d G (p, p ′ ) ≤ d * (p ′ ) -d * (p).
Similarly any simple path joining p to p ′ can be decomposed in a finite sequence of pieces e

′ 0 = [p, p ′ 1 ], e ′ 1 = [p ′ 1 , p ′ 2 ], • • • , e ′ k-1 = [p ′ k-1 , p ′ k ], e ′ k = [p ′ k , p ′
] where e ′ 0 and e ′ k are the segments joining p and p ′ to one of the endpoint of the edges that contain them, and e ′ 1 , • • • , e ′ k-1 are edges. Now, as we do not know that d * is increasing along this path, its length is thus equal to |d

* (p ′ 1 ) -d * (p)| + |d * (p ′ 2 ) -d * (p ′ 1 )| + • • • + |d * (p ′ ) -d * (p ′ n )| ≥ d * (p ′ ) -d * (p). So, d G (p, p ′ ) ≥ d * (p ′ ) -d * (p).

Recovery of Geometry

The goal of this section is to provide an upper bound of the Gromov-Hausdorff distance between X and G, and we conclude our result as the following theorem: Theorem 3.5 [Recovery of Geometry]. Let (X, d X ) be a compact connected geodesic space, let r ∈ X be a fixed base point such that the metric Reeb graph (G, d G ) of the function d = d X (r, .) : X → R is a finite graph. If for a given ε > 0 there exists a finite metric graph

(G ′ , d G ′ ) such that d GH (X, G ′ ) < ε then we have d GH (X, G) < 2(β 1 (G) + 1)(17 + 8N E,G ′ (8ε))ε
where N E,G ′ (8ε) is the number of edges of G ′ of length at most 8ε and β 1 (G) is the first Betti number of G, i.e. the number of edges to remove from G to get a spanning tree. In particular if X is at distance less than ε from a metric graph with shortest edge larger than 8ε then d GH (X, G) < 34(β 1 (G) + 1)ε.

We decompose our proofs into two parts: the first one (Theorem 3.3) asserts that the Gromov-Hausdorff distance between X and G only depends on the first Betti number β 1 (G) of G and the maximal diameter M of the level sets of π; and the rest part that estimates an upper bound of M is given in Section 3.5.1

Theorem 3.3 d GH (X, G) < (β 1 (G) + 1)M where d GH (X, G) is the Gromov-Hausdorff distance between X and G, β 1 (G)
is the first Betti number of G and M = sup p∈G {diam(π -1 (p))} is the supremum of the diameters of in the level sets of π.

Remark that as β 1 (G) ≤ β 1 (X), from the above theorem, d GH (X, G) is upper bounded by the quantities depending only on the input X. The proof of Theorem 3.3 is deduced from two propositions comparing the distances between pairs of points x, y ∈ X and their images π(x), π(y) ∈ G whose proofs rely on the notion of merging vertex. A vertex v ∈ V is called a merging vertex if it is the end point of at least two edges e 1 and e 2 that are pointing to it according to the orientation defined in Section 3.4. Geometrically this means that there are at least two distinct connected components of π -1 (d -1 * (d * (v)ε)) that accumulate to π -1 (v) as ε > 0 goes to 0. The set of merging vertices is denoted by V m . We have Lemma 3.3 The cardinality of V m is at most β 1 (G) where β 1 (G) is the rank of the first homology group of G.

Proof: The result follows from classical persistence homology theory [Edelsbrunner 2010]. First remark that, as π(r) is the only local minimum of d * , the sublevel sets of the function d * : G → R + are all path connected. Indeed if π(x), π(y) ∈ G are in the same sublevel set d -1 * ([0, α]), α > 0, then the images by π of the shortest paths in X connecting x to r and y to r are contained in d -1 * ([0, α]) and their union is a continuous path joining π(x) to π(y). As a consequence, the 0-dimensional persistence of d * is trivial. So as we increase the α value, no merging vertices serve as connecting two different connected components. Thus, each merging vertex in V m creates at least a cycle that never dies as G is one dimensional and does not contain any 2-dimensional simplex. Thus

|V m | ≤ β 1 (G).
The following lemma shows that a shortest path in G is the projection of a shortest path in X as long as it does not meet a merging vertex and allow to prove proposition 3.1 below.

Lemma 3.4 Let p, p ′ ∈ G and let δ : [d * (p), d * (p ′ )] → G be a strictly increasing path going from p to p ′ that does not contain any point of V m in its interior. Then for any

x ′ ∈ π -1 (p ′ ) cl(π -1 (δ(d * (p), d * (p ′ )))
where cl(.) denotes the closure, there exists a shortest path γ connecting a point

x of π -1 (p) to x ′ such that π(γ) = δ and d X (x, x ′ ) = d(x ′ ) -d(x) = d * (p ′ ) -d * (p) = d G (p, p ′ ).
Proof: First assume that p ′ is not a merging point. Let γ 0 : [0, d(x ′ )] → X be any shortest path between r and x ′ and let γ be the restriction of

γ 0 to [d * (p), d(x ′ )] = [d * (p), d * (p ′ )]. If the infimum t 0 of the set I = {t ∈ [d * (p), d * (p ′ )] : π(γ(t ′ )) ∈ δ, ∀t ′ ≥
t} is larger than d * (p), then there exists an increasing sequence (t n ) that converges to t 0 such that π(γ(t n )) ̸ ∈ δ. As a consequence δ(t 0 ) is a merging point; a contradiction. So t 0 = d * (p) and γ(d * (p)) intersects π -1 (p) at a point x.

Now if p ′ is a merging point, as x ′ is chosen in the closure of π -1 (δ(d * (p), d * (p ′ )), for any sufficiently large n ∈ N one can consider a sequence of points x ′ n ∈ π -1 (δ(d * (p ′ ) -1/n)) that converges to x ′ and apply the first case to get a sequence of shortest path γ n from a point x n ∈ π -1 (p) and x ′ n . Then applying Arzelà-Ascoli's theorem (see [Dieudonné 1960] 7.5) we can extract from γ n a sequence of points converging to a shortest path γ between a point x ∈ π -1 (p) and x ′ .

To conclude the proof, notice that from Lemma 3.

2 we have d G (p, p ′ ) = d * (p ′ ) -d * (p) = d(x ′ ) -d(x).
Since γ is the restriction of a shortest path from r to x we also have d

X (x, x ′ ) = d(x ′ ) -d(x).
Notice that from Lemma 3.2, δ is a shortest path and the parametrization by the interval [d * (p), d * (p ′ )] can be chosen to be an isometric embedding.

Proposition 3.1 For any x, y ∈ X we have

d X (x, y) ≤ d G (π(x), π(y)) + 2(β 1 (G) + 1)M
where M = sup p∈G {diam(π -1 (p))} and β 1 (G) is the first Betti number of G.

Proof: Let δ be a shortest path between π(x) and π(y). Remark that except at the points π(x) and π(y) the local maxima of the restriction of d * to δ are in V m . Indeed as δ is a shortest path it has to be simple, so if p ∈ δ is a local maximum then p has to be a vertex and δ has to pass through two edges having p as end point and pointing to p according to the orientation defined in 3.4. So p is a merging point.

Since δ is simple and V m is finite, δ can be decomposed in at most |V m | + 1 connected paths along the interior of which the restriction of d * does not have any local maxima. So along each of these connected paths the restriction of d * can have at most one local minimum. As a consequence, δ can be decomposed in a finite number of continuous paths

δ 1 , δ 2 , • • • , δ k with k ≤ 2(|V m | + 1)
, such that the restriction of d * to each of these path is strictly monotonic. For any i ∈ {1, • • • , k} let p i and p i+1 the end points of δ i with p 1 = π(x) and p k+1 = π(y). We can apply Lemma 3.4 to each δ i to get a shortest path γ i in X between a point x i ∈ π -1 (p i ) and a point in

y i+1 ∈ π -1 (p i+1 ) such that π(γ i ) = δ i and d X (x i , y i+1 ) = d G (p i , p i+1 ).
The sum of the lengths of the paths γ i is equal to the sum of the lengths of the path δ i which is itself equal to d G (π(x), π(y)). Now for any i ∈ {1, • • • , k}, since π(x i ) = π(y i ) we have d X (x i , y i ) ≤ M and x i and y i can be connected by a path of length at most M (x 1 is connected to x and y k+1 is connected to y. Gluing these paths to the paths γ i gives a continuous path from x to y whose length is at most

d G (π(x), π(y)) + kM ≤ d G (π(x), π(y)) + 2(|V m | + 1)M . Since from Lemma 3.3, |V m | ≤ β 1 (G), we finally get that d X (x, y) ≤ d G (π(x), π(y)) + 2(β(G) + 1)M . Proposition 3.2 The map π : X → G is 1-Lipschitz: for any x, y ∈ X we have d G (π(x), π(y)) ≤ d X (x, y).
Proof: Let x, y ∈ X and let γ : I → X be a shortest path from x to y in X where I ⊂ R is a closed interval. The path π(γ) connects π(x) and π(y) in G.

We first claim that there exists a continuous path Γ contained in π(γ) connecting π(x) and π(y) that intersects each vertex of G at most one time. The path Γ can be defined by iteration in the following way.

Let v 1 , • • • v n ∈ V be the vertices of G that are contained in π(γ) \ {π(x), π(y)} and let Γ 0 = π(γ) : J 0 → G, J 0 = I. For i = 1, • • • n let t - i = inf{t : Γ i-1 (t) = v i } and t + i = sup{t : Γ i-1 (t) = v i } and define Γ i as the restriction of Γ i-1 to J i = J i-1 \ (t - i , t + i ).
The path Γ i is a connected continuous path (although J i is a disjoint union of intervals) that intersects the vertices v 1 , v 2 , • • • , v i at most one time. We then define Γ = Γ n : J = J n → G where J ⊂ I is a finite union of closed intervals. Notice that Γ is the image by π of the restriction of γ to J and that

Γ(t) ∈ {v 1 , • • • v n } only if t is one of the endpoints of the closed intervals defining J. Now, for each connected component [t, t ′ ] of J, γ((t, t ′ )) is contained in π -1 (e)
where e is the edge of G containing Γ([t, t ′ ]). As a consequence,

d G (π(γ)(t), π(γ)(t ′ )) = |d * (π(γ)(t)) -d * (π(γ)(t ′ ))| = |d(γ(t)) -d(γ(t ′ ))|.
Recalling that d(γ(t)) = d X (r, γ(t)) and d(γ(t ′ )) = d X (r, γ(t ′ )) and using the triangle inequality we get that |d(γ(t))d(γ(t ′ ))| ≤ d X (γ(t), γ(t ′ )). To conclude the proof, since γ is a geodesic path we just need to sum up the previous inequality over all connected components of J:

d X (x, y) ≥ [t,t ′ ]∈cc(J) d X (γ(t), γ(t ′ )) ≥ [t,t ′ ]∈cc(J) d G (π(γ)(t), π(γ)(t ′ )) ≥ d G (π(x), π(y))
where cc(J) is the set of connected components of J.

The proof of Theorem 3.3 now easily follows from Propositions 3.1 and 3.2.

Proof: (of Theorem 3.3)

Consider the set C = {(x, π(x)) : x ∈ X} ⊂ X × G.
As π is surjective this is a correspondence between X and G. It follows from Propositions 3.1 and 3.2 that for any (x, π(x)), (y, π(y)) ∈ C,

|d X (x, y) -d G (π(x), π(y))| ≤ 2(β 1 (G) + 1)M So C is a 2(β 1 (G) + 1)M -correspondence and d GH (X, G) ≤ (β 1 (G) + 1)M .

Bounding the Diameter M

The two following lemmas allow to bound M , the diameter of the level sets of π.

Lemma 3.5 Let (G, d G ) be a connected finite metric graph and let r ∈ G. We denote by d r = d G (r, .) : G → [0, +∞) the distance to r. For any edge e ⊂ G, the restriction of d r to e is either strictly monotonic or has only one local maximum. Moreover the length l = l(e) of e is upper bounded by two times the difference between the maximum and the minimum of d r restricted to e.

Proof: Let l be the length of E and let t → e(t), t ∈ [0, l], be an arc length parametrization of E. Since E is an edge of G, for t ∈ [0, l] any shortest geodesic γ t joining r to e(t) must contain either x 1 = e(0) or x 2 = e(l). If it contains x 1 then for any t ′ < t the restriction of γ t between r and e(t ′ ) is a shortest geodesic containing x 1 and if it contains x 2 then for any t ′ > t the restriction of γ t between r and e(t ′ ) is a shortest geodesic containing x 2 . Moreover in both cases, the function d r is strictly monotonic along γ. As a consequence, the set I 1 = {t ∈ [0, l] : a shortest geodesic joining r to e(t) contains x 1 } is a closed interval containing 0. Similarly the set I 2 = {t ∈ [0, l] : a shortest geodesic joining r to e(t) contains x 2 } is a closed interval containing l and [0, l] = I 1 ∪ I 2 . Moreover d r is strictly monotonic on e(I 1 ) and on e(I 2 ). As a consequence I 1 ∩ I 2 is reduced to a single point t 0 that has to be the unique local maximum of d r restricted to E.

The second part of the lemma follows easily from the previous proof: the minimum of d r restricted to E is attained either at x 1 or x 2 and d r (e(t 0 )) = d r (x 1 ) + t 0 = d r (x 2 ) + lt 0 is the maximum of d r restricted to E. We thus obtain that 2t Proof: Let x, y ∈ B and let t → γ(t) ∈ B be a continuous path joining x to y in B. Let E be an edge of G that does not contain x or y and with end points x 1 , x 2 such that γ intersects the interior of E. Then γ -1 (E) is a disjoint union of closed intervals of the form I = [t, t ′ ] where γ(t) and γ(t ′ ) belong to the set {x 1 , x 2 }. If γ(t) = γ(t ′ ) we can remove the part of γ between t and t ′ and still get a continuous path between x and y. So without loss of generality we can assume that if γ intersects the interior of E, then E is contained in γ. Using the same argument as previously we can also assume that if γ goes across E, it only does it one time, i.e. γ -1 (E) is reduced to only one interval. As a consequence, γ can be decomposed in a sequence y] are pieces of edges containing x and y respectively and The example of Figure 3.3 shows that the bound of Lemma 3.3 is tight. Lemma 3.6 Let X and Y be compact geodesic metric spaces and C ⊂ X × Y be an ε ′ -correspondence between them. Assume (x 0 , y 0 ) ∈ C(X, Y ), we define functions d x 0 (.) = d X (x 0 , .) and d y 0 (.) = d Y (y 0 , .) in X and Y respectively. Then for any path γ x in X connecting x a , x b ∈ X, we can find a path γ y in Y such that its end points, y a , y b , are corresponding to x a , x b . And further more:

0 = l + (d r (x 2 ) -d r (x 1 )). As a consequence if d r (x 1 ) ≤ d r (x 2 ) then l/2 ≤ t 0 = d r (e(t 0 )) -d r (x 1 ); similarly if d r (x 1 ) ≥ d r (x 2 ) then l/2 ≤ l -t 0 = d r (e(t 0 )) -d r (x 2 ).
[x, v 0 ], E 1 , E 2 , •, E k , [v k , y] where [x, v 0 ] and [v k ,
E 1 = [v 0 , v 1 ], E 2 = [v 1 , v 2 ]•, E k = [v k-1 , v k ]
[min y∈γy d y 0 (y), max y∈γy d y 0 (y)] ⊂ [min x∈γx d x 0 (x) -2ε ′ , max x∈γx d x 0 (x) + 2ε ′ ]
Proof: Let ε > ε ′ > 0 and u, l be the maximum and minimum of d x 0 restricted to γ x . Since C is an ε ′correspondence for any x ∈ γ x there exists a point (x, y) ∈ C such that

d x 0 (x) -ε ′ ≤ d y 0 (y) ≤ d x 0 (x) + ε ′ .
As illustrated in Figure 3.4, the set of points y obtained in this way is not necessarily a continuous path from y a to y b . However one can consider a finite sequence

x 1 = x a , x 2 , • • • , x n = x b of points in γ x such that for any i = 1, • • • n -1 we have d X (x i , x i+1 ) < ε -ε ′ . If (x i , y i ) ∈ C then we have d Y (y i , y i+1 ) < ε -ε ′ + ε ′ = ε. As a consequence, since l -ε < l -ε ′ < d y 0 (y i ) < u + ε ′ < u + ε the shortest geodesic connecting y i to y i+1 in G remains in the set d -1 y 0 ([l -2ε, u + 2ε]
) and connecting these geodesics for all

X Y x 0 y 0 x a x b x 2
x 3

x 4

x 5

x 6

x 7

x 8

x 9

. i = 1, • • • , n -1 we get a continuous path from y a to y b in d -1 r ([l -2ε, u + 2ε]
). Now decreasing ε to ε ′ , we finish the construction.

As a corollary, we have the following theorem.

Theorem 3.4 Let (G, d G ) be a connected finite metric graph and let (X, d X ) be a compact geodesic metric space such that d GH (X, G) < ε for some ε > 0. Let x 0 ∈ X be a fixed point and let d x 0 = d X (x 0 , .) : X → [0, +∞) be the distance function to x 0 . Then for d ≥ α ≥ 0 the diameter of any connected component

L of d -1 x 0 ([d -α, d + α]) satisfies diam(L) ≤ 4(2 + N E (4(α + 2ε)))(α + 2ε) + ε
where N E (4(α + 2ε)) is the number of edges of G of length at most 4(α + 2ε). In particular if α = 0 and 8ε is smaller that the length of the shortest edge of G then diam(L) < 17ε.

Proof: Let ε ′ > 0 be such that d GH (X, G) < ε ′ < ε. Let C ⊂ X × G be an ε ′ -correspondence between X and G and (x 0 , r) ∈ C. we denote by d r = d G (r, .) : G → [0, +∞) the distance function to r in G. Let x a , x b ∈ L and let (x a , y a ), (x b , y b ) ∈ C.
There exists a continuous path γ ⊆ L joining x a to x b . Following lemma 3.6, we get a continuous path from

y a to y b in d -1 r ([d -α -2ε ′ , d + α + 2ε ′ ]). It then follows from Proposition 3.3 that d G (y a , y b ) ≤ 4(2 + N E (4(α + 2ε)))(α + 2ε) and since C is an ε ′ -correspondence (and so an ε-correspondence), d X (x a , x b ) < 4(2 + N E (4(α + 2ε)))(α + 2ε) + ε.
From Theorems 3.4 and 3.3 we obtain the following results for the Reeb and α-Reeb graphs.

Theorem 3.5 Let (X, d X ) be a compact connected path metric space, let r ∈ X be a fixed base point such that the metric Reeb graph (G, d G ) of the function d = d X (r, .) : X → R is a finite graph. If for a given ε > 0 there exists a finite metric graph

(G ′ , d G ′ ) such that d GH (X, G ′ ) < ε then we have d GH (X, G) < (β 1 (G) + 1)(17 + 8N E,G ′ (8ε))ε X 1 G 0 I i+1 I i I i 1 G X Figure 3
.5: X and X 1 are metric spaces that are both close to the metric graph G ′ with respect to the Gromov-Hausdorff distance. The α-Reeb graph, G, captures the stable topological features and ignores the unstable ones like the small holes in X 1 .

where N E,G ′ (8ε) is the number of edges of G ′ of length at most 8ε. In particular if X is at distance less than ε from a metric graph with shortest edge length larger than 8ε then d GH (X, G) < 17(β 1 (G) + 1)ε.

Theorem 3.6 Let (X, d X ) be a compact connected path metric space. Let r ∈ X, α > 0 and I be a finite covering of the segment [0, Diam(X)] by open intervals of length at most α such that the α-Reeb graph G α associated to I and the function

d = d X (r, .) : X → R is a finite graph. If for a given ε > 0 there exists a finite metric graph (G ′ , d G ′ ) such that d GH (X, G ′ ) < ε then we have d GH (X, G α ) < (β 1 (G α ) + 1)(4(2 + N E,G ′ (4(α + 2ε)))(α + 2ε) + ε)
where

N E,G ′ (4(α + 2ε)) is the number of edges of G ′ of length at most 4(α + 2ε).
In particular if X is at distance less than ε from a metric graph with shortest edge length larger than 4(α+2ε)

then d GH (X, G α ) < (β 1 (G α ) + 1)(8α + 17ε).

Recovery of Topology

In this section, we consider the the α-Reeb graph G of a distance function d : X → R, and show the following theorem which asserts that G recovers some topology of X.

Theorem 3.7 Let (X, d X ) be a compact connected path metric space and

(G ′ , d G ′ ) is a metric graph so that d GH (X, G ′ ) < ε. Let r ∈ X, α > 60ε and I = {[0, 2α), (iα, (i + 2)α)|1 ≤ i ≤ m} covers the segment [0, Diam(X)
] such that the α-Reeb graph G associated to I and the function d = d X (r, .) : X → R is a finite graph. If no edges of G ′ are shorter than L and no loops' lengths of G ′ are shorter than 2L with L ≥ 32α + 9ε, then we have G and G ′ are homotopy equivalent.

First we note that under the assumption of d GH (X, G ′ ) ≤ ε, X and G ′ are not necessarily homotopy equivalent. For example, in Figure 3.5 we assume that d GH (X, G ′ ) ≤ ε, then d GH (X 1 , G ′ ) ≤ ε as long as the perimeters of the two small holes in X 1 are both less than ε. Nevertheless, X and X 1 are not homotopy equivalent since X is contractible while X 1 is not.

On the other hand, we usually care less about smaller topological features such as the holes in X 1 , since they can be mixed with noisy features introduced in the real data acquisition. One of the important properties of the mapper algorithm is that it provides a way to extract robust (opposite to the small unstable) topological features in the data. Inheriting this property, the α-Reeb graph G in Figure 3.5 manages to capture the homotopy type of X and negates the two small loops of X 1 .

The above observations motivate our setting in this section, we assume that the sizes of the topological features, i.e., the length of the shortest edges and circles, of the underlying G ′ is bounded from below and that X is allowed to have unstable small features as long as d GH (X, G ′ ) ≤ ε. And then we prove that the α-Reeb graph is capable of capturing the topology of G ′ .

Our strategy of proving Theorem 3.7 is to construct some open covers for X and G ′ and relate the α-Reeb graph G and the graph G ′ to the nerves of the open covers. Specifically, we construct an initial open cover V 0 of X whose nerve N (V 0 ) is homotopy equivalent to G. Then we obtain a new open cover Ṽ of X by merging certain elements in V 0 while preserving the homotopy type of the nerve of the open cover, i.e., N (V 0 ) and N ( Ṽ) are homotopy equivalent. Based on the open cover Ṽ, we construct an open cover Ũ for G ′ whose nerve N ( Ũ ) is isomorphic to N ( Ṽ) as graphs and at the same time is homotopy equivalent to G ′ . In the following, we describe the constructions of the above open covers for X and G ′ and show the above claimed relations between them.

Since d GH (X, G ′ ) < ε, there exists an ε-correspondence between the two spaces, denoted C(X, G ′ ). For any subset V ⊂ X, denote C(V ) = {g ′ : (x, g ′ ) ∈ C(X, G ′ ), x ∈ V }, and similarly for any subset U ⊂ G ′ , denote C(U ) = {x : (x, g ′ ) ∈ C(X, G ′ ), g ′ ∈ U }. We call C(V ) and C(U ) are the correspondence set of V and U respectively under C(X, G ′ ). Recall that r ∈ X is the root point. Choose a point g r ∈ C(r) and define a distance function b : G ′ → R by b(g) = d G ′ (g r , g). Let N = {g n 1 , g n 2 , • • • , g np } be the vertices of G ′ , i.e.
, N is the set of vertices whose degree is not equal to two. From the hypotheses of the above theorem, the distance between any pair of vertices g n i , g n j with i ̸ = j is larger than L. For convenience, we also add into the vertices of G ′ the remaining local maximal/minimal points of the distance function b, which we denote using M = {g m 1 , • • • , g mq }. Note any newly added vertex g m i ∈ M is of degree two. We call the graph G ′ before adding the vertices in M the original G ′ , and the edges in the original G ′ the original edges of G ′ . An original edge of G ′ contains at most one vertex in M and thus can be split into at most two edges in G ′ .

Construction of Open Cover for X

We start with the following open cover of X. For each

I k ∈ I, denote V k = d -1 (I k ). V k may have several connected components, which can be listed in an arbitrary order. Denote V l k the l-th connected component of V k . Then V 0 = {V l k } k,l
is an open cover of X. Since at most two elements in I are overlapped, the nerve of V 0 , denoted N (V 0 ), is a graph. The following lemma states that any loop in the nerve N (V 0 ) is large, which is useful for the proof of Theorem 3.7. We say an open set

V l 1 k 1 ∈ V 0 is lower than the open set V l 2 k 2 ∈ V 0 if k 1 < k 2 and is higher than V l 2 k 2 if k 2 > k 1 .
Lemma 3.7 Let V l k and V i j are the lowest vertex and the highest vertex of a loop respectively in the nerve N (V 0 ). Then under the hypotheses of Theorem 3.7, we have jk ≥ 15.

Proof: First notice that j > k. Let x 1 ∈ V l k ∩ d -1 (kα, (k + 1)α) and x 2 ∈ V i j ∩ d -1 ((j + 1)α, (j + 2)α).
From the hypotheses of the lemma, there are two different paths

β 1 , β 2 connecting x 1 to x 2 so that β 1 ∩ d -1 ((k + 1)α, (j + 1)α) and β 2 ∩ d -1 ((k + 1)α, (j + 1)α) are in the different connected components of d -1 ((k + 1)α, (j + 1)α). Choose g i ∈ G ′ from C(x i ) for i = 1, 2. Following Lemma 3.6, the path β i in X for i = 1, 2 traces out a simple path γ i in G ′ connecting g 1 to g 2 so that γ i lies in b -1 (kα -ε, (j + 2)α + ε).
One can verify that γ 1 and γ 2 are two different paths due to the fact that β 1 and β 2 pass through different connected components of d -1 ((k + 1)α, (j + 1)α) and thus form a loop in G ′ , denoted γ. We have b(γ) ⊂ (kαε, (j + 2)α + ε).

We claim the range of the function b restricted to any loop, in particular β, covers an interval with the length at least L 2 . If the claim holds, then we have (jk + 2)α + 2ε ≥ L 2 , which implies jk ≥ 15 from the hypothesis L ≥ 32α + 9ε of Theorem 3.7. Indeed, if β contains at least two vertices in N , then it is obvious that the range of the function b restricted to γ covers an interval the length at least L 2 as any original edge of G ′ is longer than L. Now consider the case where γ contains one vertex in N , say g a . If γ does not contain g r , then there is exactly one local maximum on γ, say g b . If γ contains g r , let g b = g r . The removal of g a and g b cuts γ into two pieces. Along either piece, the function b has at most one local maximum. As the length of γ is longer than 2L. We have b(γ) covers an interval with length longer than L/2. Finally, if β contains no vertex in N , then G ′ is a single loop γ and the claim obviously holds.

In the following, we modify this open cover by merging while preserving the homotopy type of its nerve. The main purpose of the merging operation is make it easy to relate the open cover of X to the open cover of G ′ as constructed in Section 3.6.2. The merging operation is done in two steps.

For any vertex

g ∈ M ∪ N of G ′ , we construct a connected open set V (g) as the union of a subset of the open cover V 0 as follows. If b(g) ≥ α 2 , then there exists a unique positive integer k ′ s.t. k ′ α 2 ≤ b(g) < (k ′ + 1) α 2 . Let k = ⌊ k ′ +1 2 ⌋ -1 ≥ 0, and one can verify that (k + 1 2 )α ≤ b(g) ≤ (k + 3 2 )α. Therefore for all x ∈ C(g), d(x) ∈ [(k + 1 2 )α -ε, (k + 3 2 )α + ε] ⊂ I k . Moreover C(g) is contained in V l k ⊂ V k for some l. Indeed, if not, assume x 1 , x 2 ∈ C(g) with x i ∈ V i k for i ∈ {1, 2}
. By the definition of V i k , the geodesic connecting x 1 and x 2 must pass through a point

x 0 outside of V k , which means d X (x i , x 0 ) ≥ |d(x i ) -d(x 0 )| ≥ α 2 -ε. Then d X (x 1 , x 2 ) ≥ α -2ε which contradicts to the fact that d X (x 1 , x 2 ) ≤ d G ′ (g, g) + ε ≤ ε. Now we construct the open set V (g) as the union of the elements in the V V a V b l x 1 x 2 (j + 2)↵ (j + 1)↵ Figure 3.6: V with two lower neighborhoods.
open cover V 0 having non-empty intersection with V l k , i.e.,

V (g) = V ∈V 0 and V ∩V l k ̸ =∅ V.
In the case where b(g

) < α 2 , we construct the open set V (g) = V 0 ∪ V 1 = d -1 ([0, 3α)). Note in both cases, V (g) is a connected open set of X.
We abuse the notation and also denote V (g) the subset of V 0 whose union is the open set V (g). What V (g) represents will be clear from the context. For convenience, we call

V l k containing C(g) the center of V (g). Note that it is possible that V (g) = V (g ′ ) for two different vertices g, g ′ . Now we obtain an intermediate open cover of X V = {V (g) : g ∈ M ∪ N } {V ∈ V 0 : V / ∈ V (g), ∀g ∈ M ∪ N }
Note as a set, V does not have duplicated elements, i.e., if V (g) = V (g ′ ) for g ̸ = g ′ , then V only contains one copy of V (g). We call an open set V (g) ∈ V for any g ∈ M ∪ N critical and the remaining ones regular.

The following two lemmas describe the properties of the critical open sets and the regular open sets.

Lemma 3.8 Under the hypotheses of Theorem 3.7, we have for any vertex g ∈ M ∪ N , i) d(V (g)) ⊂ [sα, (s + 4)α] for some integer s ≥ 0, and ii) for any point x ∈ V ∈V 0 \V (g) V and any

g x ∈ C(x) ⊂ G ′ , d G ′ (g, g x ) ≥ α 2 -2ε. Proof:
The claim (i) is obvious from the construction of V (g). We now prove claim (ii). In the case where b

(g) < α 2 , for any x ∈ V ∈V 0 \V (g) V , we have d(x) > 3α and b(g x ) > 3α -ε. Thus d G ′ (g, g x ) ≥ |b(g x ) -b(g)| > 3α -ε -α 2 > α 2 -2ε. Now consider the case where b(g) ≥ α 2 . Assume V l k is the center of V (g). If d(x) ̸ ∈ I k , then d X (x, y) ≥ α 2 -ε for any point y ∈ C(g) from the construction of V (g), which implies d G ′ (g x , g) ≥ α 2 -2ε. Otherwise d(x) ∈ I k . Then x is not in V l k and the geodesic from x to any point y ∈ C(g) must pass x 0 / ∈ V k . This implies that d X (x, y) > d X (x 0 , y) ≥ α 2 -ε and d G ′ (g x , g) ≥ d X (x, y) -ε ≥ α 2 -2ε.
This proves the lemma.

Lemma 3.9 For any regular open set V ∈ V, V is also an open set in V 0 . Moreover, under the hypotheses of Theorem 3.7, it is of degree two in the nerve of N (V 0 ) with one neighboring vertex higher than V and one neighboring vertex lower than V .

Proof: We prove the lemma by contradiction. Assume V ∈ V 0 \ g∈M ∪N V (g) has two neighboring vertices, say V a , V b , which are lower than V . Without loss of generality, assume d(V ) ⊂ I j and d(V a ) and

d(V b ) are subsets of I j-1 . Let x a ∈ V a and x b ∈ V b such that (j -1)α < d(x a ), d(x b ) < jα.
As V a and V b both have non-empty intersection with V , there exist a path in d -1 ((j -1)α, (j + 2)α). Now let l = inf {s : there exists a path connecting

x a , x b in d -1 ((j -1)α, s] (V ∪ V a ∪ V b )}. We have (j + 1)α ≤ l < (j + 2)α as V a , V b are disconnected.
We can choose two points

x 1 , x 2 ∈ V ∪ V a ∪ V b from a path connecting x a and x b such that d(x 1 ) = d(x 2 ) = l -2ε, and x 1 , x 2 are disconnected in d -1 ([l -2ε, l)) (V ∪ V a ∪ V b ) but are connected by a path, say β, in d -1 ([l -2ε, l]) (V ∪ V a ∪ V b ). Obviously d X (x 1 , x 2 ) ≥ 2(l -(l -2ε)) = 4ε. Let g i ∈ C(x i ) for i = 1, 2. Then b(g i ) ∈ [l -3ε, l -ε] for i = 1, 2 and d G ′ (g 1 , g 2 ) ≥ d X (x 1 , x 2 ) -ε ≥ 3ε. Following Lemma 3.6, the path β traces out a simple path in G ′ denoted γ connecting g 1 and g 2 which lies in b -1 (l -4ε, l + 2ε). We claim γ must contain a vertex g c in M ∪ N . If not, d G ′ (g 1 , g 2 ) = |b(g 1 ) -b(g 2 )| ≤ |b(g 1 ) -d(x 1 ) + d(x 2 ) -b(g 2 )| ≤ |(b(g 1 ) -b(g r )) -(d(x 1 ) -d(r))| + |(d(x 2 ) -d(r)) -(b(g 2 ) -b(g r ))| = |d G ′ (g 1 , g r ) -d X (x 1 , r)| + |d G ′ (g 2 , g r ) -d X (x 2 , r)| ≤ 2ε
This contradicts to the fact that d G ′ (g 1 , g 2 ) ≥ 3ε. From the construction of γ, there exists a point g ∈ γ so that d G ′ (g, g c ) ≤ ε and there exists a point x ∈ β ∩ C(g). Then for any point

x c ∈ C(g c ), we have d X (x, x c ) ≤ 2ε. For any vertex x 0 ∈ V ∩ d -1 (l), since x 0 and x are connected in d -1 ([l -2ε, l]) (V ∪ V a ∪ V b ), d X (x 0 , x) ≤ 25ε from Theorem 3.4, and therefore d X (x 0 , x c ) ≤ 27ε. However, since x 0 ∈ V which is regular, x / ∈ V (g c ). From Lemma 3.8, we have d X (x 0 , x c ) ≥ α 2 -2ε > 27ε
. This is a contradiction. Therefore V can not have more than one neighboring vertices that are lower than V .

Using a similar argument we can also prove that V can not have more than one neighboring vertices that are higher than V .

We now perform a second step of merging. Two critical open set V (g 1 ) and V (g 2 ) in N (V) are said to be close if there is a simple path γ in the nerve N (V 0 ) connecting the center V l 1 k 1 of V (g 1 ) and the center V l 2 k 2 of V (g 2 ) so that γ consists of at most 4 edges. If there is a regular open set along the above path, we say this regular open set connects the critical open sets V (g 1 ) and V (g 2 ).

We have the following properties for two close critical open sets.

Lemma 3.10 Under the hypotheses of Theorem 3.7, we have (i) for any two vertices g n 1 , g n 2 ∈ N , V (g n 1 ) and V (g n 2 ) can not be close;

(ii) for any g m ∈ M , there exists at most one g n ∈ N such that V (g m ) and V (g n ) are close;

(iii) if V (g m 1 ) and V (g m 2 ) are close for any two vertices g m 1 , g m 2 ∈ M , then there must exist a vertex

g n ∈ N such that at least one of V (g m 1 ) and V (g m 2 ) is close to V (g n ).
Moreover, there is a path in N (V 0 ) of at most 5 edges connecting the center of V (g n ) to the center of V (g m i ) for any i = 1, 2.

Proof: Let V p j , V q k are the centers of V (g 1 ) and V (g 2 ) respectively for g 1 , g 2 ∈ N ∪ M . If V (g 1 ) and V (g 2 ) are close, then |k -j| ≤ 4. Assume j ≤ k. Then for any x 1 ∈ C(g 1 ) and x 2 ∈ C(g 2 ), there is a path in d -1 ((jα, (k + 2)α)) connecting x 1 and x 2 . Note that k + 2j ≤ 6. We claim that V (g 1 ) and V (g 2 ) are not close provided that d G ′ (g 1 , g 2 ) > 12α + 9ε. Indeed, since d G ′ (g 1 , g 2 ) > 12α + 9ε, the range of the function b restricted to any path connecting g 1 and g 2 in G ′ covers an interval of the length at least 6α + 4.5ε. This implies that the range of the function d restricted to any path in X connecting x 1 and x 2 covers an interval of the length at least 6α + 0.5ε. This means that V (g 1 ) and V (g 2 ) can not be close. Since

d G ′ (g n 1 , g n 2 ) ≥ L > 12α + 9α, V (g n 1 ) and V (g n 2 ) are not close. This proves (i). Assume V (g m ) is close to both V (g n 1 ) and V (g n 2 ) with g n 1 , g n 2 ∈ N . We have d G ′ (g n 1 , g n 2 ) ≤ d G ′ (g m , g n 1 ) + d G ′ (g m , g n 2 ) ≤ 24α + 18ε < L, which means g n 1 = g n 2 .
This proves (ii).

We now prove (iii). Since at most one vertex in M is added into an original edge of G ′ , any path in G ′ connecting g m 1 and g m 2 passes through at least one vertex from N . Furthermore, let γ be a geodesic in G ′ connecting g m 1 and g m 2 . If γ passes more than one vertices in N , d G ′ (g m 1 , g m 2 ) ≥ L > 12α + 9ε, which contradicts to the fact that V (g m 1 ) and V (g m 2 ) are close. Therefore γ contains exactly one vertex in N . Denote this vertex by g n . Let V l 1 k 1 and V l 2 k 2 be the centers of V (g m 1 ) and V (g m 2 ) respectively, and δ be the simple path from

V l 1 k 1 to V l 2 k 2 in N (V 0
) so that δ consists of at most 4 edges, or equivalently at most five elements in V 0 .

Recall V (g n ) consists of a subset of V 0 . We claim that δ must pass through an element in V (g n ).

If the claim holds, it is easy to verify that at least one of V (g m 1 ) and

V (g m 2 ) is close to V (g n ).
In addition, if we let V l k be the center of V (g n ), then there is a path in N (V 0 ) with as most 5 edges connecting V l i k i and V l k for any i = 1, 2. This proves (iii). It remains to show the above claim. We prove by contradiction. If we let

V (δ) = {V ∈ V 0 : V is on the path of δ},
then V (g n ) as a subset of V 0 does not intersect with V (δ). We have C(g m 1 ) and C(g m 2 ) are contained in V l 1 k 1 and V l 2 k 2 respectively. For any x 1 ∈ C(g m 1 ) and any x 2 ∈ C(g m 2 ), there is a path β in X connecting x 1 and x 2 so that β is contained in V ∈V (δ) V . From Lemma 3.8, for any x ∈ β and any

g x ∈ C(x), d G ′ (g x , g n ) ≥ α 2 -2ε.
From the construction in the proof of Theorem 3.4, the path β can trace out a simple path γ ′ in G ′ connecting g m 1 and g m 2 so that for any point g

∈ γ ′ , d G ′ (g, g n ) ≥ α 2 -3ε. This means that γ and γ ′ form a loop in G ′ . Since β is contained in V ∈V (δ) V , d(β)
is contained in an interval with the length of 6α, which implies that b(γ ′ ) is contained in an interval with the length of 6α + 4ε. Thus the length of γ ′ is at most 2(6α + 4ε) as it passes through at most one vertex in N . Since γ is a geodesic in G ′ connecting g m 1 and g m 2 , the length of the above loop is at most 4(6α + 4ε), which contradicts to the hypotheses of Theorem 3.7. This proves the above claim.

We say g a , g b ∈ M ∪ N are equivalent, denoted g a ∼ c g b , if there exists a finite sequence g a = g 1 , g 2 , ..., g k = g b such that V (g i ) and V (g i+1 ) are close for any i = 1, • • • , k -1. This is an equivalence relation. From Lemma 3.10 (iii), if an equivalence class contains at least two vertices in M ∪ N , it must contain a vertex in N . We have the following lemma Lemma 3.11 Under the hypotheses of Theorem 3.7, an equivalence class contains at most one vertex from N .

Proof: If not, assume g n 1 ̸ = g n 2 and g n 1 ∼ c g n 2 . Let g n 1 = g 1 , g 2 , ..., g k = g n 2 be a sequence so that V (g i ) and V (g i+1 ) are close for any i = 1, • • • , k -1. WLOG, we can further assume g i ∈ M for i = 2, • • • , k -1.
We first show that k > 5. Assume not. let andany x 2 ∈ C(g n 2 ), there is a path β in X connecting x 1 and x 2 so that d(β) is contained in an interval with the length at most 12α. The path β traces out a path γ in G ′ connecting g n 1 and g n 2 so that b(γ) is contained in an interval with the length at most 12α + 4ε, which implies

V l i k i be the center of V (g i ) for i = 1, • • • , k. From Lemma 3.10 (iii), there is a path in N (V 0 ) with at most 2 × 5 = 10 edges connecting V l 1 k 1 to V l 5 k 5 . Thus for any x 1 ∈ C(g n 1 )
d G ′ (g n 1 , g n 2 ) ≤ 2(12α + 4ε
). This contradicts to the hypothesis concerning the lengths of the edges in G ′ . Now we assume k > 5. Since V (g 3 ) and V (g 4 ) are close and g 3 , g 4 ∈ M , from Lemma 3.10 (iii), there exists a

g n ∈ N so that V (g n ) is close to at least one of V (g 3 ) and V (g 4 ). Assume V (g n ) is close to V (g 3 ). If V (g n ) ̸ = V (g n 1 ), we obtain a sequence of g ′ 1 = g n 1 , g ′ 2 = g 2 , g ′ 3 = g 3 , g ′ 4 = g n so that V (g ′ i ) and V (g ′ i+1 ) are close for any i = 1, • • • , 3. If V (g n ) ̸ = V (g n 2 ), we obtain a sequence g ′ 1 = g n , g ′ 2 = g 3 , • • • g ′ k-1 = g n 2 so that V (g ′ i ) and V (g ′ i+1 ) are close for any i = 1, • • • , k -2.
In either case, the new sequence has a length less than k. Similarly, we can obtain a shorter sequence if V (g n ) is close to V (g 4 ). One can keep shortening the sequence so that its length is no longer than 5, which however has been proven to be impossible. This proves the lemma. Now we are ready to further merge the open sets in V to obtain the final open cover Ṽ of X as follows. For any vertex and(3 • Ṽ (g 1 ) and Ṽ (g 2 ) are disjoint for two different g 1 , g 2 ∈ N .

g n ∈ N of G ′ , Let Ṽ (g n ) be the subset of V consisting of (1) V (g n ), and (2) any critical open set V (g) ⊂ V with g ∼ c g n ,
(g n ). What Ṽ (g n ) represents will be clear from the context. Let ṼN = {V ∈ V : V ∈ Ṽ (g n ) for some g n ∈ N }. The open cover Ṽ = Ṽ1 ∪ Ṽ2 ∪ Ṽ3 of X consists of three types of open sets: (1) Ṽ1 = { Ṽ (g n ) : g n ∈ N }; (2) Ṽ2 = {V (g) : g ∈ M and V (g) ̸ ⊂ ṼN } (3) Ṽ3 = {V ∈ V : V is regular and V ̸ ⊂ ṼN }.
• For any two open sets Ṽ1 , Ṽ2 ∈ Ṽ1 ∪ Ṽ2 , any path in the 1-skeleton of the nerve N ( Ṽ) connecting Ṽ1 , Ṽ2 consists of at least two elements from Ṽ3 .

• Any open set Ṽ ∈ Ṽ3 is also a regular element in V and thus an element in V 0 . Moreover any point

g ∈ C( Ṽ ) ⊂ G ′ is at least α 2 -2ε away from any vertex of G ′ .
Proposition 3.4 Under the hypotheses of Theorem 3.7, N ( Ṽ) and N (V 0 ) are homotopy equivalent.

N (V) N (V 0 ) N ( Ṽ) V (g m1 ) V (g n1 ) V (g n2 ) V (g m2 ) in Ṽ1 in Ṽ2 in Ṽ3

Proof:

We obtain the elements in the open covering Ṽ by merging a subset of open sets in V 0 . Think of any element Ṽ ∈ Ṽ as a subset of V 0 . The nerve N (V 0 ) restricted to Ṽ is a subgraph of N (V 0 ), whose vertex set is Ṽ and edge set includes the edges in N (V 0 ) with both endpoints in Ṽ . We call this subgraph the nerve of Ṽ , denoted N ( Ṽ ). The nerve of N ( Ṽ) as a topological space is the quotient space N (V 0 )/ Ṽ ∈ Ṽ N ( Ṽ ). From Proposition 0.17 in [Hatcher 2002], it is sufficient to show that N ( Ṽ ) is a tree for any Ṽ ∈ Ṽ.

For Ṽ ∈ Ṽ2 ∪ Ṽ3 , N ( Ṽ ) is obviously a tree. Consider Ṽ ∈ Ṽ1 . There exists a g n ∈ N so that V (g n ) ⊂ Ṽ (g n ) ∈ Ṽ . Let V t s be the center of V (g n ). For any g m ∈ M and V (g m ) ⊂ Ṽ (g n ), if let v t ′ s ′ be the center of V (g m ), from Lemma 3.10, |s -s ′ | < 5. Therefore, if V i
l is the element in Ṽ with the smallest sub-index and V j h is the element in Ṽ with the largest sub-index, then we have |h -l| ≤ 5 + 5 + 2 + 2 = 14. From Lemma 3.7, there is no loop in the subgraph N ( Ṽ ). This proves the proposition.

Construction of Open Cover for G ′

In this section, we construct an open cover

G ′ based on the open cover Ṽ of X. For an open set V ∈ V 0 , we construct a connected open set U V ⊂ G ′ so that C(V ) ⊂ U V as follows. Let l = min{d(V )} and u = max{d(V )}. We have u -l ≤ 2α. Let Ū = b -1 ([l -2ε, u + 2ε]), and then C(V ) ⊂ Ū . Since u -l + 4ε < 2α + 4ε < L
4 , one can verify that there is no loop in Ū and thus Ū consists of a set of trees. We claim C(V ) is contained in one of the trees. Indeed, for any two

g 1 , g 2 ∈ C(V ), we have l -ε < b(g 1 ), b(g 2 ) < u + ε. Now let x i ∈ V so that g i ∈ C(x i ) for i = 1, 2. Let β be a path in V connecting x 1 and x 2 . Following Lemma 3.6, β can trace out a path γ in Ū connecting g 1 and g 2 , which implies that C(V ) is contained in a tree in Ū . Let U V denote that tree. Let U 0 = {U V : V ∈ V 0 }. It is obvious that U 0 is an open cover of G ′ .
We now merge the elements in U 0 to construct a new open cover according to the way in which the elements in V 0 are merged to obtain Ṽ. Specifically, from our construction of Ṽ, any open set Ṽ ∈ Ṽ is the union of a subset of open sets in V 0 . We also denote this subset using Ṽ . Let

U Ṽ = {U V : V ∈ Ṽ ⊂ V 0 }. We also denote U Ṽ is the open set of the union of the open sets in U Ṽ .
Consider an open set Ṽ ∈ Ṽ3 . As it is also a regular open set in V and thus an open set in V 0 , d( Ṽ ) = (pα, (p+2)α) for some integer p > 0. From Corollary 3.1, any point in C( Ṽ ) is at least α 2 -2ε away from any vertex in M ∪N and any point in

U Ṽ is at least α 2 -4ε away from any vertex in M ∪N . Thus U Ṽ is a segment in G ′ without any branches. We shrink U Ṽ to obtain a new open set Ũ Ṽ = U Ṽ ∩ b -1 (pα +2ε, (p + 2)α -2ε), which is also a segment in G ′ . For any open set Ṽ ∈ Ṽ1 ∪ Ṽ2 , let Ũ Ṽ = U Ṽ . Thus we obtain Ũ = { Ũ Ṽ : Ṽ ∈ Ṽ}.
One can verify that Ũ is an open cover of G ′ . Moreover we have the following two lemmas which relate the nerve N ( Ṽ) to G ′ . Proposition 3.5 Under the hypotheses of Theorem 3.7, the nerve N ( Ṽ) and the nerve N ( Ũ ) are isomorphic as graphs.

Proof: It suffices to prove the following three claims.

• Claim (i): For any two Ṽi , Ṽj ∈ Ṽ1 ∪ Ṽ2 , Ũ Ṽi ∩ Ũ Ṽj = ∅.
Any path in N ( Ṽ) connecting Ṽi and Ṽj must pass through at least two open sets in Ṽ3 , which are regular open sets in V. From Lemma 3.9, any regular set has two neighbors in the nerve N (V) one lower and one higher, WLOG, assume Ṽi is higher than Ṽj . We have inf{d(x) :

x ∈ Ṽi } ≥ α + sup{d(x)|x ∈ Ṽj }, which implies inf{b(g)|g ∈ Ũ Ṽi } ≥ α+sup{b(g)|g ∈ Ũ Ṽj }-2ε > sup{b(g)|g ∈ Ũ Ṽj }. Thus Ũ Ṽi ∩ Ũ Ṽj = ∅. • Claim (ii): For any two Ṽi , Ṽj ∈ Ṽ3 , Ṽi ∩ Ṽj = ∅ if and only if Ũ Ṽi ∩ Ũ Ṽj = ∅. If Ṽi ∩ Ṽj ̸ = ∅, assume Ṽi is the only neighboring vertex in the nerve N (V) higher than Ṽj . Let d( Ṽj ) = (pα, (p + 2)α) and d( Ṽi ) = ((p + 1)α, (p + 3)α). Choose a point x from Ṽi ∩ Ṽj so that d(x) = (p + 3 2 )α. We have C(x) ∈ Ũ Ṽi ∩ Ũ Ṽj , which shows Ũ Ṽi ∩ Ũ Ṽj ̸ = ∅. If Ṽi ∩ Ṽj = ∅. Let d( Ṽi ) = (pα, (p + 2)α) and d( Ṽi ) = (qα, (q + 2)α). If |p -q| ≥ 2, it is obvious that Ũ Ṽi ∩ Ũ Ṽj = ∅. Now assume that q -p ≤ 1,
i ∈ C( Ṽi ) d G ′ (g i , g) ≥ α 2 -2ε and for any g j ∈ C( Ṽj ) d G ′ (g j , g) ≥ α 2 -2ε for any vertex g ∈ M ∪ N such that V (g) ∈ Ṽ . Thus d G ′ (g i , g j ) ≥ α -4ε, which implies Ũ Ṽi ∩ Ũ Ṽj = ∅. • Claim (iii): For any Ṽi ∈ Ṽ1 ∪ Ṽ2 and any Ṽj ∈ Ṽ3 , Ṽi ∩ Ṽj = ∅ if and only if Ũ Ṽi ∩ Ũ Ṽj = ∅.
First assume that Ṽi and Ṽj have a non-empty intersection. As Ṽj ∈ Ṽ3 , it is a regular open set in V which has one higher neighboring vertex and one lower neighboring vertex in N (V 0 ). Since Ṽj is regular, we have d( Ṽj ) = (pα, (p + 2)α) for some integer p > 0. We know Ṽi consists of a subset of open sets in V 0 and let V ∈ Ṽi be the open set in V 0 so that V ∩ Ṽj ̸ = ∅. WLOG, assume V is the higher neighboring vertex of Ṽj and we have Proof: As we have proved, Ũ is an open covering of G ′ . Since any edge on the original G ′ has a length longer than L, one can verify that any element of Ũ contains no loop and thus is a tree, and in particular is contractible. Furthermore, the union of any two elements of Ũ does not contains a loop. This means that if two elements of Ũ intersect with each other, their intersection is connected and thus contractible. Following from Nerve lemma, we have N ( Ũ ) is homotopy equivalent to G ′ . Proof of Theorem 3.7. From Proposition 3.4, Proposition 3.5, Proposition 3.6, it remains to show that the nerve N (V 0 ) is homotopy equivalent to the α-Reeb graph G. Indeed, we represent each node V l k in N (V 0 ) using a copy of the interval I k . If V l 1 k 1 and V l 2 k 2 with k 1 < k 2 are the endpoints of an edge in N (V 0 ), then we glue the upper half of I k 1 to the lower half of I k 2 . We identify any two points which are glued together directly or indirectly. By definition, the α-Reeb graph is the quotient space of the disjoint union of these intervals. From Lemma 3.7, there are more than one node between the top node and the bottom node of any loop in N (V 0 ). Thus, we have a one-to-one correspondence between the loops in N (V 0 ) and the loops in the α-Reeb graph. This proves the theorem.

d(V ∩ Ṽj ) ⊃ ((p + 1)α, (p + 2)α). We choose a point in x ∈ Ṽj ∩ V so that d(x) = (p + 2)α -4ε. Since b(C(x)) ⊂ ((p + 2)α -5ε, (p + 2)α -3ε), C(x) ∈ Ũ Ṽj ∩ Ũ Ṽi and thus Ũ Ṽi ∩ Ũ Ṽj ̸ = ∅.

Metric Reconstruction from Discrete Sampling

In the previous sections, we propose constructions of Reeb-type graphs for a given metric space (X, d X ) which is supposed to be close to an underlying metric graph (G ′ , d G ′ ) in the Gromov-Hausdorff distance. Our theoretical analysis shows that the Reeb-type graph approximation is both geometrically and topologically reliable under certain conditions. In practice, however, we usually only have access to a finite collection of isolated discrete points X sampled from X and need to build a metric d X on top of it. In general, it's unlikely to recover d X exactly with a finite sampling. Therefore it's necessary to understand the distortion between d X and its approximation d X , i.e., the Gromov-Hausdorff distance between metric spaces (X, d X ) and ( X, d X ).

Particularly, in this section we consider a case in which (X,

d X ) = (G, d G ) is a finite metric graph embedded in a d-dimensional Euclidean space R d . And we let X = {x 1 , x 2 , ..., x n } ⊂ G ⊂ R d be an ε-sampling as defined below. Definition 3.9 An ε-sampling of the embedded metric graph (G, d G ) is a finite set of points X = {x 1 , x 2 , x 3 , • • • , x n } ∈ G ⊂ R d , such that for any g ∈ G, there exists at least one x ∈ X satisfying d G (g, x) ≤ ε.
With a positive parameter σ, we construct a σ-neighborhood graph (see Definition 3.5), R σ = ( X, E σ ), on top of X with parameter σ > 0, and we denote the metric induced by R σ as in Definition 3.6 by d X,σ .

On the other hand, since X ⊂ G, restricting the graph metric d G to X give rise to another metric on X. The main goal of this section is to provide a uniform upper bound for d X,σ (x, y)-d G (x, y). Similar analysis on smooth manifolds of dimension higher than 1 is conducted in [START_REF] Bernstein | [END_REF]]. We in this section study the 1-dimensional case, which requires some different treatments.

Regularity Conditions.

The graph metric is intrinsic and independent of the ambient space. The sampling points, on the other hand, carry no information about the intrinsic metric. This fact arises difficulties in approximating graph metric with the sampling points. We illustrate a typical difficulty with the three singular cases of graphs embedded in R 2 in Figure 3.8. In each case, the red point and the blue point are distant in the graph but close in R 2 . Without prior knowledge about the graph metric, we probably end up creating a shortcut from the red points to the blue ones, therefore introduce significant distortions in metric approximation.

If the red point and the blue point are arbitrarily close to each other in the ambient space while keeping distant in the graph, then by no means can we avoid introducing the distortions. From this point of view, we propose some regularity conditions to exclude the singular situations.

To do that, we first introduce the reach of a manifold M embedded in R d . Given a point x ∈ R d and a set A ⊂ R d , the distance from x to A is defined as the infimum of the distances from x to points in A, i.e., d(x, A) = inf y∈A ∥x -y∥. The medical axis of M , denoted by med(M ), is defined as the following: Definition 3.10 The medical axis of M is the set of points which have more than one nearest neighborhoods in M , i.e.,

med(M ) = {x ∈ R d : there exist at least two points y 1 , y 2 ∈ M such that d(x, M ) = ∥x -y 1 ∥ = ∥x -y 2 ∥}
And finally we define the reach of M : Definition 3.11 The reach of M is the infimum of the distances from a point in M to the medical axis of M :

reach(M ) = inf p∈M d(p, med(M ))
It's routine to assume that reach(M ) is lower bounded by a positive constant in the study of approximating intrinsic geometric quantities with sampling points from the ambient space, however in the case of metric graph this condition is too restrictive. For example, the reach of the second and the third graphs in Figure 3.8 are zero due to the existences of sharp corners at the green points. Thus instead of bounding the global reach of G, we bound the reaches of some specified subgraphs of G from below.

Let V d be the vertex set of G (i.e. vertices of degree not equal to 2) and V l be the middle points of self-loops in G. We call V d the original vertex set of G and V a = V d ∪ V l the augmented vertex set of G. We break the self-loops into edges from the middle points of the loop, i.e., the points in V l . Formally, the augmented edge set E a is defined as

E a = i {(v i , v j ) : v i , v j ∈ V a
and there is no other vertex in the geodesic connecting v i to v j .}

We regularize a metric graph G embedded in R d with the following conditions on the augmented vertex set V a and the augmented edge set E a of G. These conditions depend on three positive parameters b, τ and β, whose roles are clear in the context. Condition 3.1 1. For any edge e ∈ E a , the reach of e is not smaller than τ , where τ is a positive constant.

For any two edges e

1 , e 2 ∈ E a , if e 1 ∩ e 2 = ∅, then reach(e 1 ∪ e 2 ) ≥ τ . 3. For any v ∈ V a , let V (v) = {v} ∪ {w ∈ V a : (v, w) ∈ E a } and E(v) = {(v, w) : w ∈ V (v), w ̸ = v}
V (v) = {w 1 , w 2 , w 3 }, E(v) = {e 1 , e 2 , e 3 , e 4 }
The first three conditions exclude the singular cases illustrated in Figure 3.8 with sufficiently large τ and β. First, thanks to condition 1, the reach of a single edge is lower bounded, meaning that if two points are distant in graph then they are distant in the ambient space. Secondly, condition 2 assures two nonintersecting edges is well-separated in Euclidean space, and therefore the case in the middle is prevented. Lastly, condition 3 takes care of the singularity between two incident edges such as the case in the right panel of Figure 3.8.

Upper-bound of |d

G -d H |.
Instead of directly comparing d X,σ and d G , we further construct two more metrics on X as auxiliaries. We start with constructing a graph H = ( X, E H ) on X with parameter δ > 0 as the following:

E H = {(x i , x j ) : d G (x i , x j ) ≤ δ, x i , x j ∈ X} (3.1)
Then two metrics d H , d S are defined on top of E H with respect to two measurements on the edges.

d H (x, y) = min Q k i=2 ∥x i-1 -x i ∥ d S (x, y) = min Q k i=2 d G (x i-1 , x i ) where Q = (x 1 , x 2 , • • • , x k ) varies over all paths along E H connecting x = x 1 to y = x k .
In the following, we present a uniform upper bound for |d G (x, y)d H (x, y)|. And then in Section 3.7.3, a uniform upper bound for |d H (x, y)d X,σ (x, y)| is given. Combining them we arrive at an upper bound for |d G (x, y)d X,σ (x, y)|.

Since

|d H (x, y) -d G (x, y)| ≤ |d S (x, y) -d G (x, y)| + |d H (x, y) -d S (x, y)|, we estimate |d H (x, y) - d G (x, y)| in two steps.
In step 1, we asserts that if the sampling density is high enough (or equivalently ε is small enough), then Let λ be the geodesic from x to y in G, we prove that λ is identical to a connected path in H. We sort points in λ ∩

d S = d G . x i+1 x i v 0 v x 0 x i+1 x i v (a) 2-vertices case (b) 1-vertex case.
X = (x, x 2 , x 3 , • • • , y) so that d G (x, x) ≤ d G (x, x 2 ) ≤ d G (x, x 3 ) ≤ • • • d G (x, y).
We claim that for any two consecutive points x i , x i+1 , d G (x i , x i+1 ) ≤ δ, i.e., x i and x i+1 are connected in H.

Denote by γ the geodesic between x i , x i+1 . If there are two vertices, say v, v ′ ∈ V a , be in the interior of γ (see Figure 3.10(a)). Then

d G (x i , x i+1 ) ≥ d G (v, v ′ ) ≥ b > δ ≥ 4ε.
As a consequence, there must be another sample x ′ between v, v ′ : otherwise let g be the middle point of the geodesic between v and v ′ , and the nearest sampling point to g is of distance larger than b 2 > 2ε, violating the assumption that X is an ε-sampling. On the other hand, according to the construction of (x, x 2 , • • • , y), there is not another sampling point between the consecutive points x i-1 , x i , and therefore there is at most one vertex in V a in the interior of γ.

If there exists exactly one vertex of V d (of degree larger than 2) between x i and x i+1 . As mentioned above, there is not a sampling point in the interior γ. Let x ′ be the sampling point that is closest to v and not in the interior of γ (as depicted in Figure 3.10(b)). With the same argument before,

d G (x i , x ′ ) ≤ 2ε and d G (x i+1 , x ′ ) ≤ 2ε, meaning that d G (x i , x i+1 ) ≤ 4ε ≤ δ.
If there is no vertex in V d between x i and x i+1 . Let g be the middle point of geodesic from

x i to x i+1 . Since X is an ε-sampling, d G (x i , x i+1 ) ≤ d G (g, x i ) + d G (g, x i+1 ) ≤ 2ε < δ. Noticing that (x, x 2 , x 3 , • • • , y) is a path in graph H connecting x to y, we have d S (x, y) ≤ d G (x 1 , x 2 )+ d G (x 2 , x 3 ) + ... + d G (x k-1 , x k ) = d G (x 1 , x k ) = d G (x, y)
, which finishes the proof.

Then we proceed to bound |d H (x, y)d S (x, y)|. A lemma by Herbert Federer [START_REF] Federer | Curvature measures[END_REF]] as the following is repeatedly used in our analysis.

Lemma 3.12 [Federer] Let K be a manifold embedded in R d , and

reach(K) = τ > 0. Let K δ be a δ-offset of K, i.e., K δ = x∈K {y ∈ R d : ∥x -y∥ ≤ δ}. For any δ < τ , the map Π δ : K δ → K is τ τ -δ -lipschitz. (Π δ maps each point in K δ to its projection in K.)
Applying Federer's lemma to a sub-graph whose reach is bounded from below by τ , we have the following corollary.

Corollary 3.2 Let K be a geodesically convex sub-graph of G, reach(K) = τ > 0. Then for any two points x, y on K, if ∥x -y∥ = l < 2τ , then d G (x, y) ≤ τ τ -l 2 ∥x -y∥.
Proof: We consider a parameterization of the straight line segment from x to y in R d such that s( 0

) = x, s(l) = y. The whole segment s([0, l]) is contained in K l 2 since that for any point p in this segment, ∥p -x∥ + ∥p -y∥ = l. Now let a partition of [0, l] be 0 = t 0 < t 1 < ... < t n = l, with Federer's lemma we have ∥Π l 2 (s(t i+1 )) -Π l 2 (s(t i ))∥ ≤ τ τ -l 2 ∥t i+1 -t i ∥.
Adding up both side of the inequality, we obtain

d G (x, y) ≤ τ τ -l 2 ∥x -y∥.
0 ≤ d S (x, y) -d H (x, y) ≤ δ 2τ diam(G) + (2 + v∈V d (deg(v) -1))δ
Proof: Let P = (x 1 , ..., x k ) be a path in H connecting two points x = x 1 and y = x k such that d H (x, y) is realized by P . We can further assume that P is a simple path in H. It follows from the definition of d H that the path realizing the minimal length between x and y must be simple, otherwise one can remove loops without increasing the path lengths.

Let γ i be the geodesic connecting x i-1 and x i in G. Following the construction of H, the length of γ i is not more than δ. We first claim that there must not be any point x j with j ≥ i + 1 in γ i (including the end points x i-1 , x i ): otherwise, let x j be on γ i , then d G (x i-1 , x j ) ≤ d G (x i-1 , x i ) ≤ δ, meaning that we can directly connect x i-1 , x j to reduce the length of P . (see Figure 3.11(a) for an illustration). That contradicts the optimality of P .

A corollary of this observation is that, a point v ∈ V d (the vertex set of G with degree larger than 2) is contained in at most deg(v) -1 pieces of geodesic γ i . Particularly we analyze the case of deg(v) = 3 as shown in Figure 3.11(b). Let γ i be the first geodesic containing v, and then we consider the position of x i+1 : it can't be in γ i as proven before, moreover it can't be in the edge containing x i-1 and v, otherwise it will be in the geodesic of some γ k , k ≤ i -1. thus it either jumps to the third branch or stays in the same edge with

x i+1 x i 2 x i v x i 1 x i 1 x i x j x y v x i x i 1 x i+1 x i+1 x i 2 (a) (b) (c)
x i and d G (x i+1 , v) > d G (x i , v).
• If it's the former case, then v is contained by 2 geodesics: γ i and γ i+1 . We claim that v can't be contained by γ j , j ≥ i + 2. That is because that neither of x j-1 , x j can be in the tripod formed by x i-1 , x i , x i+1 and v. thus γ j must contains one of the three geodesics: γ i , γ i+1 and the one connecting x i-1 , x i+1 . In either case we can show that we can connect one of x i-1 , x i , x i+1 with one of x j , x j-1 , and again that contradicts the optimality of P .

• If it's the latter case, then similarly we can prove that v can't be crossed by γ j , j ≥ i + 1.

For the case deg(v) > 3, the analysis is the same (see Figure 3.11(c) for a situation that deg(v) = 4 and v is crossed by 3 geodesics).

We classify all the segments

(x i-1 , x i ), i = 2, 3, • • • k into three types. Note that since (x i-1 , x i ) are connected in H, d G (x i-1 , x i ) ≤ δ < b.
Thus it follows the same arguments in the proof of Proposition 3.9 that γ i , ∀i can not pass more than 1 vertices in the augmented vertex set V a .

• Type A: There is no vertex in γ i .

• Type B: There is exactly one vertex of V d , the original vertex set of G, in γ i .

• Type C: There is exactly one vertex of V l , the set of middle points of self-loops in G, in γ i .

For

(x i-1 , x i ) in type A, we have d G (x i-1 , x i ) ≤ δ < τ . Following lemma 3.2, we have d G (x i-1 , x i ) ≤ τ τ -δ 2 ∥x i-1 -x i ∥, meaning that d H (x i-1 , x i ) ≤ d S (x i-1 , x i ) ≤ τ τ -δ 2 d H (x i-1 , x i ). For (x i-1 , x i ) in type B, we have δ ≥ d G (x i-1 , x i ) ≥ ∥x i-1 -x i ∥ ≥ 0, thus 0 ≤ d G (x i-1 , x i )-∥x i-1 - x i ∥ ≤ δ.
Since the simple path can't pass v more than deg(v) -1 times, the number of segments of type B is then upper-bounded by

v∈V d (deg(v) -1).
In the last case, since P is a simple path in H connecting x and y, P at most pass two vertices in V l . Notice that each vertex in V l corresponds to a self-loop of G, if there are 3 vertices in V l on the geodesic between x and y, then there must be at least one self-loop contains neither x nor y, meaning that P is not simple. To conclude, there are at most two segments falling into type C, for each of them, say (x j , x j+1 ), 0 < d G (x j , x j+1 ) -∥x jx j+1 ∥ ≤ δ Lastly, since for any path d H (P ) ≤ d S (P ), d H (x, y) ≤ d S (x, y), ∀x, y ∈ X. On the other hand, let P = P A P B P C be a decomposition of path P regarding the above classification, then d

S (P A ) -d H (P A ) ≤ d S (P A ) -τ -δ/2 τ d S (P A ) ≤ δ 2τ d S (P ) ≤ δ 2τ d G (x, y) ≤ δ 2τ diam(G). Similarly d S (P B ) -d H (P B ) ≤ v∈V d (deg(v) -1)δ, and d S (P C ) -d H (P C ) ≤ 2δ
Putting them together, we have

0 ≤ d S (x, y) -d H (x, y) ≤ δ 2τ diam(G) + (2 + v∈V d (deg(v) -1))δ
Combing Proposition 3.7 and 3.8, we prove the following:

Proposition 3.9 Let G be a metric graph embedded in R d satisfying condition 3.1, and X be an ε-sampling of G. Let H be a graph constructed with parameter δ as in Definition 3.

1. diam(G) is the diameter of G and V d is the original vertex set of G. τ, b are parameters defined in condition 3.1. If min{τ, b} > δ ≥ 4ε, then |d G (x, y) -d H (x, y)| ≤ δ 2τ diam(G) + (2 + v∈V d (deg(v) -1))δ, ∀x, y ∈ X . 3.7.3 Upper-bound of |d H -d X,σ |
Fixing δ = 4ε, it follows from Proposition 3.9 that if ε tends to 0, then d H (x, y) converges to d G (x, y) for any x, y ∈ X. However in the practice, we can't construct graph H directly since the metric d G is unknown.

In the following, we prove that if G satisfies condition 3.1, then we can approximate d G with metric of the Rips graph constructed on X with a proper parameter.

Before proceeding, we claim the following lemma which are used repeatedly in the proof of the upcoming theorem.

Lemma 3.13 Let x, y be two points lying on the same edge e ∈ E whose reach is at least

τ > δ. If ∥x -y∥ ≤ δ/2, then d G (x, y) ≤ δ.
Proof: This lemma is essentially a rephrase of corollary 3.2, e is obviously a geodesically convex sub-graph of G,

thus if ∥x -y∥ ≤ δ/2, d G (x, y) ≤ τ τ -δ/4 ∥x -y∥ ≤ 2δ 3 < δ.
Theorem 3.8 G is a metric graph satisfying conditions 3.1, X is an ε-sampling of G, δ is a positive parameter such that min{b, τ } > δ, where b, τ and β below are constants defined in conditions 3.1. Then for any two points x, y in X, if ∥x -y∥ Proof: For the sake of clarity, in the proof we let c = β 2β+2 < 1/2. We assume that x, y are on the edges e 1 , e 2 ∈ E a , where E a is the augmented edge set of G, respectively.

≤ β 2β+2 δ, then d G (x, y) ≤ δ.
The problem is analyzed on the following several cases:

1. e 1 ∩ e 2 = ∅.

If e 1 ∩ e 2 = ∅, we claim ∥x -y∥ > cδ: due to condition 3.1, reach(e 1 ∪ e 2 ) ≥ τ > δ. We let m = inf

x 1 ∈e 1 ,x 2 ∈e 2 ∥x 1 -x 2 ∥ and assume that x ′ 1 ∈ e 1 , x ′ 2 ∈ e 2 satisfy ∥x ′ 1 -x ′ 2 ∥ = m (
the existences are assured by that e 1 ∪ e 2 is compact). The interior of the Euclidean ball centered at the middle point between x ′ 1 and x ′ 2 with radius m 2 is not intersecting with e 1 ∪ e 2 , otherwise we can find another pair of points whose Euclidean distance is smaller than m, contradicting to the assumption on (x ′ 1 , x ′ 2 ). Therefore the center of the ball, r c , must be on the medial axis, meaning that ∥x -y∥

≥ ∥x ′ 1 -x ′ 2 ∥ = ∥x ′ 1 -r c ∥ + ∥x ′ 2 -r c ∥ ≥ τ + τ = 2τ > 2δ > cδ.
2. e 1 ∩ e 2 ̸ = ∅.

We first deal with the case in which e 1 and e 2 are identical, the latter two cases are illustrated in Figure 3.12.

To avoid heavy notations, we denote e Equivalently x, y are on the same edge with positive reach larger than τ and ∥x -y∥ ≤ cδ.

Directly by corollary 3.2 we know that d G (x, y)

≤ τ τ -∥x-y∥/2 ∥x -y∥ ≤ τ τ -cδ/2 cδ ≤ c 1-c/2 δ ≤ 2 3 δ < δ. (b) e 1 ∩ e 2 = {v a , v b }. First, if x, y ∈ U \(B a ∪ B b ), then it follows from condition 3.1 that ∥x -y∥ > 2βr = cδ. Secondly if x, y ∈ B a ∪ B b : either x, y are in the same geodesic ball, say B a , then d G (x, y) ≤ 2 cδ 2β = cδ β ≤ δ 2 < δ; or x ∈ B a and y ∈ B b , then given d G (v a , v b ) ≥ b > δ, ∥v a -v b ∥ > δ/2. Thus ∥x -y∥ ≥ ∥v a -v b ∥ -∥x -v a ∥ -∥y -v b ∥ > δ/2 -cδ β = cδ. Lastly, without loss of generality, we let x ∈ U \(B a ∪ B b ) and y ∈ B a . If d G (x, v a ) > δ, by lemma 3.13, ∥x-v a ∥ > δ/2. Obviously ∥v a -y∥ ≤ r = cδ 2β , so ∥x-y∥ ≥ ∥x-v a ∥-∥v a -y∥ > δ/2 -cδ 2β > cδ. Now we assume d G (x, v a ) ≤ δ and d G (x, v a ) + d G (v a , y) > δ. Following corollary 3.2 one has ∥x -v a ∥ ≥ τ -d G (x,va)/2 τ d G (x, v a ) ≥ τ -δ/2 τ (δ -d G (v a , y)) ≥ τ -δ/2 τ (δ -cδ 2β ). On the other hand ∥y -v a ∥ ≤ cδ 2β . Then it amounts to proving τ -δ/2 τ (δ - cδ 2β ) - cδ 2β > cδ.
Reorganizing, we have the following equivalent form:

τ -δ/2 τ > ( τ -δ/2 2βτ + 1 2β + 1)c.
The left-hand-side of the above inequity is not smaller than (δδ/2)/δ = 1/2 (as τ ≥ δ), while the right-hand-side is smaller than

(1 + 1/β)c = 1/2, thus the inequality is verified, meaning that ∥x -y∥ ≥ ∥x -v a ∥ -∥y -v a ∥ > cδ given d G (x, y) > δ. (c) e 1 ∩ e 2 = v a . First if x, y ∈ U \B a , then by the first regular condition, ∥x -y∥ > 2βr = cδ. Secondly if x, y ∈ B a , d G (x, y) ≤ 2r = cδ β < δ/2 < δ.
The last scenario is x ∈ U \B a and y ∈ B a , which is obviously identical to the last scenario in the case of e 1 ∩ e 2 = {v a , v b }.

Equipped with the above propositions, we conclude the main theorem of this section as the following: Theorem 3.9 If metric graph G satisfies conditions 3.1, X is an ε-sampling of G, if there exists a parameter σ such that 4ε ≤ σ < β 2β+2 min{τ, b}. Then the Rips graph with parameter σ induces a metric on X, d X,σ , as in Definition 3.5 satisfying:

0 ≤ d G (x, y) -d X,σ (x, y) ≤ (2 + 2/β)σ 2τ diam(G) + ( v∈V d (deg(v) -1) + 2)(2 + 2/β)σ
.

Proof: Let δ = (2 + 2/β)σ, and build a graph H with δ as in Definition 3.1. Then δ satisfies conditions in Proposition 3.9 thanks to δ = (2 + 2/β)σ < (2 + 2/β) β 2β+2 min{τ, b} = min{τ, b} and δ > σ ≥ 4ε.

Therefore we have |d G (x, y) -d H (x, y)| ≤ δ 2τ diam(G) + ( v∈V d (deg(v) -1) + 2)δ.
We argue that R σ is a subgraph of H. In fact, if x i , x j are connected in R σ , then ∥x ix j ∥ ≤ σ, according to Theorem 3.8, d G (x i , x j ) ≤ (2 + 2/β)σ = δ, therefore x i , x j are connected in H.

The last step is to prove that d X,σ (x, y) ≥ d H (x, y). Indeed if x, y are connected in G, then following the proof of Proposition 3.7 we have x, y are connected with a path

P = (x 1 , x 2 , • • • , x k ) in H where d G (x i , x i+1 ) ≤ 4ε, ∀1 ≤ i ≤ k -1. On the other hand, since ∥x i -x i+1 ∥ ≤ d G (x i , x i+1 ) ≤ 4ε ≤ σ, P is also a path in R σ , therefore
x, y are connected in both H and R σ . Now that R σ ⊂ H, any path connecting them in R σ must connect them in H. Following the definitions of d X,σ and d H , we have d X,σ (x, y) ≥ d H (x, y). Putting them together, we finish the proof.

Gromov-Hausdorff Distance between (G, d G ) and ( X, d X,σ ). As a corollary, we derive a bound for the Gromov-Hausdorff distance between (G, d G ) and ( X, d X,σ ). To start with, we construct a correspondence between points the two metric spaces:

C ⊂ G × X. For x ∈ X ⊂ G, we let (x, x) ∈ C. For each g ∈ G, we let x(g) = {x ∈ Xs.t.d G (x, g) = min y∈ X d G (y, g)} and add (g, x(g)) to C. For any g 1 , g 2 , since X is an ε-sampling of G, we have d G (g i , x(g i )) ≤ ε, i = 1, 2. It follows from Theorem 3.9 that |d G (g 1 , g 2 ) -d X,σ (x(g 1 ), x(g 2 ))| ≤ |d G (g 1 , g 2 ) -d G (x(g 1 ), x(g 2 ))| + |d G (x(g 1 ), x(g 2 )) -d X,σ (x(g 1 ), x(g 2 ))| ≤ |d G (g 1 , g 2 ) -d G (g 1 , x(g 2 ))| + |d G (g 1 , x(g 2 )) -d G (g(x 1 ), x(g 2 ))| + (2 + 2/β)σ 2τ diam(G) + v∈V d (deg(v) -1)(2 + 2/β)σ ≤ d G (g 2 , x(g 2 )) + d G (g 1 , x(g 1 )) + (2 + 2/β)σ 2τ diam(G) + ( v∈V d (deg(v) -1) + 2)(2 + 2/β)σ ≤ 2ε + (2 + 2/β)σ 2τ diam(G) + ( v∈V d (deg(v) -1) + 2)(2 + 2/β)σ
After all, we obtain that if the assumptions of Theorem 3.9 are satisfied, then

d GH ( X, G) ≤ ε + (1 + 1/β)σ 2τ diam(G) + ( v∈V d (deg(v) -1) + 2)(1 + 1/β)σ
This bound of the Gromov-Hausdorff distance between (G, d G ) and ( X, d X,σ ) suggests two facts that impact the quality of approximation: one is the sampling density and the other is the regularity of the underlying graph G. Higher sampling density (or equivalently smaller ε) allows for a smaller σ: as the bound suggests, if ε tends to zero, we are guaranteed to recover the underling metric d G . On the other hand, the larger β and τ are, the smaller the bound is. Stated differently, we need denser sampling for approximating the metric of a graph with less regularity.

I 4 I 2 I 1 I 3 V 1 4 V 1 0 V 1 4 V 1 3 V 2 2 V 1 2 V 1 0 I 0 V 2 1 V 1 1 V 2 1 V 1 1 I 4 I 3 I 2 I 2 I 0 I 1 I 1 X d I G disjoint union of copies of intervals ↵-Reeb graph V 1 3 V 1 2 V 2 2 Figure 3
.14: Illustration of the different steps of the algorithm for computing α-Reeb graph. In the disjoint union of copies of intervals, the subintervals marked with same labels are identified in the α-Reeb graph.

Algorithm

In this section, we describe an algorithm for computing the α-Reeb graph for some α > 0. We assume the input of the algorithm includes a neighboring graph X = (V, E), a function l : E → R + specifying the edge length and a parameter α. In the applications where the input is given as a set of points together with pairwise distances, i.e., a finite metric space, one can generate the neighboring graph X as a Rips graph of the input points with the parameter chosen as a fraction of α. We assume X is connected as one can apply the algorithm to each connected component otherwise. Our algorithm, whose different steps are illustrated in Figure 3.14, can be described as follows. In the first step, we fix a node of X as the root r and then obtain the distance function d : V → R + by computing d(v) as the graph distance from the node v to r. In the second step, we apply the Mapper algorithm [Singh 2007] to the nodes V with filter d to construct a graph G. Specifically, let I = {(iα, (i + 1)α), ((i + 0.5)α, (i + 1.5)α)|0 ≤ i ≤ m} so that I k ∈I I k covers the range of the function d. We say an interval

I k 1 ∈ I is lower than another interval I k 2 ∈ I if the midpoint of I k 1 is smaller than that of I k 2 . Now let V k = d -1 (I k ) and V l k be the lth component of V k Then of {V l k } k,l
is a cover of H and the graph G constructed by the Mapper algorithm is the 1-skeleton of the nerve of that cover. Namely, each node in G represents an element in

{V l k } k,l . Two nodes V l 1 k 1 and V l 2 k 2 are connected with an edge if V l 1 k 1 ∩ V l 2 k 2 ̸ = ∅.
In fact, when we check if

V l 1 k 1 ∩ V l 2 k 2 ̸ = ∅,
we only need to check if their vertices are overlapped or not as we assume the lengths of the edges in H are fractions of α.

In the final step, we represent each node V l k in G using a copy of the interval I k . As mentioned in the Section 3.4, α-Reeb graph is a quotient space of the disjoint union of those copies of intervals. Specifically, for an edge in G, let V l 1 k 1 and V l 2 k 2 be its endpoints. Then I k 1 and I k 2 must be partially overlapped. We identify the overlap part of these two intervals. After identifying the overlapped intervals for all edges in G, the resulting quotient space is the α-Reeb graph. Algorithmically, the identification is performed as follows.

We split each copy of internal I k into two by adding a point in the middle. Now think of it as a graph with two edges and label one of them upper and the other lower. Notice that two overlapped intervals I k 1 and I k 2 can not be exactly the same. One must be lower than the other. To identify their overlapped part, we identify the upper edge of the lower interval with the lower edge of the upper interval.

Experimental Results

In this section, we illustrate the performances of our algorithm in different applications.

Earthquake Data

The first data set was obtained from USGS Earthquake Search [EarthquakeSearch ]. It consists of earthquake epicenter locations collected, between 01/01/1970 and 01/01/2010 in the rectangular area between latitudes -75 degrees and 75 degrees and longitude -170 degrees and 10 degrees with magnitude greater than 5.0.

This raw earthquake data set contains the coordinates of the epicenters of 12790 earthquakes that are mainly located around geological faults. We follow the procedure described in [Aanjaneya 2012] to remove outliers and randomly sampled 1600 landmarks. Finally, we computed a neighboring graph from these landmarks with parameter 4. The length of an edge in this graph is the Euclidean distance between its endpoints. For each connected component, we fix a root point and compute the graph distance function d to the root point as shown in Figure 3.15(a). We also set α = 4 and apply our algorithm to the above data to obtain the α-Reeb graph. In general, the α-Reeb graph is an abstract metric graph. In this example, for the purpose of visualization, we use the coordinates of the landmarks to embed the graph into the plane as follows. Recall that for a copy of interval I k representing the node V l k in G, we split it into two by adding a point in the middle. We embed the endpoints of the interval to the landmarks of the minimum and the maximum of the function d in V l k , and the point in the middle to the landmark of the median of the function d in V l k . Figure 3.15(b) shows the embedding of the α-Reeb graph. Note this embedding may introduce metric distortion, i.e., the Euclidean length of the edge may not reflect the length of the corresponding edge in the α-Reeb graph.

GPS Data

The second data set is that of 500 GPS traces tagged "Moscow" from OpenStreetMap [Openstreetmap ]. Since cars move on roads, we expect the locations of cars to provide information about the metric graph structure of the Moscow road network. We first selected a metric ε-net on the raw GPS locations with ε = 0.0001 using furthest point sampling. Then, we computed a neighboring graph from the samples with To evaluate the quality of our α-Reeb graph for each data set, we computed both original pairwise distances, and pairwise distances approximated from the constructed α-Reeb graph. For GPS traces, we randomly select 100 points as the data set is too big to compute all pairwise distances. We also evaluated the use of α-

GPS Data with Crossing

The third data set we consider is also obtained from GPS traces. Roads are often split so that cars in different directions run in different lanes. In particular, this is the true for highways. In addition, when two roads cross in GPS coordinates, they may bypass through a tunnel or an evaluated bridge and thus the road network itself may not cross. Such directional information is contained in the GPS traces. We encode this directional information by stacking several consecutive GPS coordinates to form a point in a higher dimensional space.

In this way, we obtain a new set of points in this higher dimension space. Then we build a neighboring graph for this new set of points based on L 2 norm and apply our algorithm to recover the road network. In particular, although the paths intersect at the cross in GPS coordinates, the road network does not and this should be detected by our algorithm. To test the above strategy, we extract those GPS traces from the above "Moscow" dataset which pass through a highway crossing as shown in Figure 3.17(a). Since GPS records the position based on time, we resample the traces so that the distances between any two consecutive samples is the same among all traces. Then we apply the above algorithm to the resampled traces. Figure 3.17

Introduction

Shape comparison is a fundamental problem in geometry processing. In the most general setting, this problem consists of encoding and quantifying similarities and differences across pairs or collections of shapes. This can be especially useful for shape retrieval [START_REF] Tangelder | [END_REF][START_REF] Bai | [END_REF]], interpolation [Xu 2006[START_REF][END_REF], or visualization [Praßni 2010]. However, even when a map between shapes is given, encoding and visualizing the differences between them are still challenging. Approaches based on the point-to-point correspondences usually suffer issues such as sensitivity to noise, difficulty of selecting an appropriate scale of analysis and inconvenient visualization. The discrete nature of point correspondences is one of the major reasons of these issues. The framework of functional maps, which is introduced in [START_REF] Ovsjanikov | [END_REF]], alleviates the issues to some extent by converting point-to-point correspondences into linear operators across function spaces on each of the shapes. As demonstrated in [START_REF] Ovsjanikov | [END_REF]], the functional map is a compact, informative representation and suitable to incorporate with tools from spectral analysis.

N M T f

1 V (f 1 ) f 2 V (f 2 ) (a) (b) k = 30 k = 50 k = 100 T N Ñ M
In this chapter, we present theoretical analyses for two frameworks based on the notion of functional map. Given a pair of shapes M, N and a map T : M → N , the framework of [Rustamov 2013] encodes differences between shapes into a pair of linear operators (so-called shape difference operators) acting on function spaces on N . The other framework, proposed in [START_REF] Ovsjanikov | [END_REF]], is for map analysis and visualization. In the same setting as above, this framework generates a collection of multi-scale indicators, which are functions on N highlighting areas that deformed by T . As we will show later, the latter framework can be unified into the first one. Thus our analyses are essentially all performed on the shape difference operators.

The theoretical formulation of the shape difference operators is well-established, yet associated stability analyses remain absent. In practice, we observe robustness of these frameworks. For example, in Figure 4.1(b), we feed the framework of [START_REF] Ovsjanikov | [END_REF]] two pairs of meshed shapes: a meshed bumpy sphere M compared to two spheres N with distinct mesh structures, and generate indicators at three different scales. Two types of consistency are evidenced: horizontally, the scale k increases from 30 to 100, yet the indicators of each rows highlight nearly the same areas; vertically at each scale, though the input meshes of N are distinct, the resulting indicators are comparable.

It is then appealing to study the stability properties of the shape difference operators. As hinted by the above numerical results, two types of stability properties are worth considering: one is with respect to perturbations of input shapes in both [Rustamov 2013] and [START_REF] Ovsjanikov | [END_REF], and the other one is regarding the changes in scale, which is peculiar to the framework of [START_REF] Ovsjanikov | [END_REF]].

Though all the above frameworks assume that M and N are 3D shapes, i.e., 2-dimensional Riemannian manifolds embedded in R 3 , the formulations are actually well-defined when M and N are Riemannian manifolds of an arbitrary dimension. This fact allows for potentials of these frameworks in dealing with objects beyond shapes. From this point of view, in our analyses, we assume that M and N are smooth Riemannian manifold of dimension n.

How to exploit such potentials in practical applications, however, is another interesting yet challenging problem. For example, the above frameworks have so far only constructed on triangle meshes that are only available when dealing with 3D shapes. For manifolds that are of intrinsic dimension more than 2 or embedded in R d with d ≥ 4, the implementations are not obvious. In Chapter 5, we initiate our exploration to this problem by comparing shapes in a more primitive setting, where only discrete sampling points from the shapes are given.

Overview

We assume that M and N are two connect compact, smooth, n-dimensional Riemannian manifolds without boundary, endowed with metrics g M and g N . And let a map T : M → N be a smooth map between them.

In [Rustamov 2013], the authors introduce a pair of self-adjoint operators acting on real-value functions on N , each of which captures one type of differences between the two shapes regarding T . Particularly, for one of the shape difference operators -V as illustrated in Figure 4.1, a functional proposed in [START_REF] Ovsjanikov | [END_REF]] evaluates the deviation from a function f on N to its image V f (we will verify this connection at the beginning of Section 4.6). The maximizer of this functional is supposed to highlight the most deformed regions under map T . Instead of searching for the global maximizer, the authors maximize this functional within the space spanned by the first k eigenfunctions of the Laplace-Beltrami operator on N , and the maximizer is viewed as an indicator at scale k. By changing the scale k, a collection of indicators and a sequence of the corresponding maxima are then obtained, allowing the users to select output(s) at one or more scales for better understandings(see Figure 4.2 for an illustration). As mentioned before, we provide stability analyses regarding two type of perturbations.

• Perturbations with respect to the input manifolds: We start by introducing a model in Section 4.4 characterizing perturbations on the input manifolds. Then in Section 4.5, we discuss the stability of the two different shape difference operators under our perturbation model.

• Perturbations with respect to the changes in scale: As mentioned above, the multi-scale framework depends on a collection of subdomains (in [START_REF] Ovsjanikov | [END_REF]], that is, the spaces spanned by the first k eigenfunctions, with k = 1, 2, • • • ) in which we maximize a specified functional. However in Section 4.6.1, we demonstrate that the original subdomain construction indexed by k is not suitable for stability analysis and construct a new one which is closely related to the original but controlled by a continuous C. Then further in Section 4.6, we verify this stability with respect to this continuous scale C.

Especially, we extend the framework of [START_REF] Ovsjanikov | [END_REF]] in Section 4.6.5 by adapting the other shape difference operator to it, which enables visualizing another type of shape differences other than the areabased one. And we as well prove that this extension enjoys a similar stability property as the original one.

Finally, in Section 4.7, we demonstrate numerical experiments that reflect the stability in practice.

Related Work

The two frameworks we analyze in this chapter are based on the notion of functional maps, which has been a key ingredient of various applications in geometry processing, to name a few, analyzing maps between shapes [START_REF] Huang | [END_REF]], vector field processing [Azencot 2013[START_REF] Azencot | [END_REF]] and image segmentation [Wang 2013].

Perturbation analysis has a long and rich history. In this chapter, we perform perturbation analysis on both shape difference operators (which are linear operators, see [Kato 1995] for an introduction of perturbation analysis on them) and a spectral method based on such operators. The spectral methods have long been applied in various areas: spectral clustering [von Luxburg 2006], shape analysis [Reuter 2006] and so on. Besides demonstrating practical usefulness of the spectral methods, providing theoretical justifications is attracting more and more research interests. Theoretical guarantees for spectral clustering algorithms often stem from Cheeger's inequality, which is powerful if there exists a significant spectral gap. Assuming such a gap, several works [START_REF] Kwok | [END_REF][START_REF] James R Lee | Multi-way spectral partitioning and higher-order cheeger inequalities[END_REF], Louis 2012, Oveis Gharan 2014[START_REF] Dey | [END_REF]] present theoretical guarantees on the quality (measured by some graph conductance) of the output of the respective algorithms. It is worth noting that the above works only concern a single object, while in this paper, we study pairwise objects. From this point of view, our work has a similar flavor to the ones by Mémoli [Mémoli 2009, Mémoli 2011], who proposes metrics among shapes based on spectral invariant and discusses their robustness with respect to perturbations on the input shapes.

Beyond spectral methods, in geometric and topological data analysis, several approaches have been proposed for guaranteeing stability of the data processing and analysis techniques. In particular, stability has been theoretically proven in many works on estimating geometric quantities. In [START_REF][END_REF]], the authors analyze both effects of practical and theoretical facts on the accuracy of normal estimation process. In [Mérigot 2011], it is assured that the sharp feature detection algorithm there is stable with respect to Hausdorff noise. In the same noise model, the stability of the curvature measures is proven under certain conditions in [Chazal 2009]. Similar problems are also actively studied in the community of topological data analysis (TDA). The stability of persistence diagram is verified in [Cohen-Steiner 2007], which lays down a solid theoretical foundation for further research in TDA. Some more recent developments in TDA come with stability assurance as well, for example the notion of distance to a measure [START_REF] Chazal | [END_REF].

A rich body of research has been devoted to provide such analysis for convergence properties of various discrete Laplacian operator. In [Wardetzky 2005, Xu 2007, Dey 2010] the converging behaviors of the cotangent Laplacian operators on meshes to the underlying Laplace-Beltrami operators are investigated from diverse perspectives. While in [START_REF] Belkin | Constructing Laplace operator from point clouds in Rd[END_REF][START_REF] Liu | [END_REF], Hein 2007, Dey 2013b], similar problems are considered in a different setting, where the discrete Laplacian operators are built on point clouds.

Preliminaries

We first give a brief review of concepts about differential geometry and the (weighted) Laplace-Beltrami operators on manifolds. We refer the readers to [Grigoryan 2000] for a more detailed introduction (especially to the latter topic).

Then we briefly take a review over the functional-maps-based frameworks, which are the main focuses Riemannian Metrics A Riemannian manifold is a manifold M endowed with a Riemannian metric, g M , which is a family of {g M (x)} x∈M such that at any point x, g M (x) is a symmetrical, positive definite and bilinear form on the tangent space T x M , i.e., g M (x) :

T x M × T x M → R + .
A metric g M at x gives rise to an inner product on T x M : for η, ξ ∈ T x M , ⟨η, ξ⟩ g M (x) = g M (x)(η, ξ).

The length of η is then defined as ⟨η, η⟩ g M (x) . Particularly, in the local coordinate system g M admits a matrix representation, whose (i, j)-th element, g i,j , equals

⟨ ∂ ∂x i , ∂ ∂x j ⟩ g M (x) . Matrix [g i,j
] is invertible, and we denote by [g i,j ] its inverse matrix.

The matrix of a single metric g M is depending on the choice of the local coordinate system. Interestingly, we obtain some invariance when taking another metric into consideration. Let g M and gM be two metrics on M . Provided a local coordinate system (x 1 , x 2 , • • • , x n ), we denote by [g] and [g] the matrices of g M and gM in this system. Now assumed that another local coordinate system is obtained by letting (y

1 , y 2 , • • • , y n ) = F (x 1 , x 2 , • • • , x n ),
where F : R n → R n is a smooth map. The matrices with respect to y is denoted as [g] y and [g] y . It follows from lemma 3.12 in [Grigoryan 2000] that [g] = J(F ) T [g] y J(F ) and [g] = J(F ) T [g] y J(F ), in which J(F ) is the Jacobian matrix of F . Hence we have the following equation:

det[g] y det[g] y = det(J(F ) T [g]J(F )) det(J(F ) T [g]J(F )) = det[g] det[g] , (4.1)
which suggests the ratio of determinants of the metrics is invariant to the changes of the local coordinate systems.

Gradient Given f ∈ C ∞ (M ), its gradient at x is defined as

∇ g M f = n i,j=1 g i,j ∂f ∂x j ∂ ∂x i . Notice that ∇ g M f
is spanned by { ∂ ∂x i }, thus it is also a tangent vector in T x M . For a lighter notation, we write ∇ g M f as ∇f and by ⟨∇f, ∇g⟩ g M we always mean ⟨∇ g M f, ∇ g M g⟩ g M , i.e., the inner product of the function gradients are defined upon the metric.

Riemannian Measure A Riemannian metric g M induces a volume (also known as a Riemannian measure), ν M , on the family of all measurable subsets of M . The following theorem assures the existence and uniqueness of ν M . We refer the reader to [Grigoryan 2000] for a detailed proof (see Theorem 3.11 there).

Theorem 4.1 For any Riemannian manifold (M, g M ), there exists a unique measure ν M on Λ(M ) (the family of all measurable sets in M ) such that, in any chart U ,

dν M = det[g i,j ]dλ (4.2)
where [g i,j ] is the matrix of g M and λ is the Lebesgue measure in U .

The following proposition characterizes how a perturbation on the metric on a manifold affects the associated volume and gradient norm.

Proposition 4.1 Let g M and gM be two metrics on a n-dimensional manifold M , and ν M , νM be the volumes induced by them respectively. If there exists a constant a ≥ 1 such that for any point x ∈ M and any η ∈ T x M , the following inequality holds:

a -1 ≤ ⟨η, η⟩ g M ⟨η, η⟩ gM ≤ a (4.3) , then a -n/2 dν M ≤ dν M ≤ a n/2 dν M (4.4)
and for any f ∈ C ∞ (M ),

a -1 ≤ ⟨∇f, ∇f ⟩ g M ⟨∇f, ∇f ⟩ gM ≤ a (4.5)
Proof: To prove inequality 4.4, we show that the a -n ≤ det[g] det[g] ≤ a n , where [g] and [g] are matrices of g M and gM in a common local coordinate system.

The invariance of Equation 4.1 allows us to pick freely a local coordinate system without affecting the ratio of the determinants . We then consider {e 1 , e 2 , • • • , e n }, an orthonormal basis of T x M , in which the matrix representation of gM is an identity matrix. Thus we have det[g] = 1, and in the following we bound det[g]. Let [g] be the matrix of g M (x) in basis {e 1 , e 2 , • • • , e n }. If there exists an eigenvalue of [g], λ, that is larger than a, then we let u be the associated eigenvector and η = n i=1 u i e i . As a consequence, ⟨η, η⟩ g M = u T [g]u = λ > a. Since [g] is an identity matrix, we have ⟨η, η⟩ gM = 1 and the ratio ⟨η,η⟩g M ⟨η,η⟩ gM exceeds a, which contradicts the assumption 4.3. On the other hand, since g M (x) is positive definite, thus all its eigenvalues are positive and not larger than a. Thus the determinant of [g], which equals the product of all the eigenvalues, is not larger than a n . As a consequence of Equation 4.2, we have dν M ≤ a n/2 dν M . The other side of inequality 4.4 can be proven in the same way.

To prove inequality 4.5, we notice that the definition of ∇f depends on the matrix of metric, thus we need to do some transformation. We consider df = [g i,j ]∇f , which is independent of g. Then we have ⟨∇f, ∇f

⟩ g M = ([g i,j ]∇f, ∇f ) = ([g i,j ][g i,j ]∇f, [g i,j ]∇f ) = ([g i,j ]df, df ) = ⟨df, df ⟩ g -1 M and correspond- ingly ⟨∇f, ∇f ⟩ gM = ⟨df, df ⟩ g-1 M .
We use the same basis (e 1 , e 2 , • • • , e n ) of T x M as above. For gM , in this basis [g i,j ] and its inverse [g i,j ] are identity matrices. As we've proven that all eigenvalues of [g i,j ] do not exceed a, we have the smallest eigenvalue of [g i,j ] is at least a -1 . Therefore we derive that ⟨∇f, ∇f

⟩ g M = ⟨df, df ⟩ g -1 M ≤ a⟨df, df ⟩ g-1 M = a⟨∇f, ∇f ⟩ gM
The other side of inequality 4.5 is proven symmetrically.

Weighted Riemannian Manifold For a Riemannian manifold (M, g M ), besides the volume (or the Riemannian measure) induced by g M , we can define other measures. For example, a measure µ M on M defined by dµ M = ρdν N . where ρ is a smooth positive-value function on M . Such a measure is a weighted version of the volume ν M , and it gives rise to a weighted Riemannian manifold. Definition 4.3 A weighted Riemannian manifold is a triple (M, g M , µ M ), where g M is Riemannian metric inducing a volume ν M and there exists a positive smooth function ρ such that dµ M = ρdν M .

Particularly, if ρ is a constant function: ρ(x) = 1, ∀x ∈ M , then µ M is the volume induced by g M .

Function spaces

In this part we introduce some function spaces on (weighted) Riemannian manifolds. We denote the set of all real-value functions on M by F(M ). The set of smooth functions on M , C ∞ (M ), is obviously a subset of F(M ).

Let µ M be a measure on (M, g M ) (not necessarily a volume), by L 2 µ (M ) we denote the function space in which all functions have finite square norm with respect to µ M , i.e.,

L 2 µ (M ) = {f ∈ F(M ) : M f 2 dµ M < ∞} A special subspace of L 2 µ (M ) with smoothness constraints is H 1 0,µ (M ), which is defined on top of H 1 µ (M ): H 1 µ (M ) = {f : M f 2 + ⟨∇f, ∇f ⟩ g M dµ M < ∞} H 1 0,µ (M )
is the closure of all infinitely differentiable functions in H 1 µ that is compactly supported by M . An alternative understanding of H 1 0,µ (M ) is the following:

H 1 0,µ (M ) = {f : f ∈ H 1 µ (M ) and f ∂M = 0}.
Inner Products After introducing L 2 µ (M ) and H 1 0,µ (M ), we define two inner products on them respectively, which are crucial in formulating the shape difference operators (Section 4.3.5).

The first inner product takes in two functions of L 2 µ (M ) and returns a real-number:

⟨f, g⟩ L 2 µ = M fgdµ M , f, g ∈ L 2 µ (M ) (4.6)
And we define the second inner product on function space H 1 0,µ (M ):

⟨f, g⟩ H 1 0,µ = M ⟨∇f, ∇g⟩ g M dµ M , f, g ∈ H 1 0,µ (M ) (4.7)

The Laplace-Beltrami Operator

We denote the Laplace-Beltrami operator on manifold (M, g M ) by ∆ M .

Definition 4.4 Given a smooth Riemannian manifold (M, g M ) of dimension n, ∆ M : C ∞ (M ) → C ∞ (M ) is an operator on C ∞ (M ), which in a local coordinate system (x 1 , x 2 , • • • , x n ) is defined as: ∆ M = n i,j=1 1 √ detg ∂ ∂x i detgg i,j ∂ ∂x j
We call (ϕ, λ) an eigensolution if it satisfies the following eigen equation. Especially ϕ ∈ C ∞ (M ) is called an eigenfunction matching an eigenvalue λ ∈ R.

∆ M ϕ + λϕ = 0 (4.8)
It is well-known that the Laplace-Beltrami operator is self-adjoint and negative semidefinite. Thus in any eigensolution, λ is guaranteed to be a non-negative real number. As we assume M to be a compact manifold, the set of eigenvalues is discrete and infinite. Thus we can sort the eigenvalues with an ascending order so that 0 ≤ λ 1 ≤ λ 2 ≤ λ 3 ≤ • • • and denote the associated eigenfunctions by ϕ 1 , ϕ 2 , ϕ 3 , • • • . The Green's first identity states:

M ⟨∇u, ∇v⟩ g M dν M + M u∆ M vdν M = ∂M u∇v • nds M (4.9)
where u, v ∈ C ∞ (M ), ∂(M ) is the boundary of M , and n is the unit normal vector field to ∂M defined at each point on ∂M . Since we assume M to be boundaryless, the right hand side of Equation 4.9 is zero, resulting in the following: For any smooth functions u, v on manifold M which is compact and without boundary.

M u(-∆ M )vdν M = M ⟨∇u, ∇v⟩ g M dν M (4.10)
Considering two eigenfunctions, ϕ i , ϕ j associated to distinct eigenvalues λ i , λ j , it follows from Equation 4.10 that

M ϕ i (λ j ϕ j )dν M = M ϕ i (-∆ M )ϕ j dν M = M ⟨∇ϕ i , ∇ϕ j ⟩ g M dν M = M ϕ j (-∆ M )ϕ i dν M = M ϕ j (λ i ϕ i )dν M .
Therefore we have M ϕ i ϕ j dν M = 0 for λ i ̸ = λ j .

For each eigenvalue λ i , the functions satisfying ∆ M ϕ + λ i ϕ = 0 form a linear function space. If λ i is of multiplicity 1, then this space is of dimension 1, which is spanned by {ϕ i }. If λ i is of multiplicity more than 1, the corresponding function space might be of dimension k > 1. In this case, with Gram Schmidt method we can always find {ϕ l i } l=1,••• ,k forms a basis of this k-dimension space and

M ϕ k 1 i ϕ k 2 i dν M = 0, 1 ≤ k 1 ̸ = k 2 ≤ k.
In fact, the eigenfunctions the Laplace-Beltrami operator form a basis of L 2 ν (M ). According to above arguments, we always assume that the basis is orthogonal. Moreover, by simply rescaling ϕ i so that

M ϕ 2 i dν M = 1, we obtain an orthonormal basis of L 2 ν (M ). Proposition 4.2 Let {ϕ i } i≥1 be an orthonormal basis of L 2 ν (M ) consisting of eigenfunctions to ∆ M ϕ + λϕ = 0. For any function u ∈ L 2 ν (M ), it admits a eigendecomposition u = i≥1 a i ϕ i , a i = M uϕ i dν M . And we have M u 2 dν M = i≥1 a 2 i (4.11) If we further assume that u is differentiable, then M ⟨∇u, ∇u⟩ g M dν M = i≥1 a 2 i λ i (4.12) Proof: Assuming that u = i≥1 a i ϕ i , we have M uϕ j dν M = M i≥1
a i ϕ j ϕ j dν M . It follows the orthonormality of the basis that M uϕ j dν M = a j . 4.12 follows from Equation 4.10: M ⟨∇u, ∇u⟩

Similarly, M u 2 dν M = M ( i a i ϕ i ) 2 dν M = M a 2 i i≥1 ϕ 2 i dν M = i≥1 a 2 i as M ϕ i ϕ j dν M = 0 and M ϕ 2 i dν M = 1. The Equation
g M dν M = M u(-∆ M )udν M = M ( i≥1 a i ϕ i )( i≥1 a i λ i ϕ i )dν M = i≥1 a 2 i λ i
Weighted Laplace-Beltrami Operator We have defined a weighted Riemannian manifold (M, g M , µ M ) (see Definition 4.3), where µ M = ρν M . A weighted Laplace-Beltrami operator is defined as:

∆ M,µ = n i,j=1 1 ρ 1 √ detg ∂ ∂x i ρ detgg i,j ∂ ∂x j (4.13)
Notice that if ρ = 1, then ∆ M,µ boils down to ∆ M . Equation 4.10 is still valid with respect to the weighted Laplace-Beltrami operator, but with a change of measures.

Functional Maps

The framework of functional maps is introduced in [START_REF] Ovsjanikov | [END_REF]], which converts point correspondences between two manifolds to function correspondences between function spaces on the manifolds. Let M and N be two smooth manifolds and T be a map from M to N . A functional map, T F , is a pull-back induced by T . Namely, given a real-value function w ∈ F(N ), we define T F (w) = w • T ∈ F(M ). Therefore T F is a map from F(N ) to F(M ). Below we review two important properties of T F : informativeness and linearity.

First we claim that if T is bijective, then it can be fully recovered by T F . In fact, for any point a ∈ N , let δ a (x) be an indicator function on N such that δ a (x) equals 1 if x = a and δ a (x) = 0 otherwise. By construction, g = T F (δ a ) satisfies that g(y) equals 1 if T (y) = a and 0 otherwise. Since that T is bijective, there exist a unique b such that T (b) = a ,thus g is also an indicator function but on M . After all, for each a ∈ N , we search for a point b ∈ M such that T F (δ a )(b) equals 1. Such a point b is unique and satisfies T (b) = a, meaning that we can completely recover T with T F .

Secondly, T F is a linear map across the two function spaces, i.e., for any α, β ∈ R and f, g ∈ F(N ), we have T F (αf + βg) = αT F (f ) + βT F (g). This linearity property indicates that T F admits a (potentially infinite) matrix representation. Suppose that {ϕ N i }, {ψ M j } form a basis of F(N ) and of F(M ) respectively. Let w = i=1 a i ϕ N i and accordingly T F (w) = j=1 b j ψ M j , then there exists a unique matrix C T such that

C T a = b holds for any w, where a = (a 1 , a 2 , • • • , a n , • • • ) and b = (b 1 , b 2 , • • • , b n , • • • ).
We then call C T the matrix representation of T F with respect to basis {ϕ N i } and {ψ M j }. In practice, we usually consider {ϕ N i }, {ψ M j } as eigenbasis of the Laplace-Beltrami operator on N and M respectively. And then truncation on each basis is taken, say, we only take the first k N and k M eigenfunctions on each shape. Then C T is a k M by k N matrix instead of an infinite-dimensional one. More importantly, we can control the rank of this matrix representation to reach a trade-off between accuracy and complexity. The authors of [START_REF] Ovsjanikov | [END_REF] show that hundreds of eigenfunctions are enough to reconstruct a reasonable functional map between meshed shapes that consist of tens of thousands of vertices.

Map Analysis

Based on the functional maps, the authors of [START_REF] Ovsjanikov | [END_REF]] present a multi-scale framework to detect and visualize the area-distortions induced by maps T : M → N . Unlike the approaches which perform comparisons via point correspondences, this framework analyzes the function correspondences, which are easily converted from the point-to-point maps with the above functional map framework.

The key ingredient of this framework is a functional that measures area-distortions induced by transforming a function w ∈ L 2 (N ) to L 2 (M ) with the functional map, T F :

E(w) = M T F (w) 2 dν M N w 2 dν N (4.15)
As discussed in [START_REF] Ovsjanikov | [END_REF], if E(w) is large, then w is supposed to take high absolute value within the areas where measure are less preserved when mapped from M to N via T .

Therefore, it is natural to search for functions such that E(w) is large. Instead of maximizing E(w) in L 2 ν (N ), a multi-scale approach is taken by forcing w to reside in a subspace of L 2 ν (N ) spanned by the first k eigenfunctions of ∆ N , the Laplace-Beltrami operator on N . We define S(k) as the following

S(k) = span{ϕ i } k i=1
, where ϕ i is the i-th eigenfunction of ∆ N .

(4.16)

The constrained optimization problem considered in [START_REF] Ovsjanikov | [END_REF]] is then:

max E(w) s.t. w ∈ S(k) (4.17)
The advantage of choosing such a collection of subdomains {S(k)} k∈N + is multi-fold: (1) the span space is easy to construct and the associated constrained optimization is straightforward to solve; (2) the feasible solutions are spanned by low-frequency eigenfunctions, which are fit for the visualization purpose and stable with respect to noises of input discrete shapes. (3) the flexibility of choosing k provides a multiscale understanding of the problem. As demonstrated in Figure 4.2, the regions highlighted by w * become more and more localized as k, the dimension of spanned space, increases.

Shape Difference Operators

Another map-based approach for shape comparison is proposed in [Rustamov 2013], where a pair of Shape Difference Operators is introduced.

The area-based shape difference operator, V : L 2 ν (N ) → L 2 ν (N ), is defined as a linear operator such that for any f, g in square-integrable space

L 2 ν (N ) = {f : N f 2 dν N < +∞}, N fV (g)dν N = M T F (f )T F (g)dν M (4.18)
Note that unless T is an area-preserving map, equation N fgdν N = M T F (f )T F (g)dν M does not hold for any pair of f, g ∈ L 2 (N ), V captures the difference and compensates the discrepancy.

With respect to a different inner product, ⟨, ⟩ H 1 0,ν , the conformal-based shape difference operator, R :

H 1 0,ν (N ) → H 1 0,ν (N )
, is a linear operator such that for any f, g in Sobolev space

H 1 0,ν (N ), N ⟨∇f, ∇R(g)⟩ g N dν N = M ⟨∇T F (f ), ∇T F (g)⟩ g M dν M (4.19)
It follows from the Riesz representation theorem that V, R exist and are unique. Moreover, thanks to the commutative properties of inner products 4.6 and 4.7, V and R are as well self-adjoint operators, i.e., N fV (g)dν N = N gV (f )dν N and N ⟨∇f, ∇R(g)⟩ g N dν N = N ⟨∇g, ∇R(f )⟩ g N dν N hold for any f, g residing in the corresponding domains.

Particularly, the authors of [Rustamov 2013] prove that T is locally area-preserving if and only if V is an identity operator. That proof (of theorem 1 therein) applies to the case where M and N are Riemannian manifolds of dimension n. Regarding R, the author of [START_REF] Schumacher | Conformal maps and p-Dirichlet energies[END_REF] shows that if the dimensions of M, N are n = 2, then T is a conformal map if and only if R is an identity operator, and if n > 2 then T is an isometric map if and only if R is an identity operator. Therefore, regardless of the dimension of the input Riemannian manifolds, the shape difference operators are well-defined and both carry information about the deviation from T to an isometric map.

Perturbation Model and Bounded-distortion Condition

As mentioned in Section 4.1, one of the main results in this chapter is the stability analyses of the shape difference operators with respect to perturbations on the input manifolds. In this section, we propose our perturbation model and a bounded-distortion condition throughout our analyses.

Perturbation Model

Let M be a connected, compact, smooth, n-dimensional manifold without boundary endowed with a Riemannian metric g M . And let ν M be the volume induced by g M . We first define (a, b)-closeness between Riemannian structures on the same smooth manifold, where a, b are positive constants not smaller than 1.

We now introduce our model for characterizing perturbations on the input shapes.

Definition 4.5 A Riemannian manifold (N, gN , νN ) is a-close to another one (N, g N , ν N ) if the following holds: For any x ∈ N and any tangent vector η in T x N , the tangent plane at x: a -1 ≤ ⟨η,η⟩g ⟨η,η⟩ g ≤ a holds for a constant a ≥ 1. It is clear that the a, b-closeness characterizes perturbations on the metric and on the measure, respectively. Combining them together, a weighted Riemannian manifold, (N, gN , μN ), is said to be (a, b)-close to a Riemannian manifold (N, g N , ν N ) if

• (N, gN , μN ) is b-close to the corresponding Riemannian manifold (N, gN , νN ). • (N, gN , νN ) is a-close to (N, g N , ν N ).
Intuitively, we view (N, gN , μN ) as a perturbed version of (N, g N , ν N ). It is obvious that (1, 1)closeness implies that the two are identical.

The smooth function ρ M is usually modeled as a sampling density on manifold M . The second restriction on two Riemannian metrics might seem a bit abstract. We describe it in a special case, where M is a 2dimensional manifold and M 1 = (M, g M ), M 2 = (M, gM ) are two Riemannian manifolds embedded in the same Euclidean space R 3 . Assume that there exists Φ :

M 1 → M 2 that is a diffeomorphism. Let x ∈ M 1 ⊂ R 3 and y = Φ(x) ∈ M 2 , which are both in R 3 . If g M ,
gM satisfy the above condition, then the tangent plane of M 1 at x is parallel to the one of M 2 at y. On the other hand, the differentiation of Φ gives rise to a push-forward between the tangent planes:

a tangent vector v ∈ T x M 1 is mapped to D(Φ)v ∈ T y M 2 .
Due to the fact that M 1 , M 2 are both embedded in R 3 , the inequality is re-written as:

a -1 ≤ ⟨v, v⟩ g R 3 ⟨D(Φ)v, D(Φ)v⟩ g R 3 ≤ a,
where g R 3 is the Euclidean metric. In other words, denote the matrix representation of D(Φ) by A, then the eigenvalues of A T A both lie in the interval [a -1 , a].

The following proposition characterizes the quantitative relationships between (M, g M , ν M ) and (M, gM , μM ).

Proposition 4.3 If (M, gM , μM ) is (a, b)-close to (M, g M , ν M ), then for any smooth function w on M . a -1 ≤ ⟨∇w, ∇w⟩ g M ⟨∇w, ∇w⟩ gM ≤ a, and
a -n/2 b -1 dμ M ≤ dν M ≤ a n/2 bdμ M .
Proof: The first inequality has been proven in Proposition 4.1.

Regarding the second one, again thanks to Proposition 4.1,

a -n/2 dν M ≤ dν M ≤ a n/2 dν M . By definition of (a, b) -closeness, we have dν M ≥ a -n/2 dμ M = a -n/2 ρ -1 M dν M ≥ a -n/2 b -1 dμ M .
And similarly we obtain that dν M ≤ a n/2 bdμ M , which finishes the proof.

We always consider (M, g M , ν M ) as the original input Riemannian manifold and (M, gM , μM ) as a perturbed one. The magnitudes of the perturbations on the metric and on the measure are controlled by constants a and b separately. If (M, gM , μM ) is (1, b)-close to (M, g M , ν M ), then gM = g M . On the other hand, if the former is (a, 1)-close to the latter, then the perturbation is purely on metric. However in this case, according to Proposition 4.1, the measure is perturbed as well. Lastly, if the former are (1, 1)-close to the latter, then they are obviously isometric.

It's worth noting that our perturbation model assumes only the basic smoothness of gM and ρ M without more restrictive constraints such as bounded higher-order derivatives, thus it allows to include a large class of perturbed manifolds.

Now let (M, g M , ν M ) and (N, g N , ν N ) be a pair of original Riemannian manifolds, and we perturb them to (M, gM , μM ) and (N, gN , μN ), which are (a M , b M )-close and (a N , b N )-close to the original ones respectively. We then study how the shape difference operators constructed from the perturbed pair deviate from the ones constructed from the original pair.

Bounded-distortion Condition

Throughout our analyses in the coming sections of this chapter, we assume the input Riemannian manifolds, (M, g M , ν M ) and (N, g N , ν N ), together with the map T between them satisfy the following boundeddistortion condition.

Condition 4.1 (Bounded-distortion) Let T F be the functional map induced by T : (M, g M , ν M ) → (N, g N , ν N ), the area and conformal distortion induced by T F (or equivalently by T ) are bounded:

For any w ∈ L 2 ν (N ), M T F (w) 2 dν M ≤ B T N w 2 dν N For any w ∈ H 1 0,ν (N ), M ⟨∇T F (w), ∇T F (w)⟩ g M dν M ≤ D T N ⟨∇w, ∇w⟩ g N dν N
where B T and D T are finite positive constants.

The bounded-distortion condition rules out singular maps which can be meaningless, for example, a map T that maps M to a single point of N . Though we only regularize the distortions between the original input Riemannian manifolds, the following proposition ensures the boundness of distortions between a perturbed pair. If (M, g M , ν M ), (N, g N , ν N ) and T satisfy condition 4.1, then the distortion between the perturbed manifolds are bounded as the follows:

For any w ∈ L 2 μ(N ), M T F (w) 2 dμ M ≤ BT N w 2 dμ N For any w ∈ H 1 0,μ (N ), M ⟨∇T F (w), ∇T F (w)⟩ gM dμ M ≤ DT N ⟨∇w, ∇w⟩ gN dμ N where BT = (a M a N ) n/2 b M b N B T , DT = (a M a N ) 1+n/2 b M b N D T
and n is the dimension of the manifolds.

Proof:

According to Proposition 4.3, we have a

-n/2 M b -1 M dμ M ≤ dν M ≤ a n/2 M b M dμ M and a -n/2 N b -1 N dμ N ≤ dν N ≤ a n/2 N b N dμ N .
Since T F (w) 2 and w 2 are both non-negative on, it follows from condition 4.1 that:

M T F (w) 2 dμ M ≤ M T F (w) 2 a n/2 M b M dν M ≤ a n/2 M b M B T N w 2 dν N ≤ a n/2 M b M B T N w 2 a n/2 N b N dμ N = (a M a N ) n/2 b M b N B T N w 2 dμ N = BT N w 2 dμ N .
The second inequality is proven in the same way:

M ⟨∇T F (w), ∇T F (w)⟩ gM dμ M ≤ M a M ⟨∇T F (w), ∇T F (w)⟩ g M dμ M ≤ a M M ⟨∇T F (w), ∇T F (w)⟩ g M a n/2 M b M dν M ≤ a 1+n/2 M b M D T N ⟨∇w, ∇w⟩ g N dν N ≤ a 1+n/2 M a 1+n/2 N b M b N D T N ⟨∇w, ∇w⟩ gN dμ N = DT N ⟨∇w, ∇w⟩ gN dμ N .
As we will demonstrate soon, condition 4.1 implies that both the shape difference operators, V and R, are bounded operators, i.e., their operator norms are finite. An important fact revealed by Proposition 4.4 is that if the operators constructed from the original pair of Riemannian manifolds are bounded, then the ones constructed from a perturbed pair are still bounded as long as the perturbations are finite.

Stability of the Shape Difference Operators

As in this chapter we are concentrating on Riemannian manifolds, for the sake of simplicity, from now on we denote by M the original Riemannian manifold (M, g M , dν M ) and by M the perturbed one (M, gM , μM ), unless stated otherwise.

The map T between manifolds is preserved while Riemannian structures are perturbed, and it induces two area-based shape difference operators, V and Ṽ , such that:

N fV (g)dν N = M T F (f )T F (g)dν M , ∀f, g ∈ L 2 ν (N ) (4.20) N f Ṽ (g)dμ N = M T F (f )T F (g)dμ M , ∀f, g ∈ L 2 μ(N ) (4.21)
The stability of the area-based shape difference operator is stated in the following theorem:

lim a M ,a N ,b M ,b N →1 + N (V g -Ṽ g) 2 dν N = 0 Proof: We first prove L 2 ν (N ) = L 2 μ(N ) so that Ṽ g is well-defined for g ∈ L 2 ν (N ). According to Proposition 4.3, we have a -n/2 N b -1 N dμ N ≤ dν N ≤ a n/2 N b N dμ N . Then for any f ∈ L 2 ν (N ), N f 2 dμ N ≤ N f 2 a n/2 N b N dν N = a n/2 N b N N f 2 dν N < ∞, therefore L 2 ν (N ) ⊂ L 2 μ(N ).
On the other hand, one can similarly verify that

L 2 μ(N ) ⊂ L 2 ν (N ), which implies L 2 ν (N ) = L 2 μ(N ). For f ∈ L 2 ν (N ) = L 2 μ(N )
, it follows from the triangle inequality that

| N fV fdν N - N f Ṽ fdν N | ≤ | N fV fdν N - N f Ṽ fdμ N | + | N f Ṽ fdμ N - N f Ṽ fdν N | := |P 1 | + |P 2 |
Then we estimate P 1 and P 2 separately. According to Proposition 4.3, measures ν M (resp.ν N ) and μM (resp.μ N ) satisfy

a -n/2 M b -1 M dμ M ≤ dν M ≤ a n/2 M b M dμ M and a -n/2 N b -1 N dμ N ≤ dν N ≤ a n/2 N b N dμ N .
Thus, we have

P 1 = M T F (f ) 2 dν M - M T F (f ) 2 dμ M ( by definitions of V, Ṽ ) ≥ (1 -a n/2 M b M ) M T F (f ) 2 dν M ( since a -n/2 M b -1 M dμ M ≤ dν M )
And similarly

P 1 ≤ (1 -a -n/2 M b -1 M ) N T F (f ) 2 dν M . Noticing that 0 ≤ M T F (f ) 2 dν M ≤ B T N f 2 dν N < ∞ as f ∈ L 2 ν (N ), we have |P 1 | vanishes as a M , b M → 1 + .
Regarding P 2 , we define two complementary subsets of N with respect to f :

I + = {x ∈ N : f (x) Ṽ f(x) ≥ 0} and I -= {x ∈ N : f (x) Ṽ f(x) < 0}. P 2 = I + f Ṽ f(dμ N -dν N ) + I - f Ṽ f(dμ N -dν N ) ≤ (1 -a -n/2 N b -1 N ) I + f Ṽ fdμ N + (1 -a n/2 N b N ) I - f Ṽ fdμ N
And we obtain the lower bound similarly:

P 2 = I + f Ṽ f(dμ N -dν N ) + I - f Ṽ f(dμ N -dν N ) ≥ (1 -a n/2 N b N ) I + f Ṽ fdμ N + (1 -a -n/2 N b -1 N ) I - f Ṽ fdμ N If I + f Ṽ fdμ N , I -f Ṽ fdμ N
are both finite, then P 2 vanishes as a N , b N tend to 1. In fact, by the formulation of Ṽ , both f and Ṽ f are functions in L 2 μ(N ), thus it follows from the Cauchy-Schwarz inequality that

| I ± f Ṽ fdμ N | ≤ I ± f 2 dμ N I ± ( Ṽ f) 2 dμ N ≤ N f 2 dμ N N ( Ṽ f) 2 dμ N < ∞. To summarize, f (V f -Ṽ f)dν N → 0 for any f ∈ L 2 ν (N ) = L 2 μ(N ). Given f, g ∈ L 2 ν (N ), since both V and Ṽ are self-adjoint operators, direct computation shows that 4 N f (V g -Ṽ g)dν N = N (f + g)(V -Ṽ )(f + g)dν N - N (f -g)(V -Ṽ )(f -g)dν N As f + g, f -g ∈ L 2 ν (N ), it follows immediately that for any pair of f, g, N f (V g -Ṽ g)dν N vanishes as a M , b M , a N , and b N converge to 1 simultaneously. Especially we let f = V g -Ṽ g ∈ L 2 ν (N ), and conclude that lim a M ,a N ,b M ,b N →1 + N (V g -Ṽ g) 2 dν N = 0.
Similar stability guarantee holds for the conformal shape difference operators as well. We start with defining the conformal shape difference operators for both pairs of manifolds.

N ⟨∇f, ∇R(g)⟩ g N dν N = M ⟨∇T F (f ), ∇T F (g)⟩ g M dν M , ∀f, g ∈ H 1 0,ν (N ) (4.22) N ⟨∇f, ∇R(g)⟩ gN dμ N = M ⟨∇T F (f ), ∇T F (g)⟩ gM dμ M , ∀f, g ∈ H 1 0,μ (N ) (4.23)
The following theorem suggests that as a M , b M , a N , b N converge to 1 simultaneously, the norm of the gradient of Rf -Rf converges to zero.

Theorem 4.3 Let M, N and M, Ñ be smooth n-dimensional manifolds under the same assumptions of Theorem 4.2, then for any f ∈ H 1 0,ν (N ), Rf is well-defined. Moreover, we have

lim a M ,b M ,a N ,b N →1 + N ⟨∇(Rf -Rf ), ∇(Rf -Rf )⟩ g N dν N = 0
This theorem is proven with the same idea of proving Theorem 4.2, we refer the readers to Section 4.8 for a detailed proof.

Remark 4.1 Our proofs for theorems 4.2 and 4.3 do not require the involved manifolds to be compact or boundaryless. The stability properties proven in this section are valid for any smooth n-dimensional Riemannian manifolds and maps satisfying condition 4.1.

Stability of the Shape Difference Operators in a Multi-Scale Framework

In this section, we study the stability properties of the shape difference operators in the framework of [START_REF] Ovsjanikov | [END_REF], where they are employed in a multi-scale way. The main results of this section are Theorem 4.4 and Theorem 4.6, of which the proofs heavily depend on properties of the (weighted) Laplace-Beltrami operators introduced in Section 4.3.2, thus we follow the assumptions there and emphasize that in this section, (N, g N ) is a connected, compact, smooth, n-dimensional Riemannian manifold without boundary.

We start by pointing out the connection between the multi-scale framework and the shape difference operators. Recall that Definition 4.15 defines a functional measuring the area distortion induced by T for a given function w ∈ L 2 ν (N ). Given a pair of manifolds M, N and a map T : M → N , let V be the area-based shape difference operator. It follows directly from Definition 4.18 that:

E(w) = M T F (w) 2 dν M N w 2 dν N = N wV (w)dν N N w 2 dν N Since V is a positive-definite self-adjoint operator acting on L 2 ν (N ), the maximum of E(w) within L 2 ν (N ) is simply the L 2 -norm of V .
The framework of [START_REF] Ovsjanikov | [END_REF]] computes the constrained norm of V with respect to a special collection of subdomains of L 2 ν (N ) -the function spaces spanned by the first several eigenfunctions of the LB operator on N . In general, given a subdomain Ω of L 2 ν (N ), the maximum of E(w) constrained in Ω provides a quantitative characterization of to what extent V can distort functions in Ω. The maximizer (which we call the indicator after), w * , is a function in Ω that is the most distorted by V .

An Alternative Collection of Multi-scale Subdomains

Recall that in Section 4.3.4, we review the subdomain constructed at scale k: S

(k) = span{ϕ 1 , ϕ 2 , • • • , ϕ k },
where ϕ i is the i-th eigenfunction of the LB operator on N. Despite several advantages listed in Section 4.3.4, this subdomain construction suffers issues that are rooted in its discrete nature and is not suitable for our stability analysis.

First, we do not have a general criterion of selecting k. In general, he spectrum varies across shapes that are not isometric to each other, therefore we need to choose k according to the spectrum of the Laplace-Beltrami operator on the given input shape N .

Second, it can lead to confusing results when the truncation is not done appropriately. If a eigenvalue is of multiplicity more than 1, then as discussed in Section 4.3.2, there might be a subspace spanned by the eigenfunctions associated to this eigenvalue, which is of dimension more than 1. Though we argue that we can pick a orthogonal basis for this subspace, there does not exist a natural order for this basis functions (since the associated eigenvalues are identical) and thus we order them arbitrarily. Now if we truncate in the middle of this subspace, the resulting subdomain S(k) carries the randomness of the ordering. Therefore the constrained norm with respect to the space spanned by the first several eigenfunctions is not even welldefined, In practice, we observe instability in the more subtle case of analyzing conformal differences (see Figure 4.3).

indicators. However, it is not always the case. For example, imagine that we deform the bottom of shape M in Figure 4.1 so that the deformation from M to N is symmetrical. In this case, at every scale, the maximum of E(w) is realized by two indicators w t , w b which highlight respectively the top and the bottom of shape N , therefore we will no longer observe consistency in indicators. We then turn to study the stability of the maxima with respect to the changes in scale.

Instead of the original collection of subdomains S(k) indexed by the integer k, we consider the new collection of multi-scale subdomains A(C) that is controlled by the continuous parameter C. As mentioned before, the maximum of E(w) constrained in a subdomain of L 2 ν (N ) is a constrained norm of V . For subdomain A(C), we define:

∥V ∥ C = max E(w) s.t. w ∈ A(C)
Let C go through interval [0, +∞), then we obtain a curve (C, ∥V ∥ C ). The following theorem justifies its continuity.

Theorem 4.4 Given two connected compact smooth Riemannian manifolds M and N , and a map T between them. If M, N, T satisfy condition 4.1, then for any positive constant C > 0, C ′ = C + ε > 0, we have:

∥V ∥ C ′ -∥V ∥ C ≤ 4B T |ε|/C + 2B T |ε|/C.
We first prove an auxiliary lemma.

Lemma 4.1 Under the same assumptions of Theorem 4.4, for any w 0 , w 1 ∈ L 2 ν (N ) such that neither N w 2 0 dν N nor N (w 0 +w 1 ) 2 dν N is zero, the absolute difference between E(w 0 +w 1 ) and E(w 0 ) is bounded as below:

|E(w 0 + w 1 ) -E(w 0 )| ≤ 4B T √ r + 2B T r
where B T is the constant in condition 4.1 and r = N w 2 1 dν N N w 2 0 dν N .

Proof: We estimate the difference by two parts:

|E(w 0 + w 1 ) -E(w 0 )| ≤ E(w 0 + w 1 ) -M T F (w 0 + w 1 ) 2 dν M N w 2 0 dν N + M T F (w 0 + w 1 ) 2 dν M N w 2 0 dν N -E(w 0 ) := P 1 + P 2
For P 1 , direct computation shows:

P 1 ≤ B T N (w 0 + w 1 ) 2 dν N 1 N w 2 0 dν N - 1 N (w 0 + w 1 ) 2 dν N ≤ B T 2| N w 0 w 1 dν N | + N w 2 1 dν N N w 2 0 dν N ≤ 2B T N w 2 0 dν N N w 2 1 dν N N w 2 0 dν N + B T N w 2 1 dν N N w 2 0 dν N ≤ 2B T √ r + B T r
The third line follows from the Cauchy-Schwarz inequality applied on w 0 and w 1 which are both in L 2 ν (N ). P 2 is bounded in the same way:

P 2 ≤ M T F (2w 0 w 1 + w 2 1 )dν M N w 2 0 dν N ≤ 2| M T F (w 0 )T F (w 1 )dν M | + B T N w 2 1 dν N N w 2 0 dν N ≤ 2B T N w 2 0 dν N N w 2 1 dν N N w 2 0 dν N + B T N w 2 1 dν N N w 2 0 dν N ≤ 2B T √ r + B T r
Putting them together yields

|E(w 0 + w 1 ) -E(w 0 )| ≤ 4B T √ r + 2B T r, where r = N w 2 1 dν N N w 2 0 dν N
We then prove Theorem 4.4:

Proof: We first consider ε > 0, i.e., C ′ > C. By definition, A(C) ⊂ A(C ′ ), thus ∥V ∥ C ′ -∥V ∥ C ≥ 0.
We then estimate the upper bound for the difference. Given w ∈ A(C + ε), our strategy is to construct a function w ∈ A(C), such that |E(w) -E( w)| is bounded.

If w itself lies in A(C), then it's trivial to set w = w. We now consider w ∈ A(C + ε)\A(C). We assume that N w(-∆ N )wdν N = (C + δ) N w 2 dν N , where 0 < δ ≤ ε. Since E(w) = E(cw), ∀c ̸ = 0, without loss of generality, we further assume that w = i≥1 a i ϕ i and i≥1 a 2 i = 1, where (ϕ i , λ i ) is the i-th eigensolution to ∆ N ϕ + λϕ = 0. According to Proposition 4.2, the constraint on w can be written as:

i≥2 a 2 i λ i = C + δ.
Let b 1 be a real number satisfying b 2 1a 2 1 = δ/C and b 1 a 1 ≥ 0. The existence of b 1 is assured by the fact that a 1 is finite (in fact |a 1 | ≤ 1). Then we set w = b 1 ϕ 1 + i≥2 a i ϕ i . Direct computation shows that

N w(-∆ N ) wdν N = i≥2 a 2 i λ i = C + δ and N w2 dν N = b 2 1 + i≥2 a 2 i = b 2 1 + 1 -a 2 1 = 1 + δ C = C+δ C = 1 C N w(-D N ) wdν N , meaning that w ∈ A(C). Thanks to Lemma 4.1, |E(w) -E( w)| = |E(w) -E(w + w -w)| ≤ 4B T √ r + 2B T r, where r = N ( w -w) 2 dν N / N w 2 dν N = (b 1 -a 1 ) 2 . Without loss of generality, we assume that a 1 , b 1 ≥ 0, thus b 1 = a 2 1 + δ/C > a 1 . Moreover, b 1 -a 1 = a 2 1 + δ/C -a 1 = δ/C √ a 2 1 +δ/C+a 1 ≤ δ/C. Therefore r ≤ δ/C ≤ ε/C. Now assuming that w * is the maximizer of E(w) constrained in A(C + ε), the above derivation shows that there exists a w * ∈ A(C) such that |E(w * ) -E( w * )| ≤ 4B T ε/C + 2B T ε/C. That implies ∥V ∥ C+ε -∥V ∥ C ≤ E(w * ) -E( w * ) ≤ |E(w * ) -E( w * )| ≤ 4B T ε/C + 2B T ε/C.
Regarding the case ε < 0, i.e., C ′ < C. We simply replace C and C + ε in the previous analysis with Cε and C, respectively. With identical derivations, for each w ∈ A(C), we construct a function

w ∈ A(C -δ) such that |E(w) -E( w)| = |E(w) -E(w + w -w)| ≤ 4B T √ r + 2B T r, where r ≤ |ε|/C. Similarly, we have ∥V ∥ C -∥V ∥ C+ε ≤ 4B T |ε|/C + 2B T |ε|/C for ε < 0.
Putting them together, we finish the proof of this theorem.

Notice that B T is in fact an upper-bound for the constrained norms, i.e., ∥V ∥ C ≤ B T , ∀C > 0. Thus the inequality proven in Theorem 4.4 only makes sense when ε is close to zero. On the other hand, the inequality suggests that for a perturbation of fixed magnitude |ε|, the larger C is, the more stable ∥V ∥ C is.

Stability With Respect To Perturbed Input Manifolds

In this section, we fix the scale C and add perturbations on M and N . As in Section 4.5, we perturb M and N to M and Ñ , which are (a M , b M )-close and (a N , b N )-close to the unperturbed ones respectively. V and Ṽ are the corresponding area-based shape difference operators defined in Equation 4.20 and 4.21.

In order to define the constrained norm for Ṽ , we first construct the corresponding functional Ẽ(w), which is obtained by changing measures in E(w).

Ẽ(w) = N w Ṽ (w)dμ N N w 2 dμ N = M T F (w) 2 dμ M N w 2 dμ N (4.25)
The corresponding subdomain Ã(C) can be written with respect to the weighted Laplace-Beltrami operator (see Definition 4.13). We denote by ∆N the Laplace-Beltrami operator on the weighted manifold (N, gN , μN ) and define

Ã(C) = {w : N w(-∆N )wdμ N ≤ C N w 2 dμ N }} (4.26)
One advantage of our construction of subdomains is that it allows to associate subdomains with respect to different manifolds. Particularly we observe the following interleaved structures between A(C) and Ã(C): Lemma 4.2 Let Ñ be a n-dimensional Riemannian manifold that is (a N , b N )-close to N , the two collections of subdomains defined in 4.24 and 4.26 are interleaved:

A(C) ⊂ Ã(Ca 1+n N b 2 N ) ⊂ A(Ca 2+2n N b 4
N ) Proof: For w ∈ A(C), it follows from the Definition 4.24 that N ⟨∇w, ∇w⟩ g N dν N ≤ C N w 2 dν N . According to Proposition 4.3, we have

N ⟨∇w, ∇w⟩ gN dμ N ≤ N a N ⟨∇w, ∇w⟩ g N dμ N ≤ N a N ⟨∇w, ∇w⟩ g N a n/2 N b N dν N ≤ Ca 1+n/2 N b N N w 2 dν N ≤ Ca 1+n/2 N b N N w 2 a n/2 N b N dμ N ≤ Ca 1+n N b 2 N N w 2 dμ N . meaning that w ∈ Ã(Ca 1+n N b 2 N ), thus A(C) ⊂ Ã(Ca 1+n N b N ).
Similarly we can derive the other inclusion relationship.

On the other hand, the functionals, E(w) and Ẽ(w) are as well related to each other.

Lemma 4.3 For any w ∈ L 2 ν (N ) satisfying E(w) > 0 , the ratio of Ẽ(w) to E(w) is two-sided bounded:

a -n/2 M a -n/2 N b -1 M b -1 N ≤ Ẽ(w)/E(w) ≤ a n/2 M a n/2 N b M b N Proof: According to Proposition 4.3, we have a -n/2 M b -1 M dμ M ≤ dν M ≤ a n/2 M b M dμ M and a -n/2 N b -1 N dμ N ≤ dν N ≤ a n/2
N b N dμ N And by definition:

Ẽ(w) = M T F (w) 2 dμ M N w 2 dμ N ≤ M T F (w) 2 a n/2 M b M dν M N w 2 a -n/2 N b -1 N dν N ≤ a n/2 M a n/2 N b M b N E(w).
Symmetrically, the other side is bounded:

Ẽ(w) = M T F (w) 2 dμ M N w 2 dμ N ≥ M T F (w) 2 a -n/2 M b -1 M dν M N w 2 a n/2 N b N dν N ≥ a -n/2 M a -n/2 N b -1 M b -1 N E(w).
Based on the above constructions of Ã(C) and Ẽ(w), the constrained norm in the perturbed case is defined as

∥ Ṽ ∥ C = max Ẽ(w) s.t. w ∈ Ã(C).
The main result of this section is stated in the following theorem, which claims that at each fixed scale C, the constrained norm is stable with respect to perturbations on the manifolds. Theorem 4.5 Let M, N be two connected compact smooth n-dimensional Riemannian manifolds without boundary, and T be a map from M to N . Let M (resp. Ñ ) be a smooth manifold that is (a M , b M )-close (resp.(a N , b N )-close) to M (resp.N ). V and Ṽ are the area-based shape difference operators constructed with M, N and M, Ñ respectively. If M, N, T satisfy condition 4.1, then at any fixed scale C, the following convergence is valid:

lim a M ,b M ,a N ,b N →1 + ∥ Ṽ ∥ C = ∥V ∥ C Proof: Given a fixed scale C > 0, we denote C 1 = Ca -1-n N b -2 N and C 2 = Ca 1+n N b 2 N . Let u 1 ∈ A(C 1 ), u 2 ∈ A(C 2 ) and v ∈ Ã(C) be functions satisfying E(u 1 ) = ∥V ∥ C 1 , E(u 2 ) = ∥V ∥ C 2 and Ẽ(v) = ∥ Ṽ ∥ C .
Thanks to Lemma 4.2, we have A(C 1 ) ⊂ Ã(C), thus according to lemma 4.3:

∥ Ṽ ∥ C ∥V ∥ C 1 = Ẽ(v) E(u 1 ) ≥ Ẽ(u 1 ) E(u 1 ) ≥ a -n/2 M a -n/2 N b -1 M b -1 N .
Lemma On the other hand, as Ã(C) ⊂ A(C 2 ), we have:

∥ Ṽ ∥ C ∥V ∥ C 2 = Ẽ(v) E(u 2 ) ≤ Ẽ(v) E(v) ≤ a n/2 M a n/2 N b M b N
Lastly, Putting the above two inequalities together, we have

a -n/2 M a -n/2 N b -1 M b -1 N ∥V ∥ C 1 ≤ ∥ Ṽ ∥ C ≤ a n/2 M a n/2 N b M b N ∥V ∥ C 2 It follows from Theorem 4.2 that ∥V ∥ C 2 -∥V ∥ C 1 ≤ 4B T C 2 -C 1 C 1 + 2B T C 2 -C 1 C 1 Obviously, letting a N , b N tend to 1, we have ∥V ∥ C 2 → ∥V ∥ C 1 .
Moreover, as a M , b M tend to 1 as well, according to the squeeze theorem we have

lim a M ,b M ,a N ,b N →1 + ∥ Ṽ ∥ C = ∥V ∥ C

Approximating ∥V ∥ C

By investigating the behavior of operator within the continuously evolving subdomains A(C), we have more stable and richer understanding of V than we've got from S(k). However, in practice, calculating ∥V ∥ C is far from being obvious. Since neither E(w) nor A(C) is convex, there is no guarantee on achieving the global optimum with constraint A(C).

For the sake of consistency, we denote by ∥V ∥ k the maximum of E(w) within subdomain S(k) . As discussed in [START_REF] Ovsjanikov | [END_REF]], computing ∥V ∥ k in the case where M and N are finite discrete meshed shapes is straightforward and is solved with an eigen-decomposition of a k by k matrix.

First note that the construction of A(C) and S(k) are closely related. The following proposition quantifies this relationship.

Proposition 4.6 Let M, N and T be a pair of manifolds and a map between them, which satisfy condition 4.1. Let λ k , λ k+1 be two consecutive eigenvalues of the LB operator on N, the constrained norms with respect to A(λ k ) and S(k) satisfy the following inequality:

0 ≤ ∥V ∥ λ k -∥V ∥ k ≤ 4B T λ k /λ k+1 + 2B T λ k /λ k+1
Proof: First of all, as shown in Preposition 4.5, S(k) is a proper subset of A(λ k ), which proves the left-side inequality.

Regarding the right-side inequality, we assume that w is the maximizer realizing ∥V ∥ λ k and decompose it into the eigenbasis:w = i≥1 a i ϕ i . Now let w = i≤k a i ϕ i , obviously w ∈ S(k) and therefore E( w) ≤ ∥V ∥ k .

We then estimate the difference E(w)-E( w), according to Lemma 4.1, it's upper-bounded by

4B T √ r+ 2B T r, where r = i≥k+1 a 2 i / i≥1 a 2 i . As w ∈ A(λ k ), it follows that λ k i≥1 a 2 i ≥ i≥1 a 2 i λ i ≥ i≥k+1 a 2 i + λ i ≥ λ k+1 i≥k+1 a 2 i , thus r ≤ λ k /λ k+1 . Plugging it back, we have ∥V ∥ λ k -∥V ∥ k ≤ E(w) -E( w) ≤ 4B T λ k /λ k+1 + 2B T λ k /λ k+1 .
As a direct corollary, if λ k /λ k+1 ≈ 0, then ∥V ∥ k is a good approximation of ∥V ∥ λ k . It's worth noting that this proposition indicates a general criterion of choosing k: it's preferable to choose k such that the gap between λ k and λ k+1 is significant. And as we will discuss soon, this proposition suggests that if the spectral gap is significant, then the maximizer realizing ∥V ∥ k is a nice candidate of the initial guess for iterative algorithms maximizing E(w) constrained in A(λ k ).

Secondly, a major obstacle of optimizing within A(C) is that it is of infinite dimension. Even in the discrete case, the problem scale is still determined by the number of points, which can range in the tens or hundreds of thousands. The following proposition suggests that there is a trade-off between complexity and accuracy:

Proposition 4.7 For a fixed parameter C, let ε > 0 and λ l+1 be the smallest eigenvalue of the LB operator on N such that C ≤ ελ l+1 . Let

∥V ∥ C,l = max E(w) s.t. w ∈ A(C) ∩ S(l) (4.27) , then ∥V ∥ C -∥V ∥ C,l is of order √ ε. Proof: Let function w ∈ A(C) satisfies E(w) = ∥V ∥ C .
We assume that w = i≥1 a i ϕ i , and let w be

i≤l a i ϕ i . Then obviously E(w) -∥V ∥ C,l ≤ E(w) -E( w), because w ∈ A(C) ∩ S(l) as well.
Repeating the proof of Proposition 4.6, we have

∥V ∥ C -E( w) ≤ 4B T √ r + 2B T r, where r = i≥l+1 a 2 i / i≥1 a 2 i . Similarly we can deduce that r ≤ C/λ l+1 ≤ ε, therefore ∥V ∥ C -∥V ∥ C,,l ≤ 4B T √ ε + 2B T ε.
Optimization of problem 4.27 is then conducted in a space of limited dimension, l, which is controlled by the users.

A brief summary: Before we proceed to analyze the stability properties of the conformal shape difference operators under this multi-scale framework, we make a brief review of the following results on the area-based shape difference operators.

• Section 4.6.2: (C, ∥V ∥ C ) is a continuous curve and always above point sets {(λ k , ∥V ∥ k )} +∞ k=1 . In Figure 4.4, the black curve is continuous and always above the blue points.

• Section 4.6.3: As M converges to M and Ñ converges to N , any point (C, ∥ Ṽ ∥ C ) on the black curve in Figure 4.4 converges to point (C, ∥V ∥ C ) on the red dashed one.

• Section 4.6.4: Discussion on the optimization problem: max E(w) s.t. w ∈ A(C).

Analysis for the Conformal Shape Difference Operator

In essence, with functional E(w), the framework of [START_REF] Ovsjanikov | [END_REF]] casts the problem of extracting information on shape (manifold) differences as a series of constrained optimization problems. The functional evaluates how the area-based shape difference operator distorts functions living in L 2 ν (N ). Note that the framework of [Rustamov 2013] introduces two shape difference operators which encode different types of differences between shapes (manifolds), a natural extension of the multi-scale framework of [START_REF] Ovsjanikov | [END_REF]] is to construct parallel functionals and subdomains with respect to the conformal shape difference operators, R.

1 = 0 k k+1 ( k , kV k k ) : (C, kV k C ) : (C, k Ṽ k C ) : ## a M , b M , a N , b N ! 1 C • • • • • • • • • k • k C
Given a pair of manifolds M, N , and a map T between them, we define a functional, F , acting on H 1 0,ν (N ) (the domain of the conformal shape difference operator) as the following:

F (w) = N ⟨∇w, ∇R(w)⟩ g N dν N N ⟨∇w, ∇w⟩ g N dν N = M ⟨∇T F (w), ∇T F (w)⟩ g M dν M N ⟨∇w, ∇w⟩ g N dν N
.

(4.28)

The second inequality follows from the definition of R.

The following lemma about F (w) is a counterpart of Lemma 4.1.

Lemma 4.4 Let M, N be two connected compact smooth Riemannian manifolds, and T be a map between them. If M, N, T satisfy condition 4.1, and for any w 0 , w 1 ∈ H 1 0,ν (N ) such that neither w 0 nor w 0 + w 1 is a constant function on N , then the difference between F (w 0 + w 1 ) and F (w 0 ) is bounded as below:

|F (w 0 + w 1 ) -F (w 0 )| ≤ 4D T √ s + 2D T s
where D T is the constant in condition 4.1 and s = N ⟨∇w 1 ,∇w 1 ⟩g N dν N N ⟨∇w 0 ,∇w 0 ⟩g N dν N .

On the other hand, modifying the multi-scale subdomain construction is necessary to suit the new functional. Let ∥R∥ C = max F (w) s.t. w ∈ A(C), and in the following we prove that ∥R∥ C ′ = ∥R∥ C , ∀C ′ > C > 0. In fact, given C ′ and a function on the boundary of A(C ′ ), repeating the proof of Theorem 4.4, we find another function w lying in A(C) and ww is a constant function. According to Lemma 4.4, we have

F (w) = F ( w) since N ⟨∇(w -w), ∇(w -w)⟩ g N dν N = 0. Therefore, ∥R∥ C = ∥R∥ C ′ if we use A(C) as the multi-scale constraint.
To solve this issue, we construct a new subdomain, A conf (C), which is orthogonal to the space spanned by constant functions.

A conf (C) = A(C) ∩ {w : N wdν N = 0} (4.29) and let ∥R∥ C = max F (w) s.t. w ∈ A conf (C). It's worth noting that if C < λ 2 , the first non-zero eigenvalue of -∆ N , then A conf (C) = ∅. That follows from Proposition 4.2: let w = i≥1 a i ϕ i , since N wdν N = 0, a 1 = 0 . Therefore N w(-∆ N )wdν N = i≥2 a 2 i λ i ≥ i≥2 a 2 i λ 2 ≥ λ 2 i≥2 a 2 i = λ 2 N w 2 dν N .
In the other words, if C < λ 2 , then A(C) ∩ {w :

N wdν N = 0} is empty. Thus C must be at least λ 2 so that ∥R∥ C is well-defined. In practice, it is more convenient to maximize F (w) in the subdomains spanned by finite number of eigenfunctions. Following the same arguments above, we modify S(k) to obtain

S conf (k) = span{ϕ 2 , ϕ 3 , • • • , ϕ k },
where k is at least 2.

With the above formulations, we validate the stability of R with respect to the changes in scale.

Theorem 4.6 Let M, N be two connected compact smooth Riemannian manifolds with boundary, and T be a map between them. Let λ 2 be the second eigenvalue (as well as the first non-zero eigenvalue) of -∆ N . If M, N, T satisfy condition 4.1, then for C > λ 2 , C ′ = C + ε > λ 2 we have:

∥R∥ C ′ -∥R∥ C ≤ 4D T λ 2 |ε| (C -λ 2 )(C -|ε|) + 2D T λ 2 |ε| (C -λ 2 )(C -|ε|)
Then we consider perturbations on the input manifolds. As before, we denote by M and Ñ the perturbed version of M and N . The perturbed conformal shape difference operator, R, is defined in Equation 4.23. The associated functional, F (w), is defined as the following: Unfortunately, the strategy of proving Theorem 4.5 does not work in the case of the conformal shape difference operators. That is because the interleaved structure described in Lemma 4.2 is not guaranteed between the new subdomains A conf (•) and Ãconf (•): a function satisfying N wdν N = 0 does not necessarily fulfill N wdμ N = 0 simultaneously.

F (w) = M ⟨∇T F (

Experimental Results

In this section, we demonstrate experimental results that are related to our theoretical analyses. We conduct all the experiments on meshed shapes, i.e., discrete polygon surfaces embedded R 3 . To start with, we show in the discrete case how to approach a local maximum of the constrained optimization problem max E(w) s.t. w ∈ A(C), which is discussed in Section 4.6.4. Then the extension proposed in Section 4.6.5 is applied to detect and visualize conformal differences between a pair of shapes induced by a given map. Lastly, we demonstrate how the frameworks based on the shape difference operators react to perturbations on input shapes.

Approximating ∥V ∥ C

Now suppose that we are given a pair of meshed shapes, we demonstrate how to search for a local optimum of the constrained non-linear optimization with the barrier function method. Let M, N be two meshed shapes consisting of v M and v N vertices respectively, and let L M (resp., L N ) and A M (resp., A N ) be the stiffness matrix and the matrix of area elements for M (resp.N ) (see [Pinkall 1993]). The functional map T F induced by T is represented by a matrix P ∈ R v M ×v N in the discrete setting. Let Φ k ∈ R v N ×k be a matrix whose columns are the first k eigenvectors solved by L N f = λA N f .

Then calculating ∥V ∥ C in this setting is equivalent to maximize the following function:

max f T P T A M P f f T A N f , s.t. f T L N f f T A N f ≤ C
Based on that a barrier function is constructed

G(β, f) = - f T P T A M P f f T A N f -βlog(C - f T L N f f T A N f )
As suggested in Proposition 4.6, we take the maximizer of E(w) constrained in S(k) as the initial guess for minimizing G(1, f). After obtaining f 1 as a local minimizer, we take it as the initial guess for G( 1 2 , f). The iteration continuous until there is no more significant improvement or β is sufficient small. Optimizing 4.27 is done with an extra constraint that f = Φ T l a, where a ∈ R l . Proposition 4.7 suggests that as l tends to infinity, the optimum of problem 4.27 converges to ∥V ∥ C . We illustrate this claim by comparing a bumped sphere M with a sphere N (the same pair as in As mentioned in Section 4.6.1, both subdomains S(k) and A(C) are designed to control the Dirichlet's Energy of feasible solutions. The difference between them is that in the former case the energy is controlled by truncating high frequency components while in the latter case high frequency components are allowed but with implicit bounds on their weights. Figure 4.6 demonstrates this difference intuitively. In this experiment, we compute the local maxima/maximizers of ∥V ∥ C with different scales C range from 0.5 to 2.

Capturing Conformal Differences

In Section 4.6.5 we introduce a framework to analyze and visualize conformal differences, which can't be captured by the original framework in [START_REF] Ovsjanikov | [END_REF]. A simple example is demonstrated in Figure 4.7, where the map T from M to N is an area-preserving map. Therefore E(w) = 1 for any function w and analysis based on E(w) does not provide any information about the differences between the two surfaces.

On the other hand, we maximize F (w) in S conf (50), and plot the indicator on N (see the right plot of Figure 4.7). It is obvious that this indicator highlights the regions where changes in angles take place.

In practice, we also observe stability of the multi-scale framework based on the conformal shape difference operator. We compare again the two horses in In general, the conformal differences between shapes are less intuitive to imagine and more subtle to capture. Thus this extension based on the conformal shape difference operator helps to provide a more complete picture for users to understand the distortions between shape induced by a given map.

Stability of the Area-based Shape Difference Operators

At the beginning of this chapter, we have shown in Figure 4.1 the robustness of the multi-scale framework in different senses. Here we consider more complicated shapes than deformed spheres.

The first example is demonstrated in Figure 4.9, we compare the two horses in Besides changing the mesh structure, we perturb the input meshes by disturbing the vertices. In this example, the vertices are especially perturbed along the normal direction so that the point-to-point correspondences are roughly preserved. We first compute the mean distance of edges of each mesh: dM = 0.0144, dN = 0.0141, and the vertex normal vectors. Given a parameter σ, we perturb a point p of mesh M to p ′ = p + σ dM x p • n p , where x p is a one-dimensional random variable distributed normally with mean 0 and variance 1, and n p is the unit normal vector at vertex p. And we use the original mesh connectivity to connect perturbed points, since they are in a one-to-one correspondence to the unperturbed points.

We perturb both M and N in the same manner, and consider 4 choices of σ: 0, 0.1, 0.5, 1.0. At each level of perturbations, we generate indicators with respect to the area-based shape difference operator at 3 scales k = 20, 50, 200. The results are shown in Figure 4.10. We observe that when σ = 0.1, the indicators are well consistent with the ones of the first row. In fact, even when σ = 0.5, meaning that the standard deviation of the perturbations is half the mean distance, the indicators are still reasonable. At the end, we also notice that in the most noisy row, the indicator at k = 200 deviates from the ground-truth significantly while the first two at k = 20, 50 are still relevant. As we mentioned before, as k increases, the corresponding indicator is supposed to be more and more localized. The high-frequency indicators are more sensitive to the noises.

4.8 Proofs for Theorems in Section 4.6.5

Proof of Theorem 4.3 Proof: We first prove H 1 0,ν (N ) = H 1 0,μ (N ) so that Rg is well-defined for g ∈ H 1 0,ν (N ). Since (N, g N , ν N ) and (N, gN , μN ) are on top of the same smooth topological manifold N , they share the same boundary of N . Therefore if f is zero on the boundary of the former, then it is zero on the boundary of the latter.

We now prove if N f 2 + ⟨∇f, ∇f ⟩ g N dν N < ∞, then N f 2 + ⟨∇f, ∇f ⟩ gN dμ N < ∞ and vice versa. First according to Proposition 4.3, we have a

-n/2 N b -1 N dμ N ≤ dν N ≤ a n/2 N b N dμ N and for any smooth function f, a -1 N ⟨∇f, ∇f ⟩ gN ≤ ⟨∇f, ∇f ⟩ g N ≤ a M ⟨∇f, ∇f ⟩ gN . Then for any f ∈ H 1 0,ν (N ), N f 2 + ⟨∇f, ∇f ⟩ gN dμ N ≤ N f 2 a n/2 N b N dν N + N a N ⟨∇f, ∇f ⟩ g N a n/2 N b N dν N ≤ a n/2 N b N ( N f 2 dν N + a N N ⟨∇f, ∇f ⟩ g N dν N ) < ∞, therefore H 1 0,ν (N ) ⊂ H 1 0,μ (N )
On the other hand, one can similarly verify that H 1 0,μ (N ) ⊂ H 1 0,ν (N ), which implies H 1 0,ν (N ) = H 1 0,μ (N ). Now considering f ∈ H 1 0,ν (N ) = H 1 0,μ (N ), it follows from the triangle inequality that

| N ⟨∇f, ∇Rf ⟩ g N dν N - N ⟨∇f, ∇ Rf ⟩ g N dν N | ≤| N ⟨∇f, ∇Rf ⟩ g N dν N - N ⟨∇f, ∇ Rf ⟩ gN fdμ N | + | N ⟨∇f, ∇ Rf ⟩ gN dμ N - N ⟨∇f, ∇ Rf ⟩ gN dν N | + | N ⟨∇f, ∇ Rf ⟩ gN dν N - N ⟨∇f, ∇ Rf ⟩ g N dν N | :=|P 1 | + |P 2 | + |P 3 |
Then we estimate P 1 , P 2 and P 3 separately. According to Proposition 4.3, measures ν M (resp.ν N ) and μM (resp.μ N ) satisfy

a -n/2 M b -1 M dμ M ≤dν M ≤ a n/2 M b M dμ M a -n/2 N b -1 N dμ N ≤dν N ≤ a n/2 N b N dμ N .
and for any smooth functions f M , f N on M and N respectively,

a -1 M ⟨∇f M , ∇f M ⟩ gM ≤ ⟨∇f M , ∇f M ⟩ g M ≤ a M ⟨∇f M , ∇f M ⟩ gM a -1 N ⟨∇f N , ∇f N ⟩ gN ≤ ⟨∇f N , ∇f N ⟩ g N ≤ a N ⟨∇f N , ∇f N ⟩ gN
Thus, we have by the definitions of R, R

P 1 = M ⟨∇T F (f ), ∇T F (f )⟩ g M dν M - M ⟨∇T F (f ), ∇T F (f )⟩ gM dμ M ≤ M ⟨∇T F (f ), ∇T F (f )⟩ g M dν M -a -1 M M ⟨∇T F (f ), ∇T F (f )⟩ g M dμ M ≤ (1 -a -1-n/2 M b -1 M ) M ⟨∇T F (f ), ∇T F (f )⟩ g M dν M
The lower bound of P 1 is estimated in the same way, and we have

P 1 ≥ (1 -a 1+n/2 M b M ) M ⟨∇T F (f ), ∇T F (f )⟩ g M dν M Noticing that 0 ≤ M ⟨∇T F (f ), ∇T F (f )⟩ g M dν M ≤ D T N ⟨∇f, ∇f ⟩ g N dν N < ∞ as f ∈ H 1 0,ν (N ), we have |P 1 | vanishes as a M , b M → 1 + .
Regarding P 2 , we define two complementary subsets of N with respect to f :

I + = {x ∈ N : ⟨∇f, ∇ Rf ⟩ gN ≥ 0} and I -= {x ∈ N : ⟨∇f, ∇ Rf ⟩ gN < 0}. P 2 = I + ⟨∇f, ∇ Rf ⟩ gN (dμ N -dν N ) + I - ⟨∇f, ∇ Rf ⟩ gN (dμ N -dν N ) ≤ (1 -a -n/2 N b -1 N ) I + ⟨∇f, ∇ Rf ⟩ gN dμ N + (1 -a n/2 N b N ) I - ⟨∇f, ∇ Rf ⟩ gN dμ N
And we obtain the lower bound similarly:

P 2 = I + ⟨∇f, ∇ Rf ⟩ gN (dμ N -dν N ) + I - ⟨∇f, ∇ Rf ⟩ gN (dμ N -dν N ) ≥ (1 -a n/2 N b N ) I + ⟨∇f, ∇ Rf ⟩ gN dμ N + (1 -a -n/2 N b -1 N ) I - ⟨∇f, ∇ Rf ⟩ gN dμ N
If the integrals of ⟨∇f, ∇ Rf ⟩ gN on I + and I -are both finite, then P 2 vanishes as a N , b N tend to 1. In fact, by the formulation of R, both f and Rf are functions in H 1 0,μ (N ), thus it follows from the Cauchy-Schwarz inequality that

| I ± ⟨∇f, ∇ Rf ⟩ gN fdμ N | ≤ I ± ⟨∇f, ∇f ⟩ gN dμ N I ± ⟨∇ Rf, ∇ Rf ⟩ gN dμ N ≤ N ⟨∇f, ∇f ⟩ gN dμ N N ⟨∇ Rf, ∇ Rf ⟩ gN dμ N < ∞
Lastly, we argue that P 3 vanishes as a M , a N , b M and b N tend to 1. Since the metric defines a symmetrical inner products on the tangent spaces, we have ⟨∇f, ∇ Rf

⟩ gN = 1 4 (⟨∇(f + Rf ), ∇(f + Rf )⟩ gN -⟨∇(f - Rf ), ∇(f -Rf )⟩ gN ) and ⟨∇f, ∇ Rf ⟩ g N = 1 4 (⟨∇(f + Rf ), ∇(f + Rf )⟩ g N -⟨∇(f -Rf ), ∇(f -Rf )⟩ g N ).
On the other hand, as a N , b N → 1, the following convergence holds

lim a N ,b N →1 ⟨∇(f ± Rf ), ∇(f ± Rf )⟩ gN = ⟨∇(f ± Rf ), ∇(f ± Rf )⟩ g N Therefore we have lim a N ,b N →1 ⟨∇f, ∇ Rf ⟩ gN = ⟨∇f, ∇ Rf ⟩ g N which assures that P 3 vanishes as a N , b N → 1. To summarize, ⟨∇f, ∇(Rf -Rf )⟩ g N dν N → 0 for any f ∈ H 1 0,ν (N ) = H 1 0,μ (N ). Given f, g ∈ H 1
0,ν (N ), since both V and Ṽ are self-adjoint operators, direct computation shows that

4 N ⟨∇f, ∇(Rg -Rg)⟩ g N dν N = N ⟨∇(f + g), ∇(R -R)(f + g)⟩ g N dν N - N ⟨∇(f -g), ∇(R -R)(f -g)⟩ g N dν N
As f +g, f -g ∈ H 1 0,ν (N ), it follows immediately that for any pair of f, g, N ⟨∇f, ∇(Rg -Rg)⟩ g N dν N vanishes as a M , b M , a N , and b N converge to 1 simultaneously. Especially we let f = Rg -Rg ∈ H 1 0,ν (N ), and conclude that

lim a M ,a N ,b M ,b N →1 + N ⟨∇(V g -Ṽ g), ∇(V g -Ṽ g)⟩ g N dν N = 0.
Proof of Lemma 4.4 Proof: We estimate the difference by two parts:

|F (w 0 + w 1 ) -F (w 0 )| ≤ F (w 0 + w 1 ) -M ⟨∇T F (w 0 + w 1 ), ∇T F (w 0 + w 1 )⟩ g M dν M N ⟨∇w 0 , ∇w 0 ⟩ g N dν N + M ⟨∇T F (w 0 + w 1 ), ∇T F (w 0 + w 1 )⟩ g M dν M N ⟨∇w 0 , ∇w 0 ⟩ g N dν N -F (w 0 ) := P 1 + P 2
For P 1 , direct computation shows:

P 1 ≤ D T N ⟨∇(w 0 + w 1 ), ∇(w 0 + w 1 )⟩ g N dν N 1 N ⟨∇w 0 , ∇w 0 ⟩ g N dν N - 1 N ⟨∇(w 0 + w 1 ), ∇(w 0 + w 1 )⟩ g N dν N ≤ D T 2| N ⟨∇w 0 , ∇w 1 ⟩ g N dν N | + N ⟨∇w 1 , ∇w 1 ⟩ g N dν N N ⟨∇w 0 , ∇w 0 ⟩ g N dν N ≤ 2D T N ⟨∇w 0 , ∇w 0 ⟩ g N dν N N ⟨∇w 1 , ∇w 1 ⟩ g N dν N N ⟨∇w 0 , ∇w 0 ⟩ g N dν N + D T N ⟨∇w 1 , ∇w 1 ⟩ g N dν N N ⟨∇w 0 , ∇w 0 ⟩ g N dν N ≤ 2D T √ s + D T s
The third line follows from the Cauchy-Schwarz inequality applied on w 0 and w 1 which are both in L 2 ν (N ). P 2 is estimated in the same way, and the bound is identical:

P 2 ≤ 2 M ⟨∇T F (w 0 ), ∇T F (w 1 )⟩ g M + ⟨∇T F (w 1 ), ∇T F (w 1 )⟩ g M dν M N ⟨∇w 0 , ∇w 0 ⟩ g N dν N ≤ 2| M ⟨∇T F (w 0 ), ∇T F (w 1 )⟩ g M dν M | + D T N ⟨∇w 1 , ∇w 1 ⟩ g N dν N N ⟨∇w 0 , ∇w 0 ⟩ g N dν N ≤ 2D T N ⟨∇w 0 , ∇w 0 ⟩ g N dν N N ⟨∇w 1 , ∇w 1 ⟩ g N dν N N ⟨∇w 0 , ∇w 0 ⟩ g N dν N + D T N ⟨∇w 1 , ∇w 1 ⟩ g N dν N N ⟨∇w 0 , ∇w 0 ⟩ g N dν N ≤ 2D T √ s + D T s
Putting them together yields |F (w 0 + w 1 ) -F (w 0 )| ≤ 4D T √ s + 2D T s, where s = N ⟨∇w 1 ,∇w 1 ⟩g N dν N N ⟨∇w 0 ,∇w 0 ⟩g N dν N Proof of Theorem 4.6 Proof: The strategy of this proof is similar to the one for Theorem 4.4. We first consider the case of ε > 0, i.e., C ′ > C > λ 2 . By definition, A conf (C) ⊂ A conf (C ′ ), thus ∥R∥ C ′ -∥R∥ C ≥ 0. We then estimate the upper bound for the difference. Given w ∈ A conf (C + ε), our strategy is to construct a function w ∈ A conf (C), such that |E(w) -E( w)| is uniformly bounded.

If w itself lies in A conf (C), then it's trivial to set w = w. We now consider the case w ∈ A conf (C + ε)\A conf (C). Assume that N w(-∆ N )wdν N = (C + δ) N w 2 dν N , where 0 < δ ≤ ε. Without loss of generality, we further assume that w = i≥1 a i ϕ i and i≥1 a 2 i = 1, where (ϕ i , λ i ) is the i-th eigensolution to ∆ N ϕ + λϕ = 0. Note that N wdν N = 0, thus a 1 = 0, i.e., w = i≥2 a i ϕ i . According to the orthogonality of eigenbasis, the constraint on w can be written as

ϕ 2 + i≥3 a i ϕ i . Direct computation shows that N w(-D N ) wdν N = b 2 2 λ 2 + i≥3 a 2 i λ i = C + Cδ C-λ 2 and N w2 dν N = b 2 2 + i≥3 a 2 i = b 2 2 + 1 -a 2 2 = 1 + δ C-λ 2 = 1 C N w(-D N ) wdν N , meaning that w ∈ A conf (C). Thanks to Lemma 4.4, |F (w) -F ( w)| = |F (w) -F (w + w -w)| ≤ 4D T √ s + 2D T s, where s = N ⟨∇( w-w),∇( w-w)⟩g N dν N N ⟨∇w,∇w⟩g N dν N = λ 2 (b 2 -a 2 ) 2

C+δ

. Without loss of generality, we assume that

a 2 , b 2 ≥ 0, thus b 2 = a 2 2 + δ C-λ 2 > a 2 . Moreover, b 2 -a 2 = a 2 2 + δ C-λ 2 -a 2 = δ/(C-λ 2 ) √ a 2 2 +δ/(C-λ 2 )+a 2 ≤ δ/(C -λ 2 ). Therefore s ≤ λ 2 δ (C-λ 2 )(C+δ) ≤ λ 2 ε (C-λ 2 )C
. Now assuming that w * is the maximizer of F (w) constrained in A conf (C + ε), the above derivation shows that there exists a w *

∈ A conf (C) such that |F (w * ) -F ( w * )| ≤ 4D T √ s+2D T s. That implies ∥V ∥ C+ε -∥V ∥ C ≤ E(w * )-E( w * ) ≤ |E(w * )-E( w * )| ≤ 4D T λ 2 ε C(C-λ 2 ) + 2D T λ 2 ε C(C-λ 2 ) .
Regarding the case ε < 0, i.e., C > C ′ > λ 2 . We simply replace C and C + ε in the previous analysis with Cε and C. With identical derivations, for each w ∈ A conf (C), we construct a function

w ∈ A conf (C -δ) such that |F (w) -F ( w)| = |F (w) -F (w + w -w)| ≤ 4D T √ s + 2D T s, where s ≤ λ 2 |ε| (C-ε)(C-λ 2 ) .
Putting them together, we finish the proof of Theorem 4.6. The framework of the functional map and most of its follow-ups are proposed to deal with problems about 3D shapes, i.e., 2-dimensional Riemannian manifold embedded in R 3 . Besides the fact that 3D shapes are ubiquitous in both theory and application, the other reason might be that in practice, 3D shapes can be wellapproximated by polygon meshes, which allow people to estimate geometric quantities of the underlying shape from the discrete data. For example, the eigenbases with respect to the Laplace-Beltrami operators are heavily utilized in the frameworks related to the functional maps. Approximating this operator on a meshed shape is well-studied and solved by reliable and efficient methods (e.g., the cotangent scheme proposed in [Pinkall 1993]). However, as we mentioned in Chapter 4, the frameworks of [START_REF] Ovsjanikov | [END_REF][START_REF] Ovsjanikov | [END_REF], Rustamov 2013] can be naturally adapted to a more general setting where the intrinsic dimension of input manifolds is not necessarily 2. Stated differently, these frameworks are potentially suitable for analyzing objects of higher dimensions. Constructing polygon meshes, unfortunately, is not straightforward in this more abstract setting. Moreover, only considering 2-dimensional manifolds, we would encounter the same problem when the manifolds are embedded in an Euclidean space of dimension more than 3. Even in the very original setting, generating mesh on top of a bunch of sampling points is not trivial. For instance, noises are sort of inevitable in the data acquisition, sometimes they give rise to locally non-manifold structures and thus make mesh generation a challenging task.

Due to the above limitations, in this chapter we consider implementing the above mentioned functionalmaps-based frameworks on data in a more primitive form-point cloud data (PCD). There are nice works on approximating geometric quantities with PCD. The framework of [START_REF] Belkin | Constructing Laplace operator from point clouds in Rd[END_REF]] approximates the Laplace-Beltrami operator in PCD sampled from a general d-dimensional Euclidean space. However, this approach and a later one [START_REF] Liu | [END_REF]] both rely on reconstructing local mesh structure. Rather than falling back to build some specific structures to associate input points, we build simply k-nearest neighborhood (k-NN) graphs on the input PCD, and for all necessary ingredients for implementing the functional-maps-based frameworks that are computed with meshes, we propose a pipeline to construct their counterparts on top of the k-NN graphs.

Among the ingredients, perhaps a reliable approximation to the Laplace-Beltrami operator is the most important. The top candidate for our pipeline is the graph Laplacian, whose convergence to the Laplace-Beltrami operator on the underlying manifold has been well studied (see [START_REF] Belkin | Towards a theoretical foundation for Laplacian-based manifold methods[END_REF], Hein 2007]). And it has been wildly applied in practice, say, spectral clustering [von Luxburg 2006], dimensional reduction [START_REF] Belkin | Laplacian Eigenmaps and Spectral Techniques for Em-bedding and Clustering[END_REF], Tenenbaum 2000a]. However, it is not so popular in the area of geometry processing, one of the main reasons is the lack of robustness: for example, in practice it is not clear how a graph Laplacian varies with respect to different k for k-NN graphs, not to mention that we are unlikely to approximate the underlying Laplace-Beltrami operator with graph Laplacian based on naive graph construction such as k-NN over finite sampling points. On the other hand, the most important observation in this chapter is that despite its instablility when applied in analyzing a single object, in the context of the functional-maps-based frameworks where objects are analyzed in pair, the results gained with the graph Laplacian are surprisingly stable. Figure 5.1 shows parts of the eigenfunctions on a horse generated from the mesh setting and from our PCD setting, where two rows are distinct from each other. However, the indicators generated from the two settings are comparable.

It is worth noting that our approach is empirical: our pipeline design is inspired by results on the convergence of the graph Laplacians in [Hein 2007]. The results characterize the limit behavior of the graph Laplacian as the size of sample tends to infinity and are proven under certain technique conditions on the underlying manifolds and the sampling density, which can hardly be guaranteed in practice.

In the experimental section, we show the relevance and robustness of applying the functional-maps-based frameworks with our pipeline directly on PCD. Our pipeline is fit for point clouds sampled from Euclidean spaces of dimension more than 3, yet in this chapter we exhibit results generated from 3D shapes, so that we can compare our results to the ones from the mesh setting.

Pipeline for PCD

Before we present our pipeline, we quickly check the necessary ingredients for implementing the frameworks of functional maps [START_REF] Ovsjanikov | [END_REF]], shape difference operators [Rustamov 2013] and map analysis and visualization [START_REF] Ovsjanikov | [END_REF]].

In the mesh setting, for an input meshed shape M consisting of n M vertices, we need to compute W M and A M , which are both n M × n M matrices. W M is a stiffness matrix such as the standard cotangent weighted matrix. The formulation of W M with the cotangent scheme [Pinkall 1993] is illustrated in Figure 5.2.

W M (i, j) = 1 2 (cot↵ + cot ) v i v j ↵ v i v j ↵ W M (i, j) = 1 2 cot↵ (a) (b)
Figure 5.2: (a) W M (i, j) for an edge (v i , v j ) that is shared by two triangles. (b) W M (i, j) for an boundary edge that belongs to a single triangle. If v i , v j are not connected by an edge, then W (i, j) = 0. Notice that in a mesh, an edge is at most shared by two triangles.

A M is a diagonal matrix, whose (i, i)-th entry is assigned with the area element of the i-th vertex. To compute A M (i, i), we use the mixed Voronoi cell proposed in [Meyer 2003]. We collect all the triangles in the mesh containing vertex v i . Each triangle contributes a part to the area element of v i : if it's an acute triangle, then the contribution is one third of its area; if it is a right or an obtuse triangle, the contribution is one half of its area when v i is opposite to the longest edge and one quarter of its area otherwise.

The discrete approximation of the Laplace-Beltrami operator is then given as A -1 M W M . To compute the eigenfunctions, we solve the eigenequation W M f = λA M f . (We refer the readers to [Botsch 2010] for an introduction to general mesh processing. )

Now suppose we are given a point cloud X = {x 1 , x 2 , • • • , x n } consisting of n points, and X is sampled from the same smooth shape that is approximated by M . In the following we present how we compute the counterparts A X and W X in the PCD setting.

Building Connectivity In the mesh setting, the polygons, say, triangles, govern connections among the discrete vertices. While in the PCD case, there are two common ways to construct a graph on top of a point cloud: (1) the ε-neighborhood graph construction and (2) the k-nearest neighborhood (k-NN) graph construction. We choose the latter over the former due to two reasons: firstly we can control the sparsity of the output graphs by intuitively tuning k; Secondly, we might face the situations where the points are non-uniformly sampled from a underlying Riemannian manifold. In such a situation, the parameter ε might be difficult to pick so that the resulting graph is connected (meaning that ε is not too small) and properly sparse (meaning that ε is not too large).

Building Weighted Graphs Fixing an integer k, we first build a k-NN graph on X, and denote the graph by G X = (X, E X ), where the edge set E X = {(x i , x j ), x i ∈ N (x j , k) or x j ∈ N (x i , k)} and N (x i , k) is the set of the nearest k neighborhoods of x i among X\x i .

Then we turn this k-NN graph into a weighted graph. Typically, we use the Gaussian kernel to assign a weight to each of the connected edges in the k-NN graph above:

K X (i, j) = exp(-∥x i -x j ∥ 2 /2t 2 ) if (x i , x j ) ∈ E X 0 otherwise (5.1)
t is a parameter of K X to be determined. We in practice let t be

|X| i=1 x j ∈N (x i ,k) ∥x i -x j ∥ k|X| , which is depending on the scale of X. let D X (i) = 1 deg(i) (i,j)∈E X K X (i, j)
, where deg(i) is the degree of vertex x i in graph G X . Then we normalize K X with the D X (i) ′ s:

KX (i, j) = K X (i, j) D X (i)D X (j) (5.

2)

Building graph Laplacian There are various ways of defining a graph Laplacian, We use the unnormalized weighted graph Laplacian with respect to KX constructed above:

W X (i, j) = ⎧ ⎪ ⎨ ⎪ ⎩ -KX (i, j) if i ̸ = j j KX (i, j) if i = j (5.3)
Estimating Weighted Measures In the mesh setting, to each vertex p, we associate an area element, which is positively correlated to the sum of the areas of triangles around that point. In the point cloud setting, however, how to compute an area element is not clear because the surface is not regularly tessellated with polygons.

In our pipeline, to each sampling point x i , we assign the multiplicative inverse of the sampling density estimation at x i . There is a considerable amount of literature on density estimation from sampling, especially we chose the framework of [Biau 2011], which estimates the sampling density at point x i ∈ R d proportionally to:

ρ(x i ) ∝ ( x j ∈N (x i ,kn) ∥x i -x j ∥ 2 ) -d 2
The k n in the subscript is a parameter related to n, the total number of sampling points. The convergence of the estimator to the underlying density is assured under some conditions, one of which is that k n → ∞, k n /n → 0 as n → ∞. In practice, we simply let k n be the same as k in our k-NN graph construction. Therefore we obtain these quantities as a by-product of constructing the k-NN graph.

Since we mostly consider 2-dimensional Riemannian manifolds embedded in R 3 . A X is then obtain by letting A X (i, i) = (

x j ∈N (x i ,k) ∥x i -x j ∥ 2 ) 3 2 .
At the end, we view the graph Laplacian W X as a counterpart of the stiffness matrix in the mesh setting and A -1 X W X as an approximation to the Laplace-Beltrami operator, and compute the eigenfunctions by solving W X f = λA X f .

Rationale of Our Pipeline

According to one of the main results proposed in [Hein 2007] (see 'Main Result' at the beginning of section 3.3 there, in the case of 'unnormalized' with λ = 1), the graph Laplacian constructed in Equation 5.3 is viewed as an approximation of -ρ -1 ∆. From this point of view, the motivation of A X above is clear: since this approximation is affected by the sampling density, we let A -1 X be an approximation of ρ so that A -1 X W X is an approximation of ∆.

In fact, both A X and A M from different settings work in a similar way. To see this, we illustrate the area elements in Figure 5.3. We assume that all the points are sampled in a flat local area. First we consider the mesh case, if the sampling is uniform, then the area elements are equal to each other. However, the area element of point v 1 (the area circled by the blue segments) is larger than that of point v 2 (area circled by the red segments). The unbalance in area distribution compensates the effect of the non-uniform sampling: obviously the sampling density around point v 1 is smaller than that around v 2 . In the point cloud setting, we plot the 7-nearest neighborhoods of x 1 and x 2 and connect them with edges in different colors. It is easy to see that the A X (1, 1) > A X (2, 2), which mimics the relationship between the area element of v 1 and of v 2 .

Triangle Mesh

Point Cloud

v 1 v 2 x 2 x 1 M X Figure 5
.3: Left: the area elements in A M compensate the effect of non-uniform sampling density around points v 1 , v 2 ∈ M : larger element is assigned to the point locating in an area with lower sampling density; Right: the counterpart of A M in our pipeline, A X , that functions in a similar way.

On the other hand, as proven in [START_REF] Bernstein | [END_REF]], under certain sampling condition, the metric induced by the graph (especially by the kNN graph) converges to the intrinsic manifold metric with high probability.

Again we emphasize that our construction is experimental without rigorous theoretical guarantees. The convergence results in [Hein 2007] is proven under several technique conditions, which in our case are not satisfied.

Experimental Results

A brief review of the three functional-maps-based frameworks has been given in Section 4.3.3, Section 4.3.4 and Section 4.3.5. In this section, we discuss the how to perform these frameworks with only point clouds sampled from interested shapes. We assume that a map T is from a source shape S S to a target shape S T . Let a mesh M (resp.N ) and a point cloud X(resp.Y ) be two representations of S S (resp.S T ). The map T in the discrete case is simply a point-to-point correspondence from M to N (or from X to Y ). And the corresponding functional map T F is a pull back from the function space on N (resp.Y ) to the one on M (resp.X).

In Section 5.3.1, we compute the low-rank approximations of T F with truncated eigenbasis on both S S and S T and compare the results generated with M, N to those generated with X, Y .

After that, we turn to the frameworks of multi-scale map analysis and the shape difference operators. We induce T F with the known point-to-point correspondence T and use it in the following discrete settings. Therefore, in the mesh setting, the optimization problem of 4.17 is:

max E(w) = T F (w) T A M T F (w) w T A N w s.t. w ∈ span{φ N 1 , φ N 2 , • • • , φ N ks }
where φ N i is the i-th eigenvector of A -1 N W N . We call in this section the maximizer the area-based indicator at scale k s .

We write the functional with respect to the conformal shape difference operator as follows(see Equation 4.28)

max F (w) = T F (w) T W M T F (w) w T W N w s.t. w ∈ span{φ N 2 , • • • , φ N ks }
Note that the constraint is different. We then call the maximizer the conformal indicator at scale k s . As described in Section 5.2, we compute A X , W X as the counterparts of A M , W M in our pipeline. Therefore by replacing them with A X , W X , we are able to apply this framework on PCD.

The main target of our experiments is to demonstrate the relevance and empirical robustness of implementation of the functional-maps-based frameworks on PCD with the pipeline proposed in section 5.2. Especially, we run test on several collections of shapes. The ground truths of the differences among shapes within each of the collections are visually intuitive, and we use the results from the meshed setting as ground truths.

Functional Maps on PCD

As shown in Figure 5.4, we are given two meshed horses M and N and a list of point-to-point correspondences T from M to N . Let X and Y be the vertex set of M and N respectively. Each horse consists of 8431 vertices, thus all the function spaces related in this experiment are vector spaces of dimension 8431. Since a functional map is a linear operator across the function spaces on different shapes, it is a 8431 × 8431 matrix, which is inconvenient for analysis. In the following we compute low-rank approximations of it in different ways.

We start with meshes M and N . The first 100 eigenfunctions of W M f = λA M f and W N f = λA N f are solved as the truncated bases on M and N respectively. Based on these truncated eigenbases, we approximate and represent the functional map with a 100×100 matrix C T . With the full information of T , T F is faithfully represented as a permutation matrix. Thus C T is explicitly expressed as

C T (i, j) = ⟨T F (φ N i ), φ M j ⟩ M , 1 ≤ i, j ≤ 100
where φ N i is the i-th eigenvector on N and φ M j is the j-th eigenvector on N . On the other hand, it is also possible to compute C T with much less information. As shown in the top-left of To evaluate the qualities of approximations with respect to different settings and methods, we employ the method from [START_REF] Ovsjanikov | [END_REF]] to convert an approximated functional map, C ′ , to a point-to-point map T ′ . For each point p on M , the Euclidean distance from the ground truth T (p) to T ′ (p) (which are both on N ) is computed. The chart in Figure 5.4 shows comparison of functional maps approximated in different ways. A point (e, c) on that chart reads 100 × c% of the vertices satisfying ∥T (p) -T ′ (p)∥ ≤ e. As expected, with the ground truth T , the approximations are significantly more accurate than the ones without. And the approximation from the mesh setting is more accurate than the other two from the PCD setting. With the ground truth, the approximations in the PCD setting with different k are indistinguishable in this chart. In the more challenging situation, we observe that the using 40-NN graph results in a better approximation that is closer to the one with mesh information.

More intuitively, in the top-left of Figure 5.4, we plot a function f and T F (f ) on N and M respectively. The bottom row shows how f is transformed by functional maps computed with only landmark correspondences in different settings (for an easier comparison, we render all results with the mesh N ). Again the one regarding the 40-NN graph is closer to the mesh result, and visually the function is well transformed to Y (see the middle of the bottom row of 5.4). The result regarding the 15-NN graph, however, is not satisfying.

Selection of k for k-NN

In our PCD pipeline, the only parameter to be determined by the users is k, the number of nearest neighborhoods for constructing k-NN graphs. Usually we choose the same k for all point clouds analyzed in pair (or in collection). The previous experiment seems to suggest large k is more preferable. In the following we compare the selections from a different point of view.

For the sake of consistency, we consider the same pair of shapes as in Section 5.3.1 and follow the notations there. Recall that we analyze and extend a framework for detecting and highlighting the differences between shapes induced by a given map in Section 4.6. In this experiment we implement this method with our pipeline and especially we build k-NN graphs on the input point clouds with different k.

As before, the results obtained with meshed shapes M and N are seen as ground truths. Regarding the point clouds X and Y , we build k-NN graphs with k = 15, 25 and 40. Secondly, we take a look at the corresponding eigenfunctions in the right part of Figure 5.5. As mentioned in Section 5.1, one of the issues of applying graph Laplacian in geometry processing is its instability, which is evidenced in the plot. Judging from this short window (only four eigenfunctions from the ninth to the twelfth are plotted), we notice that when the number of nearest neighborhoods changes from 15 to 25, the eigenfunctions are distinct. The last two rows are more stable with respect to a change of k, yet they are still dissimilar to the top row. Moreover, we plot eigenfunctions associated with relatively small eigenvalues, the divergences are more prominent as we pick the ones of higher frequencies.

Lastly, we make some remarks on the selection of k for k-NN. From the implementation point of view, a smaller k results in a sparser graph, which expedites the computation. However, because of the presence of noise, small k would bring instability into the construction and therefore into the result (as discussed above and shown in Section 5.3.1). In the subsequent experiments, we select k = 40. As will be shown, this parameter works well in analyzing other collections of shapes.

Reliability of PCD Setting in Multi-scale Framework

We've observed consistency between the results from the mesh setting and the PCD setting in Figure 5.1 at a single scale, now we go for multiple scales. Two near-isometric human poses in comparison are depicted in Again we take the vertices from the meshes as the input point clouds. In Figure 5.6(c) and (e) are the corresponding results from the PCD setting. In the case of area distortion detection, the indicators in (b) and (c) match well. Though the stability of the conformal based multi-scale framework with respect to perturbed inputs is not proven in theory (see the discussion in Section 4.6.5), the outcomes regarding the conformal differences show some consistency: the highlighted areas evolve in the same pattern along rows (d) and (e)

Analyzing Shape Collections

The above experiments shows the stability of the shape difference operators for analyzing maps between a single pair of shapes in a multi-scale way. Now we perform a higher level analysis about recovering the intrinsic structure of a collection of shapes. As we prove in Section 4.5, the shape difference operators on their own are stable with respect to perturbations on input shapes. To demonstrate that, we repeat one of the experiments in [Rustamov 2013] (see figure 3 on page 7 there), but in the PCD setting. Given a collection of shapes with pairwise correspondences, our goal is to recover a reasonable layout in the intrinsic shape space. Two collections of shapes with different structures are considered-deformed spheres and galloping horses. In each collection we select the first shape as the base shape and construct both shape difference operators with maps from all shapes to the base shape. After approximating the areabased (resp. conformal) shape difference operators, we vectorize and embed them into R 2 with the PCA algorithm.

The top row of Figure 5.7 depicts the embeddings for the deformed spheres. Both layouts uncover the grid structure of the original shape collection. The results in [Rustamov 2013] suggest that in both the areabased and the conformal cases, the variances of the first two principal components are both close to 50 percents. Regarding our result: (1) Area-based case: though the sum of percentages add up to almost 100, the grid is unbalanced and stretched along the direction of the first principle component; (2) Conformal case: balance preserved, the shapes of the first and the second rows are not well differentiated, suggesting that the operators generated in the PCD setting fail to capture the small differences.

PCA layout for the area-based Shapes in dataset, indexed left to right and top to bottom shape di↵erence operators.

PCA layout for the conformal shape di↵erence operators. The bottom row shows the layouts for the galloping horse sequence. As demonstrated in the left panel, the sequence consists of two circles of continuous movements of the horse. Our results successfully capture the circular structure of the sequence: depicted in the layout, point i is close to point i+12 (i = 1, 2, • • • , 12). The group of i = 1, • • • 12 and the group of i = 13, • • • , 24 each form a circle. The result also reveals the fact that there are more conformal distortions than area distortions in this data, as the range of layout in the right panel is larger than that in the middle one.

Overall, we conclude from these experiments that although the results from the PCD setting is not as accurate as those from the mesh setting, they capture most of the desired information hidden in the data. Considering that we start from a much coarser understanding of the input shapes, these results are non-trivial and remarkable.

Non-uniformly Sampled Data

In the previous experiments, the point clouds are nearly uniformly sampled from the underlying shapes and clean without noises or outliers. In this section, we test our method with more tricky data.

In order to test the robustness of the pipeline with respect to non-uniform sampling, a collection of synthetic data is generated as in Figure 5.8. There are two pairs of shapes: S, B and S ′ , B ′ . The former pair is uniformly sampled, while in S ′ , B ′ , the densified half is triply denser than the unperturbed half. Below the shapes are area-based indicators generated with 4 possible maps between spheres S, S ′ and bumped spheres B, B ′ . Obviously the four indicators all highlight the correct area, regardless of how the sampling densities are perturbed. Further more, we take noisy point clouds into consideration. For the sake of simplicity, we only add noises to Y , therefore we can use the same map inherited from the mesh setting (though it is not surjective any more).

We first randomly select n p points in Y . Then for p from these n p points, we perturb p = (p x , p y , p z ) ∈ R 3 to (p x + dx, p y + dy, p z + dz) where dx, dy, dz are one-dimension random variables distributed normally with mean 0, and standard deviation d Y (d Y is the mean length of edges in mesh N ). Repeating the displacements r times for each p, we enlarge Y to Y ′ with n p r more points.

In this experiment, three pairs of PCD are involved: human poses (12500 points), horses (8431 points) and cats (7207 points). In Particularly, we find in the horse case that the indicator on Y ′ 2 is a bit off-it also highlights a part of the horse back, while the one on Y ′ 1 (to which 1000 more noisy points added than to Y ′ 2 ) is more consistent. We interpret this by comparing the ways noises are distributed: in Y ′ 1 , the noisy points are more decentralized, whereas in Y ′ 2 , more points are generated around each of the selected point in Y . Thus the sampling density is more distorted in Y ′ 2 , resulting in a less consistent indicator.

Limitation and Perspective

In this section, we show some limitations and perspectives of the pipeline proposed and empirically tested above.

Limitations

First, we point out some limitations of our pipeline. The first one is about scale. The way we build connectivity among points determines that the resulting graph is insensitive to local changes across objectives. On the one hand, it gives rise to robustness with respect to noisy data or even outliers (as shown in the last example in Section 5.3). On the other hand, it makes detecting differences of finer scales (i.e., larger k s ) between objectives difficult. In fact, as shown in Figure 5.7, the conformal shape difference operators from the PCD setting can not clearly differentiate the conformal changes between deformed spheres in the first and the second row depicted in the top-left panel. Plus, this limitation is the reason we fix the scale k s relatively small (mostly k s = 50).

The second limitation is that when the reach (see Definition 3.11) of the underlying shape is small, it is difficult to capture the right geometric information with the k-NN graphs. That is because the k-NN graphs are constructed with purely extrinsic information-the Euclidean distance between points. If there exist points that are distant on the shape but close in the ambient space, then the k-NN graph will connect them and introduce distortions. A typical example is illustrated in Figure 5.10. In the left-most panel of this figure, we plot a function f on N and on Y (which are the same as being rendered on the mesh). And then we repeat the procedures of the experiment in Section 5.3.1 and compute low-rank approximation of T F with the ground truth map T from M (resp.X) to N (resp.Y ). The transformed functions, T MESH F f and T P CD F f are plotted respectively on M and X. The left hind paw of X is supposed to be in cold color (like the left hind paw of M ), however, we observe that the tip of the paw close to the abdomen is rendered with a warmer color, meaning that T P CD F f is problematic in this area. Obviously, this defect is due to the huddled pose of the lion, which makes the points belonging to the abdomen and to the limbs hard to be differentiated in the k-NN graphs.

Beyond Shapes

We now demonstrate the potential of applying the functional-maps-based frameworks in a more abstract setting. In all above experiments, we consider point clouds sampled from closed shapes embedded in R 3 . Nevertheless, our theoretical analysis applies for any pair of manifolds without constraint on its intrinsic or embedding dimension. We now give an example showing one key advantage of our pipeline-we need no prior information about the object manifolds, and build everything from point sets.

We re-visit a classic example from [Tenenbaum 2000a], where a point cloud sampled from a Swissroll embedded in R 3 is processed. The goal there is to unfold the Swiss-roll into R 2 so that local distance are as preserved as possible. It is then demonstrated that by applying the Isomap algorithm proposed in [Tenenbaum 2000a] we obtain a significantly improved result than using other algorithm like MDS.

In our experiment, we obtain two embeddings in R By comparing the constrained norms of both operators with respect to the first 40 eigenvalues of the source PCD-X sw , we observe a clear distinction in the conformal-based case, as shown in the right part of Figure 5.11. The constrained norm of R MDS , constructed from the difference between X sw and Y mds is far larger than the other one (7x).

Furthermore, the corresponding conformal indicator highlights where distortions take place when mapping from X sw to Y mds .

The results indeed match the ground truth. First, we observe dramatic difference in conformal-based cases because if an algorithm preserves the intrinsic two-dimension structure of the Swiss-roll, then the edge connection should also be preserved. In this sense, we can conclude that the MDS algorithm is outperformed by the Isomap in this task as it connects pairs of points which are far away in the original embedding. Second, the red points marked out on the right of Figure 5.11 are actually mapped closely to their antipodal points, thus they are reasonably problematic points. The red points marked out in the isomap side, are close to the boundary of the original embedding. CHAPTER 6

Conclusion and Perspectives

In this thesis, we study two problems on geometric data analysis-metric reconstruction for filamentary structures and shape comparison, and contribute in both theory and practice. Our investigation also raises some interesting questions and open perspective for future research. In the following, we quickly list a few of them that might be of particular interest.

Metric Reconstruction for Filamentary Structures

In theory, we propose the α-Reeb graph and show that the (α-)Reeb graph of the distance function is an appropriate tool for approximating a filamentary structure, especially from the metric reconstruction point of view.

Still, there are many open problems following. For instance, with another base point, we obtain a new distance function and thus a new (α-)Reeb graph. Can we gain more information about the filamentary structure by comparing/combining (in some sense) these two graphs?

Stability Analysis on Functional-maps-based Frameworks We prove in the general case where M, N are Riemannian manifolds of dimension n that the shape difference operators and the associated multi-scale framework is stable in different senses.

Notice that the functional-maps-based frameworks take a trio (M, N, T ) as input. The missing part of perturbations on the input is then about T . In fact a perturbation analysis with respect to T is of practical importance, because in the discrete setting T is combinatorial and sensitive to noises.

Performing Functional-maps-based Frameworks on PCD As shown in the experiments, our pipeline is easy to implement, and enforces less constraints on the inputs. We perform the functional-maps-based frameworks with our pipeline on some point clouds from 3D shapes and obtain results that are relevant and robust. Undoubtedly this combination extends the range of potential applications. Though there exists deficiencies due to the extrinsic nature of our approach, it is still a nice complementary to the more delicate analysis methods.

It's worth noting that this approach is purely experimental and without any theoretical guarantee. Giving theoretical guarantees for the PCD setting is both appealing and challenging.

Titre : Deux contributions à l'analyse géométrique de données : approximation de structures filamentaires et stabilité des approches fonctionnelles pour la comparaison de formes. Mots clés : graphe de Reeb, rapprochement métrique, cartes fonctionnelles, perturbation (mathématiques). Résumé : En ce moment même, d'énormes quantités de données sont générées, collectées et analysées. Dans de nombreux cas, ces données sont échantillonnées sur des objets à la structure géométrique particuliére. De tels objets apparaissent fréquemment dans notre vie quotidienne. Utiliser ce genre de données pour inférer la structure géométrique de tels objets est souvent ardue. Cette tâche est rendue plus difficile encore si les objets sous-jacents sont abstraits ou encore de grande dimension. Dans cette thèse, nous nous intéressons à deux problèmes concernant l'analyse géométrique de données. Dans un premier temps, nous nous penchons sur l'inférence de la métrique de structures filamentaires. En supposant que ces structures sont des espaces métriques proches d'un graphe métrique nous proposons une méthode, combinant les graphes de Reeb et l'algorithme Mapper, pour approximer la structure filamentaire via un graphe de Reeb. Notre méthode peut de plus être facilement implémentée et permet de visualiser simplement le résultat. Nous nous concentrons ensuite sur le probléme de la comparaison de formes. Nous étudions un ensemble de méthodes récentes et prometteuses pour la comparaison de formes qui utilisent la notion d'applications fonctionnelles. Nos résultats théoriques montrent que ces approches sont stables et peuvent être utilisées dans un contexte plus général que la comparaison de formes comme la comparaison de variétés Riemanniennes de grande dimension. Enfin, en nous basant sur notre analyse théorique, nous proposons une généralisation des applications fonctionnelles aux nuages de points. Bien que cette généralisation ne bénéficie par des garanties théoriques, elle permet d'étendre le champ d'application des méthodes basées sur les applications fonctionnelles. Title: Two contributions to geometric data analysis: filamentary structures approximations, and stability properties of functional approaches for shape comparison. Keywords: Reeb graph, metric approximation, functional maps, perturbation analysis. Abstract: Massive amounts of data are being generated, collected and processed all the time. A considerable portion of them are sampled from objects with geometric structures. Such objects can be tangible and ubiquitous in our daily life. Inferring the geometric information from such data, however, is not always an obvious task. Moreover, it is not a rare case that the underlying objects are abstract and of high dimension, where the data inference is more challenging. This thesis studies two problems on geometric data analysis. The first one concerns metric reconstruction for filamentary structures. We in general consider a filamentary structure as a metric space being close to an underlying metric graph, which is not necessarily embedded in some Euclidean spaces. Particularly, by combining the Reeb graph and the Mapper algorithm, we propose a variant of the Reeb graph, which not only faithfully approximates the metric of the filamentary structure but also allows for efficient implementation and convenient visualization of the result. Then we focus on the problem of shape comparison. In this part, we study the stability properties of some recent and promising approaches for shape comparison, which are based on the notion of functional maps. Our results show that these approaches are stable in theory and potential for being used in more general setting such as comparing high-dimensional Riemannian manifolds. Lastly, we propose a pipeline for implementing the functional-maps-based frameworks under our stability analysis on unorganised point cloud data. Though our pipeline is experimental, it undoubtedly extends the range of applications of these frameworks.

Figure 2

 2 Figure 2.1: (a) scattering points with a filamentary structure, which is approximated by the red curves. (b) a sphere and a deformed one, the deformation in-between can be visually captured.

  .2. Theorem 3.2 [Homotopy Classification of graphs] Any connected graph G = (V, E) is homotopy equivalent to a wedge of |E| -|V | + 1 circles, where |V | (resp., |E|) is the cardinality of V (resp., E).

Figure 3 . 2 :

 32 Figure3.2: The black part of G 1 is contractible to the black point in G 2 . After adding two red edges, G 1 is no longer contractible and is homotopy equivalent to a wedge of 2 cycles as G 2 .

  Reeb graph. Define relation x ∼ y if and only if d(x) = d(y) and x, y are in the same path connected component of d -1 (d(x)). This relation is an equivalence relation. The quotient space G = X/ ∼ is called the Reeb graph of d. We denote by π : X → G the quotient map. Notice that π is continuous and as X is path connected, G is path connected. The function d induces a function d * : G → R + that satisfies d = d * • π. A relation on

Lemma 3. 2

 2 If δ is a path joining two points p, p ′ ∈ G such that d * • δ is strictly increasing then δ is a shortest path between p and p ′ and d G (p, p ′ ) = d * (p ′ )d * (p).

Proposition 3. 3

 3 Figure 3.3: Tightness of the bound in Lemma 3.3: there are 3 edges of length at most 4α and the diameter of B is equal to 20α. The range of the distances from r to the points on the red curve is [dα, d + α].

  are pairwise distinct edges of G contained in B. It follows from Lemma 3.5 that the lengths of the edges E 1 , • • • E k and of [x, v 0 ] and [v k , y] are upper bounded by 4α. As a consequence the length of γ is upper bounded by 4(k + 2)α which is itself upper bounded by 4(N E (4α)+2)α since the edges E 1 , • • • E k are pairwise distinct. It follows that d G (x, y) ≤ 4(N E (4α)+2)α.

Figure 3

 3 Figure 3.4: The path correspondence between two metric spaces X and Y .

  ) any regular open set V ⊂ V connecting two critical open sets which are equivalent to g n . We abuse the notation and also denote Ṽ (g n ) the open set of the union of the open sets in Ṽ

Figure 3 .

 3 Figure 3.7 shows different types of elements in Ṽ. We summarize the properties for the open cover Ṽ in the following corollary, which follows from Lemma 3.8, Lemma 3.9, Lemma 3.10, and Lemma 3.11. Corollary 3.1 Under the hypotheses of Theorem 3.7, the open sets in Ṽ satisfy the following properties.

Figure 3

 3 Figure 3.7: Illustration of the merging strategy: From top-left to top-right (First step): We select the unions of open sets from V 0 based on the critical points w.r.t. b(g) in G ′ and merge them respectively; From topright to bottom-left (Second step): We further merge the unions which are mutually close. Bottom-right: the nerve of the merging result N ( Ṽ).

  Second assume Ṽi ∩ Ṽj = ∅. If any path in the nerve N ( Ṽ) connecting Ṽi and Ṽj passes through some open set in Ṽ1 ∪ Ṽ2 , then we are done based on Claim (i). Now assume there is a path β in the nerve N ( Ṽ) connecting Ṽi and Ṽj only passing through open sets in Ṽ3 . Since any open set in Ṽ3 is a regular set in V, the worst scenery is that β contains no intermediate open sets. In this worst scenery, due to the shrinking operation on Ũ Ṽj , one can verify that Ũ Ṽi ∩ Ũ Ṽj = ∅. Proposition 3.6 Under the hypotheses of Theorem 3.7, N ( Ũ ) is homotopy equivalent to G ′ .

Figure 3

 3 Figure 3.8: 3 types of close-self-intersection.

  (see figure 3.9). Denote the geodesic ball of radius r centered at v by B(v, r) = {g ∈ G : d G (g, v) ≤ r}. For any two edges in e 1 , e 2 ∈ E(v), reach[(e 1 e 2 )\ w∈e 1 ∩e 2 B(w, r)] ≥ βr holds for some positive constant β ∈ (0, 1]. 4. No edges in E a are shorter than b, where b is a positive constant.

  Figure 3.9: V (v) = {w 1 , w 2 , w 3 }, E(v) = {e 1 , e 2 , e 3 , e 4 }

x 0 Figure 3

 03 Figure 3.10: Illustration of Lemma 3.7

  Proposition 3.8 Let G be a metric graph satisfying condition 3.1 and τ, b be the constants involved in the condition. Let diam(G) be the diameter of graph G, i.e., diam(G) = max g 1 ,g 2 ∈G d G (g 1 , g 2 ). |V d | is the cardinality of the original vertex set of G. d S and d G are metrics induced by graph H with parameter δ as in Definition 3.1. If δ < min{τ, b}, then we have

Figure 3 .

 3 Figure 3.11: Illustrations for the proof of Proposition 3.8.

Figure 3

 3 Figure 3.12: Illustration of case 2(b) and 2(c).

  1 ∪ e 2 by U and define B a = {g ∈ G : d G (g, v a ) ≤ r} where r = cδ 2β (similarly define B b and B c with respect to v b and v c as illustrated in Figure 3.12). (a) e 1 = e 2 .

  The time complexity of the above algorithm is dominated by the computation of the distance function in the first step, which is O(|E| + |V | log |V |). The computation of the connected components in the second step is O(|V | log |V |) based on union-find data structure. In the final step, there are at most O(|V |) number of the copies of the intervals. Based on union-find data structure, the identification can also be performed in O(|V | log |V |) time.

  Figure 3.15: Earthquake data -(a) The distance functions d on each connected components. The value increases from cold to warm colors. (b) The reconstructed α-Reeb graph.

  Figure 3.17: (a) GPS traces passing through a highway crossing in Moscow. (b) The distance function. (c) and (d)The reconstructed α-Reeb graph viewed from two perspectives.

Figure 4

 4 Figure 4.1: (a) Given shapes M, N and T a map between them, V is one of the shape difference operators formulated in [Rustamov 2013]. Intuitively, f 2 , which is supported in a region that undergoes deformation via T , is significantly distorted by V . Whereas f 1 , being supported in area-preserved region, remains the same after V acting on it. (b) We perturb one of the shapes, N , to Ñ and generate indicators with the multiscale framework of [Ovsjanikov 2013]. Two types of consistency are evidenced: horizontally, the scale k increases from 30 to 100, yet the indicators of each rows highlight nearly the same areas; vertically, at each scale, the indicators are stable with respect to the changes of the input shapes.

Figure 4 . 2 :

 42 Figure 4.2: Stability across ranges of scales: the indicators from k = 20 to 100 consistently highlight the hip of the horse, meanwhile the ones from k = 120 to 200 highlight the root of its front right leg. The corresponding quantitative measurements of the distortions of the indicators are marked to the top-left of each shapes.

M

  u(-∆ M,µ )vdµ M = M ⟨∇u, ∇v⟩ g M dµ M (4.14)

  Definition 4.6 A weighted Riemannian manifold (N, g N , µ N ) is b-close to a Riemannian manifold (N, g N , ν N ) if the following holds: µ N is obtained by perturbing ν N (the volume induced by g N ) with ρ N : dµ N = ρ N dν N . And b -1 ≤ ρ N ≤ b holds for a constant b ≥ 1.

  Proposition 4.4 (M, gM , μM ), (N, gN , μN ) are respectively (a M , b M )-close and (a N , b N )-close to (M, g M , ν M ) and (N, g N , ν N ), which are smooth Riemannian manifolds of dimension n.

Theorem 4. 2

 2 Let M, N be two smooth n-dimensional Riemannian manifolds, and T be a map from M to N . Let M be (a M , b M )-close to M and Ñ be (a N , b N )-close to N . If a M , a N , b M and b N are finite real numbers not smaller than 1, then L 2 ν (N ) = L 2 μ(N ). Moreover, if M, N and T satisfy condition 4.1, we have the following convergence:

  Remark 4.2 M, N, T satisfying condition 4.1 guarantees that ∥V ∥ C 1 and ∥V ∥ C 2 are finite. The finiteness of ∥ Ṽ ∥ C is assured by Proposition 4.4.

Figure 4 . 4 :

 44 Figure 4.4: Curves representing constrained norms with respect to different subdomains construction and input manifolds, due to Proposition 4.6, if C = λ k ≪ λ k+1 , then the difference between ∥V ∥ C and ∥V ∥ k is well-bounded.

  w), ∇T F (w)⟩ gM dμ M N ⟨∇w, ∇w⟩ gN dμ N . (4.30) Accordingly, we define Ãconf (C) = Ã(C) ∩ { N wdμ N = 0} and ∥ R∥ C = max F (w) s.t. w ∈ Ãconf (C).

  Figure 4.1), the indicators with respect to different constraints are depicted on N in Figure 4.5.

Figure 4 . 5 :

 45 Figure 4.5: We compare a bumped sphere to a sphere (the same as in Figure 4.1) and plot on the sphere the (local) maximizers different constraints (note that only the left-most one is a global maximizer). The middle three are locally optimized by taking the constraints as intersection of A(λ 37 ) and S(l), where l = 145, 325, 430. The eigenvalues involved are: λ 37 = 2.3857, λ 145 = 8.3682, λ 325 = 16.6667, and λ 430 = 20.0167.

Figure 4

 4 Figure 4.6: The X-axis indicates the index of eigenvalues/eigenfunctions, and the Y-axis represents the ratiok

Figure 4

 4 Figure 4.7: T : M → N is an area-preserving map. The conformal indicator at scale k = 50 captures and highlights the areas undergo conformal deformations.

Figure 4

 4 Figure 4.9: Indicators at scale k = 50 with different meshes. We densify the original shape (mesh-a) by adding points in the body of the horse (mesh-b) and simplify it by down-sampling the limbs (mesh-c). The corresponding distortion measurements are marked to the top-left.

  λ i = C + δ. Let b 2 be a real number satisfying b 2 2a 2 2 = δ/(Cλ 2 ) and b 2 a 2 ≥ 0.The existence of b 2 is assured by the fact that a 2 is finite (in fact |a 2 | ≤ 1). Then we set w = b 2

Figure 4

 4 Figure 4.10: Four pairs of meshed shapes are compared, σ indicates the strength of perturbations added on each of the shapes in the same row. At each row, three indicators are plotted on mesh N , which are obtained by maximizing E(w) within S(20), S(50) and S(200) respectively. Note that the human poses in the second column stand with their backs towards us, thus the highlighted areas are the hip and the right elbow.

Figure 5

 5 Figure 5.1: Left: We compare a pair of meshed shapes M and N and a pair of PCD X and Y , and generate indicators from the mesh setting (top) and the PCD setting (bottom) both at scale k s = 50; Right: the ninth to the twelfth eigenfunction of the Discrete LB operator on mesh (top) and those of the Graph Laplacian on PCD (bottom).

Figure 5

 5 Figure 5.4: Top-left: two shapes with different representations-M, N are meshes and X, Y are point clouds; six pair of landmarks are selected in correspondences (red balls). Top-right: quantitative evaluations of functional maps approximated with different conditions and in different settings. Bottom: transforming the function f depicted in the top-left panel with functional maps approximated with landmarks in the three settings: mesh, PCD with 40-NN graph and PCD with 15-NN graph.

First

  of all, we compare the area-based and the conformal indicators shown in the left part of Figure 5.5. All the indicators are computed at a uniform scale k s = 50. The area-based indicators along the left-most column are consistent, meaning that the results from the PCD with different parameters are relevant and consistent. Regarding the conformal case, the one of 15-NN graph fails to capture the right area, whereas the left two with larger k manage to highlight the right hind leg.

Figure 5

 5 Figure 5.5: Comparison of different k for constructing k-NN graphs. The left-most two columns depict area-based and conformal indicators generated with meshes and k-NN graphs with respect to k = 15, 25, 40, the scales of all the indicators are fixed as k s = 50. The right part shows part of the eigenfunctions on N (mesh) and on Y (PCD), which are visually distinguishable.

  Figure 5.6(a). The top row (Figure 5.6(b), (d)) shows the indicators generated with mesh inputs: part (b) shows area differences and part (d) shows conformal differences. From left to right, the corresponding scales are respectively k s = 20, 60, 180.

Figure 5

 5 Figure 5.6: (a) shapes in comparison; Top-row on the right: indicators generated from the mesh setting-(b) are the area-based ones and (c) are the conformal ones. Bottom-row on the right: indicators generated from the PCD setting-(c) are the area-based ones and (e) are the conformal ones. In each quadrant, the scales of the indicators are k s = 20, 60, 180, and increasing from left to right.

Figure 5

 5 Figure 5.7: PCA plots of the two shape difference operators.

T

  Figure 5.8: Top row: S is a uniformly sampled sphere and S ′ is generated by tripling sampling density on the east half of S. B ′ is generated by densifying B in the same fashion. Bottom row: The area-based indicators at scale k s = 50, the point clouds in comparison are marked below.

  Figure 5.9, all the area-based indicators are computed at scale k s = 50. The right-most two columns are indicators with respect to noisy point clouds regrading the target shape and the parameters (n p , r) we use to generate the random noises are marked below the corresponding point clouds. Compared to the second column, where the point clouds Y are unperturbed, the indicators from noisy data are fairly robust. Note that the number of added noisy points to each point cloud is at least 3000, which is not ignorable in any case. The noisy points are clearly visible in the figure.

  Figure 5.9: Robustness of results from the PCD setting with respect to noisy point clouds: X and Y are the original point clouds extracted from meshes. Y ′ i , i = 1, 2 are noisy versions of Y , which are generated with the parameters marked below. The functions plotted on Y, Y ′ 1 and Y ′ 2 are area-based indicator at the same scale k s = 50.

Figure 5

 5 Figure 5.10: The low-rank approximation of the functional from a cat shape to a lion shape is computed with 50 eigenvectors on the cat and 150 eigenvectors on the lion. The superscripts indicate in which setting T • F is produced. The function visualized on the cat is transformed to M and X with respective functional maps. Defects are evidenced on the tip of the left hind paw in X, where the values of T P CD F f are distorted as the left hind paw is close to the front part of the lion.

  2 , Y Isomap and Y MDS (see respectively the top-left and the bottom-left of Figure 5.11), from applying Isomap and MDS on the point cloud data consists of 1000 points X sw (in the middle of Figure 5.11) sampled from a Swiss-roll embedded in R 3 . With the ground truth of the unfolding, we can visualize and compare the two embeddings Y Isomap , Y MDS as being colored in the left of Figure 5.11. Now assume that we only have point-to-point correspondences between X sw and Y Isomap (resp.Y MDS ). Our goal is to compare the two embeddings in R 2 and figure out where the algorithms have problem in unfolding. Still we construct 40-NN graphs on the three point cloud data, and then get two pairs of shape difference operators by comparing X sw to Y Isomap (resp.Y MDS ).

  Figure 5.11: X sw ∈ R 3 is a point cloud sampled from a swiss roll. Y Isomap (resp.Y MDS ) is the embedding of X sw in R 2 generated with the Isomap (resp. MDS) algorithm. The three point clouds are colored with the ground truth. On the right are the conformal indicator computed by comparing the two embeddings with X sw respectively. The constrained norms suggest that there are more distortions introduced by the embedding regarding the MDS algorithm in this case, and the indicators manages to capture the problematic parts.

  which forces the shortest path connecting Ṽi and Ṽj in N ( Ṽ) must pass through some open set Ṽ ∈ Ṽ1 ∪ Ṽ2 . By Lemma 3.8, for any g

Table 3 .

 3 Reeb graph to speed up distance computations by showing reductions in computation time. Only 1: #OP (#OE, #N, #E) stands for the number of original points (original edges, nodes, edges in α-Reeb graph). The graph reconstruction time (GRT) is the total time of computing distance function and reconstructing the graph. The original (ODT), respectively approximate (ADT), distance computation time shows the total time of computing these distances using the original, respectively reconstructed, graph. All times are in seconds. The last two columns show the mean and median metric distortions. pairs of points in the same connected component are included because we obtain zero error for the pairs of vertices that are not. Statistics for the size of the reconstructed graph, error of approximate distances, and reduction in computation time are given in Table3.1.

	#OP	#OE	#N	#E	GRT ODT ADT Mean Median
	GPS traces 82541 313415 21644 21554 46.8 15.27 0.96	6.5%	5.3%
	Earthquake 1600	26996	147	137	0.32 1.12 0.01 14.1% 12.5%

see Section

3.3.1 for the definition

strictly speaking we should call it the α-Reeb graph associated to the covering I but we assume in the sequel that some covering I has been chosen and we omit it in notations

Acknowledgment

First of all, I would like to thank my adviser Frédéric Chazal for offering the opportunity to conduct the research presented in this thesis. Without his patient guidance and insightful suggestions along these years, this thesis would not have been possible. I would also like to extend my appreciation to Jian Sun and Maks Ovsjanikov for the fruitful collaborations, which respectively turn into the two major parts of the thesis.

I would like to thank Pierre Alliez and Yusu Wang, who took the time to carefully review my thesis. And I would also like to thank Tamy Boubekeur, Xavier Goaoc, Pooran Memari, Maks Ovsjanikov and Boris Thibert, who accepted to be the jury of this thesis.

I would like to express my gratitude to my colleagues in the Datashape team for the delightful working atmosphere. Especially I'd thank our team assistant Christine Biard. Her kind help makes my life as a non-French speaker much more easier than it could have been.

Last but not least, I owe my deepest gratitude to my parents for their consistent encouragement and support over the past 3 years.

in this chapter.

Differential Geometry

We start by defining a smooth n-dimensional manifold M . Roughly speaking, a n-dimensional manifold is a topological space that locally resembles an Euclidean space of dimension n. To give a precise definition, we introduce a chart on M : Definition 4.1 A n-dimensional chart is a couple (U, h) where U is an open subset of M and h : U → R n is homeomorphism between U and an open subset of R n .

We call M a C-manifold of dimension n if M is a Hausdorff topological space with countable bases and each point of M belongs to a n-dimensional chart. Furthermore, we define a smooth atlas which consists of smooth, compatible charts.

For any pair of i, j ∈ Γ satisfying U i ∩ U j ̸ = ∅, the transition map h ij is a smooth function (note that the domain and the image of h ij are both subsets of R n ).

Two smooth atlases are compatible if and only if their union is a smooth atlas. A smooth structure on M is the union of all compatible atlases. A smooth n-dimensional manifold is then a C-manifold endowed with a smooth structure. Now we consider a point x ∈ M and a chart (U, h) such that U contains x. This chart induces a local coordinate system (x 1 , x 2 , • • • , x n ), where x i is the i-th coordinate function of h(y) ∈ R n , where y ∈ U . The following notions from Riemannian geometry is introduced with a local coordinate system, however, most of them are independent of the selection of the local coordinate systems.

It is well-known that T x M is a n-dimensional linear space, and that, in chart, { ∂ ∂x i x } n i=1 form a basis of T x M . Lastly, since k is an integer, since λ k+1λ k can be large, small, or even 0, it is difficult to analyze the changes from S(k) to S(k + 1), though this is the minimal perturbation in discrete scales already.

To overcome these issues, we construct a new collection of multi-scale subdomains which evolves continuously. We notice that if

Moreover, thanks to the orthogonality of the eigenbasis {ϕ i } i≥1 , it's straightforward to verify that for any

It's then natural to consider the following multi-scale subdomains controlled by a continuous parameter C:

In fact, thanks to the Green formula 4.10, A(C) can be equivalently defined as:

From this point of view, this expression suggests that (the normalized) Dirichlet's energy of w ∈ A(C) is upper-bounded by C. In general, a small C prohibits large variations of w over a short distance with a global control of the magnitude of the gradient of w, therefore it forces w ∈ A(C) to be smooth.

Stability with respect to the Changes in Scale

We first investigate the stability with respect to the changes in scale, which only involves the original input manifolds M and N . As demonstrated in Figure 4.1 and 4.2, the results show consistency of the areas on N highlighted by the indicators across a range of scales. It is then attempting to validate the stability of the