The Implementation of the Colored Abstract Simplicial Complex and its Application to Mesh Generation

07/04/2018 ∙ by C. T. Lee, et al. ∙ University of California, San Diego 0

We introduce CASC: a new, modern, and header-only C++ library which provides a data structure to represent arbitrary dimension abstract simplicial complexes (ASC) with user-defined classes stored directly on the simplices at each dimension. This is accomplished by using the latest C++ language features including variadic template parameters introduced in C++11 and automatic function return type deduction from C++14. Effectively CASC decouples the representation of the topology from the interactions of user data. We present the innovations and design principles of the data structure and related algorithms. This includes a meta-data aware decimation algorithm which is general for collapsing simplices of any dimension. We also present an example application of this library to represent an orientable surface mesh.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

For problems in computational topology and geometry, it is often beneficial to use simple building blocks to represent complicated shapes. A popular block is the simplex, or the generalization of a triangle in any dimension. Due to the ease of manipulation and the coplanar property of triangles, triangulations have become commonplace in fields such as geometric modeling and visualization as well as topological analysis. Discretizations are also used for efficient solving of Partial Differential Equations (PDE). The use of meshes has become increasingly popular even in the fields of computational biology and medicine

(Zhang, 2016).

As methods in structural biology improve and new datasets become available, there is interest in integrating experimental and structural data to build new predictive computer models(Roberts, 2014). A key barrier that modelers face is the generation of multi-scale, computable, geometric models from noisy datasets such as those from Electron Tomography (ET)(Yu et al., 2008b). This is typically achieved in at least two steps: (1) segmentation of relevant features, and (2) approximation of the geometry using meshes. Subsequently, numerical techniques such as Finite Elements Modeling or Monte Carlo can be used to investigate the transport and localization of molecules of interest.

While many have studied mesh generation in the fields of engineering and animation, few methods are suitable for biological datasets. This is largely due to noise introduced by limits in image resolution or contrast. Even while using state-of-the-art segmentation algorithms for ET datasets there are often unresolved or missed features. Due to these issues, the generated meshes often have holes and other non-manifolds which must be resolved prior to mathematical modeling. Another challenge is the interpretation of a voxel valued segmentation. The conversion of zig-zag boundaries into a mesh can lead to other problems such as extremely high aspect ratio triangles or, in general, poorly conditioned elements(Yu et al., 2008b). To remedy this, various smoothing and decimation algorithms must also be applied prior to simulation.

Previous work by us and others have introduced a meshing tool for biological models, GAMer, for building 3D tetrahedral meshes which obey internal and external constraints, such as matching embedding and/or enclosing molecular surfaces. It also provides the ability to use various mesh improvement algorithms for volume and surface meshes(Gao et al., 2013; Yu et al., 2008a). GAMer uses the Tetgen library as the primary tetrahedral volume generator(Si, 2015). While the algorithms are sound, the specific implementation is prone to segmentation faults even for simple meshes. Careful analysis of the code has identified that the data structures used for the representation of the mesh is primarily at fault. This article will focus entirely on the representation of topology in very complex mesh generation codes. We note that the algorithms which handle geometric issues like shape regularity and local adaptivity are well understood(Babuška and Aziz, 1976; Liu and Joe, 1994), among others. Similarly there is a large body of literature related to local mesh refinement and decimation(Bank et al., 1983; Bank and Xu, 1996). Our innovations serve to enable the implementation of these algorithms in the most general and robust way.

GAMer currently employs a neighbor list data structure which tracks the adjacency and orientation of simplices. Neighbor lists are quick to construct, however the representation of non-manifolds often leads to code instability. Algorithms must check for aberrant cases creating substantial overhead. We note that while the need to gracefully represent 2D and 3D non-manifolds for ET applications drove our initial focus, we are also interested in mesh generation in higher dimensions with applications to: numerical general relativity (3D+1)(Holst et al., 2016; Regge, 1961), computational geometric analysis (nD)(Holst and Tiee, 2018), phase space simulations (6D), and arbitrary collective variable spaces in molecular modeling for enhanced sampling(Vanden-Eijnden and Venturoli, 2009). We therefore chose following requirements for a mesh data structure to serve as design goals:

  • General and capable of representing non-manifold, mixed dimensional, oriented and non-oriented meshes in arbitrary dimensions.

  • Support for inline and flexible data storage. In some applications, data must be associated with the topology. For example, problems in general relativity typically require the storage of metric tensors on all simplex dimensions.

  • Support for intuitive and simple manipulations and traversals.

Here we describe the development of a scalable colored abstract simplicial complex data structure called CASC. Simplices are stored as nodes on a Hasse diagram. For ease of traversal all adjacency is stored at the node level. An additional data object can be stored at each node which is typed according to the simplex dimension at compile time. This means that, for example, for a mesh the -simplices can be assigned a vertex type while the -simplices can store some material property instead. Typing of each -simplex is achieved using variadic templates introduced in C++11. CASC thus provides a natural separation between the combinatorics represented by the ASC from the underlying data types at each simplex dimension and their interactions. In §2 we briefly define an ASC and some relevant definitions followed by the introduction of the CASC data structure and it’s construction in §3. We then demonstrate the use of CASC to represent a surface mesh and compute vertex tangents in §5.

1.1. Related Work

Although many data structures to represent simplicial complexes have been developed, to the best of our knowledge there currently exists no data structure which supports meshes of arbitrary dimension with user-selected typed data stored directly on each simplex. A full review of all existing data structures is beyond the scope of this work, however we highlight several representative examples. Many data structures such as the half-edge and doubly-connected edge list among others are restricted to the representation of two-manifolds only(De Floriani and Hui, 2005). Other data structures such as SIG(De Floriani et al., 2004), IS(De Floriani et al., 2010), IA*(Canino et al., 2011), SimplexTree(Boissonnat and Maria, 2014), AHF(Dyedov et al., 2015; Ray et al., 2015), LinearCellComplex, and dD Triangulations(Boissonnat et al., 2009) support or can be extended to represent arbitrary dimensional simplicial complexes. However, their current implementations either do not consider the storage of data beyond possibly embedding, or do not support inline storage of user data. LinearCellComplex from CGAL supports only a linear geometrical embedding(CGA, [n. d.]). AHF implemented in MOAB uses separate arrays of data which are then referenced using a handle(Dyedov et al., 2015; Ray et al., 2015). In addition to the limitations of data storage, some make assumptions limiting their generality. dD Triangulations, for example, assumes that a simplicial complex is pure and therefore does not support the representation of mixed dimensional complexes(Boissonnat et al., 2009).

2. Background – Abstract Simplicial Complexes

An Abstract Simplicial Complex (ASC) is a combinatorial structure which can be used to represent the connectivity of a simplicial mesh, independent of any geometric information. More formally, the definition of an ASC is as follows.

Definition 2.1 ().

Given a vertex set , an abstract simplicial complex of is a set of subsets of with the following property: for every set , every subset is also a member of .

The sets are called a simplex or face of ; similarly a face is said to be a face of simplex if . Since is a face of , is a coface of . Each simplex has a dimension characterized by , where is the cardinality of set . A simplex of is also called a -simplex. The dimension of the complex, , is defined by the largest dimension of any member face. Simplices of the largest dimension, are referred to as the facets of the complex.

If one simplex is a face of another, they are incident. Every face of a -simplex with dimension is called a boundary face while each incident face with dimension is a coboundary face. Two -simplices, and are considered adjacent if they share a common boundary face, or coboundary face. The boundary of simplex , , is the sum of the boundary faces.

Having introduced the concept of an ASC, we can also define several operations useful when dealing with ASCs. A subcomplex is a subset that is a simplicial complex itself. The Closure () of a simplex, , or some set of simplices is the smallest simplicial subcomplex of that contains :

(1)

It is often useful to consider the local neighborhood of a simplex. The Star () of a simplex is the set of all simplices that contain :

(2)

The Link () of consists of all faces of simplices in the closed star of that do not intersect :

(3)

For some algorithms, it is often useful to iterate over the set of all vertices or edges etc. We use the following notation for the horizontal “level” of an abstract simplicial complex.

(4)

A subcomplex which contains all simplices where is the -skeleton of :

(5)
Figure 1. Hasse diagrams of several Abstract Simplicial Complexes and a geometric realization from left to right: the empty set, a vertex, an edge, a triangle, a tetrahedron.

By Definition 2.1, an ASC forms a partially ordered set, or poset. Posets are frequently represented by a Hasse diagram, a directed acyclic graph, where nodes represent sets, and edges denote set membership. Several example simplicial complexes and their corresponding Hasse diagrams are shown in Fig. 1. Colloquially we will use up and down to refer to the boundary and coboundary of a simplex respectively. In Hasse diagrams, we follow a convention that simplices shown graphically on the same horizontal level have the same simplex dimension. Furthermore, simplices of greater dimension are drawn above lesser simplices.

Figure 2. Data structure diagram of a triangle represented using CASC. Each simplex is represented as a node containing a dictionary up and/or down which maps the vertex index to a pointer to the next simplex. Data can be stored at each node with type determined at compile time. Effectively each level can contain different metadata as defined by the user, separating the interactions of user data from the representation of topology.

3. Colored Abstract Simplicial Complex

In this section we introduce the CASC data structure and its implementation. For a given simplicial complex, each simplex is represented by a node (asc_Node) in the Hasse diagram, and defined by a set of keys corresponding to the vertices which comprise the simplex. Note that we use node to refer to objects in CASC Hasse diagram and not -simplices. Instead, -simplex are referred to as the vertices of the mesh. Furthermore we refer to the -simplex or -simplex as the root simplex interchangeably. When a node is instantiated, we assign it a unique Integer Internal Identifier (iID) for use in the development of CASC algorithms. The iID is constant and never exposed to the end-user except for debugging purposes. Instead nodes can be referenced by the user using the SimplexID which acts as a convenience wrapper around an asc_Node*, providing additional support for move semantics for fast data access. All topological relations (i.e., edges of the Hasse diagram) are stored in each node as a dictionary which maps user specified keys to SimplexIDs up and down. An example data structure diagram of triangle {1,2,3} is shown in Fig. 2. Based upon this example, if a user has the SimplexID of -simplex and wishes to get -simplex , they can look in the Up dictionary of SimplexID for key which maps to a SimplexID. The vertices which constitute each simplex are not stored directly, but can be accessed by aggregating all keys in Down.

We note that while the representation of all topological relations is redundant and may not be memory optimal, it vastly simplifies the traversals across the complex. Furthermore, the associate algorithms and innovations using variable typing are general and thus compatible with other more condensed representations.

3.1. Variable Typing Per Simplex Dimension

We achieve coloring by allowing user-defined data to be stored at each node. The typical challenge for strongly typed languages such as C++ is that the types must be defined at compile time. Typical implementations would either hard code the type to be stored at each level or use a runtime generic type, such as void*. However, each of these have drawbacks. For the former, this requires writing a new node data structure for every simplicial complex we may wish to represent. For the latter, using void* adds an extra pointer dereference which defeats cache locality and may lead to code instability. Another possible implementation might be to require users of the library to derive their data types from a common class through inheritance. This solution puts an unnecessary burden on users who may have preexisting class libraries, or simply wish to store a built in type, such as an int. Furthermore, under the inheritance scheme, changes to the underlying container may require users to update their derived classes. To avoid this cumbersome step, we have employed the use of variadic templates introduced in C++11 to allow for unpacking and assignment of data types. The user specifies the types to be stored at each level in a list of templates to the object constructor, see Fig. 3.

Figure 3. Template arguments for ASC are unpacked and assigned to nodes accordingly. The first argument “int” is the key type for labeling vertices while the following arguments define node data types. Notably, Node<2> does not allocate memory for data as the corresponding template argument was “void”. By passing in additional types, simplicial complexes of higher dimensions are instantiated.

The variadic templating allows CASC to represent complexes of any user-defined dimension. To specify an -simplicial complex, CASC requires the definition of an index/key type followed by data types and edge types. The first data type provided after the key type corresponds to data stored on the -simplex which can be thought of as global metadata. For example, suppose we have a 2-simplicial complex intended for visualization and wish to store locations of vertices and colors of faces. A suitable ASC can be constructed using the following template command:

auto mesh = AbstractSimplicialComplex<int, void, Vertex, void, Color>()

If we now wish to represent a tetrahedral mesh, instead of constructing a new data structure, we can simply adjust the command:

auto mesh = AbstractSimplicialComplex<int, void, Vertex, void, Color, Density>()

In both cases, the first template argument is the key type for referring to vertices followed by the data type for each -simplex. Supposing that the user does not wish to store data on any given level, by passing “void” as the template argument, the compiler will optimize the node data type and no memory will be allocated to store data. In both cases, the 0- and -simplices will have no data.

By using variadic templates, we allow the user to specify both the dimension of the simplicial complex as well as the types stored at each level. Because the type deduction is performed at compile time, there is no runtime performance impact on user codes. There is, however, some additional code complexity introduced. If the user wishes to retrieve the data stored in a simplicial complex, they must know what level they are accessing at compile time. A consequence is that the exposed identifier object, SimplexID, is templated on the integral level, so that types can be deduced. This does not present a problem for simple use cases, such as:

// Create alias for 2-simplicial complex
using ASC = AbstractSimplicialComplex<int,int,int,int>;
ASC mesh = ASC();  // construct the mesh object
mesh.insert<2>({1,2},5); // insert edge {1,2} with data 5
ASC::SimplexID<2> s = mesh.get_simplex<2>({1,2}); // get the edge
// NOTE: the type is templated on the integral level
std::cout << *s << std::endl; // prints data "5"

However, when implementing algorithms intended to be generic on any simplicial complex, templated code must be written. We discuss the implementation of several such algorithms in the following section.

4. Implemented Algorithms

The following algorithms are provided with the CASC library:

  • Basic Operations

    • Creating and deleting simplices (insert/remove)

    • Searching and traversing topological relations (GetSimplexUp/GetSimplexDown)

  • Traversals

    • By level (get_level)

    • By adjacency (neighbors_up/neighbors_down)

    • Traversals across multiple node types (Visitor Design Pattern/double dispatch)

  • Complex Operations

    • Star/Closure/Link

    • Metadata aware decimation

Input: []: Indices of simplices to describe new simplex
: The simplex to insert relative to (most commonly )
Output: The new simplex
Function  insert(keys[n], rootSimplex)
      for (i = 0; i ¡ n; i++)
           newSimplex = createSimplex() if (i ¿ 0) /* Recurse to insert sub-simplices */
                 /* Pass only the first part of up to index */
                return insert(keys[0:i], newSimplex)
                else /* Terminal conditional */
                     return newSimplex
                    
                    
ALGORITHM 1 Insertion of a new simplex
Figure 4. Recursive insertion of tetrahedral simplex {1,2,3,4}. The order of node insertions is represented by the numbered red arrows. When each node is created, the black arrows to parent simplices are created by backfilling.

4.1. Basic Operations

4.1.1. Creating and deleting simplices

Since the CASC data structure maintains every simplex in the complex and all topological relations, inserting a -simplex, , into the complex means ensuring the existence of, and possibly creating, nodes and edges (see derivation in SI). Fortunately, the combinatorial nature of simplicial complexes allows this to be performed recursively. A generalized recursive insertion operation for any dimensional complex and user specified types, is described in Algorithm 1. The insertion algorithm defines an insertion order such that all dependent simplices exist in the complex prior to the insertion of the next simplex.

As an illustrative example of the template code used in this library, Algorithm 1 is rewritten in C++ template function-like pseudocode shown in Algorithm 6. While the templated code is more complicated, it provides many optimizations. For example, since the looping and recursion are performed at compile time, for any -simplex we wish to insert, any modern compiler should optimize the code into a series of insertNode() calls; the setupForLoop() and forLoop() function calls can be completely eliminated. As a result, the optimized templated code will exhibit superior run time performance. To illustrate the insertion operation, a graphical representation of inserting tetrahedron {1,2,3,4} by Algorithm 6 is shown in Fig. 4, and step-by-step in Fig. S6. In the example, new simplex is added sequentially to the complex and any missing topological relations are found by traversing the faces of and backfilling.

The removal of any simplex is also performed using a recursive template function. When removing simplex , in order to maintain the property of being a simplicial complex, all cofaces of or must also be removed along with any boundary/coboundary relations of . The implemented removal algorithm traverses up the complex and removes simplex and all cofaces of level by level.

Input: : the simplex to start the search from.
[]: the relative name of the desired simplex
Output: The simplex
Function  getSimplexUp(root, keys[n], index)
      if (index ¡ n)
           Find keys[index] in root.up /* Look for the edge up */ if (Edge up is found)
                return getSimplexUp(, keys, index+1)
                else /* Edge keys[index] doesn’t exist */
                     return null pointer
                    
                     else {  return root  /* Terminal condition */
ALGORITHM 2 Searching for a new simplex

4.1.2. Searching and traversing topological relations

The algorithms for retrieving a simplex as well as for basic traversals from one simplex to another across the data structure are the same. Given a starting simplex, and an array of keys up, the new simplex can be found recursively by Algorithm 2; The annotated code used is shown in §B.3. Traversals from one simplex to another require a key lookup followed by a pointer dereference and therefore occur in approximately constant time (). Since all topological relations are stored, the traversal order across the array of keys does not matter. The same algorithm can be applied going down in dimension. For the retrieval of an arbitrary simplex, we start the search up from the root node of the complex.

4.2. Traversals

Thus far, we have presented algorithms for the creation of a simplicial complex as well as the basic traversal across faces and cofaces. For many applications, other traversals, such as by adjacency, may be more useful. We present several built-in traversal algorithms as well as the visitor design pattern for complicated operations.

4.2.1. By level

It is often useful to have a traversal over all simplices of the same level. For example, iterating across all vertices to compute a center of mass. To support this in an efficient fashion, simplices of the same dimension are stored in a level specific map of iIDs to node pointers. Notably, the map for each level is instantiated with the correct user specified node type with respect to level at compile time. To achieve this, we again use variadic templates to generate a tuple of maps, where each tuple element corresponds to the map for a specific level’s node type.

Since asc_nodes are templated on the integral level, we can use a template type map to map an integral sequence to the node pointer type,

producing a tuple of integrally typed node pointers. Subsequently, we can map again to generate a tuple of maps,

By using this variadic template mapping strategy we now have the correct typenames assigned. Any level of the tuple can be accessed by getting the integral level using functions in the C++ standard library. Variations of this mapping strategy are also used to construct the SimplexSet and SimplexMap structures below.

For end users, the implementation details are entirely abstracted away. Continuing from the example above, iteration over all vertices of simplicial complex, mesh, can be performed using the provided iterator adaptors as follows.

// Deduces the type of nid = ASC::SimplexID<1>
for(auto nid : mesh.get_level_id<1>()){
  std::cout << nid << std::endl;
}
Listing 1: Example use of iterator adaptors for traversal across vertices of mesh.

The function get_level_id<k>() retrieves level from the tuple and returns an iterable range across the corresponding map.

Input: : The simplex to get the neighbors of.
Output: List of neighbors
Function  getNeighborsUp(s) /* Get the neighbors */
      Create an empty list of /* Follow all coboundary relations up from s. */
      for ()
           SimplexID id = for () /* Go down from id */
                if () {  Add to  /* Do not add self to neighbors */
               
                return neighbors
ALGORITHM 3 Get the neighbors of a simplex, , by inspecting the faces of

4.2.2. By adjacency

Many geometric algorithms operate on the local neighborhood of a given simplex. Unlike other data structures such as the halfedge, CASC does not store the notion of the next simplex. Instead, adjacency is identified by searching for simplices with shared faces or cofaces in the complex. The algorithm for finding neighbors with shared faces is shown in Algorithm 3. We note that the set of simplices with shared faces may be different than the set of simplices with shared cofaces. Both adjacency definitions have been implemented and we leave it to the end user to select the relevant function. Once a neighbor list has been aggregated, it can be traversed using standard methods. While the additional adjacency lookup step is extra in comparison to other data structures, in many cases, the generation of neighbor lists need only be done once and cached. The trade off is that CASC offers facile manipulations of the topology without having to worry about reorganizing neighbor pointers.

4.2.3. Traversals over multiple node types.

When performing more complicated traversals, such as iterating over the star of a simplex, multiple node types may be encountered. In order to avoid typename comparison based branch statements, we have implemented visitor design pattern-based breadth first searches (BFS). The visitor design pattern refers to a double dispatch strategy where a traversal function takes a visitor functor which implements an overloaded visit() function. At each node visited, the traversal function will call visit() on the current node. Since the functor overloads visit() per node type, the compiler can deduce which visit function to call. Example pseudocode is shown in Listing 2. This double dispatch strategy, eliminates the need for extensive runtime typename comparisons, and enables easy traversals over multiple node types. We provide breadth first traversals up and down the complex from a set of simplices. These visitor traversals are used extensively in the complex operations described below.

template <typename Complex>
struct Visitor{
  template <std::size_t k>
  using Simplex = typename Complex::template SimplexID<k>;
  // General template prototype
  template <std::size_t level>
  bool visit(Complex& F, Simplex<level> s){
    return true;
  }
  // Specialization for vertices
  bool visit(Complex& F, Simplex<1> s){
    s.position = s.position*2;
    return true;
  }
  // Specialization for faces
  bool visit(Complex& F, Simplex<3> s){
    s.color = ’green’;
    return true;
  }
};
void BFS(ASC& mesh, Visitor&& v){
  // ... traversal code
  v.visit(mesh, currentSimplex);
  // NOTE: visit is overloaded and called based on function prototype.
}
void main(){
  // ... define simplicial complex traits for a surface mesh
  ASC mesh = ASC(); // construct the mesh object
  // ... insert some simplices etc.
  BFS(mesh, Visitor());
  }
}
Listing 2: Example pseudocode of double dispatch to traverse the complex while scaling the mesh by 2 and coloring the faces green.

4.3. Complex Operations

4.3.1. Star/Closure/Link

The star, link, and closure can be computed using the visitor breadth first traversals to collect simplices. These operations typically produce a set of simplices spanning multiple simplex dimensions, and thus simplex typenames, which cannot be stored in a traditional C++ set. We have implemented a multi-set data structure called the SimplexSet, which is effectively a tuple of typed sets corresponding to each level. The SimplexSet is constructed using the same mapping strategy as the tuple of maps used for the iteration across levels. For convenience, we provide functions for typical set operations such as insertion, removal, search, union, intersection, and difference. Using a combination of the star and closure functions with SimplexSet difference we can get the link by Eq. 3.

4.3.2. Metadata aware decimation

(a) A geometric realization of the complex before and after decimation.
(b) Explicitly drawn out Hasse diagrams for the constructed example where the TOP is before, and BOTTOM is after decimation. Grey arrows mark the relationships between sets of simplices before and after. Because there is always a mapping, users can define strategies to manage the stored data.
Figure 5. The example decimation of edge in a constructed example.

We have implemented a general decimation algorithm which operates by collapsing simplices of higher dimensions into a vertex. While edge collapses for 2-manifolds are well studied, a general dimensional collapse is useful for decimating higher-dimensional meshes used to solve PDEs such as those encountered in general relativity. Since simplices are being removed from the complex, user data may be lost. Our implementation is metadata aware and allows the user to specify what data to keep post decimation. This is achieved by using a recursive algorithm to produce a map of removed simplices to new simplices. The user can use this mapping to define a function which maps the original stored data to the post decimation topology. This decimation strategy is implemented as an inplace operation yielding decimated mesh containing data mapped according to a user specified callback function.

Input: : the simplicial complex;
: the simplex to decimate;
: a callback function to handle the data mapping
Output: Simplicial complex with simplex collapsed according to Def. 6
Simplex np = F.newVertex() /* Create a dummy new vertex to map to */ SimplexSet N /* For the complete neighborhood */ SimplexDataSet data /* Data structure to store (simplex name, simplex data) pairs */ SimplexMap simplexMap /* Data structure to store New Simplex -> SimplexSet map */  /* Get the complete neighborhood (all simplices which are associated with ). */
for (vertex of )
      BFSup (simplex )
           if () { N.insert (j)
          
            /* Backup the complete neighborhood. These simplices will be destroyed eventually. */
           SimplexSet doomed = N  /* Generate the before-after mapping */
           BFSdown (simplex ) /* MainVisitor */
                BFSup (simplex ) /* InnerVisitor */
                      /* maps to so we need to connect to instead */
                     Name newName = SimplexSet grab BFSdown (simplex ) /* GrabVisitor */
                           /* Grab dependent simplices which have not been grabbed yet. Grabbed simplices will map to simplex */
                          N.remove(k) grab.insert(k)
                          if () {  simplexMap.insert(pair(newName, grab)) else {  simplexMap[newName].insert(grab)
                         
                          for ()
                                /* Run the user’s callback to map simplex data. */
                               DataType mappedData = (*clbk)(F, name(j), newName, grabbed) data.insert({newName, mappedData})/* Insert a pair containing new simplex name and data */
                                /* Iterate over the complete neighborhood and remove simplices */
                               performRemoval(F, doomed)  /* Iterate over data and append mapped simplices and data */
performInsertion(F, data)
ALGORITHM 4 Decimate a simplex by collapsing it into a vertex.

This decimation algorithm is a generalization of an edge collapse operation to arbitrary dimensions. It is formally defined as follows:

Definition 4.1 ().

Given simplicial complex , simplex to decimate , vertex set of , and new vertex , we define function,

(6)

where is any simplex . We define the decimation of by replacing with as .

Note that decimation under this definition is not guaranteed to preserve the topology, as can by seen by decimating any edge or face of a tetrahedron.

Decimation of a simplicial complex must result in a valid triangulation. Here we show that decimation by Definition 6 produces a valid abstract simplicial complex.

Theorem 4.2 ().

is an abstract simplicial complex.

Proof.

Given simplices and , let and let . We will show that . There are two cases for , as they can be disjoint or intersecting.

Considering the disjoint case where . This implies that and . Since is a simplicial complex, implies that . Furthermore since and , then implying that and thus .

Alternately in the intersecting case where , that is . Since there are two sub-cases where either contains or does not.

Supposing that . Then implying that and . Therefore by the disjoint case, implies that .

Supposing that . We can rewrite any such simplex as where . Furthermore we can write that where and and thus . There exists some set such that and such that Since then . Therefore and thus .

For all cases and sub-cases we have shown that therefore is an abstract simplicial complex. ∎

Order MainVisitor InnerVisitor GrabVisitor Maps to
1 {3,4} {3,4} {3,4}, {3}, {4} {6}
2 {3,4} {1,3,4} {1,3,4}, {1,3}, {1,4} {1,6}
3 {3,4} {3,4,5} {3,4,5}, {3,5}, {4,5} {3,6}
4 {3} {0,3} {0,3} {0,6}
5 {3} {0,1,3} {0,1,3} {0,1,6}
6 {3} {0,3,5} {0,3,5} {0,5,6}
7 {4} {2,4} {3,5} {5,6}
8 {4} {1,2,4} {1,2,4} {1,2,6}
9 {4} {2,4,5} {2,4,5} {2,5,6}
Table 1. Traversal order of the visitors for the decimation shown in Figure 4(a).

A pseudocode implementation for this decimation is provided in Algorithm 4. Given some simplicial complex and simplex to decimate, this algorithm works in four steps. First, we compute the complete neighborhood, , of . Simplices not in the complete neighborhood will be invariant under and are ignored. Next, we use a nested set of breadth first searches to walk over the complete neighborhood and compute for each simplex in the neighborhood. The results are inserted into a SimplexMap which maps to a SimplexSet of all which map to . Third, we iterate over the SimplexMap and run the user defined callback on each mapping to generate a list of new simplices and associated mapped data stored in SimplexDataSet. Finally, the algorithm removes all simplices in the complete neighborhood and inserts the new mapped simplices.

An example application of this decimation operation is shown in Figure 5. In Figure 4(a) we show the geometric realization of the complex before and after the decimation of simplex {3,4}. In Figure 4(b) we show the detailed Hasse diagrams for the constructed example. Note that there are two possible mapping situations. In one case, , groups of simplices are merged. In the other case, simplices only need to be reconnected to the new merged simplices. By carefully choosing the traversal order, some optimizations can be made.

We apply the decimation on the constructed example shown in Figure 5 and show the order of operations with respect to the current visited simplex for each visitor function in Table 1. Starting out, MainVisitor and InnerVisitor will visit {3,4}. At this point, GrabVisitor will search BFS down from {3,4} to grab the set {{3,4}, {3}, {4}} and remove it from the neighborhood, eliminating some future calculations at {3} and {4}. All simplices in this set will map to new simplex {6} after decimation. Continuing upwards, InnerVisitor will find simplex {1,3,4} and GrabVisitor will grab set {{1,3,4}, {1,3}, {1,4}}. Again this set of simplices will map to common simplex {1,6} post decimation. A similar case occurs with simplex {3,4,5}. At this point, all simplices in have been visited and removed from the neighborhood and MainVisitor continues BFS down and finds {3} and calls BFS up (InnerVisitor). Note that since simplex {3} has already been grabbed, InnerVisitor will continue upwards and find {0,3}. Looking down there are no simplices which are faces of {0,3} in the neighborhood. So on and so forth.

To reiterate, GrabVisitor grabs the set of simplices which will be mapped to a common simplex. We show here that the order in which simplices are grabbed by Algorithm 4 will preserve that all simplices where will map to .

When visiting any simplex where and corresponds to simplices visited by MainVisitor. We can write as where the sets and are disjoint. Looking down from all simplices fall into two cases: where or where . All simplices of form , at worst case, will been grabbed while InnerVisitor proceeded BFS up from . Remaining simplices can be grouped with and correctly mapped to .

We note that in some non-manifold cases GrabVisitor will not always grab set members in one visit. Supposing that we removed simplex {1,3,4} from the constructed example, in this case, InnerVisitor cannot visit {1,3,4} and simplices {1,3} and {1,4} will not be grouped. Instead {1,3} and {1,4} will be found individually when MainVisitor visits {3} then {4}. To catch this case and correctly map {1,3} and {1,4} to {1,6}, we use a SimplexMap to aggregate all maps prior to proceeding. We note that in all cases starting with a valid simplicial complex, this implementation of the general collapse of simplex visits each member in and maps according to Def. 6 producing a valid simplicial complex. There is no guarantee that the result will have the same topological type as the pre-decimated mesh. The preservation of the topological type under decimation is often a desirable trait. We will show how to verify the Link Condition for when edge collapse of a 2-manifold will preserve the topological type in §5.1.

5. Surface Mesh Application Example

CASC is a general simplicial complex data structure which is suitable for use in mesh manipulation and processing. For example, we can use CASC as the underlying representation for an orientable surface mesh. Using a predefined Vertex class which is wrapped around a tensor library, and a class Orientable which wraps an integer, we can easily create a surface mesh embedded in .

using Vector = tensor<double,3,1>;
struct Vertex {
  Vector position;
  // ... other helpful vertex functions;
};
struct Orientable {
  int orientation;
};
struct complex_traits
{
    using KeyType = int;
    using NodeTypes = util::type_holder<void,Vertex,void,Orientable>;
    using EdgeTypes = util::type_holder<Orientable,Orientable,Orientable>;
};
using SurfaceMesh = simplicial_complex<complex_traits>;

In this case, -simplices will store a Vertex type while faces and all edges will store Orientable. Using SurfaceMesh we can easily create functions to load or write common mesh filetypes such as OFF as shown in the included library examples.

Input: Simplices and where and is a vertex.
Output: The orientation of edge
int orient = 1 for (Vertex u : ) /* For each vertex in simplex */
      if (v¿a) {  orient *= -1;
return orient;
ALGORITHM 5 Define the orientation of a topological relation.

We can define a boundary morphism which applies on an ordered -simplex,

(7)

where . Using Algorithm 5, we can apply this morphism to assign a orientation to each topological relation in the complex. Subsequently, for orientable manifolds, we can compute orientations of faces and which share edge such that,

(8)

where and correspond to the edge up from to and respectively. Doing so, we create an oriented simplicial complex.

Supposing that we wish to compute the tangent of a vertex as defined by the weighted average tangent of incident faces. This is equivalent to computing the oriented wedge products of each incident face. This can be written generally as,

(9)
(10)

where is the number of incident faces, is incident face , and are indices of vertex members of not equal to , and . This can be easily computed using a templated recursive function.

// Terminal case
auto getTangentH(const SurfaceMesh& mesh,
    const Vector& origin,
    SurfaceMesh::SimplexID<SurfaceMesh::topLevel> curr){
    return (*curr).orientation;
}
template <std::size_t level, std::size_t dimension>
auto getTangentH(const SurfaceMesh& mesh,
    const Vector& origin,
    SurfaceMesh::SimplexID<level> curr){
    tensor<double, 3, SurfaceMesh::topLevel - level> rval;
    auto cover = mesh.get_cover(curr); // Lookup coboundary relations
    for(auto v : cover){
        auto edge = *mesh.get_edge_up(curr, v); // Get the edge object
        const auto& vtx = (*mesh.get_simplex_up({v})).position; // Vertex v
        auto next = mesh.get_simplex_up(curr,v); // Simplex curr union v
        rval += edge.orientation * (v-origin) * getTangentH(mesh, origin, next);
    }
    return rval/cover.size();
}

This demonstrates the ease using the CASC library as an underlying simplicial complex representation. Using the provided API, it is easy to traverse the complex to perform any computations.

5.1. Preservation of Topology Type of Surface Mesh Under Edge Decimation by Contraction

Supposing we wish to decimate a surface mesh by edge contraction under Def. 6 without changing the topology of the complex. This can be verified by checking the Link Condition, defined with proof from Edelsbrunner(Edelsbrunner, 2001), stating,

Lemma 5.1 ().

Let be a triangulation of a 2-manifold. The contraction of preserves the topological type if and only if .

Revisiting the example from Fig. 4(a), we can construct the topology and check the Link Condition using operations supported by the CASC library.

// Construct the topology
SurfaceMesh mesh;
mesh.insert({0,1,3});
mesh.insert({0,3,5});
mesh.insert({1,3,4});
mesh.insert({3,4,5});
mesh.insert({1,2,4});
mesh.insert({2,4,5});
// Get simplices and edge
SimplexSet<SurfaceMesh> A,B,AB, AcapB;
auto ab = mesh.get_simplex_up({3,4});
auto a = mesh.get_simplex_up({3});
auto b = mesh.get_simplex_up({4});
// Compute the links of each
getLink(mesh, a, A);
getLink(mesh, b, B);
getLink(mesh, ab, AB);
// Link(a): Simplices {1}, {4}, {5}, {1,4}, {0,1}, {0,5}, {4,5}
std::cout << A << std::endl;
// Link(b): Simplices {5}, {2}, {3}, {1}, {3,5}, {1,3}, {1,2}, {2,5}
std::cout << B << std::endl;
// Link(AB): Simplices {1}, {5}
std::cout << AB << std::endl;
set_intersection(A, B, AcapB);
// Link(a) \cap Link(b): Simplices {1}, {5}
std::cout << AcapB << std::endl;
// Check the Link Condition
std::cout << (AB == AcapB) << std::endl; // Evaluates to true

The getLink() function utilizes a series of visitor breadth first traversals down then up the complex, collecting the set of simplex into a SimplexSet object. Then using set operations provided by SimplexSet the set difference and equality comparison are performed. This example highlights the simplicity, clarity, and transparency of using the CASC library.

6. Conclusions

CASC provides a general simplicial complex data structure which allows the storage of user defined types at each simplex level. The library comes with a full-featured API providing common simplicial complex operations, as well as support for complex traversals using a visitor. We also provide a metadata aware decimation algorithm which allows users to collapse simplices of any dimension while preserving data according to a user defined mapping function. Our implementation of CASC using a strongly-typed language is only possible due to recent innovations in language tools. The CASC API abstracts away most of the complicated templating, allowing it to be both modern and easy to use. We anticipate that CASC will not only be of use for the ET community but microscopy and modeling as a whole along with other fields like applied mathematics and CAD.

One limitation is the ease of extending to other languages. The generality of CASC is reliant upon the C++ compiler. Sacrificing this, specific realizations of CASC can be wrapped using tools like SWIG for use in other languages. Another limitation of CASC is the memory efficiency. The current Hasse diagram based implementation was selected for the sake of transparency, ease of traversal, and manipulation. Optimizations to the memory efficiency of CASC can be made by employing more compact representations. The variadic template approach we use to attach user data to simplices is compatible with data structures which explicitly represent all simplices but only a subset of topological relations. This includes data structures such as SimplexTree(Boissonnat and Maria, 2014), IS(De Floriani et al., 2010), and SIG(De Floriani et al., 2004) among others. Other compressed data structures which skip levels of low importance using implicit nodes are not compatible with the current CASC implementation. Skipped levels would need to be implemented as exceptions to the combinatorial variadic rules. Similarly, although CASC in its current form is restricted to the representation of simplicial complexes, the combinatorial strategy can be easily adapted to support other regular polytopes by changing the boundary relation storage rules. In the future we hope to incorporate parallelism into the CASC library. A copy of CASC along with online documentation can found on GitHub https://github.com/ctlee/casc.

Acknowledgements.
The authors would like to thank the anonymous referees for their valuable comments and helpful suggestions. This work is supported by the Sponsor National Institutes of Health, NIGMS Rlhttp://dx.doi.org/10.13039/100000057 under grant numbers Grant #3 and Grant #3. CTL also acknowledges support from the NIH Molecular Biophysics Training Grant Grant #3. MJH was supported in part by the Sponsor National Science Foundation, Division of Mathematical Sciences Rlhttp://dx.doi.org/10.13039/100000121 under awards Grant #3 and Grant #3.

References

  • (1)
  • CGA ([n. d.]) [n. d.]. CGAL, Computational Geometry Algorithms Library. http://www.cgal.org
  • Babuška and Aziz (1976) I. Babuška and A. K. Aziz. 1976. On the Angle Condition in the Finite Element Method. SIAM J. Numer. Anal. 13, 2 (apr 1976), 214–226. https://doi.org/10.1137/0713021
  • Bank et al. (1983) Randolph E. Bank, Andrew H. Sherman, and Alan Weiser. 1983. Refinement Algorithms and Data Structures for Regular Local Mesh Refinement. In Sci. Comput. Appl. Methematics Comput. to Phys. Sci., R. S. Stepleman (Ed.). North-Holland, 3–17.
  • Bank and Xu (1996) Randolph E. Bank and Jinchao Xu. 1996. An Algorithm for Coarsening Unstructured Meshes. Numer. Math. 73, 1 (mar 1996), 1–36. https://doi.org/10.1007/s002110050181
  • Boissonnat et al. (2009) Jean-Daniel Boissonnat, Olivier Devillers, and Samuel Hornus. 2009. Incremental Construction of the Delaunay Graph in Medium Dimension. In Annu. Symp. Comput. Geom. Aarhus, Denmark, 208–216. https://hal.inria.fr/inria-00412437
  • Boissonnat and Maria (2014) Jean-Daniel Boissonnat and Clément Maria. 2014. The Simplex Tree: An Efficient Data Structure for General Simplicial Complexes. Algorithmica 70, 3 (nov 2014), 406–427. https://doi.org/10.1007/s00453-014-9887-3
  • Canino et al. (2011) David Canino, Leila De Floriani, and Kenneth Weiss. 2011. IA*: An Adjacency-Based Representation for Non-Manifold Simplicial Shapes in Arbitrary Dimensions. Comput. Graph. 35, 3 (jun 2011), 747–753. https://doi.org/10.1016/j.cag.2011.03.009
  • De Floriani et al. (2004) Leila De Floriani, David Greenfieldboyce, and Annie Hui. 2004. A Data Structure for Non-Manifold Simplicial D -Complexes. In Proc. 2004 Eurographics/ACM SIGGRAPH Symp. Geom. Process. - SGP ’04. ACM Press, New York, New York, USA, 83. https://doi.org/10.1145/1057432.1057444
  • De Floriani and Hui (2005) Leila De Floriani and Annie Hui. 2005. Data Structures for Simplicial Complexes: An Analysis and a Comparison. In Proc. Third Eurographics Symp. Geom. Process. Eurographics Association, Vienna, 119. http://dl.acm.org/citation.cfm?id=1281920.1281940
  • De Floriani et al. (2010) Leila De Floriani, Annie Hui, Daniele Panozzo, and David Canino. 2010. A Dimension-Independent Data Structure for Simplicial Complexes. In Proc. 19th Int. Meshing Roundtable. Springer Berlin Heidelberg, Berlin, Heidelberg, 403–420. https://doi.org/10.1007/978-3-642-15414-0_24
  • Dyedov et al. (2015) Vladimir Dyedov, Navamita Ray, Daniel Einstein, Xiangmin Jiao, and Timothy J. Tautges. 2015. AHF: Array-Based Half-Facet Data Structure for Mixed-Dimensional and Non-Manifold Meshes. Eng. Comput. 31, 3 (jul 2015), 389–404. https://doi.org/10.1007/s00366-014-0378-6
  • Edelsbrunner (2001) Herbert Edelsbrunner. 2001. Geometry and Topology for Mesh Generation. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511530067
  • Gao et al. (2013) Zhanheng Gao, Zeyun Yu, and Michael Holst. 2013. Feature-Preserving Surface Mesh Smoothing Via Suboptimal Delaunay Triangulation. Graph. Models 75, 1 (jan 2013), 23–38. https://doi.org/10.1016/j.gmod.2012.10.007
  • Holst et al. (2016) M. Holst, O. Sarbach, M. Tiglio, and M. Vallisneri. 2016. The Emergence of Gravitational Wave Science: 100 Years of Development of Mathematical Theory, Detectors, Numerical Algorithms, and Data Analysis Tools. Bull. Amer. Math. Soc. 53 (2016), 513–554. http://dx.doi.org/10.1090/bull/1544
  • Holst and Tiee (2018) M. Holst and C. Tiee. 2018. Finite Element Exterior Calculus for Parabolic Evolution Problems on Riemannian Hypersurfaces. Journal of Computational Mathematics 36, 6 (2018), 792–832.
  • Liu and Joe (1994) A. Liu and B. Joe. 1994. Relationship Between Tetrahedron Shape Measures. BIT 34, 2 (jun 1994), 268–287. https://doi.org/10.1007/BF01955874
  • Ray et al. (2015) Navamita Ray, Iulian Grindeanu, Xinglin Zhao, Vijay Mahadevan, and Xiangmin Jiao. 2015. Array-Based Hierarchical Mesh Generation in Parallel. Procedia Eng. 124 (2015), 291–303. https://doi.org/10.1016/j.proeng.2015.10.140
  • Regge (1961) T. Regge. 1961. General Relativity Without Coordinates. Nuovo Cim. 19, 3 (feb 1961), 558–571. https://doi.org/10.1007/BF02733251
  • Roberts (2014) Elijah Roberts. 2014. Cellular and Molecular Structure As a Unifying Framework for Whole-Cell Modeling. Curr. Opin. Struct. Biol. 25 (apr 2014), 86–91. https://doi.org/10.1016/j.sbi.2014.01.005
  • Si (2015) Hang Si. 2015. TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator. ACM Trans. Math. Softw. 41, 2 (feb 2015), 1–36. https://doi.org/10.1145/2629697
  • Vanden-Eijnden and Venturoli (2009) Eric Vanden-Eijnden and Maddalena Venturoli. 2009. Markovian Milestoning with Voronoi Tessellations. J. Chem. Phys. 130, 19 (may 2009), 194101. https://doi.org/10.1063/1.3129843
  • Yu et al. (2008a) Zeyun Yu, Michael J. Holst, and J. Andrew McCammon. 2008a. High-Fidelity Geometric Modeling for Biomedical Applications. Finite Elem. Anal. Des. 44, 11 (jul 2008), 715–723. https://doi.org/10.1016/j.finel.2008.03.004
  • Yu et al. (2008b) Zeyun Yu, Michael J Holst, Takeharu Hayashi, Chandrajit L Bajaj, Mark H Ellisman, J Andrew McCammon, and Masahiko Hoshijima. 2008b. Three-Dimensional Geometric Modeling of Membrane-Bound Organelles in Ventricular Myocytes: Bridging the Gap Between Microscopic Imaging and Mathematical Simulation. J. Struct. Biol. 164, 3 (dec 2008), 304–313. https://doi.org/10.1016/j.jsb.2008.09.004
  • Zhang (2016) Yongjie Zhang. 2016. Geometric Modeling and Mesh Generation from Scanned Images. Chapman and Hall/CRC. https://doi.org/10.1201/b19466

Appendix B Supplementary Materials

b.1. Derivation of Worse Case Insertion Performance

The number of combinations at each level for a -simplex is given by the binomial theorem, . In the worst case, where the simplicial complex is empty, the total number of nodes created for inserting a -simplex is . The number of edges to represent all topological relations is then given by .

b.2. Templated Insertion Algorithm

Pseudocode for the algorithms presented in this manuscript have been vastly simplified in order to facilitate understanding. For example Algorithm 1, while the non-templated version is appears straightforward, it is impossible to be implemented in C++ directly, due to several typing related issues. First, the function prototype for insert() requires the rootSimplex as the second argument. Simplices at different levels have different types and insert() must be overloaded. Similarly the variable newSimplex and function createSimplex() must know the type of simplex which will be created at compile time.

The actual implementation uses variadic templates to resolve the typing issues. As an example, templated pseudocode for simplex insertion (Algorithm 1) is shown in Algorithm 6. Not only does the templated code automatically build the correct overloaded functions, but it provides many optimizations.

The step-by-step insertion of tetrahedron {1,2,3,4} is shown in Figure 6. Numbered red lines correspond to newNode and root in function insertNode(). Skinny black lines are the topological relations inserted by backfill().

Input: []: Indices of simplices to describe new simplex ,
: simplicial complex
Output: The new simplex
 /* User function to insert simplex {keys} */
Function   insert<n>(keys[n]) {
      return setupForLoop<0, n>(root, keys) /* ’ is the root node */
     
     // The following are private library functions...
      /* Array slice operation. Algorithm 1: keys[0:i] */
      Function   setupForLoop<level, n>(root, keys) { /* General template */
           return forLoop<level, n, n>(root, keys) /* Setup the recursive for loop */
          
          Function   setupForLoop<level, 0>(root, keys) { /* Terminal condition */
                return root
               
               /* Templated for loop. Algorithm 1: for (i = 0; i < n; i++) */
                Function   forLoop<level, antistep, n>(root, keys) { /* For loop for antistep */
                      /* defines next key to add to */
                     insertNode<level, n-antistep>(root, keys) return forLoop<level, antistep-1, n>(root, keys)
                    
                    Function   forLoop<level, 1, n>(root, keys) { /* Stop when = 1 */
                          return insertNode<level, n-1>(root, keys)
                         
                         Function   insertNode<level, n>(root, keys) { /* Insert a new node */
                               v = keys[n] if (){
                                    newNode= root.up[v]
                                    else{ /* Add simplex */
                                         newNode = createSimplex<n>() /* Create a new node, -simplex */ newNode.down[v] = root /* Connect boundary relation */ root.up[v] = newNode /* Connect coboundary relation */ backfill (root, newNode, v)/* Backfill other topological relations */
                                        
                                        /* Recurse to insert any cofaces of newNode. Algorithm 1: insert(keys[0:i], newSimplex) */
                                         return setupForLoop<level+1, n>(newNode, keys)
                                        
                                        Function   backfill<level>(root, newNode, value) { /* Backfilling pointers to other parents */
                                              for (currentNode in root.down){
                                                   childNode = currentNode.up[value] /* Get simplex */ newNode.down[value] = child /* Connect boundary relation */ child.up[value] = newNode /* Connect coboundary relation */
                                                  
                                                  
                                                  
ALGORITHM 6 Templated pseudocode implementation of Algorithm 1.
Figure 6. The Hasse diagrams for the step-by-step insertion of tetrahedron {1,2,3,4} by Algorithm 1. Red lines represent the order of creation for each simplex. The skinny black lines represent where connections to parent simplices are backfilled.

b.3. Code for Getting Neighbors by Adjacency

The C++ code for collecting the set of -simplices sharing a common coface with simplex . A function call to neighbors_up() calls the following code which serves primarily to help the compiler deduce the dimension, , of .

template <class Complex, class SimplexID, class InsertIter>
void neighbors_up(Complex &F, SimplexID s, InsertIter iter)
{
    neighbors_up<Complex, SimplexID::level, InsertIter>(F, s, iter);
}

With the simplex dimension determined, we call an overloaded function which defines the operation for a -simplex.

template <class Complex, std::size_t level, class InsertIter>
void neighbors_up(
        Complex &F,
        typename Complex::template SimplexID<level> s,
        InsertIter iter)
{
    for (auto a : F.get_cover(s))
    {
        auto id = F.get_simplex_up(s, a);
        for (auto b : F.get_name(id))
        {
            auto nbor = F.get_simplex_down(id, b);
            if (nbor != s)
            {
                *iter++ = nbor;
            }
        }
    }
}

Neighbors of are pushed into an insert iterator provided by the user. In this fashion, depending upon the container type the iterator corresponds to, the user can specify whether or not duplicate simplices are returned (std::vector) or not (std::set).