Our goal in this paper was to develop and compare different implementations of a computational topology algorithm needed to solve efficiently and robustly the space decomposition (3D arrangement generated by polyhedral complexes) requested to evaluate Boolean formulas of Constructive Solid Geometry [2019arXiv191011848P].
Within the space decomposition pipeline to compute the arrangement [Halperin:2017] generated by a collection of cellular complexes [2017arXiv170400142P, 2019arXiv191108130P], the intermediate step is the assessment of the operator in 3D, after independent fragmentation of 2-faces in 2D, and before the topological gift wrapping (TGW) algorithm in 3D.
In particular, the (co)boundary chain is produced by merging a set of local coboundary matrices computed independently from each other on every fragmented input 2-cell. In paper [2019arXiv191108130P] this computational pipeline is discussed in full detail. In the present manuscript we provide the reduction of a chain complex to its “quotient topology” through congruence quotient computation.
We recall that the 2-skeleton of 3D arrangements, and its matrix, is computed by merging the local chain complexes produced by mutual intersections of 2-cells of input complexes. Each local chain complex , with , is generated by a two-dimensional input cell in the input collection, and by all 2-cells intersecting it. We call accumulator complex the set , . Two figures are geometrically congruent iff one can be transformed into the other by an isometry [2017arXiv170400142P]. The congruences between -cells of geometric complexes in are equivalence relations, so to compute the chain complex of quotient spaces by merging chains and operators by computing the graded quotient .
The merging is based on the discovering of all equivalence classes of the congruence relations between0-, 1-, and 2-cells of decomposed input objects, computed independently [2017arXiv170400142P] for each input 2-face. In order to discover all congruent pairs among decomposed -cells, topological tests are introduced through the computation of local topological invariants. Our Cell Congruence Enabling algorithm may be used to merge the topologies of several cellular complexes (even a large number), after having computed the intersections of their 2-cells. Two examples of actual application of the CCC algorithm are (a) the merging of plan drawings of city districts or buildings within a digital cartographic map, and (b) the mutual fragmentation of cell complexes before computing Boolean expressions [2019arXiv191011848P] between solid models.
The description of topology, geometry and physics of models using the chain of coboundary operators, i.e., of vector algebra operators, and their representation as sparse matrices, started with[PALMER1995733, Palmer1993], at least in the knowledge of the authors. While still unusual in the geometric and solid modeling areas (see [tuprints11291] for example), they seem currently to us the most promising approach to FEA methods [Desbrun:2006:DDF:1185657.1185665, Elcott:2006:BYO:1185657.1185666, arnold_falk_winther_2006, Arnold:2010, Arnold:2018, Tonti:1975, Tonti:2013, Ferretti:2014]. Also, the GraphBLAS standard, with linear algebra methods based on sparse matrices for fast algorithms on huge graphs, is currently having a big momentum [graphblas:20].
The computation of (co)boundary matrices from multiplication of sparse characteristic matrices of -complexes, was introduced in [Dicarlo:2014:TNL:2543138.2543294]. It may be useful to understand that the chain of (co)boundaries is an algebraic representation of the Hasse diagram of a cellular complex [10.1145/1236246.1236259, MILICCHIO2008172]. The main advantage of the matrix approach is that topology, geometry and physics may coexist in a unique sparse matrix framework, concurring together to define, represent and simulate the behavior of a model, without any restrictions on type, dimension, codimension, orientability, manifoldness, or connectedness [10.1145/1236246.1236259].
For the sake of reader convenience, the theoretical minimum about linear spaces of chains and their linear operators is given in Section 2. Section 3 introduces the concepts of -nearness of points and -congruence of cells in order to define the quotient topology. In Section 4.1 we introduce the block-diagonal sparse matrices which set the scene for the Cell Congruence Enabling algorithm, whose Julia code is discussed in Section 4. In Section 5 both an elegant implementation using native Julia syntax for sparse matrix and vector algebra is given, and its translation to fast GraphBLAS primitives [graphblas:20] is provided. The concluding section supplies some hints about possible uses of this approach.
2 Chain spaces and chain complex
All computer applications including geometric content may require some computation of topology, i.e., of the relations of incidence or adjacency between cells (vertices, edges, faces, volumes), i.e, between 0-, 1-,2-, and 3-cells of a cellular decomposition of either the input objects or their boundaries. Cells may be represented as basis elements in a graded vector space of chains, closed with respect to addition of cells with the same dimension and w.r.t. the product times a scalar element in a field [ieee-tase, Arnold:2018].
Given a finite collection S of geometric objects in , the arrangement is the decomposition of into connected open cells of dimensions induced by . We remark that the topology of a space partition is fully described by a chain complex, i.e. by a sequence of chain spaces of different dimension (), and by linear operators of boundary and coboundary between chain spaces. In particular, the union of bases of chain spaces provide a minimal set of generators of the induced topology. Once chosen an ordering of the cells, providing independent chain generators, such linear operators are described by their matrices, holding the scalar values which supply the coordinates with respect to the bases, that are normally either in (non-oriented chains) or in (oriented chains).
An important problem in computational geometry is to asses the arrangement [Halperin:2017, 2019arXiv191108130P] of Euclidean 3-space, i.e., the partition of produced by a collection of cellular 3-complexes. The simplest computational topology solution is provided by the computation of the chain complex
of the space partition generated by the input collection of geometric objects.
The chain of linear boundary operators is sufficient to characterize the linear chain spaces, since they contain the domain bases by columns, or better, they contain the coordinate vectors representing them as linear combinations of the target space bases. Furthermore, we have by duality, for every . The computation of the matrix, via the TGW (Topological Gift Wrapping) algorithm [2017arXiv170400142P], requires after the geometrical intersection of input 2-cells. In the present manuscript, we discuss how to compute the matrices , and , via elementary algebraic operations, implemented by using the standard GraphBLAS library for very large graphs.
3 Quotient topology
We introduce in this section a simple way to glue together a set of local geometric and topological information generated independently from each input 2-cell, by giving few definitions of nearness and congruence, which will allow to transform the union of local topologies into a single global “quotient topology”.
Let be a small positive value. Given a set of points, we can partition it into subsets such that where is a representative of , being the Euclidean distance. Such partition of a point set, characterized by a clustering of subsets of close elements, can be constructed using the inrange function of the NearestNeighbors.jl package, which implements the -d tree data structure [10.1145/361002.361007]. See the Section 4 for details.
We say that two points are -near, and write , when their Euclidean distance is .
The -nearness is an equivalence relation, since it is reflexive, symmetric, and transitive. In particular, it is transitive since any pair of points are -near, because both have a distance less than from , and hence have a distance no more than from each other. More formally, if and , then , since the distance from of every point (e.g., ) in is less than .
Let be the quotient set of points w.r.t. the -nearness relation , for some fixed . The elements of can be chosen as the points representatives of the equivalence classes of . Even better, the representative of each class can be chosen as its centroid, and the average distance from of the other points of the class certainly decreases.
3.2 Topological congruence
Two geometric figures are congruent iff one can be transformed into the other by an isometry [Coxeter:1967], i.e. by an affine transformation with unit determinant.
We say that two -cells are -congruent, and write , when there exists a bijection between their 0-faces that pairwise maps vertices to -near vertices. The identification or quotient topology gives a method of getting a topology on from a topology on . The quotient topology is exactly the one that makes the resulting space ‘look like’ the original one, with the identified elements glued together.
A subset of open sets in the topological space over the cell complex is a basis for the topology of when all other open sets can be written as union of elements of . If we have an equivalence relation on a cell complex we get a natural projection map , by mapping each cell to its equivalence class. The identification or quotient topology on is defined as follows: a set is open if and only if is open in . In particular, the topology and the quotient topology are equivalent because any non-empty open set of contains a non-empty open set of and, conversely, every non-empty open set of contains non-empty open sets of .
It is easy to see that -congruence between elementary chains (aka cells) , denoted , is a graded equivalence relation, so that a chain complex may be represented by a much smaller , where , and where
4 Chain Complex Congruence
In this section we discuss the simple operations needed to transform the initial sparse block matrices and into the final operator matrices and concerning the space partition as a whole.
As we have seen, a chain -complex is completely specified by the sequence of its boundary or coboundary maps. Hence we intend to construct the coboundary maps induced by the congruence topology on the union set of local topological spaces generated by independent decompositions of input 2-cells, in the computational pipeline [2017arXiv170400142P] for the building of the space arrangement produced by a collection of cellular complexes.
4.1 Block-diagonal accumulators of coboundaries
We fix our attention here to and in particular to the computation of and . The matrix will be computed by Topological Gift Wrapping (TGW) algorithm, not described in this paper, in a following stage of the computational pipeline [2017arXiv170400142P], starting from .
In the algorithmic specification below, a dense array W and two sparse arrays Delta_0, Delta_1 respectively provide the input vertex coordinates and the sparse block-matrices and .
After independent computation, possibly in parallel, of local chain complexes (see [2017arXiv170400142P]), both accumulator matrices and have a sparse block-diagonal structure, made by two nested levels of (sparse) diagonal blocks. Each outer block concerns one of the input geometric objects, and inner blocks, , store the matrices of each decomposed 2-cell. We call , the exterior blocks (light gray), and , the interior blocks (dark gray), where is the number of 2-cells in -th input geometric object.
4.2 Nearness of vertices
First we need to discover the -nearness of vertices on the input matrix V of 3D coordinates, by querying the kdtree data structure [10.1145/361002.361007], and getting classes of -congruent vertices, represented by the array of arrays vclasses. The output matrix W holds by columns the coordinates of identified vertices in each -congruent class, that are mapped to their centroids.
4.3 Chain Complex Congruence Algorithm
In Section 4.1 we have discussed the block diagonal marshaling and of local coboundary matrices. With the function cellcongruence we replace each subset of columns of Delta sparse matrix corresponding to -near vertices with their vector sum. In this way we produce a new matrix from an array of new vectors. Finally, equal rows of this new matrix, discovered by a dictionary, are substituted by a single representative.
The same function cellcongruence is also applied to , by summing each subset of columns corresponding to each class of congruent edges, so generating a new sparse matrix from the resulting set of columns. Then, we reduce every subset of equal rows, if any, to a single row representative of congruent faces. A stepwise description of the Chain Complex Congruence (CCC) algorithm is given below. The function cellcongruence takes Delta of type SparseMatrixCSC and the inclasses given as array of arrays, and computes the set outclasses of congruent elementary -chains, represented as arrays of -chains.
Lar.cop2lar transforms the sparse matrix into a cellarray (array of arrays of cell indices); then a vector newcell is defined to accomodate a map between old facet indices into the new ones of their quotient projection . Finally, the array cells accomodates the translated representations (arrays) of input (-1)-cells corresponding to columns of ;
the empty cells, if any, are removed from okcells. Empty cells appear when they contain a number of distinct facets less than dim, i.e., for edges, and for faces;
each -cell (elementary -chain) is represented as an array of indices of its facets (-1)-chain). The selection of single instances of output cells (from okcells, where they may be repeated) is performed using a DefaultOrderedDict for congruence classes, with key given by the cell representation (array of facets). The dictionary reading action classes[face] returns  iff the key (face) is not contained in the dictionary. In this case the key is stored, and its storage of value starts with the first index [k]; otherwise, the current index is appended to the yet incomplete value.
At the end, both the dictionary keys (column indices of congruent rows, i.e. the set of non-zero positions (for the representative row) in the output , and the dictionary values, i.e., the projected elements (class representatives) are given as output.
4.4 Top-level interface
Finally, the higher-level function chaincongruence maps the input data W, Delta_0, Delta_1, into a compact representation V, EV, FE of the chain complex , , and (see Appendix A.1). For full generality, also the 2-cell congruence was checked and taken into account. Congruent faces would in fact appear when some input 2-cells of the arrangement pipeline have inner intersection and lay on the same 3D plane. We would like to remark that decomposed 1-cells (i.e., input edges) and the representatives of their congruence classes (i.e., output edges) are associated one-to-one with the rows of and the columns of , respectively, so satisfying the topological constraints , that were checked positively in all our tests.
5 Sparse matrix implementations
Two compact matrix implementations of Chain Complex Congruence algorithm are given in this section, starting from native coding in Julia, with its concise matrix syntax, then using the Julia porting of the package SuiteSparse:GraphBLAS
, enriched by a simple user-friendly interface from Julia’s multiple dispatch. The Julia code for the three implementations of the CCC algorithm is available at open-source repositoryhttps://github.com/cvdlab/LocalCongruence.jl. The three implementations are denoted as AA (array of arrays), SP (Julia sparse array), and GB (GraphBLAS), respectively. The AA code was discussed in Section 4.
5.1 Julia’s SparseArrays.jl implementation
The vertexCongruence evaluates the Vertex Congruence for 3D-points. The function determines the points of matrix “V“ closer than “err“ to each other, and builds a new Vertex Set made of the representative of each point cluster. The method returns: the new Vertex Set; a map that, for every new vertex, identifies the subset of old vertices it is made of. The function cellCongruence evaluates the Cell Congruence for a -Cochain “cop“ of type sparse array, with equivalence classes “lo_cls“, and corresponding “lo_sign“, for -cells, both given as array of arrays of (old) lower-rank indices, to denote equivalence classes. The chainCongruence function performs the Geometry “G“ congruence and reshapes the Topology “T“.
5.2 Mapping to SuiteSparse:GraphBLAS.jl
6 Analysis of results
We have introduced here three versions of our CCC algorithm. A comparison with similar or different methods is clearly not possible, by the novelty of our approach, but a mutual comparison is interesting. The first algorithm, in Section 4, introduced the sparse block-matrix reduction method by using arrays of arrays, but with the drawback of losing the ordering information (and the sparse matrix output) of the cells, which were available in the input. The second version is given in Section 5.1 by directly using the Julia’s native sparse matrices. The third implementation, in Section 5.2, makes use of the GraphBLAS primitives within a tiny Julia’s wrapping.
Below you may see, using mean times, and assuming AA = 1x (Native Julia’s array of array — last column), that we have: native Julia (first column) SparseArrays = 16.6x; Julia wrapping (second column) of SuiteSparse:GraphBLAS = 6.12x.
All implementations are quite naive. No optimizations were performed. We believe that major benefits may be ported by optimization to the two versions that use sparse matrices. A great caveat for the array of arrays version is that it looses track of cell signs, that must be recovered later with further computations, so losing all time benefits. Anyway, it let us suppose that the GraphBLAS approach can be forther optimized.
The two implementations with sparse matrices allow for maintaining the 2-cell orientation within the global chain complex. We have seen that the GraphBLAS implementation has two important benefits: (1) with parity of storage occupation, it allows to compute the largest chain complexes; (2) it also allows for block decomposed matrix computations, so opening the way to hierarchical computation of huge topological operators. Finally, we remark that the CCC algorithm is multidimensional, and its implementation may be extended, as is, to manage higher dimensional complexes, just by adding more cellcongruence call instances to chaincongruence function on page 4.4.
In this paper we have discussed an efficient algebraic algorithm to create the and sparse coboundary matrices, encoding the topological congruences between a set of chain complexes. In the algorithmic pipeline to construct the arrangement of Euclidean 3-space [2019arXiv191108130P], local chain complexes are generated independently from single fragmented input 2-cells, starting from a collection of cellular 2- or 3-complexes in 3D. The correctness of results, that might depend on numeric approximations of floating-point arithmetics, is checked by testing the matrix constraint . Therefore, we have given one possible topological solution to the robustness problem in geometric computations, in the line discussed by Christoph Hoffmann in [Hoffmann:2001], using epsilon geometry [10.1145/73833.73857]. The chain complex congruence (CCC) enabling algorithm introduced here was inplemented in Julia using the package SuiteSparseGraphBLAS.jl, and it is a key component of a computational pipeline to produce solid models of complex geometric scenes, using robust Boolean algebra methods for next-generation image understanding.
Appendix A Appendix
a.1 A simple example
For the sake of reproducibility, as well as of reader convenience and understanding, a very simple example of the congruence reduction for a 3D cube is given here, by using the initial implementation of Sections 4.2 and 4.3. The open-source realization with Julia sparse matrices and/or with GraphBLAS may by downloaded from GitHub at https://github.com/cvdlab/LocalCongruence.jl.
The cube was generated with random size and random attitude produced by a random 3D rotation, in order to start from close but unequal coordinates for each single instance of face vertices. All faces were generated independently, in order to get columns (vertex instances) in the input matrix W. The Julia’s matrix input follows (as matrix, in order to show the full Float64 digital resolution).
The two corresponding sparse block-diagonal matrices Delta_0 and Delta_1 are given below, as array triples (I,J,X) of row and column indices, and non-zero values.
The execution command from Julia terminal is given below, followed by output showing, with V output by columns, and both EV (edges by vertex indices) and FE (faces by edge indices) as array of arrays: