Computational Homology Applied to Discrete Objects

HAL Id: tel-01477399 https://hal.archives-ouvertes.fr/tel-01477399 Submitted on 27 Feb 2017 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. Computational Homology Applied to Discrete Objects Aldo Gonzalez-Lorenzo To cite this version: Aldo Gonzalez-Lorenzo. Computational Homology Applied to Discrete Objects. Discrete Mathematics [cs.DM]. Aix-Marseille Universite; Universidad de Sevilla, 2016. English. �tel-01477399� AIX-MARSEILLE UNIVERSITÉ École Doctorale en Mathématiques et Informatique de Marseille THÈSE DE DOCTORAT en Informatique UNIVERSIDAD DE SEVILLA Instituto de Matemáticas Universidad de Sevilla TESIS DOCTORAL en Matemáticas Thèse présentée pour obtenir le grade universitaire de docteur Memoria presentada para optar al título de Doctor Aldo GONZALEZ LORENZO Computational Homology Applied to Discrete Objects Soutenue le 24/11/2016 devant le jury : Massimo FERRI Università di Bologna Rapporteur Jacques-Olivier LACHAUD Université de Savoie Rapporteur Pascal LIENHARDT Université de Poitiers Examinateur Aniceto MURILLO Universidad de Málaga Examinateur Jean-Luc MARI Aix-Marseille Université Directeur de thèse Alexandra BAC Aix-Marseille Université Directeur de thèse Pedro REAL Universidad de Sevilla Directeur de thèse iii “The scientists of today think deeply instead of clearly. One must be sane to think clearly, but one can think deeply and be quite insane.” Nikola Tesla v Abstract Computational Homology Applied to Discrete Objects Homology theory formalizes the concept of hole in a space. For a given subset of the Euclidean space, we define a sequence of homology groups, whose ranks are considered as the number of holes of each dimension. Hence, β0, the rank of the 0-dimensional homology group, is the number of connected components, β1is the number of tunnels or handles and β2 is the number of cavities. These groups are computable when the space is described in a combinatorial way, as simplicial or cubical complexes are. Given a discrete object (a set of pixels, voxels or their analog in higher dimension) we can build a cubical complex and thus compute its homology groups. This thesis studies three approaches regarding the homology computation of discrete objects. First, we introduce the homological discrete vector field, a combinatorial structure which generalizes the discrete gradient vector field and allows us to compute the homology groups. This notion allows us to see the relation between different existing methods for computing homology. Next, we present a linear algorithm for computing the Betti numbers of a 3D cubical complex, which can be used for binary volumes. Finally, we introduce two measures (the thickness and the breadth) associated to the holes in a discrete object, which provide a topological and geometric signature more interesting than only the Betti numbers. This approach provides also some heuristics for localizing holes, obtaining minimal homology or cohomology generators, opening and closing holes. ∗ Homologie algorithmique pour les objets discrets La théorie de l’homologie formalise la notion de trou dans un espace. Pour un sous-ensemble de l’espace Euclidien, on définit une séquence de groupes d’homologie, dont leurs rangs sont interprétés comme le nombre de trous de chaque dimension. Ainsi, β0, le rang du groupe d’homologie de dimension zéro, est le nombre de composantes connexes, β1est le nombre de tunnels ou anses et β2est le nombre de cavités. Ces groupes sont calculables quand l’espace est décrit d’une façon combinatoire, comme c’est le cas pour les complexes simpliciaux ou cubiques. À partir d’un objet discret (un ensemble de pixels, voxels ou leur analogue en dimension supérieure) nous pouvons construire un complexe cubique et donc calculer ses groupes d’homologie. Cette thèse étudie trois approches relatives au calcul de l’homologie sur des objets discrets. En premier lieu, nous introduisons le champ de vecteurs discret homologique, une structure combinatoire généralisant les champs de vecteurs gradients discrets, qui permet de calculer les groupes d’homologie. Cette notion permet de voir la relation entre plusieurs méthodes existantes pour le calcul de l’homologie et révèle également des notions subtiles associées. Nous présentons ensuite un algorithme linéaire pour calculer les vi nombres de Betti dans un complexe cubique 3D, ce qui peut être utilisé pour les volumes binaires. Enfin, nous présentons deux mesures (l’épaisseur et la largeur) associées aux trous d’un objet discret, ce qui permet d’obtenir une signature topologique et géométrique plus intéressante que les simples nombres de Betti. Cette approche fournit aussi quelques heuristiques permettant de localiser les trous, d’obtenir des générateurs d’homologie ou de cohomologie minimaux, d’ouvrir et de fermer les trous. ∗ Homología computacional para los objetos discretos La teoría de la homología formaliza la noción de agujero en un espacio. Dado un subconjunto del espacio Euclídeo, se define una secuencia de grupos de homología, cuyos rangos se consideran el número de agujeros de cada dimensión. Así, β0, el rango del grupo de homología de dimensión 0, es el número de componentes conexas, β1es el número de túneles o asas y β2es el número de cavidades. Estos grupos son calculables cuando el espacio es descrito de manera combinatoria, como ocurre con los complejos simpliciales o cúbicos. También, dado un objeto discreto (un conjunto de píxeles, vóxeles o elementos de dimensión superior), podemos construir un complejo cúbico y así calcular sus grupos de homología. Esta tesis estudia tres enfoques relativos al cálculo de la homología en los objetos discretos. En primer lugar, introducimos el campo de vectores discreto homológico, una estructura combinatoria que generaliza el campo de vectores gradiente discreto y que permite calcular los grupos de homología. Este concepto permite ver la relación entre varios métodos existentes para el cálculo de la homología. Posteriormente presentamos un algoritmo lineal para calcular los números de Betti de un complejo cúbico 3D, y por tanto, de un volumen binario. Por último introducimos dos medidas (el espesor y la amplitud) asociadas a los agujeros de un objeto discreto, las cuales proporcionan una firma topológica y geométrica más interesante que simplemente los números de Betti. El cálculo de estas medidas además también aporta unas heurísticas para localizar los agujeros, obtener generadores de homología o cohomología mínimos, abrir o cerrar agujeros. vii French Extended Abstract – Résumé étendu La théorie de l’homologie formalise la notion de trou dans un espace. Pour un sous-ensemble de l’espace Euclidien, on définit une séquence de groupes d’homologie, dont leurs rangs sont interprétés comme le nombre de trous de chaque dimension. Ainsi, β0, le rang du groupe d’homologie de dimension zéro, est le nombre de composantes connexes, β1est le nombre de tunnels ou anses et β2est le nombre de cavités. Ces notions peuvent aussi être définies pour des dimensions supérieures, mais il n’y a plus d’intuition géométrique pour elles. Les groupes d’homologie sont calculables quand l’espace est décrit d’une façon combinatoire, comme c’est le cas pour les complexes simpliciaux ou cubiques. Puisqu’un objet discret (un ensemble de pixels, voxels ou leur analogue en dimension supérieure) peut être transformé en complexe cubique, nous pouvons aussi calculer ses groupes d’homologie. Cette thèse étudie trois approches relatives au calcul de l’homologie sur des objets discrets. Le champ de vecteurs discret homologique Le champ de vecteurs discret homologique (abrégé HDVF en anglais) est introduit au Chapitre 3. Le HDVF est une structure combinatoire définie sur un CW-complexe fini (notion généralisant les complexes simpliciaux et cubiques entre autres) qui induit une réduction (cf. Section 2.3.4), donc nous pouvons déduire ses nombres de Betti, un ensemble de générateurs de homologie ou cohomologie, etc. Étant donné un CW-complexe K, un HDVF est un pair d’ensembles disjoints de ses cellules X= (P, S)tels que la restriction de la matrice du bord sous ces deux ensembles (c’est-à-dire, la sous-matrice avec les colonnes correspondantes aux cellules de Set les lignes correspondantes aux cellules de P) est inversible. Nous démontrons (cf. Theorem 3.9) qu’un HDVF induit une réduction. Nous étudions ensuite comment calculer efficacement un HDVF avec sa réduction. En nous appuyant sur des formules connues du calcul matriciel, nous trouvons que la meilleure option consiste à ajouter les cellules dans le HDVF par couples (dont une cellule est ajoutée à Pet l’autre à S) et mettre à jour la réduction à chaque étape. Nous évitons ainsi toute inversion de matrice ou calcul du déterminant. Nous déduisons alors que le calcul d’un HDVF et sa réduction a une complexité O(n3). Nous introduisons après cinq opérations basiques pour transformer un HDVF. Bien que la notion de HDVF soit inspirée de la théorie discrète de Morse et que certains concepts sont purement des généralisations de cette théorie, nous remarquons que ces opérations sont nouvelles puisqu’elles sont basées sur le formalisme du HDVF. viii La section suivante est dédiée à l’étude de la relation entre le HDVF et d’autres méthodes d’homologie algorithmique. Nous déduisons que le HDVF généralise le champ de vecteur gradient discret (discrete gradient vector field) et le champ de vecteur gradient discret itéré (iterated discrete gradient vector field). Nous démontrons aussi comment le calcul de la forme normale de Smith (la méthode classique pour calculer les nombres de Betti) est équivalent au calcul d’un HDVF. Par conséquence, on peut calculer l’homologie persistante avec un HDVF en utilisant les opérations basiques. Nous finissons cette partie avec une étude expérimentale (et non théorique) de la complexité du calcul du HDVF. Bien que nous estimons sa complexité comme O(n3), nous apprécions qu’elle est en moyenne O(n2) ou même O(n1.4)si nous calculons seulement le HDVF sans sa réduction. Calcul rapide des nombres de Betti sur un complexe cubique 3D Cette partie est le fruit d’une étroite collaboration avec Mateusz Juda. Nous introduisons un algorithme efficient pour calculer les nombres de Betti sur un complexe cubique 3D. Cet algorithme est essentiellement basé sur le calcul du nombre de composantes connexes dans deux graphes, donc sa complexité est linéaire. Calculer les groupes d’homologie d’un CW-complexe de dimension quelconque requiert des techniques générales, telles que la méthode basée sur le calcul de la forme normale de Smith ou le HDVF. Par contre, si le CWcomplexe est un complexe cubique 3D, il existe une astuce pour calculer ses nombres de Betti s’il est un complexe cubique 3D. Soit Kun tel complexe, nous démontrons (cf. Proposition 4.3) que son nombre de Betti de dimension 0, β0(K), est le nombre de composantes connexes d’un certain graphe défini sur les cellules de dimension 0 et 1 de K. Nous prouvons ensuite qu’il existe une relation entre les nombres de Betti de Ket les nombres de Betti de son complémentaire dans un sur-complexe acyclique (cf. Proposition 4.4). Ceci permet de démontrer que β2(K)est le nombre de composantes connexes d’un certain graphe défini sur les cellules de dimension 2 et 3 de L−K, où Lest un sur-complexe de Kacyclique. Ces résultats formalisent l’idée intuitive que β0(K)est le nombre composantes connexes de Ket β2(K), le nombre de cavités ou composantes connexes bornées de son complémentaire. Ces propositions sont toutes démontrées en s’appuyant sur le formalisme des HDVFs. Étant donné que β0(K)et β2(K)peuvent être obtenus en comptant des composantes connexes, on peut déduire β1(K)grâce à la formule d’EulerPoincaré. Nous proposons une approche simple et itérative pour calculer ses quantités en utilisant l’algorithme classique pour compter des composantes connexes avec un parcours en profondeur. Ensuite nous introduisons un algorithme récursif avec une technique diviser pour régner qui permet de paralléliser partiellement le calcul du nombre de composantes connexes. Cette approche est spécialement conçue pour les complexes cubiques et n’est donc pas valable pour des complexes simpliciaux. Nous comparons notre implémentation avec la bibliothèque CAPD : :RedHom et nous montrons que nous obtenons de meilleurs temps d’exécution ainsi qu’elle permet de traiter des complexes plus grands grâce à ses moindres besoins de mémoire. xv Acknowledgements I would like to thank my advisors Jean-Luc Mari, Alexandra Bac and Pedro Real for their guidance, encouragement, and for providing me the freedom to work on the topics that interested me most. I am also indebted to many colleagues in the G-Mod team for their help in so many technical problems: Arnaud Polette, Joris Ravaglia, Jules Morel, Eric Remy and Romain Raffin to cite a few. I also want to thank all those colleagues that I have met—personally or virtually—during these three years with whom I have had so many interesting discussions: Erin Chambers, Frédéric Chazal, Guillaume Damiand, Paweł Dłotko, Herbert Edelsbrunner, Laurent Fuchs, Tao Ju, Lu Liu, Clément Maria, Patrick Min, Vidit Nanda, Steve Oudot, Sophie Viseur and the DGtal team. I am particularly grateful to Mateusz Juda, whose collaboration led to successful results that I could otherwise never have imagined. I thank my friends, flatmates, family and my love, Apolline, who always supported me and, in exchange, now know a little more about holes. xvii Contents Abstract v French Extended Abstract – Résumé étendu vii Spanish Extended Abstract – Resumen extendido xi Acknowledgements xv 1 Introduction 1 1.1 Overview .............................. 1 1.2 How to read this dissertation .................. 3 2 Common Background 5 2.1 General Topology ......................... 5 2.2 Algebraic Topology ........................ 7 2.2.1 Homotopy ......................... 8 2.2.2 Homology ......................... 8 2.3 Computational topology ..................... 10 2.3.1 Homotopy ......................... 11 2.3.2 Simplicial homology ................... 11 2.3.3 Cubical homology .................... 14 2.3.4 Effective Homology ................... 15 2.3.5 Persistent homology ................... 17 2.4 Digital Geometry ......................... 18 3 The Homological Discrete Vector Field 21 3.1 Introduction ............................ 21 3.2 Previous Works .......................... 22 3.3 Preliminaries ............................ 22 3.3.1 CW Complex ....................... 22 3.3.2 Homology of a CW complex . . . . . . . . . . . . . . 23 3.3.3 Homological Information . . . . . . . . . . . . . . . . 23 3.3.4 Discrete Morse Theory .................. 24 3.3.5 Some Matrix Properties ................. 26 3.4 Motivation ............................. 28 3.5 Introducing the HDVF ...................... 31 3.6 Computing a HDVF ....................... 35 3.6.1 Computing the Reduced Complex . . . . . . . . . . . 35 3.6.2 Computing also the reduction . . . . . . . . . . . . . 41 3.6.3 Some questions about the algorithm . . . . . . . . . . 41 3.6.4 Another algorithm for computing a HDVF ...... 43 3.7 Deforming a HDVF ........................ 43 3.7.1 Basic Operations ..................... 43 3.7.2 Delineating (co)homology generators ......... 46 xviii 3.7.3 Connectivity between HDVFs .............. 49 3.8 Relation with other Methods in Computational Homology . 49 3.8.1 Iterated Morse decomposition ............. 49 3.8.2 The Smith normal form ................. 50 3.8.3 Persistent homology ................... 51 3.9 Experimental Complexity .................... 51 3.10 Conclusion and future work ................... 54 4 Fast Computation of Betti Numbers on Three-Dimensional Cubical Complexes 59 4.1 Introduction ............................ 59 4.2 The Iterative Algorithm ..................... 61 4.3 The Recursive Algorithm ..................... 66 4.4 Results ............................... 68 4.5 Conclusion ............................. 69 5 Measuring Holes 71 5.1 Introduction ............................ 71 5.2 The measures ........................... 74 5.2.1 On the computation of the measures . . . . . . . . . . 77 5.2.2 The robustness of the measures ............. 78 5.3 Thickness and breadth balls ................... 80 5.4 Small generators .......................... 81 5.4.1 Homology generators .................. 93 5.4.2 Cohomology generators ................. 94 5.5 Opening or closing holes ..................... 95 5.5.1 Opening a hole ...................... 98 5.5.2 Closing a hole ....................... 100 5.5.3 We want voxels, not cubes ................ 102 5.6 Conclusion and future works .................. 105 6 Conclusion 109 6.1 General conclusion ........................ 109 6.1.1 Homological Discrete Vector Field . . . . . . . . . . . 109 6.1.2 Fast Computation of Betti Numbers on Three-Dimensional Cubical Complexes .................... 111 6.1.3 Measuring Holes ..................... 112 6.2 Other works ............................ 113 6.2.1 Cellular skeletons ..................... 113 6.2.2 Opening holes in discrete objects ............ 114 xix For David Robert Jones. 1 Chapter 1 Introduction 1.1 Overview Understanding an object can be addressed by determining its volume, its convexity, its curvature, its medial axis or any other geometric descriptor. A higher level analysis can be made through topology, which tolerates continuous deformations. This could be seen as a less interesting approach, since we could not distinguish a coffee mug from a donut, but it actually provides a more essential information of the object. Homology is a powerful tool as its formalizes the concept of hole in algebraic terms. One of the clearest concepts in topology is that of connectivity. It tells us how many disjoint parts there are in an object. Let us introduce this notion for a simple family of spaces: graphs. Consider a graph G= (V, E). We recall that a path is a sequence of incident edges [e1, e2, . . . , er], that is, each pair of consecutive edges share a vertex. Then, two vertices u, v ∈Vare said to be connected if there is a path connecting them, namely [{u, w1}, {w1, w2},··· ,{wr, v}]. Being connected is an equivalence relation, so we can define classes and the quotient set. Given a vertex u∈V, the class [u]is the set of all vertices which are connected to u, that is, the connected component containing u. Thus, the quotient set under this relation is the collection of the connected components of the graph, and its cardinal is the number of connected components. Homology theory extrapolates this concept to higher dimensions. Instead of graphs, we consider simplicial complexes (or any similar structure such as cubical complexes or CW complexes). Next, we define a homology group Hqfor each dimension q≥0. Let us assume that our simplicial complex is embedded in R3. Then its homology groups are isomorphic to free groups of the form Zβ, so their rank is β. For each dimension q≥0, the rank of Hq—called βq, the q-th Betti number—is the number of q-holes. We can interpret the 0-holes as connected components, so β0(as the cardinal of the previously defined quotient set) is the number of connected components. The 1-holes correspond to tunnels or handles and the 2-holes, to cavities or voids. It is straightforward to generalize these notions to higher dimensions, though we lose geometric intuition. For instance, a hollow square has one 1-hole (a handle) and a hollow cube has one 2-hole (a void), so a hollow four-dimensional cube contains a 3-hole, even if we cannot conceive what it is. Also, if the simplicial complex is embedded in a higher-dimensional space, like the Klein bottle, its homology groups can contain a torsion subgroup, which reveals the presence of “strange” holes. ∗ 2Chapter 1. Introduction Homology groups are computable, so given a space with a finite and combinatorial description (such as a simplicial complex), we can figure out how many holes of each dimension it has. Moreover, we can even “draw” the holes, although this presents many disadvantages. Consequently, homology allows us to compare and understand objects with regard to their topology, that is, their holes. This theory can be considered to have begun with the Euler-Poincaré characteristic during the 18th century. However, its practical applications have not been exploited until the last twenty years due to its computational complexity. There are applications in dynamical systems [89,93], material science [27,112], electromagnetism [62,40], geometric modeling [43], image understanding [2,52,96,30] and sensor networks [38]. The general idea is to use homology to analyze and understand high dimensional structures in a rigorous way. Persistent homology has revolutionized these applications, and more than likely a second revolution will come with the zigzag persistent homology. ∗ The content of this thesis is organized into three main chapters: The homological discrete vector field This research was motivated by the works of Helena Molina-Abril and Pedro Real about homological spanning forests (see [90] for a general picture), which were the starting point of this PhD thesis. From this, we define a combinatorial structure, namely the homological discrete vector field (HDVF), which encodes a reduction on a CW complex. This concept passed through different formalisms until we found a clear definition that allows a deep understanding of its nature. Roughly speaking, we defined it firstly as a discrete vector field (possibly with cycles) iteratively built and we ended up defining it as a collection of cells satisfying an algebraic condition regarding the boundary of the complex. This concept, which we prove to be equivalent to different methods in computational homology, however shows an interesting combinatorial relation between different computed homology groups for a same complex. The HDVF thus improves the homological spanning forest in different ways: it works for any dimension, it always computes (if the ground ring is a field) the homology groups without needing a later diagonalization and its clear definition allows us to prove theorems using its formalism. Nevertheless, there are still many interesting open questions. Fast computation of Betti numbers on three-dimensional cubical complexes This research, which was conducted in collaboration with Mateusz Juda, develops an algorithm for efficiently computing only the Betti numbers of a 3D cubical complex (that is, without homology generators nor a reduction). Its description is simple and follows from a constructive proof involving the HDVF framework. Also, the regular structure of the cubical complexes allows us to efficiently parallelize the algorithm. We show that this algorithm outperforms the existing software specialized in cubical homology. 1.2. How to read this dissertation 3 Measuring holes Let us introduce the problem that originated this research by an example. Consider a cubic portion of cheese with edge length of 5 cm. If we find out that there are 10 holes—that is, β2= 10—we cannot really tell if the cheese is full of holes or if there are just some small bubbles inside. A direct approach is to weigh the cheese to figure out the proportion of void in the portion, but this does not tell us if all holes have similar size or not. Moreover, this will not tell us anything about the 1-holes—there may be some tunnels going through the portion or torus-shaped holes inside the cheese. The idea we had is to gradually expand the cheese and see when the holes disappear. Moreover, because we are topologists and we always think about duality, we can also shrink or erode the cheese and see when the holes disappear, which gives us an idea of the fragility of the holes. This example gives a good intuition about the two measures that we introduce in Chapter 5. By using persistent homology and the signed distance transform of a discrete object we obtain a pair of values—the thickness and the breadth— for each hole, even though holes cannot be canonically located. These measures have many good properties: their definition is general for objects of any dimension and for any type of holes and they are stable under small perturbations on the boundary of the objects. More surprisingly, they seem useful to visualize holes, find small generators of homology or cohomology and close or open holes. Let us point out that we can do this for all the holes or just for some of them—which we can choose regarding their measures. These last results are quite visual, so we illustrate them by presenting some examples. 1.2 How to read this dissertation The common background for the three main chapters is presented in Chapter 2, while the specific preliminaries for each chapter are included in it. The three chapters are ordered chronologically, but this order also reveals an increasing interest in discrete objects. Roughly speaking, Chapter 3presents a framework for computing the homology of CW-complexes (including cubical complexes); Chapter 4shows how to count the number of holes in a 3D discrete object and Chapter 5studies further the geometry of its holes, even in higher dimensions. Consequently, the chapters are not completely independent: •Chapter 3can be read without the other two. •Chapter 4uses concepts from Chapter 3in the proofs, but they can be omitted if one is only interested in the algorithm. •Chapter 5can almost be read without referring to Chapter 3, except for Section 5.4 and 5.5. Note that each chapter contains its own conclusion. The general conclusion in Chapter 6recalls the main results of each topic and explains some future works not strictly related to the research developed herein. 2.3. Computational topology 11 2.3.1 Homotopy The fundamental group is classically computed with the Seifert–van Kampen theorem. Unfortunately, it provides a group representation (a set of letters with relations between them) of the fundamental group, which is computationally useless. The problem of telling if the fundamental group of a space is the trivial group is undecidable, as it reduces to the halt problem. Thus, there is no algorithm that extracts any information from the group representation given by the Seifert–van Kampen theorem and hence the fundamental group cannot be used (computationally) as a topological invariant. Nevertheless, let us point out that [11] succeeds in comparing different fundamental groups by extracting some computable invariants. 2.3.2 Simplicial homology The easiest way of understanding homology in computational terms is through simplicial complexes. Simplicial complexes Given a collection of points {p0, . . . , pm} ⊂ Rn, their convex hull is hp0, . . . , pmi=(m X i=0 λipi|λi≥0, m X i=0 λi= 1) We recall also that a collection of points {p0, . . . , pm} ⊂ Rnis in general position if the set of vectors {−−→ p0pi|i≥1}is linearly independent. In other words, no (m−1)-dimensional flat contains all the points. Definition 2.13. Aq-simplex is the convex hull of q+1 points in general position. Let σ=hp0, . . . , pqibe a simplex, its faces are the simplices τ=hIifor every subset I⊂ {p0, . . . , pq}. Aq-simplex σis said to be of dimension qand we denote it by σ(q) when this is not clear. Note that this notation is also used for other types of complexes. Definition 2.14. A (finite) simplicial complex Kis a collection of simplices such that (1) for every σ∈K, its faces are also contained in Kand (2) for every σ, τ ∈ K, its intersection is empty or a common face. Its dimension is the maximal dimension of its simplices. For each q≥0, we denote by Kqthe set of the qsimplices of K. We now define the (simplicial) chain complex (C, d)associated to the simplicial complex K: •For each q≥0,Cqis the free R-module generated by the q-simplices of K Cq=(X i λi·σi|λi∈R, σi∈Kq) •dqis defined over the q-simplices as the alternating sum of its (q−1)- faces dq(hx0, . . . , xqi) = q X i=0 (−1)i+1 ·hx0,...,ˆxi,...xqi 12 Chapter 2. Common Background where ˆximeans that the point xihas been removed. For the 0-simplices we define d0= 0. It is easy to check that for every q > 0,dq−1dq= 0. Thus (C, d)is a chain complex and the (simplicial) homology groups of Kare well defined. Note that Kis a topological space and thus we can define its singular homology groups. Fortunately, they are isomorphic to the simplicial homology groups. More about the homology groups Let us now look more closely at the homology groups. They are quotient spaces, where each element is a class under the equivalence relation ∀x, y ∈Cq, x ∼y⇔x−y∈im(dq+1). They are finitely generated R-modules, so there exists a generating set (a basis if it is a vector space). By the fundamental theorem of finitely generated abelian groups [37, §5.2], there are two different “normalizations” of this generating set: 1. The R-module is isomorphic to Rβq×R/λ1R×R/λ2R×. . ., where each λidivides λi+1. This is called the invariant factor decomposition. 2. The R-module is isomorphic to Rβq×R/λ1R×R/λ2R×. . ., where each λiis a power of some prime number. This is called the primary decomposition. As most of the literature about computational homology, we use the first decomposition. The number βqis called the q-th Betti number and λ1, . . . , λt are the torsion coefficients of dimension q. Let us recall that if the ambient space is R3there are no torsion coefficients. The homology groups depend on the ground ring R. Most of the works in computational homology choose R=Z2since the operations between chains are simpler, the homology groups are vector spaces and thus there are no torsion coefficients. However, the homology groups with any ground ring Rcan be deduced from the homology groups with coefficients in Zby the universal coefficient theorem [67, §3.A]. Can we see the homology groups? If our simplicial complex is in R3 then there are no torsion coefficients and we can consider the ground ring as Z2. Chains are just sets of simplices and the elements of the homology groups (which are equivalence classes since they are quotient groups) are collections of sets of simplices. Figure 2.3 shows three representatives for the same element of the homology group. This means that the difference between any two of them (actually their symmetric difference) is a boundary, that is, it belongs to im(d2). Instead of seeing all the elements in the homology groups it may be more interesting to visualize only a basis of these groups and choose a representative for each generator. Therefore, there is a double choice. Figure 2.4 shows two generators for the one-dimensional homology group. The first set of cycles gives the image that we expect from a set of generators: they correspond to the holes of the simplicial complex. We have also 2.3. Computational topology 13 FIGURE 2.3: Three cycles of C1belonging to the same class in H1(C). added the second set of cycles which is a basis for H1(C)but where the holes are not so well located. FIGURE 2.4: Two possible representations for the basis of H1(C). This provides a good picture about homology. On the one hand, we consider the Betti numbers of a simplicial complex as the number of holes of each dimension, which usually corresponds with the intuition. Note that a wire-frame cube and tetrahedron (see Figure 2.5) have five and three 1holes respectively, instead of six and four. One usually sees an extra hole, which is the sum of the other holes, but this depends on the point of view. Thus, the Betti numbers let us formally define the number of holes in a space regardless of its embedding. On the other hand, we cannot say where these holes are, as there is a large number of possible sets of generators and there is no canonical choice. FIGURE 2.5: A wire-frame cube and tetrahedron seen from two different points of view. How to compute the homology groups The simplicial homology groups are computable. Since the boundary operators are linear, they can be encoded as a sequence {Dq}n q=1 of matrices, called boundary matrices. The classical method for computing the homology groups was introduced in [99] and consists in computing the Smith normal form of the boundary matrices. 14 Chapter 2. Common Background Definition 2.15. The Smith normal form of a (not necessarily square) matrix A∈ Mn×m(R)with entries in a ring Ris the matrix N=∆ 0 0 0 ∈Mn×m(R) such that ∆is a diagonal square matrix ∆ =      α10 0 0 0α20 0 0 0 ...0 0 0 0 αr      with αidividing αi+1 for 1≤i < r and there are two invertible matrices Pand Q such that PAQ =N. The method present in [99] obtains the invariant factor decomposition of the homology groups. A variant of this method introduced in [102] obtains also a basis of the homology groups. A very clear description of this method can be found in [8]. Let us point out that computing the Smith normal form of a matrix is similar to perform a Gaussian elimination, except that every pivot must divide all the remaining entries. If Ris a field then there is no difficulty. If not, one has to perform elementary operations on the rows and columns of the matrix until an entry dividing all the others appears. Then this entry is chosen as a pivot and we make all the other entries in its row and column into zeros. Storjohann introduced in [111] an algorithm with super-cubical complexity for computing the Smith normal form of a matrix over the integers or over the integers modulo d. Let us also point out that the computation of the Smith normal form can produce huge integers [65]. 2.3.3 Cubical homology In this section we introduce the cubical complexes and their homology groups. These complexes are very similar to the simplicial complexes, except that they are built with q-dimensional squares instead of q-dimensional triangles. Let us fix our ambient space as Rn. An elementary interval is an interval of the form [k, k + 1] or a degenerate interval [k, k], where k∈Z. Definition 2.16. An elementary cube is the Cartesian product of nelementary intervals σ= [x1, x1+δ1]×···×[xn, xn+δn]xi∈Z, δi∈ {0,1} =: [x, δ]x∈Zn, δ ∈ {0,1}n Its Khalimsky coordinates is the vector σK= 2x+δ∈Zn, the sum of the intervals endpoints. The number of non-degenerate intervals in this product σ(or the number of odd entries in its Khalimsky coordinates) is the dimension of σ. An elementary cube of dimension qwill be called a q-cube. Given two elementary cubes σand τ, we say that σis a face of τif σ⊂τ. 2.3. Computational topology 15 For instance, the Khalimsky coordinates of the elementary cube σ= [1,1] ×[2,3] ×[1,2] are (2,5,3) and hence it is a 2-cube. Definition 2.17. A (finite) nD cubical complex Kis a collection of elementary cubes such that for every σ∈K, its faces are also contained in K. Its dimension is the maximal dimension of its elementary cubes. For each q≥0, we denote by Kqthe set of the q-cubes of K. Let us point out that we do not demand any condition about the intersection of different elementary cubes due to the regular structure of the cubical complex. We now define the (cubical) chain complex (C, d)associated to the cubical complex K: •For each q≥0,Cqis the free R-module generated by the q-cubes of K Cq=(X i λi·σi|λi∈R, σi∈Kq) •dqis defined over the q-cubes as the alternating sum of its (q−1)-faces along each axis dq([x, δ]) = n X i=1 (−1)o(i)·([x+δi·ei, δ −δi·ei]−[x, δ −δi·ei]) where o(i)denotes the number of ones in (δ1, . . . , δi)(or equivalently, the number of non-degenerate intervals among the ifirst elementary intervals of [x, δ]), x+δi·ei= (x1, . . . , xi+δi, . . . , xn)and δ−δi·ei= (δ1,...,0, . . . , δn). For the 0-cubes we define d0= 0. Again, it is easy to check that for every q > 0,dq−1dq= 0. Thus (C, d)is a chain complex and the (cubical) homology groups of Kare well defined. Note again that Kis a topological space whose singular homology groups are isomorphic to its cubical homology groups. 2.3.4 Effective Homology Computing the homology groups using the Smith normal form is practically impossible for large complexes due to its high complexity. A solution to reduce the amount of information to compute is the notion of reduction. It is a strong relation between two chain complexes that guarantees that they have isomorphic homology groups. This is the main tool in effective homology theory [106]. We typically reduce the initial chain complex to another one much smaller (called reduced complex). In the following we omit the subscripts whenever it is clear from the context. Definition 2.18. Areduction between two chain complexes (C, d)and (C′, d′)is a triplet of graded homomorphisms ρ= (h, f, g)such that: •hq:Cq→Cq+1 for every q≥0 •fq:Cq→C′ qis a chain map: fq−1dq=d′ qfq •gq:C′ q→Cqis also a chain map: gq−1d′ q=dqgq 16 Chapter 2. Common Background •gf = 1C−dh −hd •fg = 1C′ •hh, fh, hg = 0 A reduction is usually represented with the following diagram: For instance, consider the following chain complexes (C, d)and (C′, d′), whose chain groups are freely generated with R=Z: C0=hσ1, σ2, σ3i, C1=hσ4, σ5, σ6i, C2=hσ7i, C3= 0,··· d1= −1 0 −1 0−1 1 110 , d2= −1 1 1  C′ 0=hτ1, τ2i, C1=hτ3i, C2= 0,··· d′ 1=−1 1, d′ 2= 0 Hence (h, f, g)is a reduction, where h0=  0 0 0 0 0 0 −100 , h1=−1 0 0  f0=110 001, f1=110 g0=  0 0 1 0 0 1  , g1=  0 1 0  We have followed Sergeraert’s terminology. There are equivalent or similar definitions in the literature: contraction [45, §12], strong deformation retraction [82, §2], Eilenberg-Zilber data [63, §4] or trivialized extension [98, §2]. If we remove the last line of conditions, the resulting relation is called achain homotopy equivalence. It ensures that both chain complexes have isomorphic homology groups, where the isomorphisms are the induced maps fand gin the homology groups. However, these last conditions provide more information: they decompose the chain complex into two subcomplexes: ker(f)and im(g), where the former is acyclic (its homology groups are all trivial) and the latter is isomorphic to the reduced chain complex (C′, d′). This can be thought as a homotopical thinning of the complex, 2.3. Computational topology 17 where we remove parts from the original chain complex without modifying its homology. A reduction is perfect if d′= 0. In such case, H(C)∼ =H(C′) = C′and thus the homology groups are directly obtained. Moreover, g(C′)is a basis for H(C). Also, let x∈Cqbe a cycle. If it is a boundary then f(x) = fd(y) = d′f(y) = 0 since d′= 0. Hence, g f(x) |{z} 0 =x−dh(x)−hd(x) = x−dh(x)⇒x=dh(x) That is, xis a boundary if and only if f(x) = 0 and in that case x=d(y) for the chain y=h(x). These facts should justify the interest of having a perfect reduction. Let us point out that if the homology groups of a chain complex have a torsion subgroup then there is no perfect reduction, since a perfect reduction involves homology groups freely generated, and hence of the form Zβ. Also, a reduction can always be obtained via the Smith normal form computation as described in [8, p. 48]. This reduction is perfect if the homology groups are torsion-free. Otherwise, its reduced boundary matrices are in the Smith normal form. If ρ= (h, f, g)is a reduction from (C, d)to (C′, d′), it is easy to prove that ρ∗= (h∗, g∗, f∗)is a reduction between the cochain complexes (C, d∗)and (C′,(d′)∗). Consequently, a perfect reduction also provides a basis for the cohomology groups, namely f∗(C). 2.3.5 Persistent homology Persistent homology studies the global behavior of the homology groups of a complex that changes along time. Formally, we consider a nested sequence of simplicial complexes. Definition 2.19. Afiltration of a simplicial complex Kis a sequence of subcomplexes ∅=K0⊂K1⊂ ··· ⊂ Km=K It can also be described by a function f:K→[1, m]mapping each cell to the index of the first subcomplex containing it. Note that f(σ)< f(τ)for all σ < τ since a cell must appear after its faces. We can thus build a filtration from any function on a complex by taking the maximum of the images of all the faces of a cell (even itself). Persistent homology formalizes the idea that, in a filtration, some holes last longer (are more persistent) than others. The inclusion between the subcomplexes of the filtration, which induces a chain map between their chain groups and a homomorphism between their homology groups, plays a central role. Let (Ci, di)be the chain complex associated to the subcomplex Kiand ιi,p :Ci→Ci+pthe inclusion from Kito Ki+p. 18 Chapter 2. Common Background Definition 2.20. The p-persistent q-th homology group of Kiis Hi,p q= ker(di q)/(im(di+p q+1)∩ker(di q)) ∼ =im(ιi,p q) Intuitively, it is the set of cycles in Kithat remain non-bounding for the pfollowing steps. The rank of its free subgroup is called the p-persistent q-th Betti number of Ki. This definition has three parameters: i,pand q. If we were able to define the birth and death time of each homology class in the filtration, then we could easily deduce the p-persistent Betti numbers (see the k-triangle Lemma in [44]). Zomorodian and Carlsson introduced a correspondence in [119] that associates each filtration Fto its persistent diagram, a set of intervals P D(F) = ∪q≥0P Dq(F) = {(i, j)|0≤i≤j≤ ∞}. Unfortunately, this is only possible if the ground ring is a field. A simple algorithm is given, which amounts to compute the Smith normal form by choosing pivots in the order determined by the filtration. The set PDq(F)can be interpreted as points in the plane or as intervals in the line. The points of the form (i, j)with j < ∞correspond to q-holes that are born in Kiand die in Kj, while the points of the form (i, ∞)conform the q-holes that are present in Kand appear in Ki. If the filtration is described as a function f, we denote its persistent diagram by PD(f). Cohen-Steiner et al. proved in [19] that PD(f)is stable under small perturbations of f. 2.4 Digital Geometry Adiscrete object is a finite subset of Zn. It is also called binary image (if n= 2) or binary volume (if n= 3) in order to make the difference against a grayscale image or a color image. Its elements are called pixels when n= 2, voxels when n= 3 or points in general. We endow a discrete object with a connectivity relation. Let us recall some usual connectivity relations. Let be x= (x1, . . . , xn)∈Zn, ||x||1=|x1|+···+|xn|and ||x||∞= max {|x1|,··· ,|xn|} Thus, if x, y ∈Z2, •xand yare 4-connected if ||x−y||∞≤1and ||x−y||1≤1 •xand yare 8-connected if ||x−y||∞≤1and ||x−y||1≤2 Also, if x, y ∈Z3, •xand yare 6-connected if ||x−y||∞≤1and ||x−y||1≤1 •xand yare 18-connected if ||x−y||∞≤1and ||x−y||1≤2 •xand yare 26-connected if ||x−y||∞≤1and ||x−y||1≤3 With this notation, the number accompanying the “connected” word tells the number of points connected to a point in Zn. Note that these definitions can be extended to any dimension. The reflexive and transitive closure of this connectivity relation allows us to define connected components of a discrete object. 2.4. Digital Geometry 19 We can actually build cubical complexes from discrete objects in order to obtain higher-dimensional topological information through the homology groups. We introduce two kinds of cubical complexes associated to a discrete object regarding the 2n-connectivity and the (3n−1)-connectivity. Primal associated cubical complex Let Xbe a discrete object, we denote by Kp[X]its primal associated cubical complex. Let us give a constructive definition: for each point x= (x1,··· , xn)of Xwe add to the cubical complex the n-cube [x1, x1+ 1] ×···×[xn, xn+ 1] together with its faces. This construction can be found in [16]. Dual associated cubical complex We denote the dual associated cubical complex by Kd[X]. Let us first adapt the notion of clique to our context: a d-clique is a maximal (in the sense of inclusion) set of points of Znsuch that the intersection of their correponding n-cubes is a d-cube. First, for every point (in fact n-clique) x= (x1,··· , xn)of the discrete object, we add the 0-cube σ= [x1, x1]×···× [xn, xn]. Then, for every d-clique (d < n) in the discrete object, we add to the cubical complex a (n−d)-cube such that its vertices are the points of the d-clique. This approach was used in [84]. We can also define the dual associated cubical complex in a different fashion. Consider Kthe full nD cubical complex and, for each point x= (x1,··· , xn)not in X, remove from Kthe 0-cube σ= [x1, x1]×···×[xn, xn] and its cofaces. The resulting cubical complex coincides with Kd[X]. Figure 2.6 illustrates a binary volume and its two associated cubical complexes. FIGURE 2.6: Left: a binary volume. Center: its primal associated cubical complex. Right: its dual associated cubical complex 3.3. Preliminaries 27 Lemma 3.3. (Schur determinant formula) Let be M=A B C D  a block matrix (A∈Mn×n(Z), B ∈Mn×k(Z), C ∈Mk×n(Z)and D∈Mk×k(Z)). If Ais invertible then det(M) = det(A)·det(D−CA−1B). Lemma 3.4. (The Banachiewicz identity) Let be M=A B C D  a block matrix (A∈Mn×n(Z), B ∈Mn×k(Z), C ∈Mk×n(Z)and D∈Mk×k(Z)). If Aand D−CA−1Bare invertible, then Mis invertible and M−1=A−1+A−1B(D−CA−1B)−1CA−1−A−1B(D−CA−1B)−1 −(D−CA−1B)−1CA−1(D−CA−1B)−1 We recall that the transpose of a matrix Ais denoted A⊤. Lemma 3.5. (Sherman-Morrison formula) Let A∈Mn×n(Z)and u, v ∈Zn. If Ais invertible and 1 + v⊤Au 6= 0 then A+uv⊤−1=A−1−A−1uv⊤A−1 1 + v⊤Au . Lemma 3.6. Let A∈Mm×n(Z), we say that it is an [r, c]-matrix if each row contains at most rnon-zero entries and each column contains at most cnon-zero entries. Thus, •If A∈Mm×n(Z)is an [r, c]-matrix and B∈Mm×n(Z)is an [r′, c′]-matrix, then A+Bis an [r+r′, c +c′]-matrix and it can be computed within O(min(m·(r+r′),(c+c′)·n)) operations. •If A∈Mm×n(Z)is an [r, c]-matrix and B∈Mn×p(Z)is an [r′, c′]- matrix, then A·Bis an [rr′, cc′]-matrix and it can be computed within O(m·min(r, c′)·p)operations. Lemma 3.7. (Matrix inversion lemma) Let A∈Mn×n(Z)and u, v ∈Zn. If Ais invertible then det(A+uv⊤) = det(A)·1 + v⊤A−1u. Proof. We write M=1−v u A  Since det(M) = det(Mt), by the Schur determinant formula (see Lemma 3.3), det(M) = det(Mt) det(1) ·det(A+uv) = det(A)·det(1 + vA−1u) det(A+uv) = (1 + vA−1u)·det(A) 28 Chapter 3. The Homological Discrete Vector Field FIGURE 3.2: Left: an iterated Morse decomposition, where the red arrow belongs to the first DGVF and the purple one, to the second DGVF. Right: a (standard) DGVF inducing the same reduction. 3.4 Motivation The discrete Morse theory approach has a strong interest as it addresses the computation of homology as a purely combinatorial problem rather than an algebraic one. The associated reduction can be encoded just as a list of pairs of cells. It also provides an approximation of the Betti numbers that can sometimes be accurate (depending on the choice of the integral arrows) but that is always wrong for some well known spaces as, for instance, the Bing’s house [6] (also called house with two rooms) or the dunce hat [117]. We can increase a DGVF (and thus improve the approximation) by canceling pairs of critical cells: find two critical cells τ(q+1) and σ(q)connected by only one V-path and exchange the integral and differential arrows in this path. This can be seen as reversing the direction of the V-path. Note that, even though this transformation is expressed in combinatorial terms, computing the number of V-paths is equivalent to compute the associated reduction. Another approach for reducing the number of critical cells is to compute the Morse complex and to establish a new DGVF V′on it, which is useful when there is no unique V-path between the critical cells. This is known as iterated Morse decomposition [42]. Regarding the associated reduction, reversing the only V-path between two critical cells is equivalent to adding an integral arrow between them in the Morse complex. Figure 3.2 illustrates this. Thus, reversing a V-path can be seen as pushing an integral arrow from the Morse complex back to the original one. However, not all the integral arrows on the Morse complex are equivalent to reverse a V-path: this is the case when there are several V-paths between two critical cells. Figure 3.3 shows an example where there are three V-paths between two critical cells. However, the 1-cell is a face of the 2-cell in the associated Morse complex, so we can add an integral arrow which does not correspond to a unique V-path. The motivation for our work was to push all the integral arrows in the Morse complex back to the original one. There is a different (but equivalent) point of view which is more surprising. Finding an optimal DGVF, with the minimal number of critical cells, is an NP problem. Canceling pairs of critical cells by reversing V-paths could 3.4. Motivation 29 FIGURE 3.3: The same DGVF depicted in Figure 3.1. Some differential arrows are shown in purple. seem to be a solution to this problem, but we cannot do it in general because of the conditions in the definition of a DGVF. Thus, one could think of removing one of them: 1. It must be a matching: if there are several V-paths between two critical cells, we could think of reversing all of them. Both critical cells would disappear and no cycles would thus appear. Sadly enough, this idea does not seem to give any homological information. We cannot affirm that this approach is impossible, since extra conditions could be added, but we can show a very discouraging example at Figure 3.4. On the left, there is a DGVF with 3 V-paths between the two leftmost critical cells of dimension 1 and 2. If we reverse all of them, there is just one V-path between the two rightmost critical cells. If we cancel them, we would finish with just one critical cell of dimension 0, while the complex has β1=β2= 1. A more detailed description of this example would take too much space, and we only wish to show that this does not seem a good idea. 2. There cannot be closed V-paths: miraculously, this has been a successful idea. Only by adding one condition that we introduce in Section 3.6, we obtain a generalization of the DGVF. The methods for constructing such an object and the equations for computing its associated reduction are valid for a standard DGVF. We call this kind of DVF a homological discrete vector field (HDVF). Integrating all the integral arrows in the Hasse diagram of the original complex is not a simple challenge. Nonetheless, it has already been noted that a reduction from a CW complex can benefit from its geometric realization and, in our opinion, this is the real advantage of discrete Morse theory. Let us point out a few examples supporting this rather informal affirmation: 30 Chapter 3. The Homological Discrete Vector Field (A) (B) FIGURE 3.4: (a) A DGVF. (b) The result after reversing all three V-paths between the two leftmost critical cells. •We have some a priori information about the boundary matrices of a simplicial complex: the columns of the matrix dqhave exactly qnonzero entries. Moreover, the boundary matrix dqof a cubical complex has 2qnon-zero entries in its columns and less than or equal to 2(n−q) in its rows, where nis the dimension in which the cubical complex is embedded. •In [91, §6] a parallel method for establishing a DGVF was introduced for cubical complexes. This method seems impossible to extend to other kinds of CW complexes, so it is really based on the geometry of the complex. In terms of the reduction, it can be seen as doing a partial parallel diagonalization of the boundary matrices. Extending this approach to general chain complexes is not at all clear. •Given an n-dimensional cubical complex (that is, a cubical complex embedded in Rn), it is not difficult to set a DGVF such that the homology generators of dimension n−1lie on the boundary of the complex. We can identify the (n−1)-holes of the complex by considering its complement. Choose a (n−1)-cell on the boundary of the complex next to one of those holes, and add integral arrows starting from its boundary, covering all that part of the boundary. Repeat this step for every hole and then cancel the remaining critical cells without modifying these integral arrows. This is not easily generalizable to other classes of CW complexes. Such an idea, that we could name as “modeling” or “shaping” the homology generators makes no sense when we establish a reduction from a general chain complex. 3.5. Introducing the HDVF 31 3.5 Introducing the HDVF In the context of discrete Morse theory, we always try to set a DGVF with the maximum number of integral arrows (or equivalently, with the minimum number of critical cells) in order to obtain the best possible approximation of the Betti numbers. In the language of effective homology theory, the induced reduction greatly “reduces” the original chain complex. Given a DGVF, we can improve it by incrementing the number of integral arrows. If we find two critical cells σ < τ, such that inserting an integral arrow between them does not create a cycle, adding this integral arrow reduces by two the number of critical cells. More generally, if there is only one V-path between one cell σ′belonging to the boundary of a critical cell τand another critical cell σ, we can reverse it and add the arrow (σ′, τ). This means that the integral and differential arrows in the V-path are exchanged. This can be considered as the general method for improving a DGVF (actually, in the previous case, the V-path has length zero so there is no reversing). However, depending on the order in which we cancel the critical cells and on the CW complex itself, we can create several V-paths between the other pairs of critical cells, so that we cannot cancel them anymore. This gives an intuition on why this optimization problem is NP [83]. In order to avoid this situation, we propose to allow cycles in the DVF, provided that we create them “smartly”, so a reduction can still be defined. We cancel pairs of critical cells independently of the number of V-paths, but considering the information given by the associated reduction. This means that the reduction must be known at every step, but do not panic: finding a V-path amounts also to compute a reduction. We recall that a DVF induces a partition K=P⊔S⊔Cof a CW complex. Definition 3.1. Ahomological discrete vector field (HDVF) X= (P, S)on a CW complex Kis a partition K=P⊔S⊔Csuch that d(Sq+1)|Pqis an invertible matrix (in R) for every q≥0, where Pqand Sqdenote the restrictions of Pand Sto the q-cells and d(Sq+1)|Pqis the submatrix of the boundary matrix dq+1 consisting in the columns associated to the secondary (q+ 1)-cells and the rows associated to the primary q-cells. Note that the DVF is not explicit in the definition of the HDVF. When Xis a DGVF, there is a unique DVF inducing its partition, but this is not the case for a HDVF. For instance, Figure 3.5 depicts three different DVFs inducing the same HDVF, since the primary and secondary cells in each complex are the same. Deducing a DVF requires to find a perfect matching in a bipartite graph. The existence of this perfect matching, when the partition is a HDVF, follows from Proposition 3.8. Proposition 3.8. Let Kbe a CW complex endowed with a HDVF X= (P, S). Then there exists a discrete vector field Vthat induces the partition K=P⊔S⊔C. Proof. In this proof we do not use the fact that d(Sq+1)|Pqis invertible, but that det d(Sq+1)|Pq6= 0. Let us fix a dimension q. By the Laplace expansion formula, there is a pair of cells (σ, τ)such that hd(τ), σi 6= 0 and det d(Sq+1 \τ)|Pq\σ6= 0. Thus, the discrete vector field Vcan be found recursively. 32 Chapter 3. The Homological Discrete Vector Field FIGURE 3.5: Three different matchings inducing the same HDVF. The DVF can be computed using the Hopcroft-Karp algorithm [69] in O(m√n)time, where nand mdenote the number of vertices and edges in the Hasse diagram. It is interesting as it allows us to visualize the HDVF and its computation. Let us now present the reduction induced by a HDVF. We showed in Section 3.3.4 a reduction induced by a DGVF. Since a DGVF has no cycles, the chain (1 −dV )is nilpotent and hence the sum Pk≥0V(1 −dV )kis well defined. This does not hold for the HDVF, and therefore we must consider an appropriate reduction. Note that all the operators of a reduction are linear, so they can be represented by matrices. An appropriate choice of bases can provide nice matrices and we have found a very good one: the basis B=hPq, Sq, Cqifor every chain group Cq. In the following we omit the subscripts to facilitate readability. Theorem 3.9. Let Kbe a CW complex endowed with a HDVF X. Then Xinduces the reduction (h, f, g) : (C, d)⇒(R[C], d′), where the operators h,f,gand the reduced boundary d′are given by H 0 0 0 0 0 0 00 PS C P S h=f=F0 PS Ig=G 0 P S I d′=D C C C C C C C where H= (d(S)|P)−1 F=−d(S)|C·(d(S)|P)−1 G=−(d(S)|P)−1·d(C)|P D=d(C)|C+F·d(C)|P=d(C)|C+d(S)|C·G Proof. Let us see that these linear operators satisfy the conditions of a reduction. By developing the matrix products by blocks, we can easily check that hh = 0,fh = 0,hg = 0 and fg = 1C. The rest of the conditions precise more detail. 3.5. Introducing the HDVF 33 gf = 1C−dh −hd: By developing the matrix product, we obtain   00 0 GF 0G F0I =  I−d(S)|PH0 0 −d(S)|SH−Hd(P)|PI−Hd(S)|P−Hd(C)|P −d(S)|CH0I  All the equalities can be deduced directly from the definition of H,Fand G. The equality GF =−d(S)|SH−Hd(P)|Pis more difficult to see. Let us call X=GF +d(S)|SH+Hd(P)|P =Hd(C)|Pd(S)|CH+d(S)|SH+Hd(P)|P Then, d(S)|PXd(S)|P=d(S)|PHd(C)|Pd(S)|CHd(S)|P +d(S)|Pd(S)|SHd(S)|P+d(S)|PHd(P)|Pd(S)|P =d(C)|Pd(S)|C+d(S)|Pd(S)|S+d(P)|Pd(S)|P Then, the reader can check that d(S)|PXd(S)|P= (dd)(S)|P= 0, so X= 0. We need now some properties whose proof is direct by developing the matrix product: d′=fdg =fd   0 0 I =00Idg (3.1) f=00I·(1C−dh)(3.2) g= (1C−hd)·  0 0 I (3.3) d′f=fd. Using (3.1) and (3.2), d′f=0 0 Idgf =0 0 Id(1C−dh −hd) =00I(1C−dh)d=fd 34 Chapter 3. The Homological Discrete Vector Field gd′=dg. Symmetrically, gd′=gfd   0 0 I  = (1C−dh −hd)d  0 0 I  =d(1C−hd)  0 0 I =dg d′d′= 0. Using (3.2) and (3.3), d′d′= (fdg)(fdg) = fdg(d′f)g=f(dd)gfg = 0 We say that a HDVF is perfect if its associated reduction is perfect. The previous theorem allows us to prove the desired property that the number of critical cells approximates the Betti numbers also in the HDVF. Theorem 3.10. Let Kbe a CW complex endowed with a HDVF X. Then, for every q≥0, the number of q-critical cells is greater than its q-th Betti number. Proof. A HDVF induces a reduction to a chain complex C′with isomorphic homology groups, whose rank in each dimension qis the number of critical q-cells. This proves the theorem. Let us point out that this reduction is not directly a generalization of the reduction introduced in [91]. Though, it has a similar form if we consider the reduction ρ′= (h′, f′, g′) = (h, 1−dh −hd, ι)between (C, d)and (f′(C), d). Using the same language as discrete Morse theory, this class extending the DGVF allows us to find the correct number of critical cells in complexes which do not admit a perfect DGVF, such as the Bing’s house or the dunce hat [3]. Instead of providing the explicit (and enormous) description of each complex and its HDVF, we prefer to show illustrations and to comment the construction of those HDVFs. The cubical complex version of the Bing’s house has been created by the authors. It contains 60 0-cubes, 129 1-cubes and 70 2-cubes. The first DGVF defined on it contains 13 critical cubes (see Figure 3.6-(a)): 1 0-cube, 6 1-cubes and 6 2-cubes. Let us comment that it is not the best DGVF possible. Starting from this DGVF, and after canceling pairs of critical cells by reversing V-paths, it remains only 1 critical 0-cube, which corresponds to the Betti numbers of the complex. Obviously, these V-paths were chosen to preserve the HDVF structure. Consequently, the Morse graph contains two cycles. For the dunce hat we used a simplicial complex from [64] consisting of 8 0-simplices, 24 1-simplices and 17 2-simplices. We can set a DGVF 3.6. Computing a HDVF 35 (A) (B) FIGURE 3.6: (a) A DGVF over the Bing’s house. (b) A perfect HDVF obtained on the Bing’s house. There is only one critical 0-cell (in blue). containing 3 critical cells (see Figure 3.7-(a)): one of each dimension. After reversing one V-path between the critical cells of dimension 1 and 2, we obtain a HDVF with only 1 critical 0-cell, which is in accordance with the Betti numbers of the complex. The two cycles created in the homological DVF are shown in green. (A) (B) FIGURE 3.7: (a) A DGVF over the dunce hat with three critical cells in blue. (b) The HDVF obtained after improving the DGVF. The only critical cell is the 0-cell denoted by 1. The two cycles in the Morse graph are displayed in green. 3.6 Computing a HDVF We explain in this section how we can compute a HDVF and its reduction efficiently. We do it in terms of the partition K=P⊔S⊔Cinstead of the DVF V, but we briefly describe how to obtain the DVF. 3.6.1 Computing the Reduced Complex Our first proposition states when we can add a pair of cells to a HDVF so that the matrix d(S)|Pis still invertible. 36 Chapter 3. The Homological Discrete Vector Field Proposition 3.11. Let Kbe a CW complex endowed with a HDVF X= (P, S). Let σ(q)and τ(q+1) be two critical cells. If hd′(τ), σiis a unit then X′= (P∪ {σ}, S ∪{τ})is a HDVF. Proof. We only need to prove that the matrix d(S′)|P′is invertible, where S′=S∪{τ}and P′=P∪{σ}. This matrix has the form d(S′)|P′=d(S)|P Sτ P σwu v where u=d(S)|σ,v=d(τ)|Pand w=d(τ)|σ. We know that d(S)|Pis invertible. Let us prove that w−u(d(S)|P)−1vis also invertible. By hypothesis, hd′(τ), σi=±1. Since D=d(C)|C−d(S)|C·(d(S)|P)−1·d(C)|P then, by Lemma 3.1 (without specifying the indices), hd′(τ), σi=d(τ)|σ−d(S)|σ·(d(S)|P)−1·d(τ)|P =w−u·(d(S)|P)−1·v Consequently, by the Schur determinant formula (c.f. Lemma 3.3), det(d(S′)|P′) = det(d(S′)|P′)·det(w−u(d(S)|P)−1v) is a unit, so d(S)|Pis invertible. Once we have added two critical cells to a HDVF, we do not need to comute a new DVF inducing the expanded HDVF. Instead of this, we can deduce the corresponding DVF by inverting one of the V-paths connecting both critical cells. The following proposition proves that such V-path exists. Proposition 3.12. Let Kbe a CW complex endowed with a HDVF X. Let σ(q) and τ(q+1) be two critical cells. If hd′(τ), σiis a unit then there is a V-path between them. Proof. Let Vdenote the matrix associated with the DVF introduced in Section 3.3.4. Thus, d′(C) = d(C)|C−d(S)|C·(d(S)|P)−1·d(C)|P =d(C)|C−d(S)|C·V(P)|S·(V(P)|S)−1·(d(S)|P)−1·d(C)|P =d(C)|C−dV (P)|C·(dV (P)|P)−1·d(C)|P Hence hd′(τ), σi=d′(τ)|σ=−dV (P)|σ·(dV (P)|P)−1·d(τ)|P+d(τ)|σ. If σ < τ then it is obvious. Otherwise, hd′(τ), σi=−dV (P)|σ·(dV (P)|P)−1·d(τ)|P. 3.7. Deforming a HDVF 43 second question has a positive answer, then Algorithm 1can find a perfect HDVF, though it may not find it always. We have already seen that a CW complex whose homology groups have torsion coefficients does not admit a perfect HDVF. In addition, as a consequence of Proposition 3.17, every CW complex admits a perfect HDVF whenever Ris a field. Nevertheless, we ignore what happens when R=Z and the homology groups are torsion-free. In order to find a counterexample, we executed Algorithm 1for all the torsion-free simplicial complexes in Benedetti and Lutz’s library of triangulations [4] and we always found a perfect HDVF. Moreover, the HDVFs returned for the simplicial complexes with just one torsion coefficient per dimension (i.e., Hom_C5_K4, RP4,RP4#K3_17,RP4#11S2xS2 and RP5_24) had their reduced boundary matrix already in SNF. Hence, even if they are not perfect HDVFs, the homology groups can be directly read from them. We point out that the simplicial complex hyperbolic_dodecahedral_space presented an interesting behavior. Its 1-dimensional homology group is H1= (Z5)3. Due to its small size (718 simplices), we executed Algorithm 1500 000 times with random choices of pairs of cells and we only found 72 HDVFs whose reduced boundary matrices were in SNF. The other simplicial complex with more than one torsion coefficient, PG128_PG128P7, is much larger (13 462 simplices) and we still have not found any HDVF whose reduced boundary matrix is in SNF. 3.6.4 Another algorithm for computing a HDVF Algorithm 1consists in iteratively adding a pair of critical cells to the HDVF. Nevertheless, we can also add several pairs of cells to a HDVF at the same time. Let Xbe a HDVF and Σ = {σ1, . . . , σr}and T={τ1, . . . , τr}be two sets of critical cells of codimension 1 (that is, dim(σi) = dim(τi)−1). If the matrix d′(T)|Σis invertible in Rthen X′= (P∪Σ, S ∪T)is a HDVF. The proof is similar to that of Proposition 3.11. As a consequence, Algorithm 1is not the unique way of computing a HDVF. However, we prefer it for its simplicity and we do not study in this work the above alternative approach. Let us point out that, if we can add several pairs of cells at the same time, then the second question of the previous section is true since we can add all the pairs of cells in a HDVF at once. 3.7 Deforming a HDVF In Section 3.6 we described how the reduction changes after adding a pair of critical cells to the HDVF. This can be seen as a basic operation on a HDVF, in which two critical cells γand γ′are transformed into a primary and a secondary cell respectively. In this section we extend this idea to define five basic operations that allow us to modify a HDVF. 3.7.1 Basic Operations Let Kbe a CW complex endowed with a HDVF X= (P, S). Let σ∈P, τ∈Sand γ, γ′∈C. Thus, 44 Chapter 3. The Homological Discrete Vector Field •X′=A(X, γ, γ′) = (P∪{γ}, S ∪{γ′})is a HDVF identical to Xexcept for γ, which is a primary cell, and γ′, which is a secondary cell •X′=R(X, σ, τ) = (P\{σ}, S \{τ})is a HDVF identical to Xexcept for σand τ, which are critical cells •X′=M(X, σ, γ) = ((P\{σ})∪{γ}, S)is a HDVF identical to Xexcept for σ, which is a critical cell, and γ, which is a primary cell •X′=W(X, τ, γ) = (P, (S\{τ})∪{γ})is a HDVF identical to Xexcept for τ, which is a critical cell, and γ, which is a secondary cell •X′=MW(X, σ, τ) = ((P\{σ})∪{τ},(S\{τ})∪{σ})is a HDVF identical to Xexcept for τ, which is a primary cell, and σ, which is a secondary cell The operation Ahas been largely explained in Section 3.6 and Rconsists in removing a pair of cells from the HDVF. Mexchanges a primary cell with a critical one, while Wexchanges a secondary cell with a critical one. MW is like combining Mand Wexcept that no critical cell is needed. Let us see the conditions under which we can perform each operation. Proposition 3.19. Let Kbe a CW complex endowed with a HDVF X. Let σ∈P, τ∈Sand γ, γ′∈C. Thus, 1. A(X, γ, γ′)is a HDVF if hd′(γ′), γiis a unit 2. R(X, σ, τ)is a HDVF if hh(σ), τiis a unit 3. M(X, σ, γ)is a HDVF if hf(σ), γiis a unit 4. W(X, τ, γ)is a HDVF if hg(γ), τiis a unit 5. MW(X, σ, τ)is a HDVF if hdh(σ), τiand hhd(τ), σiare units Proof. The first statement only rephrases Proposition 3.11. For the second statement we need to prove that dq(S′ q+1)|P′ qis invertible after removing the two cells. In the following we omit the subscripts. We write d(S)|P=AB C d(S′)|P′, M =10 C d(S′)|P′ where A=d(τ)|σ,B=d(S′)|σand C=d(τ)|P′. Note that det(M) = det d(S′)|P′. Then det d(S′)|P′= det(M) = det d(S)|P+1−A−B 0 0  = det d(S)|P+1 0·1−A−B = det d(S)|P·1 + 1−A−B·H·1 0 (cf. Lemma 3.7) = det d(S)|P·1 + 10·H·1 0−AB·H·1 0 = det d(S)|P·(1 + H11 −1) = det d(S)|P·H11 3.7. Deforming a HDVF 45 where H11 denotes h(σ)|τ=hh(σ), τi. The third statement is also proved using Lemma 3.7. We write dq(S)|P=a M, dq(S)|P′=b M where a=d(S)|σand b=d(S)|γ. We note that F=−b N·H where N=d(S)|C\γand thus hf(σ), γi=−b·h where h=h(σ)|S. Then det dq(S)|P′= det dq(S)|P+1 0·(b−a) = det dq(S)|P′·1 + (b−a)·H·1 0 = det dq(S)|P′·(1 + (b−a)·h) = det dq(S)|P′·(1 + b·h−a·h) = det dq(S)|P′·(1 −hf(σ), γi−1) =−det dq(S)|P′·hf(σ), γi dq(S)|P′is thus invertible. We omit the proof of the last two statements since they are similar to the third one. All these operations can be applied in terms of the DVF by reversing a V-path between the two cells considered. This V-path is not unique, but its existence can be proved using the same argument present in Proposition 3.12. In the case of MW, there are two V-paths to reverse. Figure 3.8 illustrates the operations M,Wand MW on a cubical complex. Some of these operations were introduced in [91]. Namely, the arrow reversing is the operation Mbetween 0-cells, the edge rotation is MW between 1-cells and the face rotation is MW between 2-cells. These three operations were announced as local deformations of a DGVF but, since no condition was given, this is in general false: applying an edge rotation or a face rotation to a DGVF can produce a non-gradient discrete vector field. Only the arrow reversing preserved the structure of DGVF since, in dimension 0, the existence of a V-path between a primary cell σand a critical cell γimplies that hf(σ), γiis a unit. 46 Chapter 3. The Homological Discrete Vector Field M WMW FIGURE 3.8: A HDVF on a cubical complex and the result after applying M,Wand MW. Blue cells are those which are exchanged. 3.7.2 Delineating (co)homology generators Given a perfect HDVF, the operations M,Wand MW are interesting since they change the reduction and thus the generators of the homology and the cohomology groups. The next proposition specifies how a generator changes when the operations Mor Waffect its associated critical cell. Proposition 3.20. Let Kbe a CW complex endowed with a perfect HDVF X= (P, S)with R=Z2. Let σ∈P,τ∈Sand γ∈C. Then, 1. If hf(σ), γiis a unit, the cohomology generators associated to γin Xand σ in M(X, σ, γ)are the same. 2. If hg(γ), τiis a unit, the homology generators associated to γin Xand τin W(X, σ, γ)are the same. Proof. The proof of these statements is quite lengthy, but it provides a partial description of the reduction after perturbing the HDVF. For the first statement we write Hq=u A−1 =: H1H2 Fq=−v B·Hq=−vH1vH2 BH1BH2=: F11 F12 F21 F22  where u=d(S)|σ,A=d(S)|P\σ,v=d(S)|τand B=d(S)|C\γ 3.7. Deforming a HDVF 47 Then, H′=d(S)|P′−1 =d(S)P+1 0·(v−u)−1 (cf. Lemma 3.5) =H+F−1 11 H·1 0·(v−u)·H =H·I−F−1 11 F11 F12 0 0 −F−1 11 10 0 0  =H·−F−1 11 −F−1 11 F12 0I =−H1F−1 11 H2−H1F−1 11 F12  Consequently, F′=−u B·−H1F−1 11 H2−H1F−1 11 F12  =F−1 11 F−1 11 F12 F21F−1 11 F22 −F21F−1 11 F12  The proof of the second statement is similar. We write Hq=uA−1=: H1 H2 Gq+1 =−Hq·vB=−H1vH1B H2v H2B=: G11 G12 G21 G22  where u=d(τ)|P,A=d(S\τ)|P,v=d(γ)|Pand B=d(C\γ)|P. Then it is easy to prove that H′=−G−1 11 H1 H2−G−1 11 G21H1 G′=G−1 11 G−1 11 G12 G−1 11 G21 G22 −G−1 11 G21G12  We can thus use these operations to change the shape of the generators. Figure 3.9 shows a cubical complex endowed with a HDVF. We want to have a one-dimensional homology generator around the hole. For doing this, it suffices that all the 1-cells are secondary except for one which is critical. Thus, we use Mon the top 1-cell to put there the critical cell. Then, for the other three 1-cells, we use MW to make them secondary. At the end, the homology generator induced by the HDVF stands at the desired location. It is unclear whether this application is computationally feasible. The problem is: given a perfect HDVF Xand a set of cycles S, can we find a perfect HDVF X′whose homology generators contain this set? We may 48 Chapter 3. The Homological Discrete Vector Field FIGURE 3.9: Example of multiple applications of the operations on a cubical complex. Blue cells are those that have changed. The one-dimensional homology generator is depicted in green at the beginning and at the end. first check that the cycles are linearly independent. This can be done by computing the rank of the matrix f(S). If the rank is maximal, then the cycles can be part of a homology basis. But even if the cycles are linearly independent, the HDVF X′does not exist in general, and Figure 3.10 provides a counterexample. Thus, this problem must be studied further in order to find conditions under which such HDVF exists. A possible hint to follow is that every cycle must have a cell not included in any other cycle, which is the intuition that led to our counterexample. Assuming that the HDVF X′exists, it is possible to find a sequence of operations that transform one HDVF into the other: it suffices to successively apply Rto Xfor removing all the pairing in the DVF, and then build the other HDVF using A(this is guaranteed if Ris a field by Proposition 3.18). Thus, the interesting question is to find a minimal sequence of operations that transform Xinto X′. 3.8. Relation with other Methods in Computational Homology 49 FIGURE 3.10: There exists no HDVF on this simplicial complex whose homology generators are the four triangles consisting of three 1-simplices. 3.7.3 Connectivity between HDVFs The new definitions let us state that Algorithm 1can compute any HDVF which is the result of applying only the operation Ato an empty HDVF. We explained in Section 3.6.3 that we ignore if every HDVF on a simplicial or cubical complex can be found through Algorithm 1. Thus it is natural to wonder if any HDVF on a simplicial or cubical complex can be obtained by a sequence of operations (not only A) on an empty HDVF. This can also be formulated as follows: are all the HDVF on a simplicial or cubical complex connected via a sequence of operations? This is still an open question. 3.8 Relation with other Methods in Computational Homology There are several methods for computing homology in the literature which seem to be equivalent. The simple formulation of the HDVFs allows us to clearly see these equivalences. 3.8.1 Iterated Morse decomposition First, let us prove that the HDVF generalizes the notion of DGVF. Proposition 3.21. Every DGVF is a HDVF. Proof. We need to prove that the matrices d(Sq+1)|Pqare invertible for each q≥0. In the following we omit the subscripts since the proof is the same for every dimension q. Let V={(σi, τi)}m i=1 be a DGVF. Consider the weighted digraph whose vertices are the primary cells and where arrows connect two vertices whenever there is a V-path of length 1 between them. Formally, (σi, σj)is an arrow if hd(τi), σji 6= 0 and σi6=σj. Its weight is the value −hd(τi), σii · hd(τi), σji. It is immediate to see that the matrix associated to this graph is I−d(S)|PV, where S={τi}m i=1,P={σi}m i=1 and the diagonal matrix V= (vi,j)is such that for every i,vi,i =hd(τi), σiiand zero elsewhere. Note that Vis invertible. Since there are no closed V-paths in the DGVF, the matrix I−d(S)|PVis nilpotent, and thus d(S)|PVis invertible. As Vis invertible, we deduce that d(S)|Pis invertible. A DGVF is a limited tool for computing homology. A more elaborate tool is the iterated Morse decomposition [42], which consists in iteratively: 50 Chapter 3. The Homological Discrete Vector Field (1) computing a DGVF and (2) considering the resulting Morse complex for a next DGVF. We prove now that every iterated Morse decomposition is also a HDVF. Proposition 3.22. Every iterated Morse decomposition is a HDVF. Proof. For clarity, we assume that the iterated Morse decomposition consists of only two DGVFs V1and V2. We recall that the second DGVF is defined on the chain complex consisting of the critical cells of V1and the boundary operator d′=d(C1)|C1−d(S1)|C1·H·d(C1)|P1. If we write d(S1∪S2)|P1∪P2=AB C D  then d(S1)|P1=Aand d′(S2)|P2=D−CA−1B Given the previous proposition, these two matrices are invertible. Thus, by Lemma 3.3, det d(S1∪S2)|P1∪P2= det d(S1)|P1·det d′(S2)|P2 which is a unit. It is easy to see that if a HDVF has been created using Algorithm 1, then the list of pairs of cells [(σi, τi)]m i=1 is an iterated Morse decomposition. As we showed in Section 3.6, it is not true in general that every HDVF can be computed with Algorithm 1, so we cannot deduce that every HDVF is an iterated Morse decomposition. However, this does not mean that it is false. This question remains open. 3.8.2 The Smith normal form The classic algorithm for computing homology groups computes the Smith normal form (SNF) [99]. We prove in this section that the reduced boundary matrices obtained in Algorithm 1are similar to the diagonalization performed for the computation of the SNF. Let Kbe a CW complex and Xa trivial HDVF (P=S=∅). Let us choose some pivot in some boundary matrix Dq. For simplicity, we assume that the pivot is the element Dq(1,1) = hdq(τ), σiand we omit the subscript. In order to make all the other entries in its row and column into zeros we perform ∀j6= 1, D(·, j)←D(·, j)−D(1, j)D(1,1)−1D(·,1) ∀i6= 1, D(i, ·)←D(i, ·)−D(i, 1)D(1,1)−1D(i, ·) 3.9. Experimental Complexity 51 Using the notation of Proposition 3.13, this is equivalent to D′=D−D11 D21 D−1 11 0D12  D′′ =D′−0 D′ 21 D−1 11 D′ 11 D′ 12  By developing both equations we obtain that the pseudo-diagonalized boundary matrix is D11 0 0D22 −D21D−1 11 D12 , where the bottom-right block is the reduced boundary computed in Algorithm 1after inserting the pair of cells (σ, τ). Proposition 3.23. Let Kbe a CW complex. Then, Algorithm 1performs a partial diagonalization of the boundary matrices of K. Proof. The proof is direct from the previous argument. We have just seen that computing a HDVF is equivalent to compute the SNF of the boundary matrices using only the pivot operation, that is, given an invertible entry in the matrix, we make all the entries in its row and column into zeros. Computing the SNF needs also another type of operation: if there is no entry dividing all the others, we make elementary operations on the rows and columns so such an entry appears. For this reason, Algorithm 1cannot always return a perfect HDVF if R=Z, since in the computation of the SNF we can arrive to a matrix without units even if the SNF contains only units in its diagonal (see Section 3.6.3 for an example). 3.8.3 Persistent homology Proposition 3.23 also implies that persistent homology can be computed with a variation of Algorithm 1. The classical algorithm for persistent homology [119] is based on the Smith normal form. The main difference with a standard homology computation is that cells are considered in the order given by the filtration. Therefore, Algorithm 2computes the persistence intervals of a filtration using the same calculations as Algorithm 1. Algorithm 2is a mere translation of the algorithm described in [119] into the HDVF framework. The purpose of doing so is to show that we can obtain a reduction for every step of the filtration and that we can apply the conclusions of Section 3.9 to the persistent homology theory. 3.9 Experimental Complexity We fix in this section R=Z2, so we are sure that we obtain a perfect HDVF and thus we compute the homology of the CW complex. Computing the homology groups of a CW complex is considered in general a problem with O(n3)time complexity. Only [88] proves that it can be computed in matrix multiplication time, but there is no implementation of 52 Chapter 3. The Homological Discrete Vector Field Algorithm 2: Compute a HDVF associated to a filtration Input: A CW complex Kand a filtration F=σjn j=1 Output: The persistence intervals of F for k= 0 to dim(K)do Lk← ∅; X←(∅,∅); for j= 1 to ndo if d′(σj)6= 0 then i←max nj′:hd′(σj), σj′i= 1o; X←A(X, σi, σj)(and update the boundary matrices D); Ldim(σj)←Ldim(σj)∪(deg σi,deg σj); for j= 0 to ndo if σjis critical then Ldim(σj)←Ldim(σj)∪(deg σj,∞); this algorithm. Nevertheless, it has been noticed that in practice the execution time is linear for homology [39, §4] and persistent homology [119, §4]. We estimated in Theorem 3.15 the complexity of our algorithm by bounding the number on non-zero entries in rows and columns by n, obtaining that Algorithm 1can find a HDVF within (n/2) ·n2=n3/2operations or (n/2) ·n2+n2+n2+n2= 2n3if we also want to obtain the associated reduction. Since these bounds are not tight, it should not be surprising that the complexity in practice is lower than cubic. One advantage of the HDVF framework is that we can easily count the number of operations that we perform along its computation. At each step of Algorithm 1, updating the matrices H, F, G and Drequires |F11||G11|, |D21||F11|,|G11||D12|and |D21||D12|operations respectively, where |v|denotes the number of non-zero entries in the vector v. Thus, updating the reduced boundary requires |D21||D12| operations (plus some operations to remove rows and columns). Moreover, updating all the reduction requires |F11||G11|+|D21||F11|+|G11||D12|+|D21||D12|= (|F11|+|D12|)(|G11|+|D21|) operations. Let us study now the average complexity for two random models. Random cubical complexes We introduce a random model for constructing cubical complexes. We denote it by K(p, m)and it is similar to the closed faces model introduced in [116]. Let m∈Z+and p∈R,0≤p≤1. A cubical complex in K(p, m)is built by adding each cubical cell σ∈[0, m]3(with its faces) to the complex with probability p. Note that each cell σbelongs to the cubical complex with probability 1−(1 −p)c, where cdenotes the number of cofaces (in the full cubical complex [0, m]3) of σ, including itself. Thus, lower-dimensional cells are more frequent. 59 Chapter 4 Fast Computation of Betti Numbers on Three-Dimensional Cubical Complexes THIS chapter is based on the conference paper [61], which was co-written with Mateusz Juda. We explain how to efficiently compute (only) the Betti numbers of a 3D cubical complex without manipulating any boundary matrix. 4.1 Introduction Computing homology usually needs algebraic methods. It seems that they all are based on the Smith normal form as shown in Section 3.8. However, there are Betti numbers that are easier than others. Consider a simple shape in the real plane R2. It is well known that β0 is the number of connected components and β1is the number of bounded connected components of the complement. Thus, if we can count the connected components and detect which ones are bounded, we can obtain both Betti numbers without resorting any algebraic method. Let us mention two simple scenarios where this is possible: 1. A subcomplex of a simplicial complex triangulating a rectangle. We can compute the number of connected components with the usual algorithm on the connectivity graph of the subcomplex and its complement. The unbounded connected components of the complement correspond to those connected components which contain a simplex from the boundary of the rectangle. Figure 4.1 illustrates this. Two simplicial complex K⊂Lare shown in Figure 4.1a, while the connectivity graphs of Kand L−Kare depicted in Figure 4.1b. Note that L−Kis not a simplicial complex since some simplices in L−Kdo not have their faces in L−K. There are two connected components in the connectivity graph of K(in red), so β0= 2. On the other hand, there are four connected components in the connectivity graph of L−K(in green), one of them containing all the simplices in the boundary of the square, so β1= 4 −1. 60 Chapter 4. Fast Comp. of Betti Numbers on 3D Cubical Complexes (A) (B) FIGURE 4.1: Left: a simplicial complex L(in gray) and a subcomplex K(in blue). Right: the connectivity graphs of Kand L−K. 2. A 2D cubical complex. This case is particularly simple because a cubical complex is always a subcomplex of a bigger cubical complex. The same ideas apply to this setting. We now focus only on the second scenario, since the first one is very restrictive. Can we adapt this idea to the three-dimensional space? Let Kbe a 3D cubical complex. Then β0is still the number of connected components, but now β2is the number of bounded connected components of the complement. Alexander duality generalizes this fact to any dimension and any Betti number: Theorem 4.1 (Alexander duality).Let Kbe an nD cubical complex. Then Hq(K)and Hn−1−q(Sn−K)are isomorphic for reduced homology and cohomology. This is why homology starts being interesting starting from three dimensions. The first homology group H1describes the handles and tunnels of an object, which cannot be easily deduced as connected components or voids. Nevertheless, since β0and β2are easy to compute, we can compute β1via the simplest topological invariant: the Euler-Poincaré characteristic. The Euler-Poincaré characteristic of a 3D cubical complex Kis the alternating sum of its cubes. Formally, χ(K) = k0−k1+k2−k3, where kqdenotes the number of cubes of dimension qin K. Theorem 4.2 (Euler-Poincaré formula).Let Kbe a 3D cubical complex. Then χ(K) = β0(K)−β1(K) + β2(K). Therefore, β1(K) = β0(K) + β2(K)−χ(K). 4.2. The Iterative Algorithm 61 Consequently, we can obtain the Betti numbers of a 3D cubical complex only by counting connected components and computing the Euler-Poincaré characteristic. Delfinado and Edelsbrunner introduced in [29] an algorithm with almost linear time complexity that computes the Betti numbers of a filtered simplicial complex which is a subcomplex of a triangulation of S3. They also sketched the cases where no filtration is given or where the simplicial complex is embedded just in R3. This last algorithm was developed further by Dey and Guha in [32]. Its distinctive feature is that β2is found by recognizing closed surfaces on the boundary of the complex using a graph approach, so the simplicial complex does not need to be a subcomplex of a triangulation of S3. The algorithm is first defined for three-manifolds, and if the simplicial complex is not a three-manifold, a technique for converting it is described. Sadly, there is no available implementation of this brilliant algorithm. Our work shares several ideas with these articles, though we focus our research on cubical complexes and exploit their structure, which provides a simpler algorithm. Juda and Mrozek presented in [75] an optimal algorithm which computes Z2Betti numbers and homology generators of a special class of pseudomanifolds. This is an extension of Delfinado and Edelsbrunner’s work for cubical and simplicial complexes. Our work shares several ideas with these articles, though we focus our research on cubical complexes and exploit their structure, which provides a simpler algorithm. 4.2 The Iterative Algorithm In the following, Kdenotes a 3D cubical complex. We first prove that β0(K)can be computed by counting connected components. This is a well known fact, but we include this proof to introduce this kind of reasoning which we also use in Proposition 4.5. Proposition 4.3. Let Kbe a 3D cubical complex and G0(K) = (V, E)denote the graph such that •V=K0, the 0-cells of K •E={{u, v} | uand vhave a common coface}. Thus, β0(K)is the number of connected components in the graph G0(K). Proof. Observe that Ecorresponds to the 1-cells of K, which connect its two faces. Let Fbe a spanning forest of G0(K). Assume that F={T1, . . . , Tr} where each Tiis a connected component. For each connected component Ti of F, choose a 0-cell rias root and put arrows on the other 0-cells pointing to the coface following the tree towards its root. This provides a discrete vector field Von Kwhere the only critical 0-cells are r1, . . . , rt. We now prove that it is a DGVF and that the number of critical 0-cells is minimal, so the statement of the proposition follows. As the arrows of Vare induced by trees, it is clear that there are no closed V-paths. Hence, Vis a DGVF and thus a HDVF (see Proposition 3.21 ). 62 Chapter 4. Fast Comp. of Betti Numbers on 3D Cubical Complexes The number of critical 0-cells is minimal if we cannot cancel a critical 0-cell with a critical 1-cell. This is true if d′(γ) = 0 for each γ∈C1, which we prove now. Following Theorem 3.9, we recall that d′=fd and F=−d(S)|CH =−d(S)|C  H+I−d(S)−1 |P |{z } H d(S)|P   =−d(S)|CI+H(I−d(S)|P)=−d(S)|C+F(I−d(S)|P) Thus, if (σ, τ)∈ V,f(σ) = −d|C(τ) + f(σ−d|P(τ)). By induction, fof a primary 0-cell gives the root riof its connected component. Consider any γ∈C1. Since Fis a spanning forest, the two faces of γmust be contained in the same connected component of G0(since adding this edge to Fmust create a cycle). Thus, d′(γ) = fd(γ) = ri−ri= 0 Figure 4.2 illustrates the construction done in this proof. On the left there is a 2D cubical complex together with its graph G0(K). In the middle there is a spanning forest of G0(K)(in red). On the right we can appreciate the induced DGVF after choosing as root the top leftmost 0-cell of the complex. b b bb b b b b (A) b b bb b b b b (B) (C) FIGURE 4.2: Illustration of the construction done in the proof of Proposition 4.3 We now present a proposition similar to Alexander duality which tells how the homology of a 3D cubical complex and its complement in an acyclic supercomplex are linked. Let Lbe a 3D cubical complex such that K⊂Land β(L) = (1,0,0,0). We typically consider L= [0, m]3for some m > 0, assuming that the coordinates of the cells of Kare all positive. Proposition 4.4. Let Kand Lbe two 3D cubical complexes such that K⊂Land β(L) = (1,0,0,0). Then, βq(K) = (β1(L−K) + 1 if q= 0 βq+1(L−K)else 4.2. The Iterative Algorithm 63 Proof. Note that, since K⊂L, the boundary matrix (for all dimensions together) of Lis of the form D=D1· 0D2 where D1=d(K)|Kand D2=d(L−K)|L−K. As D·D= 0 and D1·D1= 0, D2·D2= 0 too. Let X1= (P1, S1)and X2= (P2, S2)be two perfect HDVFs for Kand L−Krespectively. Let X12 = (P12, S12) := (P1∪P2, S1∪S2)be their union. It is a HDVF for L: det(d(S12)|P12 ) = det A· 0B = det(A)·det(B)∈R∗ Moreover, as Aand Bare invertible, d(S12)−1 |P12 =A−1· 0B−1 Let Y= (P, S)be a perfect HDVF for Lthat extends X12. Thus d(S)|P=    A·Y1· 0B0Y2 X1· · · 0X20·     where the last two columns (resp. rows) are T1and T2(resp. Σ1and Σ2), the new secondary (resp. primary) cells belonging to Kand L−Krespectively. By the Schur determinant formula, det(d(S)|P) = det A· 0B· det  ·· 0·−X1· 0X2A−1· 0B−1Y1· 0Y2 so the rightmost determinant must be a unit. However, developing the equation we obtain ·· 0·−X1H1Y1· 0X2H2Y2=0· 0 0  since X1and X2are perfect HDVFs. Thus, Σ2and T1must be empty. This means that all the new couples of cells in Yare from Kto L−K. Let q > 0. Since βq(L) = 0 and βq+1(L) = 0, all the critical q-cells of K must cancel with the critical (q+ 1)-cells of L−Kand vice versa. Hence, βq(K) = βq+1(L−K). For q= 0, since β0(L) = 1 and β1(L) = 0, all the critical 1-cells of L−Kmust cancel with all the critical 0-cells of Kbut one and vice versa. Therefore, β0(K) = β1(L−K) + 1. Thus, in order to compute β2(K)we can compute β3(L−K). The following proposition tells that this can be achieved also via counting connected 64 Chapter 4. Fast Comp. of Betti Numbers on 3D Cubical Complexes components Proposition 4.5. Let K⊂Lbe two 3D cubical complexes. Consider the graph G3(L−K) = (V, E)such that •V= (L−K)3∪{ǫ}, the 3-cells of L−Kplus an extra (abstract) vertex •E={{u, v} | uand vhave a common face}∪{{u, ǫ} | ucontains a free face}. Thus, β3(L−K)is the number of connected components in the graph G3(L−K) minus one. Proof. Observe that Ecorresponds to the 2-cells of L−K, which can have two cofaces (type {u, v}) or just one (type {u, ǫ}). Let Fbe a spanning forest of G3(L−K). Assume that F={T1, . . . , Tr} where each Tiis a connected component and ǫ∈T1. Fix ǫas the root of T1 and, for the other connected components Tiof F, choose a 3-cell rias root. Put arrows on the 2-cells pointing to the coface in the opposite direction towards its root. This provides a discrete vector field Von Kwhere the only critical 3-cells are r2, . . . , rt. We now prove that it is a DGVF and that the number of critical 3-cells is minimal. As the arrows of Vare induced by trees, it is clear that there are no closed V-paths. Hence, Vis a DGVF and thus a HDVF. The number of critical 3-cells is minimal if we cannot cancel a critical 3-cell with a critical 2-cell. This is true if d′(γ) = 0 for each γ∈C3, or equivalently, if (d′)∗(γ) = 0 for each γ∈C2, where (d′)∗denotes the dual of d′. We recall that d′=dg, so (d′)∗=g∗d∗and G=−Hd(C)|P =−  H+I−d(S)|Pd(S)−1 |P |{z } H   d(C)|P =−I+ (I−d(S)|P)Hd(C)|P=−d(C)|P+ (I−d(S)|P)G Thus, if (σ, τ)∈V,g∗(τ) = d∗ |C(σ) + g∗(τ−d∗ |P(σ)). By induction, g∗of a secondary 2-cell gives the root riof its connected component (if i6= 1) or the empty chain (if i= 1). Consider any γ∈C2. 1. If γhas two cofaces, they must be contained in the same connected component of G3. Thus, (d′)∗(γ) = g∗d∗(γ) = ri−ri= 0 2. If γhas only one coface then it belongs to T1. Thus (d′)∗(γ) = g∗d∗(γ) = 0 Therefore, Vcontains r−1critical 3-cells and this number is minimal, which completes the proof. Once again, Figure 4.3 illustrates the construction done in this proof in the two-dimensional space. On the left there is a 2D cubical complex with the cells of L−Kcolored in light blue, with the graph G2(L−K) superimposed. In the middle there is a spanning forest of G2(L−K)(in red). On the right we can appreciate the induced DGVF after choosing as root the top leftmost 2-cell of the complex. 4.2. The Iterative Algorithm 65 b b b b (A) b b b b (B) (C) FIGURE 4.3: Illustration of the construction done in the proof of Proposition 4.5 The previous propositions can directly be extended to higher dimensions. Also, Ldoes not need to be a 3D cubical complex of the form [0, m]3. We have considered the bounding box of K, but any acyclic supercomplex can be used. We recall that these statements are true for simplicial complexes, but obtaining an acyclic supercomplex of a simplicial complex of the same dimension does not seem to be a trivial task. We can count connected components using breadth first search. Algorithm 3explicitly illustrates this. Algorithm 3: Count connected components in a graph with BFS Input: A graph G= (V, E) Output: Number of connected components of G n←0; foreach u∈Vdo if unot marked then Q.push(u); mark u; while Qnot empty do u←Q.pop(); foreach vsuch that {u, v} ∈ Edo Q.push(v); mark v; n←n+ 1; return n We actually do not need to build the graphs G0(K)nor G3(L−K)to count its connected components, since they are included in Kand L−K. Algorithm 3has complexity O(|V|+|E|). In G0(K),|V|(resp. |E|) is the number of 0-cubes (resp. 1-cube) of K. Also, in G3(L−K)|V|(resp. |E|) is the number of 3-cubes (resp. 2-cube) of L−K. In addition, computing χ(K) only needs counting the cubes of Kwith the appropriate sign. Hence, Algorithm 4—called ViteBetti from the French word “vite” (fast)—computes the Betti numbers of a 3D cubical complex in O(n)time, where ndenotes the number of cubes in L. 66 Chapter 4. Fast Comp. of Betti Numbers on 3D Cubical Complexes Algorithm 4: ViteBetti Input: A 3D cubical complex K Output: Its Betti numbers: β0,β1,β2 χ←χ(K); β0←number of connected components of G0(K); β2←number of connected components of G3(L−K)−1; β1←β0+β2−χ; return β0,β1,β2 4.3 The Recursive Algorithm We showed in the previous section that the computation of the Betti numbers of a 3D cubical complex reduces to (1) compute the Euler-Poincaré characteristic and (2) to find the number of connected components of two graphs. We present in this section a way to parallelize Algorithm 4with a divide-and-conquer approach. Computing the Euler-Poincaré characteristic can be achieved in O(log(n)) time on a parallel machine. Assuming that the 3D cubical complex Kis encoded as a binary 3D array AK(called CubeMap in [115]), it suffices to split the array into two parts, recursively sum both of them and then sum the two values. The recursion stops whenever only two elements remain, in which case we sum them with their corresponding signs. The recursion depth of this method is ⌈log2(n)⌉, so the problem can be solved in O(log(n)) with n/2 processors. A simpler method consists in dividing the binary array into p parts, summing each of them in parallel and then summing all the partial results. This approach needs O(n/p +p)steps, so it can be done in O(√n) time with p=√nprocessors. There has been an extensive research about computing connected components in parallel, see [68,107,114] for some examples. We present here a simple recursive method that works well for commodity computers with few processors. Algorithm 3described how to count connected components by traversing the graph. This approach is not suited for parallelization since it uses a queue data structure. Another well known approach to compute connected components is to use the disjoint-set data structure (see [20, Chap. 21]). This data structure maintains a collection S={S1, . . . , Sk}of disjoint sets. Each set in Sis identified by a representative, which is a member of the set. The following operations can be performed on a disjoint-set data structure: •MakeSet(u)- creates a new set whose only member (and thus representative) is u. •Find(u)- returns a pointer to the representative of the (unique) set containing u. •Union(u, v)- merges the sets containing uand vinto a new set which is the union of these two sets. 4.3. The Recursive Algorithm 67 To compute connected components of a graph it is enough to call Union(u, v) for each pair (u, v)of adjacent vertices. A parallel version of such algorithm requires synchronization, so in practice it cannot be implemented efficiently. However, the regular structure of the cubical complex allows us to propose a different approach where synchronization is not needed. The idea is to recursively cut the graph in two halves, find the connected components in each half and then merge them. We recall that both G0(K)and G3(L−K)(without the special vertex ǫ) are grid graphs and that their vertices are identified to points of Z3via the Khalimsky coordinates. We define the left slice,right slice and middle slice of a subset of vertices Win dimension dby xrespectively as S(W, x−, d) := {u∈W|u= (u1, u2, u3), ud< x} S(W, x+, d) := {u∈W|u= (u1, u2, u3), x ≤ud} S(W, x0, d) := {u∈W|u= (u1, u2, u3), x −1≤ud≤x}. These three operations allow us to divide a set of vertices into two parts (the left and the right slice) plus a small subset that intersects both. Algorithm 5recursively computes the connected components of a grid graph. Observe that, at each step of the recursion, the set Wis divided into two parts. Each of them can be treated independently since there are no edges between both sets. We then consider the middle slice to combine both parts, which is not subdivided (since ǫ=∞). Algorithm 5: RecursiveCC Input: G= (V, E)a grid graph, ǫ > 0 S←disjoint-set data structure; RecursiveCC(V, 0, ǫ); return |S|; Procedure RecursiveCC(W, d, ǫ) Input: W⊂V,d∈Z,ǫ > 0 1if |W|> ǫ then 2d←((d+ 1) mod 3) + 1; x←middle point among the dth Khalimsky coordinates of W; 3RecursiveCC(S(W, x−, d), d, ǫ); 4RecursiveCC(S(W, x+, d), d, ǫ); 5RecursiveCC(S(W, x0, d), d, ∞); 6else 7foreach u∈Wdo 8MakeSet(u) 9foreach u∈Wdo 10 foreach v∈Wincident to udo 11 Union(u, v); At the end of Algorithm 5all the edges have been treated, so the disjointset data structure Scontains the connected components of the graph. In order to use this method for counting connected components in Algorithm 4 we need to clarify what happens with the vertex ǫof G3(L−K). Instead of directly counting the connected components of G3(L−K), we do it for 68 Chapter 4. Fast Comp. of Betti Numbers on 3D Cubical Complexes the induced subgraph without ǫ, which is a grid graph. Then, by merging all the sets in Scontaining a vertex incident to ǫin G3(L−K), we obtain the number of connected components of G3(L−K). However, it is more convenient to add a flag to the sets containing such vertices (at line 8) and then consider only one of those sets when counting |S|. 4.4 Results We compare in this section our algorithm with the library CAPD::RedHom [77], specialized in Betti numbers computation on cubical complexes. We have used the default settings, which executes shaving, coreduction algorithm, discrete Morse theory reduction and finally algebraic reductions. Three versions of our algorithm are tested: VB-i for the iterative version introduced in Section 4.2,VB-r for the recursive version described in Section 4.3 and VB-rp for the same algorithm using parallel computing. Our algorithms are implemented in C++ and compiled by the GNU compiler g++ (version 5.2.1) with option -O3. The parallel algorithm VBrp uses the Threading Building Blocks (TBB) library [70] (version 4.4). We have tested the algorithm with random 3D cubical complexes. A (2·m+1)3cubical complex consists in a cube of m33-cubes with their faces. Each cubical cell (and its faces) is added with a fixed probability p(see the random model K(p, m)in Section 3.9). We have made 10 random cubical complexes of size 513,1013,2013,3013,4013and 5013and probability 0.25, 0.5and 0.75. We use ǫ= 105as threshold for VB-r and VB-rp, which seems to be the best one for this implementation. We computed these 180 cubical complexes on a Dell PC with 3.70GHz ×8 Intel Xeon E5-1630 v3 CPU and 31.3 GB RAM. The experimental results are obtained by averaging the execution time of the algorithms, while the reading time (which exceeds the execution time for the bigger complexes) is omitted. Table 4.1 shows the results obtained. The Betti numbers calculated by each of the algorithms are obviously the same. Size RedHom VB-i VB-r VB-rp 5130.1842 0.0026 0.0026 0.0023 10131.268 0.0142 0.0148 0.0091 201310.78 0.1309 0.1232 0.0552 301340.89 0.4303 0.4176 0.1583 4013101.26 1.436 0.983 0.3092 5013— 3.609 1.977 0.5494 TABLE 4.1: Execution time (in seconds) versus the size of the cubical complex. We can appreciate that our algorithm clearly improves the execution time of RedHom. Moreover, RedHom uses a memory consuming data structure for its computation which does not allow us to process complexes bigger than 4013, which can be achieved by ViteBetti. The recursive version of our algorithm, VB-r, is clearly faster than VB-i. The TBB library automatically chooses the number of parallel threads, but we can appreciate that VB-rp is more than three times faster than VB-r for cubical complexes of size bigger than 4013. 5.2. The measures 75 Our measures are only defined for discrete objects and cannot consider individual homology classes, but they can be computed faster (there is no optimization problem) and they have very good geometric properties as it will be shown in this chapter. Let O⊂Zdbe a discrete object (see Section 2.4). We fix d= 3 for simplicity, but the generalization to any dimension is direct. We have two ways of building the cubical complex Kassociated to Odepending on which connectivity relation we choose: the 6 or the 26-connectivity relation. The distance transform dtOof Ois the map that sends every voxel x∈O to dtO(x) = d(x, O) = min {d(x, y)|y /∈O} where d(x, y) = v u u t 3 X i=1 (xi−yi)2 is the Euclidean distance. However, we can also consider other distances such as the Manhattan distance (L1), the chessboard distance (L∞), distances based on chamfer masks [92,9] or sequences of chamfer masks [97, 101]. The signed distance transform sdtOof Ois the map sdtO:Z3→Rdefined as follows: sdtO(x) = (−dtO(x) = −min {d(x, y)|y /∈O}if x∈O dtZ3\O(x) = min {d(x, y)|y∈O}if x /∈O Let us point out some simple properties about the sublevel sets L− t(sdtO) := sdt−1 O(]−∞, t]) of the signed distance transform. 1. O=sdt−1 O(]−∞,0]) 2. Z3=sdt−1 O(]−∞,∞[) 3. sdt−1 O(]−∞, a]) ⊂sdt−1 O(]−∞, b]) whenever a < b Figure 5.5 shows five sublevel sets at different values. Observe that the sequence of objects sdt−1 O(]−∞, t[)0 t=−∞ looks like an erosion of O, while sdt−1 O(]−∞, t[)∞ t=0 seems a dilation. We now define the filtration associated to the signed distance transform. A simple formulation of this filtration is F=K[L− t(sdtO)]t∈R where K[Y]denotes the primal or the dual associated cubical complex of Y. PDq(F)⊂R2denotes the persistence diagram in dimension qof this filtration. We denote TBq={(x, y)∈PDq(F)|x < 0, y > 0}. It is obvious that TBqcontains βq(K)pairs, that is, there are as many pairs in TBqas q-holes in O. Definition 5.1. Let O⊂Z3be a discrete object. Let us fix a distance function d:Z3×Z3→Rand a connectivity relation. Let q≥0, 76 Chapter 5. Measuring Holes FIGURE 5.5: Sublevel sets of the signed distance transform at values -10 (top-left), -5 (top-right), 0 (middle), 5 (bottomright) and 10 (bottom-right) •The thickness of the q-holes of Oare the values {−x|(x, y)∈TBq} •The breadth of the q-holes of Oare the values {y|(x, y)∈TBq} Observe that the thickness and the breadth of the holes appear in pairs. We can thus represent them as points in R2in the thickness-breadth diagram, just like the persistence diagrams. We call these points thickness-breadth pairs. The interpretation of the thickness-breadth diagram is similar to that of the persistence diagram: points close to the axes are holes with one small measure which may be originated by the presence of noise in the discrete object. Section 5.3 contains many examples of thickness-breadth diagrams. We speak about measures of the holes since we obtain as many values as holes in the object. However, we cannot measure a given hole (that is, a cycle xsuch that [x]6= 0). This is why we have preferred to talk about measures for the homology groups in [59], which seems a more correct formulation. Nevertheless, we think that speaking about measures for the holes sounds clearer. 5.2. The measures 77 5.2.1 On the computation of the measures In this section we give a more detailed description of how the measures are computed. Let us assume that the discrete object is contained in a bounding box BB := [0, w1]×[0, w2]×[0, w3]⊂Z3. The signed distance transform can be obtained by computing the distance transform of Oand BB \O in O(w1w2w3)[74]. Next, the filtration induced by sdtOdepends on the associated cubical complex that we consider. •Primal associated cubical complex: let K:= Kp[BB]be the primal cubical complex associated to the bounding box BB. We define the filtration induced by sdtOin terms of a function fOdefined on K. Recall that each 3-cube is identified with a voxel of BB. Thus, for every 3-cube σ∈K,fO(σ)takes the value of sdtOon its associated voxel. For the rest of the cubes, the value of fOis assigned so its sublevel sets are complexes. Namely, fO(σ) = min fO(τ(3))|τ > σ. In other words, fOmaps each cube σ∈Kto the first value tsuch that σ∈Kp[L− t(sdtO)]. •Dual associated cubical complex: the description is similar. Let K:= Kd[BB]be the dual cubical complex associated to the bounding box BB. Since every 0-cube is identified with a voxel of BB,fO(σ)takes the value of sdtOon its associated voxel for each 0-cube σ∈K. For the rest of the cubes, fO(σ) = max fO(τ(0))|τ < σ. Let a0< a1<··· <−1<1<··· < anbe the different values of fOover K. Thus, we can consider the filtration F: K0=f−1 O(]−∞, a0]) ⊂ ··· ⊂ K[O]⊂ ··· ⊂ Kn=f−1 O(]−∞, an]) = K Persistent homology is a very active field of research. The computation of the persistence intervals of a filtration has cubical worst-case complexity. An algorithm in matrix multiplication time was introduced in [88]. However, the most recent algorithms [10,7] are observed to have near linear complexity. An algorithm adapted for cubical complexes was developed in [115]. Algorithms for persistent homology consider a kind of elementary filtration where each step consists in adding only one cell to the previous complex. Thus, we need to decompose the filtration F. Some heuristics are given in [7] for this decomposition in order to accelerate the computation of the persistence intervals. We thus obtain some sets Pqof pairs of cells (σ(q), τ(q))for q≥0. The q-dimensional persistence diagram is PDq(f) = {(f(σ), f(τ)) |(σ, τ)∈Pq}∪{(f(γ),∞)|γnot paired} Observe that this set does not depend on how the filtration Fis refined. However, the pair of cells associated with each point does depend. This is relevant for the following sections. There are typically many points of the form (x, x)in the persistence diagram which are usually ignored. 78 Chapter 5. Measuring Holes Regarding the thickness-breadth diagram, there exists only one point (x, ∞)corresponding to the first connected component that appears in the filtration. This point can be plotted as (x, −1). 5.2.2 The robustness of the measures We prove in this section the robustness of the measures. The idea is that if we slightly deform an object, its measures suffer small changes. The celebrated article [19] introduced a theorem for the stability of the persistence diagrams. We recall that ||x||∞= max {|x1|,|x2|}for x∈R2and ||f||∞= max {|f(a)| | a∈A}for a function f:A→R. We can compare two persistence diagrams via the Hausdorff distance. Let X, Y be two multisets (sets with repetitions) of points in R2. Their Hausdorff distance is dH(X, Y ) = max max xmin y||x−y||∞,max ymin x||x−y||∞ where x∈Xand y∈Y. Thus, if dH(X, Y ) = ǫthen for every x∈Xthere is a y∈Ysuch that ||x−y||∞≤ǫand vice versa. Let us adapt one of the theorems of [19] to our context. Let fand gbe two functions on a cubical complex defining a filtration, that is, their sublevel sets are cubical complexes (each cube contains all its faces in the sublevel set). Let PD(f)and PD(g)be their respective persistence diagrams (all dimensions taken together). Therefore, dH(PD(f), PD(g)) ≤ ||f−g||∞ Consequently, if the two functions are similar, their associated persistence diagrams are also similar. Let now Xand Ybe two discrete objects. We call fXand fYthe filtration functions induced by sdtXand sdtYrespectively. Lemma 5.1. Let x, y ∈Z3be two 26-neighbors and A⊂Z3. Then |sdtA(x)−sdtA(y)| ≤ 2√3 Proof. We recall that, as xand yare 26-neighbors, |x−y|:= ||x−y||2≤√3. We have to consider four cases: •x, y /∈A. Let us call px, py∈Atheir closest points in A. Thus sdtA(x) = |x−px| ≤ |x−py| ≤ |x−y|+|y−py|=|x−y|+sdtA(y) Thus, sdtA(x)−sdtA(y)≤ |x−y| By symmetry, |sdtA(x)−sdtA(y)| ≤ |x−y| ≤ √3≤2√3 5.2. The measures 79 •x /∈A, y ∈A. Then sdtA(x)≤ |x−y| −sdtA(y)≤ |y−x| Thus, sdtA(x)−sdtA(y)≤2·|x−y| ≤ 2√3 By symmetry, |sdtA(x)−sdtA(y)| ≤ 2√3 •The other two cases follow the same arguments. Theorem 5.2. Let Xand Ybe two discrete objects in Z3. Let us call δ=dH(X, Y ) + dH(Z3\X, Z3\Y) + 2√3. Then, for every thickness-breadth pair pX= (x, y)of Xsuch that x, y > δ, there exists another thickness-breadth pair pY= (x′, y′)of Ysuch that ||pX−pY||∞≤δ Proof. Let σbe an elementary cube. Let us consider the two possible cubical complexes associated to the discrete objects. •Primal associated cubical complex: fX(σ) = fX(τX)for a 3-dimensional coface τX. Similarly, fY(σ) = fY(τY). Let qXand qYdenote the voxels associated to the 3-cubes τXand τYrespectively. •Dual associated cubical complex: fX(σ) = fX(τX)for a 0-dimensional face τX. Similarly, fY(σ) = fY(τY). Let qXand qYdenote the voxels associated to the 0-cubes τXand τYrespectively. In both cases qXand qYare 26-neighbors, so ||qX−qY|| ≤ √3. Therefore, |fX(σ)−fY(σ)|=|fX(τX)−fY(τY)| =|sdtX(qX)−sdtY(qY)| ≤ |sdtX(qX)−sdtY(qX)|+|sdtY(qX)−sdtY(qY)| By Theorem 2 of [81], ||sdtX−sdtY||∞≤dH(X, Y ) + dH(Z3\X, Z3\Y) Thus, |fX(σ)−fY(σ)| ≤ |sdtX(qX)−sdtY(qX)|+|sdtY(qX)−sdtY(qY)| ≤dH(X, Y ) + dH(Z3\X, Z3\Y) + 2√3 Consequently, dH(PD(fX), PD(fY)) ≤ ||fX−fY||∞≤δ 80 Chapter 5. Measuring Holes As the thickness-breadth diagram is the intersection of the persistence diagram with the quadrant {(x, y)∈R2:x, y ≥0}, the theorem follows from this. Observe that, if a thickness-breadth pair is close to the axes, a small perturbation in the dicrete object can make it disappear. In other words, a hole being not thick or broad enough can easily disappear after a small perturbation. Hence, in order to compare the thickness-breadth pairs of two discrete objects, they must be far enough from the axes. This is why we ask their values to be bigger than δin Theorem 5.2. Thus, we can bound the distance between two thickness-breadth diagrams via the Hausdorff distance of the two objects and their complements. 5.3 Thickness and breadth balls In the previous section we represented the thickness and the breadth of the holes as points in the thickness-breadth diagram. There is an alternative way, in terms of balls. Each thickness-breadth pair (t, b)has an associated pair of cubes (σ, τ). Therefore, •its thickness ball is the ball centered at the barycenter of σwith radius t; •its breadth ball is the ball centered at the barycenter of τwith radius b. Observe that the thickness balls are contained in the object, while the breadth balls are outside. Note however that the cubes σand τare not unique since they depend on how we decompose the filtration induced by the signed distance transform. The thickness and breadth balls allow us to represent both measures directly on the object. Moreover, the first impression we have when we encounter the breadth balls is that they are in the center of the holes. It is well accepted to visualize holes as representatives for a set of homology generators. For each representative, which is a chain, we mark those cells with a non-negative coefficient. Nevertheless, these representatives can be visually unpleasant. In order to better formalize this aspect, some authors suggest that the best representatives are those which are minimal in terms of their length, area, volume, etc [47,17,31]. Breadth balls emerge as an interesting alternative to the representation of holes in terms of homology generators. Symmetrically, thickness balls look like cohomology generators. In the following we show several examples of thickness-breadth diagrams and balls. We have considered several meshes from the AimAtShape repository1(except for Buddha2) and we have converted them to binary volumes using the software binvox3. The thickness-diagram and balls were computed with a specific software which does not take advantage on the 1http://visionair.ge.imati.cnr.it/ 2Courtesy of the Stanford Computer Graphics Laboratory 3http://www.patrickmin.com/binvox/ 5.4. Small generators 81 latest results in persistent homology computation [10,7], so we have omitted the time spent in these calculations. Figure 5.6 illustrates a voxelized version of Buddha. The thickness balls are shown in red, while the breadth balls are in green. We only display the balls of the 1-holes, since the other ones are less visually interesting. The thickness-breadth diagram shows the 0-holes (red circle), 1-holes (green triangles) and 2-holes (blue square). Observe that the thickness of the only connected components, which is ∞, is represented as −1. Also, the small (in both thickness and breadth) 2-hole is due to an error in the voxelization process of the mesh. Figures 5.7–5.16 show other models. 5.4 Small generators In this section we introduce a heuristic for obtaining well-shaped generators for the homology and cohomology groups based on the thickness and breadth balls. In the previous section we claimed that the thickness and breadth balls seem to be a good alternative for localizing holes, instead of displaying the homology generators. Therefore, it is a natural question to wonder if we can use these balls to find well-shaped homology generators. Moreover, the duality thickness/breadth of the measures provides also results for the cohomology groups. The localization problem [18] consists in finding the smallest representative cycle of a homology class with regard to a geometric measure. It seems natural that such cycles are good representatives for the holes. We explain some of these measures in a nutshell. Let Kbe a CW complex endowed with a weight function on its cells (its q-dimensional volume or just a constant). Some measures are: Volume The volume of a chain is the sum of the weights of its cells. Diameter The diameter of a chain is the maximal discrete geodesic distance (see Section 5.2) between the 0-dimensional faces of the cells in the chain. If Kis embedded in a metric space we can also consider the maximal distance between the 0-dimensional faces. Radius The radius of a chain is the radius of the smallest geodesic ball containing the chain. Again, if Kis embedded in a metric space we can consider the radius of the smallest ball containing the chain. Chen and Freedman [18] proved that finding a cycle minimizing the volume is an NP-hard problem, even if we look for an approximation. Considering the diameter is also an NP-hard problem, but we can compute a 2-approximation considering the radius, for which there is a polynomial time algorithm. Sadly, they also showed that considering the diameter or the radius does not always provide visually pleasant generators as they can wiggle. Our heuristic provides a cycle and a cocycle associated to each thicknessbreadth pair of a discrete object. The intuition is that a minimal homology generator must be around a breadth ball, while a minimal cohomology generator must traverse a thickness ball. However, given the complexity results exposed previously, we cannot expect to prove that these generators are optimal for the volume. 82 Chapter 5. Measuring Holes FIGURE 5.6: Buddha: thickness balls (in red), breadth balls (green) and thickness-breadth diagram. 5.4. Small generators 83 FIGURE 5.7: Casting: there are two types of 1-holes according to the breadth. We also observe that the narrow (less broad) holes do not have the same thickness, as they are not equally close to the border. 84 Chapter 5. Measuring Holes FIGURE 5.8: Dancing: the breadth balls clearly localize the 1-holes of the object. 5.4. Small generators 91 FIGURE 5.15: Neptune: there are three notable holes in this object. The rest, located in the beard and in the hand, can be considered as noise. 92 Chapter 5. Measuring Holes FIGURE 5.16: Pegasus: there are five significant holes and a small one near the right fron paw (see its thickness ball). 5.4. Small generators 93 5.4.1 Homology generators Let Obe a discrete object for which we have computed the thickness and the breadth of its holes. Let (σ, τ)be a pair of cubes associated to a thicknessbreadth pair and let ρ−= (h, f, g)be the reduction associated to the persistent homology computation performed for obtaining the measures of O before adding the cells with positive signed distance transform. Note that ρ−is a reduction for K[O]. Actually, we only need the chain f∗(σ). Algorithm 6provides a cycle associated to the homology generator g(σ). Algorithm 6: Homology generator Input: A discrete object O, its reduction ρ−, a pair of cubes (σ, τ) associated to a thickness-breadth pair Output: A cycle xsuch that hf(x), σi 6= 0 ~p ←some voxel associated to τ; F←filtration induced by d~p :O→R; Compute the persistent homology on F. When a cycle xis found, check if hf(x), σi 6= 0. If true, return this cycle; Let ~p be a voxel associated to τ: a 3-dimensional coface if we considered the primal associated cubical complex or a 0-dimensional face if we considered the dual associated cubical complex. It is not unique, so we choose one arbitrarily. Let d~p :O→Rbe the function that maps every voxel of O to its distance to ~p. Consider the filtration Fassociated to this function (as we did for the signed distance transform in Section 5.2). When we compute the persistent homology of this filtration we find a boundary associated to a negative cell or a cycle xassociated to a positive cell. Among these cycles, we take the first one for which hf(x), σi 6= 0. The condition hf(x), σi 6= 0 means that if we write the class [x]in terms of the homology base g(C)associated to ρ−, the coefficient of the generator g(σ)is not zero. This is necessary to capture the hole we want. Figure 5.17 shows a binary image with two holes and its breadth balls (in blue). While searching for a small homology generator for the broader hole, Algorithm 6 finds the small cycle xin the middle before the bigger cycle yon the right. The condition hf(x), σi 6= 0 avoids to return x, which does not correspond to the hole we chose. FIGURE 5.17: Left: a binary image with its breath balls in blue. Center and right: two cycles found during the computation of Algorithm 6. It seems that we could obtain an approximation of a minimal homology base by computing the chains associated to each thickness-breadth pair, but this is unclear since we could obtain a set which is not linearly independent. 94 Chapter 5. Measuring Holes 5.4.2 Cohomology generators Representing holes by the generators of the cohomology groups seems an unused approach. This is possibly due to the fact that they are not manifoldlike (as homology generators) since they are cocycles instead of cycles. Nevertheless, they can be visually interesting if we do not only display the cells in the cohomology generators but also all its cofaces. Note that when we display the homology generators we also display their faces, which is the dual statement of the previous sentence. We are interested in computing small cohomology generators since they display holes as thickness balls do. The heuristic for obtaining a cochain associated to a thickness-breadth pair is very similar. Algorithm 7provides a cocycle associated to the cohomology generator f∗(σ). Algorithm 7: Cohomology generator Input: A discrete object O, its reduction ρ−, a pair (σ, τ)of cubes associated to a thickness-breadth pair Output: A cocycle xsuch that hg∗(x), σi 6= 0 ~p ←some voxel associated to σ; F←cofiltration induced by d~p :O→R; Compute the persistent cohomology on F. When a cocycle xis found, check if hg∗(x), σi 6= 0. If true, return this cocycle; The input is the same, unless we only need the cochain g(σ). The main difference is that we do not consider a filtration but a cofiltration. In Algorithm 6we compute the persistent homology of the filtration induced by d~p because we want to consider all the cycles that appear in this filtration. Using the algorithm for persistent homology avoids to consider cycles that are boundaries, but we actually do not need the persistence intervals of this filtration. In this context we want to find cocycles, so we take a dual approach. Let Kbe a cubical complex. L⊂Kis a sub-cocomplex if for each cube of L, all its cofaces in Kare also included in L. A cofiltration of Kis a sequence of nested sub-cocomplexes ∅=K0⊂ ··· ⊂ Km=K. Thus, a cube enters the cofiltration before its faces. We can adapt the persistent homology algorithm by considering coboundaries instead of boundaries, which gives a coboundary or a cocycle for each cube. As in Algorithm 6, we take the first cocycle such that hg∗(x), σi 6= 0. Let us present a few results of these two algorithms. Figure 5.18 depicts a thickened wire-frame cube. We can appreciate its homology generators (in green) and cohomology generators (in red) of dimension 1. Observe that two of its cohomology generators are too close to be distinguished. A double torus is shown in Figure 5.19, whose generators are tight. The object in Figure 5.20 does not have four holes, but three. They are located by the generators. Let us point out that the generators obtained for these three objects are linearly independent and hence they conform a (co)homology base. However, as mentioned above, this is not guaranteed for our two algorithms. 5.5. Opening or closing holes 95 FIGURE 5.18: An object with its homology generators (top, in green) and homology generators (bottom, in red). 5.5 Opening or closing holes If [x]is a non-trivial homology class, it means that xis a cycle, but it is not a boundary. Thus, if we add a “coface” so xbecomes a boundary, it will be a trivial class. This means that a homology generators is the boundary of a chain that is missing. Hence, in order to remove a hole, we can add such chain. Let us see now why this reasoning does not work for cohomology. A non-trivial cohomology class [x]is a cocycle that it is not a coboundary. Thus, if we want it to become a coboundary, we must add a “face” whose coboundary is x. However, this does not make sense since a complex already contains all its faces. We define now what is to open and to close a homology class. Given 96 Chapter 5. Measuring Holes FIGURE 5.19: A double torus. Model from the AimAtShape repository. two cubical complexes K⊂L, the inclusion map ι:K→Linduces a chain map ι#:C(K)→C(L)between their corresponding chain groups. Hence, ι#induces a homomorphism ι∗:H(K)→H(L)between their homology groups, that is, ι∗([x]) is the homology class of the chain xin H(L). Definition 5.2. Let Kbe a cubical complex, xa cycle and Sa set of cubes. Then, •Sopens the chain xif K−Sis a cubical complex, ι∗:H(K−S)→H(K) is injective and [x]/∈im(ι∗). We say that Sis an opening set for x. •Scloses the chain xif K∪Sis a cubical complex, ι∗:H(K)→H(K∪S) is surjective and [x]∈ker(ι∗). We say that Sis a closing set for x. Let [x]be a non-trivial homology class of H(K). The injectivity of ι∗ means that K−Sdoes not contain new holes, and [x]/∈im(ι∗)implies that we have removed the hole [x]. A similar interpretation follows from the definition of closing a chain. Note that, by closing a q-hole, several other holes may disappear. For instance, closing a 1-hole can merge two connected components (think of two chained circles) and closing the 2-hole of a torus 5.5. Opening or closing holes 97 FIGURE 5.20: An object with three 1-holes. Model from the AimAtShape repository. removes one of its 1-holes. Also, opening a 0-hole removes all the higherdimensional holes in the connected component and opening a 1-hole in a torus removes its 2-hole. The terms open and close are inspired by the work of Aktouf et al. [1]. They introduced the concept of topological hull for a three-dimensional discrete object O: it is a minimal (for inclusion) superset TH ⊃Owhich has no holes or cavities in the sense of digital topology [79]. Later, Janaszewski et al. [71] made the distinction between closing a hole (adding a minimal set of voxels to remove the hole) and filling a hole (the set is not minimal since it tries to fit the local geometry of the object). Observe, however, that these notions are defined for discrete homotopy and not for homology. These definitions evoke the following problem: what is the minimal opening/closing set (under any geometric criterion) for a given chain, or more generally, what is the minimal set that opens/closes all the holes of 98 Chapter 5. Measuring Holes a cubical complex? This seems to be a hard problem. Intuitively, the intersection of a minimal closing (resp. opening) set with the object gives a small homology (resp. cohomology) generator. Thus, we can even expect this problem to be NP. Consequently, as in Section 5.4, we only provide algorithms that seem to work well, without any proof of optimality. 5.5.1 Opening a hole The previous section, where several cohomology generators were shown, should have suggested an idea: removing a cohomology generator erases a hole. We cannot prove this fact because it is not true in general. Figure 5.21 illustrates a counter-example. It shows a small simplicial complex with one 1-hole. It has a cohomology generator involving three 1-simplices (marked in red). After removing them, and their cofaces, there is still a 1-hole in the complex. Note that we could have removed other cohomology generators which do open the hole. FIGURE 5.21: Left: a simplicial complex Kwith a cohomology generator x(in red). Right: after removing x(and its cofaces) from K, there is still a 1-hole. Let Kbe a cubical complex endowed with a perfect HDVF X= (PX, SX) whose critical cells are CX={γ1, . . . , γr}. If we have a subcomplex L⊂K endowed with a perfect HDVF Y= (PY, SY)such that PY⊂PX,SY⊂SX and CY⊂CX−{γ1}, then we have opened the homology generator gX(γ1) associated to the critical cell γ1in K: •Lis a cubical complex •Since both HDVFs are perfect, [gX(γ1)] = fXgX(γ1) = γ1, which does not belong to CY, so [gX(γ1)] /∈im(ι∗) •CY⊂CXimplies that ι∗is injective. We can obtain such cubical complex by removing a critical cell and pairs of cells until we obtain a cubical complex. Algorithm 8describes this procedure. At the end we obtain a perfect HDVF for the subcomplex K−Owhich does not contain γ. Thus, K−Ohas at least one hole less, and no holes are created. Before proving it, we need a lemma that ensures that we can find the cells τand σat lines 6and 9respectively. 5.5. Opening or closing holes 99 Algorithm 8: Hole opening Input: A CW complex Kendowed with a perfect HDVF X= (P, S), a critical cell γ Output: An opening set Ofor the homology generator g(γ) 1Q.push(d∗(γ)); 2O← {γ}; 3while Qnot empty do 4a←Q.pop(); 5if a∈Pthen 6choose τ > a st. hh(a), τi 6= 0;X←R(X, a, τ); 7O←O∪{a};Q.push(d∗(a)); 8else if a∈Sthen 9choose σ < a st. hh(σ), ai 6= 0;X←R(X, σ, a); 10 O←O∪{σ};Q.push(d∗(σ)); 11 else if a∈Cthen 12 O←O∪{a};Q.push(d∗(a)); 13 return O Lemma 5.3. Let Kbe a CW complex endowed with a perfect HDVF X= (P, S). Then, •If σ∈Pthen there exists τ∈S,τ > σ such that hh(σ), τi 6= 0. •If τ∈Sthen there exists σ∈P,σ < τ such that hh(σ), τi 6= 0. Proof. Let σ∈Kbe a primary cell. Since d(S)|P·H=I, then d(S)|σ·h(σ)|S= 1. Thus, there exist some τ∈Ssuch that hd(τ), σi 6= 0 and hh(σ), τi 6= 0 The second statement follows from H·d(S)|P=I. Next, we need yet another lemma. Lemma 5.4. Let K⊂Lbe two CW complexes endowed with two HDVFs X⊂Y respectively. If Yis perfect then Xis perfect too. Proof. The boundary matrix of Lis of the form D=D1· 0D2 where D1=d(K)|Kand D2=d(L−K)|L−K. Thus, d(S)P·H=uv 0w·xy z t =I0 0I and H·d(S)P=xy z t ·uv 0w·=I0 0I 100 Chapter 5. Measuring Holes Note that uis invertible (since Xis a HDVF). Hence, xu +y0 = I⇒x=u−1 zu +t0 = 0 ⇒z= 0 Therefore, H=u−1· 0· Consequently, as Yis perfect, d(C)|C= 0 ⇒A· 0·=B· 0··u−1· 0··D· 0· so A=Bu−1Dand thus Xis perfect. We can now prove the correctness of Algorithm 8. Proposition 5.5. Algorithm 8returns an opening set for the chain g(γ). Proof. Let Kbe a CW complex endowed with a perfect HDVF X. Let γbe a critical cell. We have to prove that Oopens the chain g(γ). First, let us prove that K−Ois a cubical complex. This is equivalent to prove that for any cube in O, its cofaces (in K) are also included in O. Indeed, whenever we add a cell to O, we add its cofaces to the queue Q. Let abe one of those cofaces. If ais primary or critical, it will be added to Owhen it is taken from Q. If it is secondary, it will become critical and will be added again to Q, so it will eventually be added to O. We denote by X′the resulting HDVF at the end of the algorithm. Let us prove now that X′on K−Ois perfect. By construction, X′is a HDVF for K−Oand, since Xis perfect for K, the result follows from Lemma 5.4. As γdoes not belong to K−Oand there are no new critical cells, it follows that Oopens the chain g(γ). Thus, given a discrete object Oand a thickness-breadth pair (σ, τ), we can remove the homology generator associated to σusing Algorithm 8with the HDVF on K[O]obtained while computing the measures. However, in practice, it seems that it suffices to remove the cohomology generator associated to σin the HDVF and its cofaces. We have performed several experiments on objects with complex geometry and we have never found an example as Figure 5.21. Note that, in that example, there are other cohomology generators which do open the hole. We are not able to prove why this works, or in what cases it does. 5.5.2 Closing a hole Closing a hole has also an intuitive answer which is false. Let xbe a homology generator in a complex K. If we add a set of cells such that its boundary is x, then it becomes a trivial class in the homology group. However, this does not necessarily close the hole. Let Kbe the boundary of a Möbius strip, which is homotopy equivalent to S1, and xthe chain with all its 1-cells. By adding the interior of the strip, xbecomes a boundary but there is still a 1-hole. This is illustrated in Figure 5.22. 5.6. Conclusion and future works 107 the space (or at least the convex hull of the complex) and assign the signed distance transform values to the simplices. The main difficulty is thus how to define/compute a triangulation Tof the space with a parameter, as we are going to compute an approximation of the measures in the continuous space. 109 Chapter 6 Conclusion APART from the individual conclusions in each chapter, we summarize here the main results of this essay and, more importantly, we discuss the future perspectives of this work. Later we describe two works that have not been developed in this dissertation. 6.1 General conclusion The main goal of this thesis is to study discrete objects—that is, binary images, volumes or their equivalent notions in higher dimensions—from a homological point of view. The initial objective was to extend the work on the homological spanning forest developed by H. Molina-Abril and P. Real [14,91,90], considering geometric features of the discrete objects such as the curvature or the medial axis. We will now consider the history of each chapter and the relation between them. 6.1.1 Homological Discrete Vector Field The homological spanning forest naturally led to the concept of the homological discrete vector field (HDVF), which has a richer structure and better properties. However, finding the correct definition was far from being a trivial task. When studying the homological spanning forest we were disappointed that it did not work well for classical examples in discrete Morse theory such as the Bing’s house or the dunce hat. The homological discrete vector field was born while studying the former. One cannot find a perfect DGVF on it since at some point we find three V-paths (instead of one) between two critical cells that should be paired and we cannot reverse any of them since this would create a closed V-path in the DGVF. Despite this, we decided to reverse one. The operator hin a reduction is the most important since the others are defined through it. In the reduction associated to a DGVF, his defined recursively and thus it is not well defined if there are closed V-paths. We thought that we could compute hin terms of the chains h(σ) where σis a confluence cell, that is a cell where two or more closed V-paths merge. Hence we obtain a system of linear equations, which we found that always had a solution. This approach was published in [57], but we did not understand why this works and we could not prove it. We later found the correct definition inspired by the formula 1 + x+x2+··· =1 1−x 110 Chapter 6. Conclusion The classical definition of the operator his similar to the left side of the formula. As in our case the geometric sequence never vanishes, we must use the right side of the formula. Thus, we need something (actually a submatrix of the boundary matrices) to be invertible. We also changed the formulas for the reduction in order to have a clean definition. This is the formalism we used in [60]. Given a CW complex K, a homological discrete vector field is a pair of disjoint subsets (P, S)of the cells of Ksuch that the boundary matrix of Krestricted to these sets is invertible. Given this property, we can define a reduction on the chain complex associated to K, so we can compute its homology groups. This definition generalizes several concepts such as the discrete gradient vector field, the iterated discrete gradient vector field and the reduction induced by the Smith normal form. This structure reveals an idea which is not evident in other methods for computing homology. There are many possible different HDVFs for the same CW complex. They are actually not completely independent, since we can define basic relations between them. Thus, a global structure containing all the possible HDVFs for a CW complex appears, which has not been completely understood. There is a natural algorithm for building a HDVF. We have studied how to efficiently compute the associated reduction and we have estimated its complexity both in theory and in practice, by considering random cubical complexes. The results show that the worst-case complexity of building a HDVF is O(n3)but significantly less in practice. Open questions We recall that a HDVF is perfect if its associated reduction is perfect, that is, if it reduces the chain complex associated to the CW complex to its homology groups. We do not know if a perfect HDVF exists for every CW complex. We have found chain complexes (actually just matrices) for which this is false, so we think that there may be exceptions. Nevertheless, the boundary matrices of simplicial, cubical or regular CW complex have certain properties that may guarantee the existence of a perfect HDVF. We have tried to find a counter-example by brute-force and we have not succeeded. Observe, however, that this problem is solved if the ring of coefficients is a field. The previous problem assumes that the considered CW complexes have torsion-free homology groups, because there exists no perfect reduction (and thus, HDVF) otherwise. However, we have found that there exist HDVFs whose reduced boundary matrix is already in the Smith normal form. Let us call such kind of HDVF a pseudo-perfect HDVF. Thus, like in the previous paragraph, we cannot conclude that any CW complex admits a pseudo-perfect HDVF. Given the operations introduced in Section 3.7, we could define a kind of edit distance [20, § 15] between HDVFs. Given two HDVFs, their distance can be the minimal number of operations A,R,M,W,MW (or a subset of them) needed to transform one into the other. It is not clear that such a number exists, since one HDVF may be impossible to transform into another, so we should consider an extended metric (with infinite value for non-connected HDVFs). We find this problem fascinating, though we cannot see any practical application for this. 6.1. General conclusion 111 Possibly, the most useful perspective to this work is to compute zigzag persistence homology with the HDVF framework. We have explained in Section 3.8 how to compute (standard) persistent homology with HDVFs, which gives a very clear idea of what persistent homology represents. Zigzag persistent homology [13] is a more recent and complex theory which lacks a clear intuition. We are working on an algorithm for computing zigzag persistent homology using HDVFs and their operations. The advantage of doing this, rather than reducing its complexity, is to make this theory clearer and easier to understand. 6.1.2 Fast Computation of Betti Numbers on Three-Dimensional Cubical Complexes In order to compute the homology groups of discrete objects, one needs to first transform them into cubical complexes. The algorithm introduced in this chapter computes (only) the Betti numbers of 3D cubical complexes. We thus can compute the Betti numbers of a discrete object in linear time with regard to the size of the bounding box. This algorithm reduces the problem of computing the Betti numbers to count connected components in two graphs defined on the cubical complex. This is proved using the HDVF framework. The advantage of the cubical complex is not only that we know its complement (which is used), but that its regular structure allows us to easily subdivide the graphs into subgraphs with a reasonable number of edges between them. Three versions of this algorithm have been implemented depending on how we count the number of connected components. Even the slowest one (the sequential version) outperforms the library RedHom, specialized in homological operations on cubical complexes. Open questions Delfinado and Edelsbrunner sketched an algorithm in [29] for computing the Betti numbers of a simplicial complex embedded in R3which is not necessarily a subcomplex of a triangulation of S3, though it is not clearly proved. Later, Dey and Guha proved in [32] this algorithm, though complexes not being a three-manifold require a pre-processing step. The idea is similar to our algorithm, except that β2is obtained by recognizing closed surfaces on the boundary of the complex. We are trying to find a simpler algorithm and proof using the HDVF framework for both simplicial and cubical complexes. We also aim at computing the Betti numbers directly on a 3D discrete object for the 6 and the 26-connectivity relation. We can do this by building the associated primal (or dual) cubical complex and using our algorithm, but it seems easy to do this directly on the discrete object. For instance, given a discrete object Xwith the 6-connectivity relation, β0is the number of connected components of X,β2is the number of connected components in its complement with the 26-connectivity relation (minus the unbounded component) and the Euler characteristic can be computed locally in each voxel of X. This idea is already present in [100], but the theoretical foundations of this work are not clear. In order to assign a topological space to a discrete object, we can use the primal associated cubical complex 112 Chapter 6. Conclusion for the 6-connectivity and the dual associated cubical complex for the 26connectivity. As these complexes are locally built, it is easy to transfer the calculations of our algorithm in the cubical complex to the original object. Finally, we want to explore other methods for counting connected components. In the image processing literature, connected components labeling is generally performed by traversing the volume in raster scan order (first incrementing the first coordinate, then the second and later the third one) and making an equivalence table between the voxels that are adjacent. This strategy avoids accessing the voxels in a random order (such as while performing a breadth first search), which provokes page faults. We intend to try these other methods to make our algorithm run faster, and also to treat huge complexes (bigger than 10013) which can be read by slices in order to fit in the virtual memory. 6.1.3 Measuring Holes It took us a long time to combine digital geometry and homology. Motivated by a problem in geostatistics [23], where the Betti numbers do not give any information about the shape or size of the holes, we developed a simple definition of the size of the holes using the signed distance transform and persistent homology. It was a great surprise to find that holes do not have one natural measure but two, since we can erase them by cracking or filling them. Fortunately, these measures are stable under small perturbations, which makes them suitable for applications where objects come from acquisition devices or computer simulations. Surprisingly, the computation of these measures provides a novel representation of holes in terms of balls (instead of homology generators) which we show to be useful through examples. We also include some recent research about obtaining homology and cohomology generators and about opening and closing holes using these measures. Open questions We want to emphasize three perspectives about this research. The definition of the breadth and the thickness is not only valid for cubical complexes, but for any subspace of Rn. We want to compute the measures for simplicial complexes embedded in R3. The challenge is how to triangulate—tetrahedrize, actually—the space according to the signed distance transform. We think that we can only compute an approximation of the measures, since we can only compute the signed distance transform in a finite number of points. We are convinced that computing small generators for (co)homology and opening or closing holes in this context will give better results since the output is not restricted to fit in the grid. While the problem of finding small homology generators has been largely studied from a computational point of view, that of closing and opening holes seems to be overlooked. We do not aim at proving that finding a minimal closing or opening set is a NP problem, but we believe it is. Even though the measures were conceived to solve a practical problem (studying the holes of different types of soils), we really do not have any application for them. Specific needs may give new ideas, such as using different distances. In particular, some recent works [108,109,113] study 6.2. Other works 113 the shape of the universe through persistent homology, and we think that the breadth and the thickness may be useful for this. 6.2 Other works We present here two works that are not described in this thesis for different reasons. 6.2.1 Cellular skeletons One of the main theorems in discrete Morse theory states that a simplicial complex endowed with a DGVF is homotopy equivalent to a CW complex with exactly one q-cell for each critical q-simplex. From a topological point of view, the skeleton of a discrete object is a subset which is homotopy equivalent. Thus, it is clear that a DGVF—moreover, a reduction—provides a kind of skeleton for a discrete object. In [58] we introduced the concept of cellular skeleton. Given a reduction from the primal (or dual) associated cubical complex of a discrete object, it is the set of chains g(C). Taking Z2as ring of coefficients, each chain is a set of cells comprehending a manifold-like part of the skeleton. Thus, this is not just a subset of voxels but a chain complex. FIGURE 6.1: Cellular skeleton of Fertility. These skeletons were computed using known algorithms for homotopy thinning for cubical complexes such as [84,41,16,21,22]. Then we used a cell clustering algorithm to reduce the number of cells without changing the shape of the skeleton. Figure 6.1 shows an example of the visualization of a cellular skeleton. The cellular skeleton has several advantages. On one hand, the resulting skeleton preserves the topology of the original object since it is defined through a reduction. On the other hand, since the previously mentioned algorithms give (hopefully) centered skeletons, thus a subsequent homology computation on the cellular skeleton should give centered homology generators which, though not minimal, are visually pleasant. The cell clustering algorithm was not fully understood in [58]. It takes a 3D cubical complex as input and it returns a reduction. We later conceived a general version for any dimension and we proved that every maximal 114 Chapter 6. Conclusion (without cofaces) cell of the complex belongs to one and only one chain g(γ). This implies that the shape of the skeleton does not change. On the other hand, the resulting skeleton is not unique (since there are multiple choices in the algorithm) and it does not always return a minimal (in the number of cells) skeleton. These results have not been published, but they seem to be equivalent to the works of Damiand et al. [25,24] in the context of generalized maps. Sadly, we did not advance more in this direction due to the lack of practical applications for this theory. 6.2.2 Opening holes in discrete objects There is a simple way of closing the holes of a discrete object. Let Xbe a 3D discrete object with the 26-connectivity. Consider an acyclic volume Y⊃X containing it (its bounding box, for instance) and remove simple points— that is, voxels that can be removed without changing the homotopy type of the object—in Y−Xuntil idempotency. The final volume contains Xand has no holes since it is homotopy equivalent to Y, so it closes the holes of X. Giving priority to voxels which are further from the object usually gives minimal fillings. This has been studied in [1,71,72]. Inspired by the duality in the measures (see Chapter 5), we would like to open holes in discrete objects. A simple algorithm consists in considering a point xinside the object and then adding simple points in the object until idempotency. The object obtained is a subset of the original object without holes. Again, considering the distance transform for the order in which the points are added should give a minimal result. Simple examples in 2D such as those in Figure 6.2 show that this natural approach does not give optimal fractures. Instead of obtaining the minimal cuts opening the holes, the fractures seem to follow some angles. FIGURE 6.2: Result for the same image with different resolution. The cuts in the figure on the right have been thickened for visibility. This idea is very recent and we have not had time enough to develop it. We decided to include it here for its simplicity and possible future perspectives. We plan to explore more advanced techniques for opening holes that converge to the minimal cuts. 115 Bibliography [1] Zouina Aktouf, Gilles Bertrand, and Laurent Perroton. A threedimensional holes closing algorithm. Pattern Recogn. Lett., 23(5):523– 531, March 2002. [2] Madjid Allili and David Corriveau. Topological analysis of shapes using Morse theory. Computer Vision and Image Understanding, 105(3):188–199, 2007. [3] Rafael Ayala, Desamparados Fernández-Ternero, and José Antonio Vilches. Perfect discrete Morse functions on 2-complexes. Pattern Recogn. Lett., 33(11):1495–1500, August 2012. [4] Bruno Benedetti and Frank H. Lutz. Random discrete Morse theory and a new library of triangulations. Experimental Mathematics, 23(1):66–94, 2014. [5] Ainhoa Berciano, Helena Molina-Abril, and Pedro Real. Searching high order invariants in computer imagery. Appl. Algebra Eng. Commun. Comput., 23(1-2):17–28, 2012. [6] R. H. Bing. Some aspects of the topology of 3-manifolds related to the Poincaré conjecture. Lectures on Modern Mathematics, 2, 1964. [7] Jean-Daniel Boissonnat, Tamal K. Dey, and Clément Maria. The compressed annotation matrix: An efficient data structure for computing persistent cohomology. Algorithmica, 73(3):607–619, 2015. [8] Dobrina Boltcheva, Sara Merino Aceitunos, Jean-Claude Léon, and Franck Hétroy. Constructive Mayer-Vietoris algorithm: Computing the homology of unions of simplicial complexes. Research Report RR-7471, INRIA, December 2010. [9] Gunilla Borgefors. Distance transformations in arbitrary dimensions. Computer Vision, Graphics, and Image Processing, 27(3):321–345, 1984. [10] Peer-Timo Bremer, Ingrid Hotz, Valerio Pascucci, and Ronald Peikert, editors. Topological Methods in Data Analysis and Visualization III, Theory, Algorithms, and Applications. Springer, 2014. [11] Piotr Brendel, Paweł Dłotko, Graham Ellis, Mateusz Juda, and Marian Mrozek. Computing fundamental groups from point clouds. Applicable Algebra in Engineering, Communication and Computing, 26(1):27–48, 2015. [12] Kenneth S. Brown and Ross Geoghegan. An infinite-dimensional torsion-free FP∞group. Inventiones mathematicae, 77(2):367–381, 1984. [13] Gunnar Carlsson and Vin Silva. Zigzag persistence. Foundations of Computational Mathematics, 10(4):367–405, 2010. 116 BIBLIOGRAPHY [14] Javier Carnero, Helena Molina-Abril, and Pedro Real. Triangle mesh compression and homological spanning forests. In Computational Topology in Image Context - 4th International Workshop, CTIC 2012, Bertinoro, Italy, May 28-30, 2012. Proceedings, pages 108–116, 2012. [15] Manoj K. Chari. On discrete Morse functions and combinatorial decompositions. Discrete Mathematics, 217(1–3):101–113, 2000. [16] John Chaussard and Michel Couprie. Surface thinning in 3D cubical complexes. In Petra Wiederhold and Reneta P. Barneva, editors, Combinatorial Image Analysis, volume 5852 of Lecture Notes in Computer Science, pages 135–148. Springer Berlin Heidelberg, 2009. [17] Chao Chen and Daniel Freedman. Measuring and computing natural generators for homology groups. Comput. Geom., 43(2):169–181, 2010. [18] Chao Chen and Daniel Freedman. Hardness results for homology localization. Discrete & Computational Geometry, 45(3):425–448, 2011. [19] David Cohen-Steiner, Herbert Edelsbrunner, and John Harer. Stability of persistence diagrams. Discrete & Computational Geometry, 37(1):103–120, 2006. [20] Thomas H. Cormen, Clifford Stein, Ronald L. Rivest, and Charles E. Leiserson. Introduction to Algorithms. McGraw-Hill Higher Education, 2nd edition, 2001. [21] Michel Couprie. Hierarchic Euclidean skeletons in cubical complexes. In Discrete Geometry for Computer Imagery - 16th IAPR International Conference, DGCI 2011, Nancy, France, April 6-8, 2011. Proceedings, pages 141–152. 2011. [22] Michel Couprie. Topological maps and robust hierarchical Euclidean skeletons in cubical complexes. Computer Vision and Image Understanding, 117(4):355–369, 2013. [23] Asmae Dahrabou, Sophie Viseur, Aldo Gonzalez-Lorenzo, Jeremy Rohmer, Alexandra Bac, Pedro Real, Jean-Luc Mari, and Pascal Audigane. Topological comparisons of fluvial reservoir rock volumes using Betti numbers: Application to CO2storage uncertainty analysis. In 6th International Workshop on Computational Topology in Image Context (CTIC 2016), Lecture Notes in Computer Science (LNCS 9667), pages 101–112. Springer International Publishing, 2016. DOI:10.1007/9783-319-39441-1_10. [24] Guillaume Damiand, Rocío González-Díaz, and Samuel Peltier. Removal operations in nD generalized maps for efficient homology computation. In Computational Topology in Image Context - 4th International Workshop, CTIC 2012, Bertinoro, Italy, May 28-30, 2012. Proceedings, pages 20–29, 2012. [25] Guillaume Damiand, Samuel Peltier, and Laurent Fuchs. Computing homology generators for volumes using minimal generalized maps. In Combinatorial Image Analysis, 12th International Workshop, IWCIA 2008, Buffalo, NY, USA, April 7-9, 2008. Proceedings, pages 63–74, 2008. BIBLIOGRAPHY 123 [110] J. R. Stallings. Lectures on polyhedral topology. Tata Institute of Fundamental Research, Bombay, 1968. [111] Arne Storjohann. Near optimal algorithms for computing Smith normal forms of integer matrices. In Proceedings of the 1996 International Symposium on Symbolic and Algebraic Computation, ISSAC ’96, pages 267–274, New York, NY, USA, 1996. ACM. [112] Takashi Teramoto and Yasumasa Nishiura. Morphological characterization of the diblock copolymer problem with topological computation. Japan Journal of Industrial and Applied Mathematics, 27(2):175–190, 2010. [113] Rien van de Weygaert, Gert Vegter, Herbert Edelsbrunner, Bernard J. T. Jones, Pratyush Pranav, Changbom Park, Wojciech A. Hellwing, Bob Eldering, Nico Kruithof, E. G. P. (Patrick) Bos, Johan Hidding, Job Feldbrugge, Eline ten Have, Matti van Engelen, Manuel Caroli, and Monique Teillaud. Transactions on computational science XIV. chapter Alpha, Betti and the Megaparsec Universe: On the Topology of the Cosmic Web, pages 60–101. Springer-Verlag, Berlin, Heidelberg, 2011. [114] Uzi Vishkin. An optimal parallel connectivity algorithm. Discrete Applied Mathematics, 9(2):197–207, 1984. [115] Hubert Wagner, Chao Chen, and Erald Vuçini. Efficient computation of persistent homology for cubical data. In Ronald Peikert, Helwig Hauser, Hamish Carr, and Raphael Fuchs, editors, Topological Methods in Data Analysis and Visualization II, Mathematics and Visualization, pages 91–106. Springer Berlin Heidelberg, 2012. [116] Michael Werman and Matthew L. Wright. Intrinsic volumes of random cubical complexes. Discrete & Computational Geometry, 56(1):93– 113, 2016. [117] E.C. Zeeman. On the dunce hat. Topology, 2(4):341–358, 1963. [118] Afra Zomorodian and Gunnar Carlsson. Localized homology. Comput. Geom. Theory Appl., 41(3):126–148, November 2008. [119] Afra Zomorodian and Gunnar E. Carlsson. Computing persistent homology. Discrete & Computational Geometry, 33(2):249–274, 2005.