Persistent Homology Computation Using Combinatorial Map Simplification

HAL Id: hal-02002963 https://hal.archives-ouvertes.fr/hal-02002963 Submitted on 1 Feb 2019 HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Persistent Homology Computation Using Combinatorial Map Simplification Guillaume Damiand, Rocio Gonzalez-Diaz To cite this version: Guillaume Damiand, Rocio Gonzalez-Diaz. Persistent Homology Computation Using Combinatorial Map Simplification. International Workshop on Computational Topology in Image Context, Jan 2019, Malaga, Spain. pp.26-39, �10.1007/978-3-030-10828-1_3�. �hal-02002963� Persistent Homology Computation Using Combinatorial Map Simplification Guillaume Damiand and Rocio Gonzalez-Diaz 1Univ. Lyon, CNRS, LIRIS, UMR5205, F-69622 France [email protected] https://liris.cnrs.fr/guillaume.damiand/ 2Universidad de Sevilla, Dpto. de Matem´atica Aplicada I, S-41012, Spain rogo[email protected] http://personal.us.es/rogodi/ Abstract. We propose an algorithm for persistence homology computation of orientable 2-dimensional (2D) manifolds with or without boundary (meshes) represented by 2D combinatorial maps. Having as an input a real function hon the vertices of the mesh, we first compute persistent homology of filtrations obtained by adding cells incident to each vertex of the mesh, The cells to add are controlled by both the function hand a parameter δ. The parameter δis used to control the number of cells added to each level of the filtration. Bigger δproduces less levels in the filtration and consequently more cells in each level. We then simplify each level (cluster) by merging faces of the same cluster. Our experiments demonstrate that our method allows fast computation of persistent homology of big meshes and it is persistent-homology aware in the sense that persistent homology does not change in the simplification process when fixing δ. Keywords: Persistent homology computation; 2D combinatorial map; mesh simplification 1 Introduction Topological data analysis (TDA) is a relatively new field in computer science. One of the most useful concept in TDA is the one of persistent homology which is an algebraic method for measuring topological features (connected components, voids, cavities, etc) of shapes and functions. Two of the crucial ingredients of persistence are: (1) a cell complex to structure the data; and (2) a filtration which is a nested sequence of subcomplexes that starts with the empty complex and ends with the whole complex. See [1, 2] for initial reports and [3, 4] for a modern exposition of the field. In [5], the authors proposed an efficient algorithm that computes persistent homology for 3D gray-scale images using the Morse-Smale complex previously obtained, which is much smaller than the input data, but with all necessary information. The authors first computed a combinatorial gradient vector field 2 Guillaume Damiand and Rocio Gonzalez-Diaz (GVF) by a process presented in [6]. To do this, the cell complex is decomposed into the lower star of its vertices. The authors then computed persistent homology from the boundary map of the chain complex associated to the Morse-Smale complex induced by GVF. In [7], we proposed an efficient algorithm for computing the homology of meshes (orientable manifolds with or without boundary), represented by 2D combinatorial maps (which are models of representation of meshes composed by vertices, edges linking two vertices, and 2D faces bounded by a closed path of edges), avoiding the time-consuming step of constructing and modifying boundaries and coboundaries of cells. The process consists of merging faces if they share a common edge, guaranteeing that the structure of combinatorial map and the homology information of the mesh is preserved until the end of the process. In this paper we extend our work to compute persistent homology of meshes. First, as in [7], a simplification process is made to improve computation time. Now, faces as dispatched in clusters depending on a parameter δand only faces of the same cluster are merged. For constructing the cluster the following rule is used: two faces are in a same cluster if there is a path of vertices of these two faces of length smaller than δ. At the end of the process, a smaller than the input 2D combinatorial map is obtained. To obtain persistent homology of the simplified mesh, lower-start filtration induced by a function hon its vertices (in our case, his the height function) is computed. Varying the parameter δ, the filtration varies and also its persistent homology. The paper is organized as follows. Section 2 recalls the background of the paper regarding combinatorial maps and persistent homology. Section 3 is the main section of the paper and presents our method to compute persistent homology starting from a particular filtration constructed from the height function and a parameter δ. Several experimental and computational results are presented in Section 4. Finally, we summarize the paper with a brief discussion about future work in Section 5. 2 Preliminary Notions In this section we recall the needed background of the paper regarding combinatorial maps and persistent homology. 2.1 2D combinatorial maps A 2D combinatorial map [8, 9], called 2-map, is a model of representation of a mesh, which is composed by i-cells:vertices or 0-cells associated with points, edges or 1-cells which link two vertices, and faces or 2-cells which are bounded by a closed path. Two cells are incident if one cell belongs to the boundary of the other one; while two i-cells c1and c2are adjacent if it exists one (i−1)-cell incident to both c1and c2. An edge eis dangling if it is incident to one vertex vsuch that no other edge than eis incident to v. An edge is isolated if it has no adjacent edge. Persistent Homology Computation Using Combinatorial Map Simplification 3 An edge incident to two different faces is called inner. Such an edge is necessarily not dangling nor isolated. Lastly, an edge is called border if it is incident to only one face and if it touches the boundary of the mesh. See Figure 1(a). f3 1 e 4 e 5 e f2 2 e3 e f4 v2 v1 f1 v3 6 e 7 e (a) 4 5 1 2 3 6 11 12 9 8 10 7 13 16 17 14 15 18 20 19 (b) Fig. 1. (a) Example of a mesh having 5 faces (the four faces incident to vertex v1, and the “degenerated one” bounded twice by edge e7), 14 edges (e6is dangling, e7 is isolated, {e1, e2, e3, e4}are inner and the rest are border) and 12 vertices. (b) The corresponding 2-map has 20 darts. Images taken from [7]. The different elements of a mesh are encoded in a 2-map by darts and two mappings between these dart: β1and β2: β2: A dart is an orientation of an edge. If an edge separates two faces, it is described by two darts d1,d2in the 2-map linked by β2(i.e., β2(d1) = d2and β2(d2) = d1). These two darts represent the two possible orientations of the edge (for example β2(8) = 11 and β2(11) = 8 in Figure 1(b)). Each border edge is described by only one dart din the 2-map, linked by β2with a special element ∅ (cf. for example dart 10 in Figure 1(b) which describes border edge e5).β1: For each dart d,β1(d) is the dart following dart dand belonging to the same face than d(for example β1(1) = 2 in Figure 1(b)). Note that a 2-map is oriented and thus described a given orientation of the mesh. A dart belongs exactly to one vertex, one edge and one face, and thus each cell of the mesh is described by a set of darts in the 2-map. For example, in Figure 1(b), vertex v1is described by the set of darts {2,5,8,12}. Note that this is a very important property of 2-map. Even an isolated edge (like e7in Figure 1(a)) belongs to one face (which explain why we have 5 faces and not 4 in Figure 1(a)). The different type of edges can be detected in a 2-map thanks to particular configurations of darts and βlinks (for example an edge is isolated if β1(β1(d)) = d,dbeing of of the dart of the edge). The algorithm presented in this paper for computing persistent homology on meshes used a modified version of Algorithm 1 detailed below which was presented in [7] to compute the minimal 2-map (i.e. with minimal number of 4 Guillaume Damiand and Rocio Gonzalez-Diaz cells) describing a given mesh. The algorithm uses two operations on 2-maps: edge removal and edge contraction. It simplifies a given combinatorial map in its minimal form while preserving all the homology information. The proof that Algorithm 1 preserves homology information is given in [10]. Algorithm 1: Simplification of a mesh (modified version of Algorithm 1 of [7]). Input: A 2-map Mrepresenting the mesh. Output: The simplified 2-map corresponding to M. foreach edge eof Mdo if eis an inner edge then remove e; foreach edge eof Mdo while eis dangling do e0←one edge adjacent to e; remove e;e←e0; foreach edge eof Mdo if eis not a loop then contract e; 2.2 Persistent Homology In this subsection we give elementary notions from topology needed to understand the rest of the paper. In particular, we introduce the notion of homology and persistent homology. Precise definitions of homology can be found for example in [11], and definition of persistent homology for example in [4]. Homology can be thought as a method for defining k-dimensional holes (connected components, tunnels, voids) in a given mesh. For example. a 1-cycle is a closed path and a 1-boundary is the boundary of a 2D manifold. Then, 1- homology classes (which represent tunnels) are equivalence classes of 1-cycles modulo 1-boundaries. This concept can be generalized to k-homology classes. Finally, k-homology groups are the groups of k-homology classes. Persistent homology captures the topological changes occurring in a growing sequence of meshes, called filtration. During the growth of a mesh, homology classes of different dimension may appear (be born) and disappear (die). Filtrations are frequently constructed using a real-valued function hon the vertices of the mesh M. For example, the lower-start filtration is computed as follows: –First, order the vertices of Min a non-decreasing way, h(v1)≤h(v2)≤ · · · ≤ h(vn). –Second, compute the lower-star of a vertex vin M, which is the set of cells of Mincident to vwhose vertices all have function values at most h(v). Persistent Homology Computation Using Combinatorial Map Simplification 5 –Define Mias the union of the lower-star of all vertices of Mwhose function value is at most h(vi). This way, if h(vi−1)< h(vi) then Mi\Mi−1is the set of cells of Mhaving a vertex with function value exactly h(vi). And if h(vi−1) = h(vi) then Mi−1=Mi. The lower-star filtration of the mesh Minduced by the function his the sequence of nested meshes: ∅=M0⊆M1⊆M2⊆ · · · ⊆ Mn−1⊆Mn=M. Intuitively, imagine we sweep the mesh Min increasing values of the function h. At any real-value α, we consider the set of cells whose function value on their vertices is below or equal to α. As αincreases, this gives us a sequence of subsets of M, growing larger and larger. The topological evolution along the filtration is expressed by the corresponding sequence of homology groups. When adding the cells in order according to the filtration, new homology classes may born and some of them may later die when they become trivial or merge with another class. If a homology class γis born at Miand dies entering Mjthen h(vj)−h(vi) is the persistence of γ. If γis born at Mibut never dies then its persistence is set to infinity. Homology classes with low persistence are considered noise and the ones that persist are considered features of the mesh. The information obtained when computing persistent homology can be visualized as a persistence barcode which consists of the set of (birth, death) intervals, each interval recording a persistent homology event. The bottleneck distance is used to compare two persistence barcodes corresponding to two different filtrations of the same mesh. Given a bijection ηbetween two persistence barcodes, we take the supremum L∞-distance3between matched points and define the bottleneck distance by taking the infimum over all supremums. In order to compute persistent homology, in this paper we have implemented a simplified version of the incremental algorithm for computing AT-models given in [12]. Given an ordering of the cells of the mesh, Algorithm 2 computes a triplet (M, H, f) where: –Mis the given mesh (decomposed in cells obtained from the combinatorial map). If σis a k-cell, then ∂(σ) is the set of (k−1)-cells in its boundary. –His a subset of cells of Mcalled surviving cells. Fixed k, the set of all the surviving k-cells together with the addition operation + (here + means the disjoint union of sets) form the group Ck(H) which is isomorphic to the k-dimensional homology group Hk of M. –f:C(M)→C(H) maps each k-cell in Mto a sum of surviving cells, satisfying that if a, b ∈Ck(M) are two homologous k-cycles then fk(a) = fk(b). Let Mσibe the set of cells {σ1, . . . , σi}. Then, in the ith step of 3The L∞-distance between points u= (u1, u2) and v= (v1, v2) in the extended plane is max{|u1−v1|,|u2−v2|}, 6 Guillaume Damiand and Rocio Gonzalez-Diaz Algorithm 2: Computing persistent homology (Algorithm 2 of [12]). Input: An ordering of the cells of M:{σ1,...,σm}. Output: Persistent homology. Initialize H:= ∅and f(σi) := 0, for 1 ≤i≤m. for i= 1 to mdo if f∂(σi) = 0 then f(σi) := σi,H:= H∪ {σi}(a new homology class was born). if f∂(σi)6= 0 then Let σj∈f∂(σi) s t. j= max{index(µ) : µ∈f∂(σi)} H:= H\ {σj}(an old homology class died). foreach x∈Msuch that σj∈f(x)do f(x) := f(x) + f∂(σi). Algorithm 2, σibelongs to a k-cycle cin C(Mσi) if and only if f∂(σi) = 0. This is why if f∂(σi) = 0 then a new homology class was born (the one represented by the k-cycle c) and σienters H. Otherwise, if f∂(σi)6= 0, then a homology class died, which is equivalent to say that an element of f∂(σ)⊆His removed from H. The element to be removed from Hwill be the youngest one: max{index(µ) : µ∈f∂(σi)}, being index(µ) the position of the cell µin the given ordered list of cells {σ1, . . . , σm}. In [13] the authors establish a correspondence between the incremental algorithm for computing AT-models given in [12] and the one for computing persistent homology [4]. Since we are only interested in computing the persistence events, we only compute the set Hand the map f. See Algorithm 2. 3 Computing Persistence Our starting point is a subdivision of a mesh M(with or without boundaries) into vertices, edges and faces, and a real-valued function hon the vertices of the mesh. Our method is based on three steps: 1. Simplification of the 2-map according to a parameter δ; 2. Computation of the lower-star filtration of the simplified mesh; 3. Computation of persistent homology of the given filtration. Our goal in step 1 is to simplify the 2-map decreasing the number of faces in each level of the filtration in order to improve the computation time in Step 3 which is the more time-consuming step. Observe that persistent homology varies when δvaries since the filtration computed is different. Nevertheless, we have observe in the experiments that our simplification can be seen as a filtering of small persistent homology events. Persistent Homology Computation Using Combinatorial Map Simplification 7 3.1 2-Map Simplification In this step, the 2-map is simplified by dispatching the faces into clusters and applying Algorithm 1 with constraints. First, faces are dispatched into clusters according to the parameter δ. To compute such clusters, vertices of the mesh are ordered in a non-decreasing way by their height values h(v). We assign a height value to each face with is the maximum value of the height of its vertices. Then in the first cluster we add the first face fin the ordering and all the faces “at distance” less than δ. which means that there exists a path of vertices of these two faces of length smaller than δ. For example, if δ= 0, only one face per cluster is added. If δ= 1 all the faces sharing an edge with fare added. For any δ > 1 all the faces at distance less than or equal to δto fare added to the cluster. We repeat the process with the next face provided by the ordering that was not included in any cluster. We repeat the process until all faces are in a cluster. After dispatching the faces in clusters, we apply Algorithm 1 with the following constraints: –Faces merge (i.e, the inner shared edge eis removed) only if they belong to the same cluster. –Besides, contrary to Algorithm 1, critical edges (separating faces belonging to two different clusters) are not removed here. Merging faces belonging to two different clusters could lead to loose a persistent event, and this is why we do not merge such faces. –We do not use the contraction step (last foreach in Algorithm 1). Indeed, the simplified 2-map obtained here has several faces, contrary to Algorithm 1 computed without constraints that always produces one face per connected component. For this reason, the number of possible edges to contract is here smaller and thus we have observed no gain (and even sometimes a loss) when using the contraction step comparing to not use it. 3.2 Filtration The second step in our algorithm for computing persistent homology is to compute the lower-star filtration (see Section 2.2) of the simplified mesh SM. Observed that increasing the value of δin Step 1 will decrease the different number of SMjsets (i.e., the number of levels in the filtration), which increases the average number of cells belonging to a same SMj, as illustrated in Figure 3 for the Neptune mesh and three different δvalues. Note that bigger δincreases the number of simplifications done and thus decreases the size of the simplified combinatorial map. In this case, the persistent homology computed is not the same than the one obtained by the lower-star filtration on the original mesh (they only coincides when δ= 0). Nevertheless, we have seen in our tests that the effect of the parameter δ > 0 is to remove small persistent homology events. However this new possibility gives to users a way to choose a level to analyze a given mesh, while allowing to speed-up the method. 8 Guillaume Damiand and Rocio Gonzalez-Diaz 3.3 Computation of persistent homology The last step of our method is the computation of persistent homology of the simplied mesh SM. We order the cells in SM according to the given filtration and obtain the ordered set of cells {σ1, . . . , σm}such that if i < j then there exist i0, j0such that i0≤j0,σi∈SMi0,σj∈SMj0and σjis not in the boundary of σi. We then apply Algorithm 2 to compute persistent homology. The persistence barcode is stored in a list Lwith the (birth, death) events as follows: if σ∈M`\M`−1is born and dies entering µ∈Mm\Mm−1, then store (birth, death) in Lbeing birth= h(vi`) and death= h(vim). Finally, bottleneck distance between different filtrations of the same mesh obtained from different values of δcan be computed to measure the effect of the parameter δin the persistent homology information obtained. 4 Experiments We have implemented our algorithm for persistent homology computation by using the CGAL implementation of combinatorial maps [14] and the additional layer, called linear cell complex, which additionally represents the geometry [15]. All our experiments were run on an Intel R i7-4790 CPU, 4 cores @ 3.60GHz with 32 Go RAM. All the computation times given here are averages of 10 consecutive runs of the same method. In our tests, we used the six meshes shown in Figure 2, having between 703,512 and 10,000,000 faces. All these meshes have only one connected component, except Blade which has 295 connected components because it contains many small isolated closed meshes inside the blade. In our experiment, we compared the persistent homology computation of the six meshes for the following values of δ: 0, 1, 2, 4, 8, 16, 32 and 64. For δ= 0, the persistent homology computed is the one of the lower-star filtration induced by the height function on the vertices of the original mesh. When δincreases, the number of faces in a same cluster increases also and thus the combinatorial map becomes more and more simplified. Nevertheless, persistent homology varies since the filtration varies, although differences are “small”. We can see an illustration of the effect of the δparameter on the size of the different clusters in Figure 3. The number of cells of the different simplified 2-maps for each value of δis given in Figure 4 (average values for the six meshes). The effect of δon the computation time is analyzed in Figure 5 where the six meshes shown in Figure 2 are used, and our method of persistent homology computation based on the 2-map simplification is ran by using different values of δ. Obviously, computation time decreases while δincreases, since more faces belong to the same cluster, and thus the combinatorial map becomes more and more simplified. We can see that the computation time decreases a lot even for small value of δwhich is very interesting. For example, for δ= 2, computation time is divided by 2.75 in average.