The Stellar decomposition: A compact representation for simplicial complexes and beyond

We introduce the Stellar decomposition, a model for efficient topological data structures over a broad range of simplicial and cell complexes. A Stellar decomposition of a complex is a collection of regions indexing the complex's vertices and cells such that each region has sufficient information to locally reconstruct the star of its vertices, i.e., the cells incident in the region's vertices. Stellar decompositions are general in that they can compactly represent and efficiently traverse arbitrary complexes with a manifold or non-manifold domain. They are scalable to complexes in high dimension and of large size, and they enable users to easily construct tailored application-dependent data structures using a fraction of the memory required by a corresponding global topological data structure on the complex. As a concrete realization of this model for spatially embedded complexes, we introduce the Stellar tree, which combines a nested spatial tree with a simple tuning parameter to control the number of vertices in a region. Stellar trees exploit the complex's spatial locality by reordering vertex and cell indices according to the spatial decomposition and by compressing sequential ranges of indices. Stellar trees are competitive with state-ofthe-art topological data structures for manifold simplicial complexes and offer significant improvements for cell complexes and non-manifold simplicial complexes. We conclude with a highlevel description of several mesh processing and analysis applications that utilize Stellar trees to process large datasets.


Introduction
Efficient mesh data structures play a fundamental role in a broad range of mesh processing applications in computer graphics, geometric modeling, scientific visualization, geospatial data science and finite element analysis. Although simple problems can be easily modeled on small low dimensional meshes, phenomena of interest might occur only on much larger meshes and in higher dimensions. Thus, we often require flexibility to deal with increasingly complex meshes including those defined by irregularly connected heterogeneous and/or multidimensional cell types discretizing spaces with complicated topology. Moreover, as advances in computing capabilities continue to outpace those in memory, it becomes increasingly important to optimize and exploit locality within the mesh as we process and locally query it. Such queries are the primary means of interacting with the mesh and have traditionally been posed in terms of a few spatial and topological primitives. However, while there are simple, intuitive models for representing polygonal surfaces, there are numerous challenges in generalizing these structures to higher dimensions and in scaling to very large meshes.
In this paper, we introduce the Stellar decomposition, a model for topological data structures that supports efficient navigation of the topological connectivity of simplicial complexes and of certain classes of cell complexes, e.g., those composed of quadrilaterals, polygons, hexahedra, prisms and pyramids. We refer to this class of complexes as Canonical Polytope complexes (CP complexes). The defining property of a Stellar decomposition is that the complex is broken up into regions indexing a collection of vertices such that each vertex within a region has sufficient information to locally reconstruct its star, i.e., the set of cells from the complex incident in that vertex.
A Stellar decomposition is general, in that it can easily represent arbitrary complexes with a manifold or non-manifold domain, it is scalable to complexes both in high dimensions and with a large number of cells, and it is flexible, in that it enables users to defer decisions about which topological connectivity relations to encode. It, therefore, supports the generation of optimal application-dependent local data structures at runtime. Due to the locality of successive queries in typical mesh processing applications, the construction costs of these local topological data structures are amortized over multiple mesh operations while processing a local region.
We also formally define and analyze the Stellar tree as a concrete instance of the Stellar decomposition model for spatially embedded complexes. Stellar trees utilize a hierarchical ndimensional quadtree, or kD-tree, as vertex decomposition, and are easily tunable using a single parameter that defines the maximum number of vertices allowed in each local region.
While Stellar trees have been previously utilized in several mesh processing applications ranging from mesh simplification [Fellegara et al., 2020b] to morphological feature extraction [Weiss et al., 2013, Fellegara et al., 2017, they have not been formally defined and their performance has not yet been characterized in relation to existing topological data structures for simplicial and cell complexes. This paper presents a careful study of the storage requirements, generation algorithms and timings and query performance for Stellar trees over a wide range of CP complexes. As we demonstrate in Section 8, Stellar trees are competitive with dimension-specific state-of-the-art topological data structures for (pseudo)-manifold triangle and tetrahedral complexes and offer significant improvements for other CP complexes, especially over data structures for general simplicial complexes in 3D and higher dimensions. The source code for our Stellar tree implementation will be released in the public domain.
Contributions The contributions of this work include: • The formal theoretical definition of a Stellar decomposition over Canonical Polytope (CP) complexes, a class of cell complexes that includes simplicial and cubical complexes of arbitrary dimension, as well as cells in the finite element 'zoo', such as polygons, pyramids and prisms. • The definition of the Stellar tree as a concrete realization of the Stellar decomposition for spatially embedded complexes. The decomposition in a Stellar tree is based on a hierarchical spatial index with a simple tuning parameter to facilitate balancing storage and performance needs. • The definition of Sequential Range Encoding (SRE), a compact encoding for the entities indexed by each region of the decomposition. When applied to CP complexes reindexed by the spatial decomposition of a Stellar tree, SRE yields compressed Stellar trees with only a small overhead relative to the original CP complex's cells.
Outline The remainder of this paper is organized as follows.
In Sections 2 and 3, we review background notions and related work, respectively. In Section 4, we define Stellar decompositions, describe our compact encoding, and provide a high-level description of the procedure for generating a Stellar decomposition. In Section 5, we define the Stellar tree, a spatio-topological realization of the Stellar decomposition. In Section 6, we describe a general mesh processing paradigm that can be followed by applications defined on a Stellar tree. In Section 7, we discuss our experimental setup and evaluate how our tuning parameter affects the quality of a Stellar tree's decomposition and its performance in extracting topological features. We then compare Stellar trees to several state-of-the-art topological data structures in Section 8. In Section 9, we describe how to extract local connectivity information from the Stellar tree and evaluate the performance of these algorithms. We provide a high-level overview of several mesh processing and analysis applications that have benefited from Stellar trees to process large datasets in Section 10 and conclude in Section 11 with some remarks and directions for future work.

Background notions
In this section, we review notions related to cell and simplicial complexes, which are the basic combinatorial structures for representing discretized shapes. Throughout the paper, we use n to denote the dimension of the ambient space in which the complex is embedded, d to represent the dimension of the complex and k to denote the dimension of a cell from the complex, where 0 ≤ k ≤ d.
A k-dimensional cell in the n-dimensional Euclidean space E n is a subset of E n homeomorphic to a closed k-dimensional ball B k = {x ∈ E k : x ≤ 1}. A d-dimensional cell complex Γ in E n is a finite set of cells with disjoint interiors and of dimension at most d such that the boundary of each k-cell γ in Γ consists of the union of other cells of Γ with dimension less than k. Such cells are referred to as the faces of γ. A cell which does not belong to the boundary of any other cell in Γ is called a top cell. Γ is a pure cell complex when all top cells have dimension d. The subset of E n spanned by the cells of Γ is called the domain of Γ. An example of a pure cell 3-complex is shown in Figure 1(a): all its top cells are 3-cells (tetrahedra).
Throughout this paper, we are concerned with a restricted class of cell complexes whose cells can be fully reconstructed by their set of vertices, e.g., via a canonical ordering [Schoof and Yarberry, 1994, Poirier et al., 1998, Remacle and Shephard, 2003, Celes et al., 2005. We refer to this class of complexes as Canonical Polytope complexes (CP complexes), and note that it includes simplicial complexes, cubical complexes, polygonal cell complexes and heterogeneous meshes with cells from the finite element 'zoo ' (e.g., simplices, hexahedra, pyramids, and prisms). In what follows, we denote a CP complex as Σ. An example of a CP complex is shown in Figure 1(b), which contains top edges, triangles, quadrilaterals, and tetrahedra.
A pair of cells in a CP complex Σ are mutually incident if one is a face of the other. They are h-adjacent if they have the same dimension k > h and are incident in a common h-face. We informally refer to vertices (0-cells) as adjacent if they are both incident in a common edge (1-cell) and, similarly, for k-cells that are incident in a common (k−1)-cell (i.e., they are (k−1)adjacent). The (combinatorial) boundary of a CP cell σ is defined by the set of its faces. The star of a CP cell σ is the set of its cofaces, i.e., CP cells in Σ that have σ as a face. An example of star for a 0-cell (vertex) is shown in Figure 2(a). In this example, the star of vertex v 0 is formed by five edges, four triangles, a quad, and a tetrahedron. Of these CP cells, tetrahedron σ 5 , quad σ 1 and triangle σ 4 are top cells. The link of a CP cell σ is the set of all the faces of cells in the star that are not incident in σ. An example of link for a 0-cell (vertex) is shown in Figure 2(b). In this example, the link of v 0 is composed of six vertices, six edges, and a triangle.
Two h-cells σ and σ in Σ are (h−1)-connected if there is a sequence, called an h-path, of (h−1)-adjacent h-cells in Σ from σ to σ . A complex Σ is h-connected, if for every pair of h-cells σ 1 and σ 2 , there is an h-path in Σ joining σ 1 and σ 2 .
We can now define a d-dimensional CP complex Σ as a set of CP-cells in E n of dimension at most d such that: (1) Σ contains all CP-cells in the boundary of the CP-cells in Σ; (2) the intersection of any two CP-cells in Σ is conforming, i.e., it is either empty, or it consists of faces shared by both CP-cells. Simplicial complexes are an important subset of CP complexes whose cells are simplices. Let k be a non-negative integer. A k-simplex σ is the convex hull of k + 1 independent points in E n (with k ≤ n), called vertices of σ. A face of a k-simplex σ is an h-simplex (0 ≤ h ≤ k) generated by h + 1 vertices of σ.
Other important notions are those of manifolds and pseudomanifolds. A subset M of the Euclidean space E n is called a dmanifold, with d ≤ n, if and only if every point of M has a neighborhood homeomorphic to the open d-dimensional ball. A more practical concept for the purpose of representing CP complexes is that of pseudo-manifold. A pure d-dimensional CP complex Σ is said to be a pseudo-manifold when it is (d−1)-connected and its (d−1)-cells are incident in at most two d-cells. Informally, we refer to the connected and compact subspace of E n not satisfying the manifold conditions as non-manifold.
Queries on a cell complex are often posed in terms of topological relations, which are defined by the adjacencies and incidences of its cells. Let us consider a d-dimensional CP complex Σ and a k-cell σ ∈ Σ, with 0 ≤ k ≤ d: • a boundary relation R k,p (σ), with 0 ≤ p < k, consists of the p-cells of Σ in the boundary of σ; • a co-boundary relation R k,q (σ), with k < q ≤ d, consists of the q-cells of Σ in the star of σ; • an adjacency relation R k,k (σ) consists of the set of k-cells of Σ that are (k−1)-adjacent to σ.

Related work
In this section, we review the state of the art on topological mesh data structures, hierarchical spatial indexes, data layouts and distributed mesh data structures.

Topological mesh data structures
There has been much research on efficient representations for manifold cell and simplicial complexes, especially for the 2D case. A comprehensive survey of topological data structures for manifold and non-manifold shapes can be found in [De Floriani and Hui, 2005]. A topological data structure over a cell complex encodes a subset of its topological relations and supports the efficient reconstruction of local topological connectivity over its cells. Topological data structures can be classified according to: (i) the dimension of the cell complex, (ii) the domain to be approximated, i.e., manifolds versus non-manifold shapes, (iii) the subset of topological information directly encoded, and (iv) the organization of topological information directly encoded, i.e., explicit or implicit data structures.
The explicit cells and topological relations can either be allocated on demand using small local data structures, such as linked lists, or contiguously, e.g. using arrays. In the former case, pointers are used to reference the elements, which can be useful when the data structure needs to support frequent updates to the underlying cells or their connectivity. In the latter case, indexes of the cells within the array can be used to efficiently reference the elements. Recently, an approach has been proposed in [Nguyen et al., 2017] to reconstruct topological relations on demand and to cache them for later reuse.
Broadly speaking, topological data structures can be categorized as incidence-based or adjacency-based. Whereas incidence-based data structures primarily encode their topological connectivity through incidence relations over all the complex's cells, adjacency-based data structures primarily encode their connectivity through adjacency relations over its top cells.
The Incidence Graph (IG) [Edelsbrunner, 1987] is the prototypical incidence-based data structure for cell complexes in arbitrary dimension. The IG explicitly encodes all cells of a given cell complex Γ, and for each p-cell γ, its immediate boundary and co-boundary relations (i.e., R p,p−1 and R p,p+1 ). Several compact representations with the same expressive power as the IG have been developed for simplicial complexes [De Floriani et al., 2004, which typically require less than half the storage space as the IG [Canino and De Floriani, 2014].
Several incidence-based data structures have been developed for manifold 2-complexes, which encode the incidences among edges. The half-edge data structure [Mantyla, 1988] is the most widely data structure of this type [CGAL, 2020, OML, 2015. Design tradeoffs for data structures based on halfedges are discussed in [Sieger and Botsch, 2011].
Indexed data structures [Lawson, 1977] provide a more compact alternative by explicitly encoding only vertices, top cells and the boundary relations from top cells to their vertices. Since the cells of a CP complex are entirely determined by their ordered list of vertices, this provides sufficient information to efficiently extract all boundary relations among the cells, but not the co-boundary or adjacency relations. The Indexed data structure with Adjacencies (IA) [Paoluzzi et al., 1993, Nielson, 1997 extends the indexed representation to manifold simplicial complexes of arbitrary dimension by explicitly encoding adjacency relation R d,d , giving rise to an adjacency-based representation. All remaining topological relations can be efficiently recovered if we also encode a top simplex in the star of each vertex (i.e., a subset of relation R 0,d ).
The Corner- Table (CoT) data structure [Rossignac et al., 2001] is also adjacency-based. It is defined only for triangle meshes, where it has the same representational power as the IA data structure. It uses corners as a conceptual abstraction to represent individual vertices of a triangle and encodes topological relations among corners and their incident vertices and triangles. Several efficient extensions of the Corner-Table data structure have been proposed that exploit properties of manifold triangle meshes [Gurung et al., 2011, Luffel et al., 2014. The Sorted Opposite Table (SOT) data structure [Gurung and Rossignac, 2009] extends the Corner-Table data structure to tetrahedral meshes and introduces several storage optimizations. Notably, it supports the reconstruction of boundary relation R d,0 from co-boundary relations R 0,d (implicitly encoded) and R d,d relations (explicitly encoded), reducing its topological overhead by nearly a factor of two. Since modifications to the mesh require non-local reconstructions of the associated data structures, this representation is suitable for applications on static meshes.
The Generalized Indexed data structure with Adjacencies (IA * data structure) [Canino et al., 2011] extends the representational domain of the IA data structure to arbitrary nonmanifold and mixed dimensional simplicial complexes. The IA * data structure is compact, in the sense that it gracefully degrades to the IA data structure in locally manifold neighborhoods of the mesh, and has been shown to be more compact than incidence-based data structures, especially as the dimension increases [Canino and De Floriani, 2014]. A similar data structure for non-manifold complexes was presented in [Dyedov et al., 2015]. A detailed description can be found in Section 8.2.
The Simplex tree [Boissonnat and Maria, 2014] also encodes general simplicial complexes of arbitrary dimension. It explicitly stores all simplices of the complex within a trie [Fredkin, 1960] whose nodes are in bijection with the simplices of the complex. A public domain implementation is available in the GUDHI library [GUDHI, 2018]. We provide a detailed description of this data structure in Section 8.2.
Boissonnat et al. [Boissonnat et al., 2017] also propose two top-based data structures targeting a compact Simplex tree representation. The Maximal Simplex Tree (M ST ) is an induced subgraph of the Simplex tree, in which only the paths corresponding to top simplices are encoded, but most operations require processing the entire complex. The Simplex Array List (SAL) is a hybrid data structure computed from the top simplices of a simplicial complex Σ that improves processing efficiency by increasing the storage overhead. Both the M ST and the SAL are interesting structures from a theoretical point-of-view, but, as described in [Boissonnat et al., 2017], the model does not currently scale to large datasets and results are limited to complexes with only a few thousand vertices. Moreover, to the best of our knowledge, there is no public domain implementation currently available.
The Skeleton-Blocker data structure [Attali et al., 2012] encodes simplicial complexes that are close to flag complexes, simplicial complexes whose top simplices are entirely determined from the structure of their 1-skeleton, i.e., the vertices and edges of the complex, and has been successfully employed for executing edge contractions on such complexes. It encodes the 1skeleton and the blockers, simplices that are not in Σ, but whose faces are. Its generation procedure is computationally intensive for general simplicial complexes since identifying the blockers requires inserting simplices of all dimensions.
We compare the Stellar tree representation with the IA, CoT, and SOT data structures as well as with the Simplex tree, and IA * data structure in Section 8.2.

Hierarchical spatial indexes
A spatial index is a data structure used for indexing spatial information, such as points, lines or surfaces in the Euclidean space. Spatial indexes form a decomposition of the embedding space into regions. This can be driven by: (i) an object-based or a space-based criterion for generating the decomposition; and (ii) a hierarchical or a non-hierarchical (flat) organization of the regions. These properties are independent, and, thus, we can have hierarchical object-based decompositions as well as flat spacebased ones.
We now consider how the regions of a decomposition can intersect. In an overlapping decomposition the intersection between the regions can be non-empty on both the interiors and on the boundary of their domain, while, in a non-overlapping decomposition, intersections can only occur on region boundaries. We say that a region is nested within another region if it is entirely contained within that region. In the remainder of this section, we focus primarily on hierarchical spatial indexes, which can be classified by the dimensionality of the underlying ambient space and by the types of entities indexed.
Hierarchical spatial indexes for point data are provided by Point Region (PR) quadtrees/octrees and kD-trees [Samet, 2006]. In these indexes, the shape of the tree is independent of the order in which the points are inserted, and the points are only indexed by leaf blocks. The storage requirements of these data structures can be reduced by allowing leaf blocks to index multiple points, as in the bucket PR quadtree/octree [Samet, 2006], whose buck-eting threshold determines the number of points that a leaf block can index before it is refined.
Several data structures have been proposed for spatial indexing of polygonal maps (PM), including graphs and planar triangle meshes. PM quadtrees [Samet and Webber, 1985] extend the PR quadtrees to represent polygonal maps considered as a structured collection of edges. While there are several variants (PM 1 , PM 2 , PM 3 and the randomized PMR), which differ in the criterion used to refine leaf blocks, all maintain within the leaf blocks a list of intersecting edges from the mesh. The PM 2 -Triangle quadtree [De Floriani et al., 2008] specializes PM quadtrees over triangle meshes and has been applied to terrain models. The PM index family has also been extended to PMoctrees encoding polyhedral objects in 3D [Carlbom et al., 1985, Navazo, 1989, Samet, 2006, where the subdivision rules have been adjusted to handle edges and polygonal faces of the mesh elements. Another proposal for triangulated terrain models are Terrain trees [Fellegara et al., 2017], that are a spatial index family for the efficient representation and analysis of large-scale triangulated terrains generated from LiDAR (Light Detection and Ranging) point clouds. A collection of spatial indexes for tetrahedral meshes called Tetrahedral trees was developed in , Fellegara et al., 2020a.
We note that data structures in the PM family are spatial data structures optimized for efficient spatial queries on a complex (e.g., point location, containment and proximity queries) and are not equipped to reconstruct the connectivity of the complex. In contrast, the PR-star octree [Weiss et al., 2011] is a topological data structure for tetrahedral meshes embedded in 3D space. It augments the bucket PR octree with a list of tetrahedra incident in the vertices of its leaf blocks, i.e., those in the star of its vertices.
In this paper, we have generalized the PR-star data structure to handle a broader class of complexes (CP complexes) in arbitrary dimensions and with an arbitrary domain (i.e., non-manifold and non-pure complexes). At the same time, our new leaf block encoding exploits the spatial coherence of the mesh, yielding a significant storage saving compared to PR-star trees (see Section 8.1). As we discuss in Section 10, Stellar trees have been shown to be effective in several geometrical and topological applications including local curvature estimation, mesh validation and simplification [Weiss et al., 2011], morphological feature extraction [Weiss et al., 2013] and morphological simplification [Fellegara et al., 2014], among others.

Optimized data layouts
Considerable effort has been devoted to reindexing meshes to better exploit their underlying spatial locality, for example, to support streamed processing [Isenburg and Lindstrom, 2005], better cache locality [Yoon et al., 2005] or compression [Yoon and Lindstrom, 2007].
Cignoni et al. [Cignoni et al., 2003b] introduce an external memory spatial data structure for triangle meshes embedded in E 3 . Whereas our aim is to enable efficient topological operations on the elements of general simplicial and CP complexes, the objective of [Cignoni et al., 2003b] is to support compact out-of-core processing of massive triangle meshes. Since the data structure in [Cignoni et al., 2003b] is dimension-specific, by exploiting geometric and topological properties of triangle meshes in E 3 , it would be difficult to generalize to CP complexes or to higher dimensions. Dey et al. [Dey et al., 2010] use an octree to index a large triangle mesh for localized Delaunay remeshing. Due to the significant overhead associated with their computations, the octrees in [Dey et al., 2010] are typically shallow, containing very few octree blocks. In the context of interactive rendering and visualization of large triangulated terrains and polygonal models, Cignoni et al. [Cignoni et al., 2003a, Cignoni et al., 2004b associate patches of triangles with the simplices of a multiresolution diamond hierarchy [Weiss and De Floriani, 2011].

Distributed mesh data structures
Stellar decompositions and Stellar trees are also related to distributed mesh data structures [Devine et al., 2009, Ibanez et al., 2016, which partition large meshes across multiple processors for parallel processing e.g. in numerical simulations [Anderson et al., 2021, Kirk et al., 2006, Edwards et al., 2010. In the latter, each computational domain maintains a mapping between its boundary elements and their counterparts on neighboring domains. To reduce inter-process communication during computation, each domain might also include one or more layers of elements from other domains surrounding its elements, typically referred to as ghost, rind or halo layers [Poirier et al., 2000, Lawlor et al., 2006, Ollivier-Gooch et al., 2010. Although each region of a Stellar decomposition (or tree) can be seen as a computational domain in a distributed data structure with a single ghost layer (i.e., the elements in the star of its boundary vertices), Stellar trees are aimed at providing efficient processing on coherent subsets of the mesh (regions), where users can generate optimized local topological data structures. In a distributed regime, we envision Stellar trees helping more with fine-grained (intra-domain) parallelism than with coarse-grained multi-domain partitions.

Stellar decomposition
The Stellar decomposition is a model for data structures representing Canonical Polytope (CP) complexes. We denote a CP complex as Σ, its ordered lists of vertices as Σ V and of top CP cells as Σ T . We provide a definition of the Stellar decomposition in Section 4.1, and describe its encoding in Section 4.2.

Definition
Given a CP complex Σ, a decomposition ∆ of its vertices Σ V is a collection of subsets of Σ V such that every vertex v ∈ Σ V belongs to one of these subsets. We will refer to the elements of decomposition ∆ as regions, which we will denote as r.
A Stellar decomposition S D defines a map from the regions of a decomposition ∆ of its vertex set Σ V to the vertices and top CP cells of complex Σ. Formally, a Stellar decomposition is defined by three components: 1. a CP complex Σ; 2. a decomposition ∆ whose regions cover the vertices of Σ; 3. a map Φ from regions of ∆ to entities of Σ. Thus, a Stellar decomposition is a triple S D = (Σ, ∆, Φ). Since Σ is entirely characterized by its vertices, and top CP cells, we define map Φ in terms of the two components: Φ V ERT , the mapping to vertices and Φ T OP , the mapping to top CP cells.
For the vertices, we have a map from ∆ to Σ V based on an application-dependent 'belongs to' property. Formally, Φ V ERT : Figure 3 illustrates an example decomposition ∆ over a point set where mapping function Φ V ERT associates points with regions of ∆. In this paper, we will assume that each vertex in Σ V is uniquely associated with a single region r in ∆.
The Stellar decomposition gets its name from the properties of its top cell map Φ T OP . For each region r of ∆, Φ T OP (r) is the set of all top CP cells of Σ T incident in one or more vertices of Φ V ERT (r). In other words, Φ T OP (r) is defined by the union of cells in the star of the vertices in Φ V ERT (r). Formally, Φ T OP : ∆ → P(Σ T ) is a function from ∆ to the powerset of Σ T , where (1) Figure 4 illustrates mapping Φ T OP for two regions of the decomposition of Figure 3(b) on a triangle mesh defined over its vertices. Note that Φ T OP is based on a topological rather than a spatial property: A top CP cell σ is only associated with a region r when one (or more) of its vertices is associated with r under Φ V ERT .
To characterize this representation, we define the spanning number χ σ of a top CP cell in a Stellar decomposition as the number of regions to which it is associated.
Definition 4.1 Given Stellar decomposition S D = (Σ, ∆, Φ), the spanning number χ σ of a top CP cell σ ∈ Σ T is the number of regions in ∆ that map to σ. Formally, ∀σ ∈ Σ T , χ σ = |{r ∈ ∆|σ ∈ Φ T OP (r)}|. (2) A consequence of the unique mapping of each vertex in Φ V ERT is that it provides an upper bound on the spanning number of a top CP cell in a Stellar decomposition. Specifically, the spanning number χ σ of a top CP cell σ is bounded by the cardinality of its vertex incidence relation R k,0 : It is also interesting to consider the average spanning number χ as a global characteristic of the efficiency of a Stellar decomposition over a complex, measuring the average number of times each top CP cell is represented.
Definition 4.2 The average spanning number χ of a Stellar decomposition S D is the average number of regions indexing its top cells. Formally, (3)

Encoding
In this section, we describe how we represent the two components of a Stellar decomposition, providing a detailed description of the data structures for representing a CP complex (subsection 4.2.1), and a compressed encoding for the regions of the partitioning (subsection 4.2.2). We do not describe how the decomposition ∆ is represented, as this is specific to each concrete realization of the Stellar decomposition model.

Indexed representation of the CP complex
We represent the underlying CP complex as an indexed complex [Lawson, 1977], which encodes the vertices, top CP k-cells and the boundary relation R k,0 of each top k-cell in Σ. In the following, we assume a d-dimensional CP complex Σ embedded in E n . We use an array-based representation for the vertices and for the top cells of Σ. Since the arrays are stored contiguously, each vertex v has a unique position index i v in the vertex array, that we denote as Σ V . Similarly, each top CP cell σ has a unique position index i σ . The top CP cells of Σ are encoded using separate arrays Σ T k for each dimension k ≤ d that has top CP cells in Σ. Σ T k encodes the boundary connectivity from the top CP k-cells of Σ to their vertices, i.e., relation R k,0 in terms of the indices i v of the vertices of its cells within Σ V . This requires |R k,0 (σ)| references for a top k-cell σ, e.g., k+1 vertex indices for a k-simplex and 2 k references for a k-cube. Thus, the total storage cost of Σ T is: We note that when Σ is pure (i.e., its top CP cells all have the same dimension d), the encoding of Σ requires only two arrays: one for the vertices and one for the top cells. For simplicity, we refer to the top cell arrays collectively as Σ T .

A compressed region representation
In this subsection, we discuss two encoding strategies for the data maps in each region of the partition ∆. We begin with a simple strategy that explicitly encodes the arrays of vertices and top CP cells associated with each region and work our way to a compressed representation of these arrays. Coupling this compressed representation with a reorganization of the vertices and cells of the CP complex (as we will describe in Section 4.3) yields a significant reduction in storage requirements. We will demonstrate  this claim in Section 8.1 on a data structure instantiating the Stellar decomposition.
Recall that under Φ, each region r in ∆ maps to an array of vertices and an array of top CP cells from the complex Σ which we denote as r V and r T , respectively. A straightforward strategy would be to encode arrays of vertices and top CP cells that explicitly enumerate the associated elements for each region r. We refer to this as the EXPLICIT Stellar decomposition encoding. An example of this encoding for a single region with six vertices in r V and twenty triangles in r T is shown in Figure 5.
It is apparent that the above encoding can be very expensive due to the redundant encoding of top CP cells with vertices in multiple regions. A less obvious redundancy is that it does not account for the ordering of the elements.
We now consider a COMPRESSED Stellar decomposition encoding that compacts the vertex and top CP cells arrays in each region r by exploiting the locality of the elements within r. The COMPRESSED encoding reduces the storage requirements within region arrays by replacing runs of incrementing consecutive sequences of indices using a generalization of run-length encoding (RLE) [Held and Marshall, 1991]. RLE is a form of data compression in which runs of consecutive identical values are encoded as pairs of integers representing the value and repetition count, rather than as multiple copies of the original value. For example, in Figure 6a, the four entries with value '2' are com- pacted into a pair of entries [-2, 3], where a negative first number indicates the start of a run and its value, while the second number indicates the remaining elements of the run in the range. While we do not have such duplicated runs in our indexed representation, we often have increasing sequences of indexes, such as {40,41,42,43,44}, within a local vertex array r V or top CP cells array r T . We therefore use a generalized RLE scheme to compress such sequences, which we refer to as Sequential Range Encoding (SRE). SRE encodes a run of consecutive non-negative indexes using a pair of integers, representing the starting index, and the number of remaining elements in the range. As with RLE, we can intersperse runs (sequences) with non-runs in the same array by negating the starting index of a run (e.g. [-40, 4] for the above example). Thus, it is easy to determine whether or not we are in a run while we iterate through a sequential range encoded array. A nice feature of this scheme is that it allows us to dynamically append individual elements or runs to an SRE array without any storage overhead (other than occasional array reallocations). Furthermore, we can easily expand a compacted range by replacing its entries with the first two values of the range and appending the remaining values to the end of the array. After the updates are finished, we can sort the array and reapply SRE compaction to recover space. Figure 6b shows an example SRE array over an array, where, e.g., sequence {1,2,3,4} is represented as [-1, 3].
To facilitate comparisons between the EXPLICIT and COM-PRESSED representations of a Stellar decomposition, we introduce a global characteristic that measures the average storage requirements to represent a top CP cell.
where |r T | is the size of the top CP cells array in a region r.
In contrast to the average spanning number χ, which is a property of the decomposition, the average reference number µ is a property of how the decomposition is encoded. An EXPLICIT representation is equivalent to a COMPRESSED representation without any compressed runs, and, thus, it is always the case that µ ≤ χ.
In the EXPLICIT representation (i.e., without any sequence-based compression), µ = χ, while in the COMPRESSED representation, µ decreases as the compression of the r V and r T arrays becomes more effective. Figure 7 illustrates a COMPRESSED representation of the mesh from Figure 5 after its vertex and triangle arrays have been reordered (in an external process) and highlights its sequential ranges, where r V is encoded by a single run and r T is encoded by four sequential runs as well as several non-run indices.

Generating a Stellar decomposition
We now describe how to generate a COMPRESSED Stellar decomposition from an indexed CP complex Σ and a given partition ∆ on its vertices Σ V . This process consists of three phases: 1. reindex the vertices of Σ following a traversal of the regions of ∆ and SRE-compress the r V arrays; 2. insert the top CP cells of Σ into ∆; 3. reindex the top CP cells of Σ based on locality within common regions of ∆ and SRE-compress the regions r T arrays. As it can be noted, the generation process ignores how the partitioning on the vertices is obtained, since this step is defined by the data structure instantiating a Stellar decomposition. The reindexing of the vertices follows a traversal of the regions of ∆ in such a way that all vertices associated with a region have a contiguous range of indices in the reindexed global vertex array Σ V (as detailed in the Appendix A.1).
The second phase inserts each top CP k-cell σ, with index i σ in Σ T k , into all the regions of ∆ that index its vertices. This is done by iterating through the vertices of σ and inserting i σ into the r T array of each region r whose vertex map Φ V ERT (r) contains at least one of these vertices. As such, each top CP k-cell σ appears in at least one and at most |R k,0 (σ)| regions of ∆. Due to the vertex reindexing of step 1, this operation is extremely efficient.
Determining if a vertex of a given cell lies in a block requires only a range comparison on its index i v . Finally, we reindex the top CP cell arrays Σ T to better exploit the locality induced by the vertex-based partitioning and compress the local r T arrays using a sequential range encoding over this new index. The reindexing and the compression of the top CP cells is obtained following a traversal of the regions of ∆ in such a way that all top CP cells associated with the same set of regions have a contiguous range of indices in the reindexed arrays Σ T . This last step is detailed in the Appendix A.2 and B. As we demonstrate in Section 8, this compression yields significant storage savings in a wide range of mesh datasets.

Stellar trees
The Stellar decomposition is a general model that is agnostic about how the partitioning is attained and about its relationship with the underlying CP complex. Thus, for example, we can define a Stellar decomposition using Voronoi diagrams, or based on a nearest neighbor clustering of the vertices of a given CP complex. In this section, we introduce Stellar trees as a class of Stellar decompositions defined over nested spatial decompositions of the CP complex and discuss some of our design decisions. Before defining a Stellar tree (Section 5.1), its encoding (Section 5.2) and its generation procedure (Section 5.3), we review some underlying notions.
The ambient space A is the subset of E n in which the data is embedded. We consider the region bounding the ambient space to be a hyper-rectangular axis-aligned bounding block, which we refer to simply as a block. A k-dimensional closed block r in E n , with k ≤ n, is the Cartesian product of closed intervals [l i , u i ] where exactly k of them are non-degenerate, i.e., Given two blocks r := [l i , u i ] and r := [l i , u i ], r is a face of r if, for each dimension i, either their intervals overlap (i.e., l i = l i and u i = u i ) or the i-th interval of r is degenerate (i.e., Given a block r, we refer to its 0-dimensional face of degenerate intervals x i = l i as its lower corner and to its 0-dimensional face where x i = u i as its upper corner. The above block definition describes closed blocks. It can be useful to allow some faces of r to be open, especially on faces of neighboring blocks that overlap only on their boundaries. We now focus on nested decompositions, hierarchical spacebased decompositions whose overlapping blocks are nested and whose leaf blocks ∆ L (i.e., those without any nested blocks) form a non-overlapping cover of the ambient space A. The nesting relationship defines a containment hierarchy H, which can be described using a rooted tree. The tree's root H ROOT covers the ambient space A; the tree's leaves H L encode the regions of the decomposition ∆; and its internal nodes H I provide access to the regions of the decomposition.
Nested decompositions can adopt different hierarchical refinement strategies. Among the most popular are those based on regular refinement and bisection refinement of simple primitives (e.g., simplices and cubes). An n-dimensional block r is regularly refined by adding vertices at all edge and face midpoints of r and replacing r with 2 n disjoint blocks covering r. This generates quadtrees in 2D, and octrees in 3D [Samet, 2006]. In bisection refinement, a block is bisected along an axis-aligned hyperplane into two blocks, generating kD-trees [Bentley, 1975].

Definition
Since a Stellar tree S T is a type of Stellar decomposition, it consists of three components: (1) a CP complex Σ embedded in an ambient space A; (2) a nested decomposition ∆ covering the domain of Σ; and (3) a map Φ from blocks of ∆ to entities of Σ. The nested decomposition is described by a containment hierarchy H, represented by a tree whose blocks use the half-open boundary convention to ensure that every point in the domain is covered by exactly one leaf block.
Since Stellar trees are defined over nested spatial decompositions that cover the ambient space, we customize the vertex mapping function Φ V ERT to partition the vertices of Σ according to spatial containment: each vertex is associated with its single containing leaf block. Formally, A two-dimensional example is shown in Figure 8, where a set of points are associated with the leaf blocks of ∆ through Φ V ERT .
The top CP cells mapping function Φ T OP for a Stellar tree has the same definition as for the Stellar decomposition (see Equation 1). Figure 9 shows the mapping Φ T OP for two blocks of the nested kD-tree decomposition of Figure 8(b) over the triangle mesh from Figure 4.
Since the nested decomposition ∆, and, consequently, the tree H describing it, are determined by the number of vertices indexed by a block, we utilize a bucket PR tree [Samet, 2006] to drive our decomposition. This provides a single tuning parameter, the bucketing threshold k V , that uniquely determines the decomposition for a given complex Σ.
Recall that a (leaf) block r in a bucket PR-tree is considered full when it indexes more than k V vertices (in our case, when |Φ V ERT (r)| > k V ). Insertion of a vertex into a full block causes the block to refine and to redistribute its indexed vertices among its children. As such, the domain decomposition of a Stellar tree depends only on the bucketing threshold k V . Smaller values of k V yield deeper hierarchies whose leaf blocks index relatively few vertices and top CP cells, while larger values of k V yield shallower hierarchies with leaf blocks that index more vertices and top CP cells. Thus, k V and the average spanning number χ of a Stellar tree are inversely correlated.
In practice, we use different spatial indexes to represent H based on the dimension n of the ambient space A. In lower dimensions, we use a quadtree-like subdivision, i.e., a quadtree in 2D, and an octree in 3D, while in higher-dimensions, we switch to a kD-tree subdivision. As discussed in [Samet, 2006], while quadtree-like subdivisions are quite efficient in low dimensions, the data becomes sparser in higher dimensions (due to the curse of dimensionality [Bellman, 1966]), and tends to be better encoded by kD-trees.

Encoding
We represent the containment hierarchy H using an explicit pointer-based data structure, in which the blocks of H use a type of Node structure that changes state from leaf to internal block during the generation process of a Stellar tree.
We use a brood-based encoding [Hunter and Willis, 1991], where each block in H encodes a pointer to its parent block and a single pointer to its brood of children. This reduces the overall storage since leaves do not need to encode pointers to their children, and also allows us to use the same representation for n-dimensional quadtrees and kD-trees. We explicitly encode all internal blocks, but only represent leaf blocks r in H with nonempty maps Φ(r).
The mapped entities of the CP complex Σ are encoded in the leaf blocks H L using the mapping arrays Φ. Note that each leaf block r encodes the arrays of vertices r V and of top CP cells r T in terms of the indices i v and i σ , respectively, that identify v and σ in the Σ V and Σ T arrays. For each block r, we have: (1)  pointers for the hierarchy: one to its parent, another to its list of children and it is pointed to by one parent; (2) a pointer to an array of vertices r V and the size of this array; (3) a pointer to an array of top CP cells r T and the size of this array. Thus, the hierarchy H of a Stellar tree requires 7|H| storage. By considering the encodings, defined in Section 4.2.2, for the CP complex Σ, and for the vertices and top cp-cells associated with the regions of H, we can estimate the storage requirements for the EXPLICIT and COMPRESSED Stellar trees. An EXPLICIT Stellar tree requires a total of |Σ V | references for its vertex arrays, since each vertex is indexed by a single leaf block, and a total of χ|Σ T | references for all top CP cells arrays. Thus, the total cost of the EXPLICIT Stellar tree, including the hierarchy (but excluding the cost of the indexed mesh) is: Conversely, in a COMPRESSED Stellar tree, we can reindex the vertex array Σ V in such a way that all vertices associated with the same leaf block are indexed consecutively (see Section A.1 in the Appendix for additional details). Thus, we can encode the r V arrays using only two integers per leaf block for a total cost of 2|H L | rather than |Σ V |. Moreover, since leaf blocks no longer need to reference an arbitrary array, these two references can be folded into the block's hierarchical representation for r V : instead of a pointer to a array and a size of the array, we simply encode the range of vertices in the same space. As the cost of representing the r T arrays is µ|Σ T |, the total cost for encoding a COMPRESSED Stellar tree (excluding the cost of the indexed mesh representation) is: 7|H| + µ|Σ T |.

Generating a Stellar tree
In this section, we describe how to generate a COMPRESSED Stellar tree from an indexed CP complex Σ and a given bucketing threshold k V . We can also deal with input complexes that are not already indexed. For example, if our input is a "soup" of CP cells in which each CP cell is specified by a list of coordinates, we can generate an indexed representation of the complex as we insert the vertices and generate the decomposition.
First, given a user-defined bucketing threshold k V , we generate a bucket PR-tree over the vertices of Σ. The procedure for inserting a vertex v with index i v in Σ V into H is recursive. We use the geometric position of v to traverse the internal blocks to reach the unique leaf block r containing v. After adding v to r (i.e., appending i v into the r V array of r), we check if this causes an overflow in r. If it does, we refine r and reinsert its indexed vertices into its children. Once all the vertices in Σ have been inserted, the decomposition is fixed.
The rest of the Stellar tree generation process follows the strategy described in Section 4.3 and detailed in the Appendix A.
One key optimization between a generic partitioning on the vertices and a nested hierarchical decomposition relates to extracting the vertex index ranges. In a Stellar tree, this step is performed through a depth-first traversal of the tree, which, for each leaf block r, generates a contiguous range of indices for the vertices in r, and, for each internal block, provides a single contiguous index range for the vertices in all descendant blocks. For example, in Figure 10, after executing this step on leaf block b, we have v s = 4 and v e = 7. Similarly, at the end of this step the root H ROOT has v s = 1 and v e = 13.
We provide an experimental evaluation of the timings for generating a Stellar tree in Section 8.3.

Processing paradigm for Stellar trees
Mesh processing applications rarely process individual mesh elements. Rather, they typically operate on the entire complex, or on large regions of interest within the complex. The structure of a Stellar tree naturally supports a batched processing strategy, i.e., a strategy in which portions of the complex are reconstructed and processed within each block of the tree. As these local blocks are processed, their representation and extracted topological relations can be customized to suit the needs of the application. This helps in amortizing the reconstruction costs and, thus, processing the entire complex efficiently.
The general paradigm for executing application algorithms on a Stellar tree is to iterate through the leaf blocks of hierarchy H, locally processing the encoded complex in a streaming manner. For each leaf block r in H, a local topological data structure catered to the application's needs is constructed and used to process the indexed subcomplex. We refer to this local data structure in a block r as an expanded leaf-block representation, and we denote it as r E . Once we finish processing leaf block r, we discard r E and begin processing the next block.
For efficiency and at relatively low storage overhead, we can cache the expanded leaf block representation r E , using a Least-Recent-Used (LRU) cache. This is especially advantageous in applications that require processing portions of the complex in neighboring leaf blocks. Adopting a fixed-size cache allows us to amortize the extraction costs of the local data structures, with a controllable storage overhead.
Algorithm 1 outlines the general strategy for processing a Stellar tree. The algorithm recursively visits all the blocks of the hierarchy H. For each leaf block r, we either recover r E from the LRU cache (rows 5-8), or construct the desired applicationdependent local topological data structure r E . After using this local data structure to process the local geometry in r (row 9), we either cache or discard r E (rows 10-13).
Within this general processing paradigm, we can have two different approaches, that we call local and global, depending on how auxiliary data structures are encoded and maintained. In for all blocks r C in CHILDREN(r) do discard r E a local approach, the scope of these auxiliary data structures is limited to that of a single leaf block r, or to a restricted subset of its neighbors. In general, a local approach is preferred for applications that extract, or analyze local features, such as those that depend only on the link or star of cells. These includes, for instance, the extraction of geometric features, like the curvature at a vertex, or the extraction of morphological features, like critical points, when the complex is a discretization of the domain of a scalar field. In these examples, the auxiliary data structures are just needed within the scope of a leaf block r, and thus, immediately discarded after extracting the corresponding feature in r. Conversely, in a global approach, data structures are maintained over the entire complex and updated during the visit of the tree. A global approach can be preferable for applications that require the analysis or the processing of the entire complex, like geometric simplification, morphological segmentation, or validation of geometric and topological properties. In these examples, auxiliary data structures are used to represent partial results over the complex.
The decision between using a local and global approach can be driven by the needs of the application or as a tradeoff balancing memory usage and execution times. Due to the limited scope of auxiliary data structures in the local approach, the storage overhead is typically proportional to the complexity of the local complex but requires an increased number of memory allocations compared to a global approach since each leaf block expansion requires memory allocations. Conversely, while auxiliary data structures in the global approach are allocated only once, these structures can require significantly more storage space compared to the local approach.
In Section 10, we present applications on mesh processing and analysis, based on Stellar trees, on which these two paradigms have been extensively applied.

Experimental setup
In this section, we describe our experimental setup, including the datasets used in our evaluation (Section 7.1). We also evaluate how the bucketing threshold k V affects the quality of a Stellar tree's decomposition and its performance in extracting topological queries (Section 7.2).

Experimental datasets
We have performed experiments on a range of CP complexes consisting of triangle, quadrilateral, tetrahedral and hexahedral meshes in E 3 as well as pure non-manifold simplicial complexes in higher dimensions and higher dimensional non-manifold simplicial complexes (embedded in E 3 ). Table 1 summarizes the datasets used in our experiments and their numbers of vertices and top cells.
Our triangle and tetrahedral meshes are native models ranging from 4 to 28 million triangles and from 24 to 29 million tetrahedra, where we use the term native to refer to models from public domain repositories discretizing objects in space. Since we only had access to relatively small native quadrilateral and hexahedral meshes (with tens to hundreds of thousand elements), we have generated some larger models ranging from 12 to 125 million elements from our triangle and tetrahedral models. The generation procedure refines each triangle into three quadrilaterals and each tetrahedron into four hexahedra by adding vertices at the face centroids.
To experiment with pure non-manifold models in higher dimensions, we have generated some models based on a process that we call probabilistic Sierpinski filtering, where we regularly refine all simplices in the complex and randomly remove a fixed proportion of the generated simplices in each iteration. For our experiments, we have created 5-, 7-and 40-dimensional models using different levels of refinement and a filtering threshold of 65%, yielding pure simplicial complexes with 16.5 million to 258 million top simplices.
Finally, to experiment with general simplicial complexes in higher dimensions, we have generated several (non-pure) Vietoris-Rips complexes, which we embed in a lower dimensional space. A Vietoris-Rips (V-Rips) complex is the flag complex defined by a neighborhood graph over a point cloud whose arcs connect pairs of points with distance less than a user-provided parameter . Given the neighborhood graph, the simplices of the V-Rips complexes are defined by its cliques, subsets of the graph vertices that form a complete subgraph. We refer to [Zomorodian, 2010] for further details. For our experiments, we have generated V-Rips complexes over the vertices of a trian- gle model (LUCY) and of two tetrahedral models (VISMALE and FOOT) from our manifold datasets and set our distance threshold to {0.1%, 0.5%, 0.4%} of the bounding box diagonal, respectively. The range of top simplices in the generated complexes goes from 6.4 million to 64 million and their dimension from 7 to 34. Although the generated complexes are synthetic, they provide a good starting point to demonstrate the efficiency of the Stellar tree in higher dimensions.
All tests have been performed on a PC equipped with a 3.2 gigahertz Intel i7-3930K CPU with 64 gigabytes of RAM. The source code will be made available at [Fellegara, 2021].

Calibrating Stellar tree bucket thresholds
Spatial indexes typically involve a careful balance among index generation times, storage costs and query performances. Stellar trees provide users with a single tuning parameter k V to control the maximum number of vertices indexed by each block of the tree. In the following, we calibrate k V on a characteristic subset of three of our experimental datasets: NEPTUNE triangle mesh, BONSAI tetrahedral mesh, and VISMALE Vietoris-Rips complex. For each dataset, we generated 195 Stellar trees using k V values ranging from 1 to 1500 and compared Stellar tree generation and query times as well as the number of blocks as a proxy for the complexity of the generated tree. Within this range, we increment k V by 1 for values between 1 and 50, and by 10 for values between 60 and 1500. This allows us to evaluate the decomposition quality and the extraction performance for a fundamental topological query at different scales. For the latter, we use the vertex co-boundary extraction, i.e., the top cells incident in each vertex (which we describe in Section 9.2).
The results are summarized in the charts of Figures 11 and 12, which compare the complexity of the generated Stellar tree (in terms of number of blocks), its generation and query times, and the average spanning (χ) and reference (µ) numbers as a function of the threshold value k V .
To better compare different units (i.e., number of blocks and timings), each chart in Figure 11 has two logarithmic y-axes, showing the time scale (blue curves using the left y-axes) and the number of blocks (red curves using the right y-axes), respectively. In this way, we can directly compare, for each dataset, how the k V value influences the decomposition and the timing performances. After an initial rapid decrease in the generation time and block number, the curves begin to level off for increasingly large k V values. While there are more than a million blocks when k V is less than 10, the number of blocks rapidly decreases to hundreds of thousands for k V 's between 50 and 200, and grows even smaller for large k V values (e.g., above 500), where the number of blocks remains steadily in the thousands to tens of thousands. This trend appears to be related to the point distribution within each dataset, which induces finer decompositions for k V values between 1 and 50, and coarser decompositions for larger k V values. This trend can also be observed for the generation times, which reduce by a factor of two for k V values between 1 and 100, and by another factor of two for large bucketing thresholds. While the topological extraction query is largely unaffected by k V size, it gets slightly faster for larger k V values. When comparing the influence of k V on χ and µ (shown in Figure 12), we observe that the behavior of these two variables is very similar to that of the number of blocks. This is expected, since the top cells distribution is directly linked to the number of blocks in the tree. As mentioned in Section 4.1, the number of leaf blocks indexing a top cell is bounded from above by the number of its vertices, and this defines a topological upper bound that reduces the overall storage requirements. We note that the SRE compression is able to reduce the number of references per top cell (µ), even for very small k V values.
Our calibration experiments indicate that, while there are slight differences in timing and storage costs, Stellar tree performance is relatively stable over a wide range of k V values. However, threshold values that are either too small or too large should be avoided, since in the first case the storage requirements and time performances are heavily affected, while in the latter case the benefit of having a hierarchical decomposition is limited, as both storage requirements and time performance are not clearly influenced by it. In the rest of this paper, for every model, we build two Stellar trees to compare how their performances depend on parameter k V . These two k V values are chosen in order to: (i) have a hierarchical decomposition that still plays a critical role in the storage requirements and time performances; and (ii) obtain trees with different characteristics: one deeper and another relatively shallower. In the following, we use k S to refer to the smaller k V value and k L to the larger one. Since there is a direct correlation between the decomposition quality and χ, these calibration choices are also reflected in the χ values across our experimental datasets. Table 2 summarizes statistics on the Stellar trees obtained from each input data set by considering two values of the vertex threshold k V , namely k S and k L . Figure 13 illustrates the k S octree decomposition for the 4M triangle NEPTUNE dataset.

Evaluation of storage costs and generation times
In this section, we evaluate the storage costs and generation times of Stellar trees. First, we compare the cost of different Stellar tree encodings (Section 8.1), then we compare the Stellar tree against several state-of-the-art topological mesh data structures (Section 8.2), and, finally, we evaluate the generation times of the Stellar tree (Section 8.3).

Storage comparison among Stellar tree encodings
We begin by comparing the EXPLICIT and COMPRESSED Stellar tree encodings as well as a VERTEX-COMPRESSED encoding, similar to the PR-star encoding for tetrahedral meshes [Weiss et al., 2011], that compresses the vertex array but not the top cells arrays. Table 3 lists the storage costs for the indexed representation of the complex ('Base Complex') as well as the additional costs required for the three Stellar tree encodings, in terms of megabytes (M Bs). In the following, we assume that pointers require 64 bits and indices 32 bits, the de-facto standard in modern computing hardware. Stellar trees based on the COM-PRESSED encoding are always the most compact. We first consider the storage requirements of the hierarchical structures with respect to our tuning parameter k V and observe that higher values of k V always yield reductions in memory requirements. As expected, this effect is more pronounced for the COMPRESSED encoding than for the other two encodings. Specifically, the EXPLICIT and VERTEX-COMPRESSED k L trees achieve a 20-50% reduction in storage requirements compared to their k S counterparts, while the COMPRESSED k L trees are a factor of 3-10 smaller than their k S counterparts. For example, on the triangular NEPTUNE dataset, storage requirements for the EXPLICIT Stellar tree reduces from 32.0 MB (k S ) to 26.2 MB (k L ), while the COMPRESSED Stellar trees reduces by more than a factor of 4 from 5.76 MB (k S ) to 1.24 MB (k L ). When comparing the three encodings, we see that compressing the vertices alone, as in the VERTEX-COMPRESSED representation, achieves only 10-20% reduction in storage requirements compared to the EXPLICIT representation, in most cases. In contrast, compressing the vertices and top cells, as in our COM-PRESSED representation, yields an order of magnitude improvement, requiring a factor of 10-20 less storage than their EXPLICIT counterparts. This trend is nicely tracked for each dataset by the differences between its average references number µ and its average spanning number χ. This is particularly evident on our probabilistic datasets, for which it is difficult to calibrate k V in order to reduce χ values. However, after SRE compression, µ values are always very small, leading to significant storage reductions in the COMPRESSED representation.
Considering the hierarchical storage requirements against those of the original indexed base complex, we observe that EX-PLICIT Stellar trees require about 50% to 80% the storage of the base complex, while COMPRESSED Stellar trees require only around 10% (k S ) and 1% (k L ) the storage of the base complex. Thus, for reasonable k V values, COMPRESSED Stellar trees impose only a negligible storage overhead with respect to the underlying indexed complex, which the Stellar tree representation does not modify. In the remainder of this paper, we restrict our attention to the COMPRESSED Stellar Tree, which we refer to as the Stellar tree, for simplicity.

Storage comparison with respect to other data structures
We compare the Stellar tree with several dimension-independent topological data structures as well as dimension-dependent topological data structures for 2D and 3D simplicial complexes. Figures 14, 15 and 16 compare the storage requirements for the different data structures normalized against the storage costs of the indexed base complex. The analysis compares the topological overhead of the data structures, and thus, we omit the cost of the geometry of the underlying complex, which is common to all the data structures. Based on our analysis of the literature (see Section 3.1),  [Edelsbrunner, 1987], the Incidence Simplicial (IS) , the Simplex tree [Boissonnat and Maria, 2014], and the Generalized Indexed data structure with Adjacencies (IA * ) [Canino et al., 2011]. Since Canino et al. [Canino and De Floriani, 2014] demonstrated that the IA * data structure is more compact than the IG and IS data structures for low and high-dimensional datasets, we restrict our comparisons to the IA * and Simplex tree data structures. The IA * data structure has been defined for dimensionindependent simplicial complexes, and for our experiments, we have extended it to dimension-independent CP complexes. It explicitly encodes all vertices and top CP k-cells in Σ, with 0 < k ≤ d, as well as the following topological relations: (i) boundary relation R k,0 (σ), for each top CP k-cell σ; (ii) adjacency relation R k,k (σ), for each top CP k-cell σ; (iii) co-boundary relation R k−1,k (τ ), for each non-manifold (k−1)-cell τ bounding a top CP k-cell; (iv) partial co-boundary relation R * 0,k (v), for each vertex v, consisting of an arbitrary top CP k-cell σ from each k-cluster in the star of v. A k-cluster is a (k−1)-connected component of the star of v restricted to its top CP k-cells. Note that for pure CP complexes, the non-manifold co-boundary relation R k−1,k is empty. Further, for pseudo-manifold complexes, the partial vertex co-boundary relation R * 0,k has car-   Costs (labels to right of each bar) are normalized to the indexed mesh representation (listed along yaxis). Note that: (1) the x-axis is truncated to a factor of 3; (2) datasets marked with or ⊗ could not be directly generated on our test machine for the Simplex tree or IA * (respectively); and (3) the Simplex tree results for the PROB.40D and LUCY34D dataset are partial (>). dinality 1, and the IA * is identical to the IA data structure [Paoluzzi et al., 1993].
The Simplex tree encodes all j-simplices in Σ, with 0 ≤ j ≤ d, like the IG, while storing a subset of the incidence relations encoded by the IG. The Simplex tree is defined over a total order on the vertices of Σ, and thus, each simplex σ is uniquely represented as an ordered path in a trie whose nodes correspond to the boundary vertices of σ. Thus, the nodes are in bijection with the simplices of the complex, and a Simplex tree over a simplicial complex with |Σ| simplices (of any dimension) contains exactly |Σ| nodes. This provides an efficient representation for extracting all boundary relations of simplices in Σ. We compare the Stellar tree to the implementation of the Simplex tree provided in [GUDHI, 2018], where each node of a Simplex tree requires a reference to the label of the vertex and three references to the tree structure (pointers to the parent node, to the first child and to the next sibling node) for a total of 4|Σ| references.
Note that the Stellar tree and our extended IA * data structure can both represent CP complexes in arbitrary dimension and, thus, have the same expressive power, while the Simplex tree can represent only simplicial complexes. Another difference is that Stellar trees require the complex to be embedded in an ambient space A, while the other data structures are purely topological and do not require a spatial embedding. We note, however, that while this is a requirement for Stellar trees, it is not a requirement for the more general Stellar decomposition.
In terms of storage requirements, the Stellar tree is always   Datasets marked with ⊗ could not be directly generated on our test machine using the standalone IA * . more compact than the IA * data structure, requiring approximately half of the storage, nearly all of which is used for encoding boundary relation R k,0 for the top cells (i.e., the indexed representation that they share in common). It is worth noting that we were unable to directly generate the IA * data structure for several of our larger datasets on our 64 GB test machine. We generated the IA * data structure on these datasets indirectly using our Stellar tree representation (see Section C in the Appendix) and we have marked these datasets with an ⊗ in Figures 14 and 15. When comparing the Stellar tree to the Simplex tree, we observe that the Stellar tree is significantly more compact: by an order of magnitude on manifold and pure models, and by two orders of magnitude or more on non-manifold models. Here too, we were unable to generate Simplex trees for several of the higher dimensional models on our test machine. For these datasets (marked with in Figure 14), we estimated the storage requirements based on the number of simplices of each dimension in the model. On two of these datasets, PROB 40D and LUCY 34D, we were unable to extract all simplices in all dimensions (even indirectly, see Section 9.1), and thus, the storage shown in Figure 14 is a lower bound of the real storage requirements. In contrast, we had no difficulty generating Stellar Trees for any of our test datasets.
For our dimension-dependent comparisons on manifold simplicial complexes, we also considered the Corner Table  (CoT ) [Rossignac et al., 2001] and the Sorted Opposite Table  (SOT ) [Gurung and Rossignac, 2009] data structures, both defined only for manifold triangle and tetrahedral complexes. The CoT data structure is similar to the IA data structure and explicitly encodes boundary relation R d,0 (σ) and adjacency relation R d,d (σ) of each top d-simplex σ. The SOT extends the CoT by implicitly encoding boundary relation R d,0 (σ). It only explicitly encodes adjacency relation R d,d (σ).
When comparing the Stellar tree to corner-based data struc-   tures, we observe that the CoT data structure has similar storage requirements as the IA and is roughly twice as large as the Stellar tree, while the SOT has similar storage requirements as the Stellar tree, requiring about 1% to 10% less space.
Finally, we consider the effect of different bucketing threshold on the size and efficiency of the Stellar tree representation. For our experimental datasets, there was only about a 10% difference in storage requirements between the large (k L ) and small (k S ) bucketing factors. Clearly, this is not always true, especially in the limit cases, i.e., with k V = 1 and k V = ∞. Very low bucketing thresholds (with k V near 1) yield deeper trees whose leaf blocks index only a few entities, leading to a high topological overhead but more efficient execution for individual mesh processing operations. Conversely, really large bucketing threshold values lead to lower storage overhead at the expense of increased query and execution times for individual operations. At the limit, when k V = ∞, the Stellar tree is effectively identical to the indexed representation.
These results confirm that the Stellar tree can efficiently represent both low-and high-dimensional complexes with only a slight storage overhead relative to that of the indexed base complex. This is largely due to the Stellar tree's exploitation of the complex's spatial locality via SRE compression.

Evaluation of Stellar tree generation times
In this section, we evaluate the generation times for the Stellar tree. Table 4 shows the timings of the four generation phases and the overall total timings. The two insert columns show the time for creating the base indexing structure H over the vertices Σ V of the complex Σ, or the time for inserting the top cells Σ T into H, while reindex columns show the timings for reordering and SRE compressing the indexed lists and arrays in H and Σ.
We first consider the relative cost of each of the generation phases. In general, the vertex reindexing phase consumes less than 10% of the overall timings. For the triangle, quadrilateral, hexahedral complexes, and the lower dimensional Vietoris-Rips complex, generating H is the most expensive phase, while for the tetrahedral, probabilistic-refinement and the two higher dimensional Vietoris-Rips models, reindexing the top cells is the most expensive phase. These results can be understood by considering the relative sizes of Σ V and Σ T . When the number of vertices is greater than or equal to the number of top cells, it is more expensive to generate the spatial hierarchy H. Otherwise, reindexing and compressing the top cells arrays dominates.
Finally, considering the effect of the bucketing thresholds (k V ) on generation times, we find that Stellar trees with higher bucketing thresholds (k L ) can be generated in less time than those with lower bucketing thresholds (k S ). This is expected since high values of k V tend to produce coarser spatial subdivisions with lower average spanning numbers χ.
Algorithm 2 EXTRACT P CELLS (p,r,Σ) Input: p is the cell dimension to extract Input: r is a leaf block in H Input: Σ is the CP complex indexed by H Variable: m p maps a p-cell vertex tuple to its local index Require: Extract boundary p-cells of top k-cells, for all p-faces τ in R k,p (σ) (with face index i τ in σ) do // Rearrange τ 's vertices into a canonical order 3: v tuple ← CANONICAL TUPLE(R p,0 (τ )) // If τ is indexed by r, add it to the local p-faces map 4: if there exists v ∈ R p,0 (τ ) such that v ∈ Φ V ERT (r) then // Insert τ as a new p-cell, if not already present

5:
if v tuple is not in m p then 6: Topological queries on a Stellar tree In this section, we describe how to perform topological queries on a CP complex Σ in the Stellar tree representation. These queries are the fundamental building blocks for locally traversing and processing the underlying complex.
Since these queries often depend on all cells in the complex, rather than just the explicitly represented top cells, we first describe how we obtain and represent the implicitly encoded boundary cells of the complex from the Stellar tree representation (Section 9.1). We then present algorithms for extracting the coboundary (Section 9.2). For brevity, we omit a description of how to extract adjacency relations, but in the Appendix C we describe how to extract the R d,d adjacency relations to generate the IA * data structure from a Stellar tree.

Extracting boundary relations
The Stellar tree's underlying indexed representation of a CP complex Σ explicitly encodes only the vertices and top CP k-cells of Σ for k ≤ d (see Section 4.2.1). However, many applications require access to non-top cells within the complex. Since they are implicitly encoded within the Stellar tree representation, we must create a local (explicit) representation to support algorithms for processing and attaching data to such cells.
Our strategy for extracting all p-cells is to iterate through the top k-cells of a leaf block for each dimension k, 0 < p ≤ k ≤ d and to extract an ordered set of p-cells (see Algorithm 2). We use an associative array m p to track the unique set of encountered p-cells with at least one vertex indexed by r (row 4). Array m p maps the tuple of vertices for a p-cell τ to an integer index id τ in the set, accounting for changes in ordering and orientation through the CANONICAL TUPLE routine (row 3). In some applications, it is useful to also explicitly maintain the boundary relation R p,0 for the p-cells and/or the incidence relations R k,p or R p,k for the top k-cells. These are encoded using the local indices within the ordered set of extracted p-cells.
We note that, for truly high-dimensional datasets, it is not feasible to extract p-cells in all cases. For example, there are 41 21 20simplices within a single 40-simplex. Encoding these 269 billion simplices would require more than 40TB of storage. However, Table 5: Summed timings (seconds) and additional storage requirements (number of references) to extract boundary p-cells from Stellar tree, IA * and Simplex tree data structures. Datasets marked with an ⊗ could not be directly generated on our test machine by the IA * . even on these datasets, we can still extract the lowest and highest dimensional p-cells. This highlights an advantage of only encoding the top cells of the complex (as in the Stellar tree and IA * data structure) compared to representations that encode all cells of the complex (as in the IG or Simplex tree). Stellar trees have no difficulty encoding and processing such high-dimensional complexes, despite the combinatorial explosion in the number of overall cells.

Experimental results
We now analyze the effectiveness of the Stellar tree representation for (batched) p-cell extractions against our implementation of the IA * data structure and the Simplex tree (as implemented in the GUDHI framework [GUDHI, 2018]). Table 5 lists the aggregate times and storage requirements for extracting all non-top p-cells from our experimental datasets. Notice that we do not consider the higher dimensional probabilistic dataset and the LUCY 34D V-Rips complex, as extracting all pcells on these datasets is unfeasible due to its computational and storage requirements. First, we analyze the influence of the bucketing threshold k V for Stellar trees. Smaller k V values lead to faster extractions on all our experimental datasets. This speedup increases with the dimension of the complex since the auxiliary data structure encoding a p-face type becomes smaller, and thus, checking for the presence of duplicates has a lower computational cost.
The IA * data structure follows a similar strategy as the Stellar trees for extracting its implicit p-cells since both data structures use an indexed representation for encoding the boundary relations of a CP complex. Table 5 demonstrates the computational and storage advantages of the Stellar tree over the IA * for this task. Namely, Stellar trees require from 20% to 55% less time for the two-dimensional datasets and approximately 10% less time on the higher dimensional ones. Notice, however, that the IA * data structure is a global data structure over the entire complex and runs out of memory (OOM) on our hexahedral datasets and on the 7D probabilistic and FOOT 10D V-Rips datasets. In addition, the Stellar tree's auxiliary storage requirements are negligible compared to those of the IA * data structure.
The Simplex tree explicitly encodes all simplices of a simplicial complex, thus, its p-cells can be enumerated by traversing all simplices at the p-th level of the tree. Explicitly encoding boundary relation R p,0 would require the same auxiliary storage as the IA * data structure, since both data structures require global structures. Table 5 demonstrates that Stellar trees are slower than Simplex trees at boundary cell extraction, but, still, competitive with respect to a representation that explicitly encodes all cells. This is possible thanks to the smaller local auxiliary data structures used by Stellar trees. Note that the Simplex tree runs out of memory (OOM) on our workstation for the 7D probabilistic dataset and the FOOT 10D V-Rips complex. Since a Simplex tree can only represent simplicial complexes, it does not support pcell extraction on quad and hexahedral datasets.

Extracting co-boundary relations
Co-boundary queries arise in a variety of mesh processing applications, including those requiring mesh simplification and refinement [Garland and Heckbert, 1997, Natarajan and Edelsbrunner, 2004, Zorin, 2000, or the dual of a complex [Hirani, 2003, Mullen et al., 2011, Weiss et al., 2013.
Co-boundary queries are naturally supported by the Stellar decomposition model. By definition, all regions of the decomposition that contain at least one vertex of a CP cell τ must index all CP cells in the star of τ (see Equation 1). Since the top cells are explicitly represented in Σ, we first describe how to extract the vertex co-boundary relation R 0,k restricted to the top k-cells of Σ, which we will refer to as the restricted co-boundary relation R 0,k . We will then discuss how to extend this to extract vertex co-boundary relation R 0,p over all p-cells in Σ, and the general co-boundary relation R p,q with 0 ≤ p < q ≤ d.
The restricted vertex co-boundary relation R 0,k in a leaf block r is generated by inverting boundary relation R k,0 on the top CP k-cells in Φ T OP (r). Since the indexed vertices in the leaf blocks of a COMPRESSED Stellar tree are contiguous, with indices in the range [v s , v e ), we encode our local data structure using an array of size |Φ V ERT (r)| = v e − v s . Each position in the array corresponds to a vertex indexed by r and points to an (initially empty) list of indexes from Σ T . As shown in Algorithm 3, we populate these arrays by iterating through relation R k,0 of the top CP k-cells in Φ T OP (r). For each cell σ such that relation R k,0 (σ) contains a vertex v with index i v ∈ [v s , v e ), the index of σ is added to vertex v's list.
Extending the vertex co-boundary relation to all p-cells in r is complicated by the fact that we only have an explicit representation for the top cells in Σ. A simple strategy we have developed for extracting R 0,p on all p-cells in r is to first extract the explicit Algorithm 3 EXTRACT RESTRICTED VERTEX CBDRY(r,Σ) Input: r is a leaf block in H Input: Σ is the mesh indexed by H Variable: r 0 k encodes R 0,k relation for the vertices in r Ensure: set of all p-cells in r, as in Algorithm 2 (see Section 9.1). We then invert R p.0 to obtain the complete relation R 0,p for the vertices in r.
In some applications, we prefer to express R 0,p entirely in terms of top cells from Σ. Thus, another strategy we have developed is to extract the restricted co-boundary relation R 0,k for all top k-cells in r, with p ≤ k ≤ d. This redundant representation is thus used as an intermediate representation for R 0,p (v) since each k-cell in R 0,k (v) contains one (or more) p-face in the co-boundary of v. For example, this provides a convenient representation for the star of a vertex v as a union of restricted coboundary relations Similarly, we have defined and implemented a strategy for generating the general co-boundary relation R p,q , where p < q. First, the sets of all q-cells, which is expressed as R q,0 , is extracted. This implicitly provides also boundary relation R q,p . Then, coboundary relation R p,q is extracted by inverting R q,p .

Experimental results
We now analyze the effectiveness of the Stellar tree representation for co-boundary extractions. Specifically, since the main co-boundary extraction in our applications (see Section 10) is the restricted vertex co-boundary relation and most of the other co-boundary extractions can be posed in terms of this primitive extraction, we compare the performance of the Stellar tree against our implementation of the IA * data structure for this query and against the Simplex tree. Table 6 lists the extraction times and storage requirements for the vertex coboundary relation R 0,d on our manifold (triangular, quad, tetrahedral and hex) and pure (probabilistic) complexes and the sum of extraction times for the restricted vertex co-boundary relations R 0,k for each dimension k with top cells on our non-manifold (V-RIPS) complexes.
We first consider the influence of the bucketing threshold k V for Stellar trees. While there is not much difference in extraction times for the two-dimensional complexes, larger k V values lead to faster extractions for three-dimensional and non-manifold datasets in most cases. While this comes with a slight increase in storage requirements for encoding the relation (see right column in Table 6), the overall storage cost per block is pretty low, requiring at most a few megabytes for the probabilistic models, and a few kilobytes in all other cases.
The IA * data structure extracts co-boundary relations through a traversal along the face adjacencies of its top cells (encoded in the R k,k adjacency relation). The traversal for a given vertex v is seeded by one top k-cell per k-cluster (encoded by partial relation R * 0,k (v), see Section 8.2; we refer to [Canino et al., 2011] for more details). Since each such traversal is run on demand, there is a negligible memory impact for this query. Table 6 demonstrates that Stellar trees are significantly faster at extracting R 0,k relations, which can be performed in about one tenth of the time in most cases. However, it is important to note that the Stellar tree extraction is batch-based (by leaf blocks of H), and individual co-boundary extractions would likely be faster on the IA * data structure. The Simplex tree extracts co-boundary relations through a traversal of the underlying trie. Given a vertex v, the procedure for extracting its restricted co-boundary first identifies the simplices incident in v (i.e., its star), and then extracts just the top simplices from the star. The former requires a trie traversal, with a worst-case complexity linear in the number of nodes in the trie, since, as stated in the GUDHI documentation [GUDHI, 2018], this corresponds to a depth-first search of the trie starting from the node with value v. Identifying the top simplices in the star of a vertex has a negligible cost on low dimensional meshes, while it becomes a costly operation on higher-dimensional ones, where it accounts for nearly 50% of the overall extraction time.  Figure 17: Extraction times (in seconds) for the restricted vertex co-boundary relations. The top dataset is the triangle mesh used in our main comparison, the second is a tetrahedral mesh with 256 thousand vertices and 1.4 million tetrahedra, and the last dataset is a probabilistic-refinement CP complex with 7-dimensional top simplices.
the IA * , since this traversal is done on demand, this query imposes negligible memory impact. On our experimental datasets, the Simplex tree was able to complete the extraction of restricted vertex co-boundary relations only on the smaller triangle mesh NEPTUNE, for which it requires nearly 72 hours. To provide a comprehensive performance comparison against the Stellar tree, we consider two additional smaller datasets for this query: a tetrahedral mesh (FIGHTER2) with 256 thousand vertices and 1.4 million tetrahedra, and a probabilistic-refinement CP complex with six thousand vertices and two million top 7-simplices. The results, shown in Figure 17, highlight the Stellar tree's significant advantage over the Simplex tree for restricted vertex co-boundary extraction (i.e., less than a second vs hours).
10 A brief tour of applications in the Stellar universe Stellar decompositions and Stellar trees have been successfully applied in several mesh processing applications. In this section, we provide a high-level overview of several such applications over large CP complexes with a focus on how Stellar trees uniquely benefit the application. As we will describe, each such application utilizes local topological data structures designed for the underlying application. Due to the streamed processing approach discussed in Section 6, the storage requirements for these data structures are proportional to the geometry indexed within a leaf block of the tree and the generation costs are amortized over all processed cells in the block.

Validation of geometric and topological properties
Many popular topological mesh data structures are valid only for a restricted class of complexes due to assumptions they exploit in their encodings, such as the cardinality of the adjacency relation among top cells. For example, popular edge-based and adjacency-based data structures, such as the half-edge [Mantyla, 1988, CGAL, 2020, OML, 2015, Corner-  SOT [Gurung and Rossignac, 2009] and IA [Paoluzzi et al., 1993, Nielson, 1997, require the underlying complex to be pseudo-manifold.
While one can verify such topological conditions using local checks on the star or link of the vertices of the complex, it can be infeasible to reconstruct such relations on large meshes without the aid of an efficient topological data structure. On the other hand, global approaches that directly build the required relations do not scale to larger complexes.
In contrast, Stellar trees are ideally suited to verify topological properties of large CP complexes, even in memory-limited environments, since each leaf block of the Stellar tree only requires a list of vertices and their incident top cells (i.e., those in the star of the vertices). A simple local topological verification operation was utilized in [Weiss et al., 2011] to mark boundary vertices of a tetrahedral mesh by checking properties of its link, such as its Euler number. This was extended in [Fellegara et al., 2016b] to a full suite of topological checks on a CP complex, implemented using global Stellar tree traversals. In particular, a graph traversal of the 1-skeleton was used, in conjunction with a Union-Find data structure [Tarjan, 1975], to count the connected components of a pure simplicial complex. A similar traversal of the 1-skeleton of its dual complex (i.e., the graph of the d-adjacency relation) was used to verify the d-connectedness of the complex and whether it is pseudo-manifold. Simplified checks for combinatorial manifolds (when applicable) were then performed on the links of each vertex to check that they were locally homeomorphic to (d−1)spheres (for internal vertices) or to (d−1)-balls (for boundary vertices).

Topology-preserving simplification
One of the earliest applications of the Stellar tree (actually, its predecessor, the PR-star octree [Weiss et al., 2011]) was to accelerate a mesh simplification algorithm for tetrahedral meshes based on edge collapses. An edge collapse is a local topological operation defined in terms of the stars of an edge's two vertices. This operation identifies the pair of vertices along an edge, removes all tetrahedra in the star of that edge and updates the mesh connectivity within this local region [Cignoni et al., 2004a]. Edge collapses are valid when they satisfy the so-called link con-dition [Dey et al., 1999], which consists of local checks on the links of the edge and its vertices.
The simplification procedure in [Weiss et al., 2011] was implemented as an iterative process that alternated between: (i) streaming through the Stellar tree blocks, where it collapsed eligible candidates, and (ii) rebuilding the Stellar tree's index over the simplified mesh.
It used discrete distortion [Mesmoudi et al., 2008] and a quadric error metric [Garland and Heckbert, 1997] to organize eligible edges into a priority queue. Applying this simplification algorithm to the leaf blocks of a Stellar tree rather than to the entire mesh provides a significant space savings due to the reduced sizes of the edge queues. In many cases, the simplification was 10-50% faster and required only 0.1% of the additional memory for auxiliary data structures. Moreover, this speedup was more pronounced as the mesh size increased.
More recently, Stellar trees have been used to perform homology-preserving edge-contractions on general simplicial complexes of arbitrary dimension [Fellegara et al., 2020b]. The core edge contraction algorithm was implemented using a custom local topological data structure over the top star of the vertices and edges within each leaf block of the tree, similar to the restricted co-boundary relations R 0,k and R 1,k (c.f. Section 9.2). To avoid regenerating these topological relations, it utilized an LRU cache for the expanded leaf blocks as it traversed the tree. In this mesh simplification application, Stellar trees were applied to datasets with dimensions up to 70 and were able to remove more than 90% of the simplices of the mesh, significantly reducing the dimensionality of the complex while preserving important topological invariants of the dataset. Compared to existing state-of-the-art data structures for edge contractions [Attali et al., 2012], Stellar trees were able to simplify complexes using comparable or less runtime and memory in all cases, and requiring significantly less memory and/or processing time in several cases. Notably, in one case, the Stellar tree was able to successfully complete the simplification process in less than 30 minutes, while [Attali et al., 2012] did not complete after more than 24 hours. Figure 18 shows simplified versions of the genus-3 NEPTUNE triangle mesh. Each simplified mesh preserves the homology of the complex, such as its number of connected components, loops and cavities. We note that the above application incorporates only topological considerations into its simplification error metric. Incorporating mesh quality considerations into the error metric could significantly improve the mesh quality [Garland and Heckbert, 1997].

Shape analysis and morphological feature extraction
While topological validation and simplification operations can be implemented in terms of local operations on the star of a vertex, shape analysis applications, such as watershed analysis [Roerdink and Meijster, 2000] and visibility queries on terrain datasets [Bittner andWonka, 2003, De Floriani andMagillo, 2003] often require algorithms that are seeded locally and span vast interwoven regions of the complex. This section discusses how Stellar trees have aided in the generation and simplification of the discrete Morse gradient field and of the associated Morse and Morse-Smale complexes of triangulated terrains [Fellegara et al., 2014] and of tetrahedralized volumetric data [Weiss et al., 2013]. The discrete Morse gradient field [Forman, 1998] is composed of arrows (ordered pairs) between incident cells of the complex and can be computed locally using scalar values associated with cells incident in the star of a vertex [Robins et al., 2011]. Since the encoding of [Weiss et al., 2013] compactly encodes the discrete Morse gradient field as a scalar field on the top simplices of the complex, the latter can be computed using a local traversal of a Stellar tree's blocks. Compared to an IA implementation, Stellar trees were able to extract the discrete Morse gradient of scalar fields defined over tetrahedral meshes in about half the time (see [Weiss et al., 2013] for details). While the experiments in [Weiss et al., 2013] were performed on VERTEX-COMPRESSED Stellar trees, yielding a 30% storage savings over the IA data structure, a COMPRESSED Stellar tree encoding would likely provide a 50% total memory savings while maintaining similar performance improvements.
Extracting the Morse complex from a Stellar tree-based encoding is more complicated since it involves traversing the directed acyclic graphs (DAGs) induced by the discrete Morse gradient field's arrows. Specifically, in the encoding of [Weiss et al., 2013], each k-dimensional critical point of the discrete Morse gradient field corresponds to a k-cell of the ddimensional Morse complex. The disjoint regions of influence of each such critical point, referred to as the k-manifolds of the Morse complex, are extracted by traversing the DAG rooted at a given critical point of the discrete Morse gradient field. Since each such graph traversal can visit the blocks of a Stellar tree multiple times, an LRU cache was used in [Weiss et al., 2013] to support global extraction algorithms for each k-manifold of the Morse complex.
Further, since the extraction algorithm for each dimension's manifolds depends on different topological connectivity relations, Morse complex extraction benefits from the Stellar tree's ability to generate customized local topological data structures. For example, since extracting the 2-manifolds from a tetrahedral complex requires only the R 2,2 adjacency relation, its extraction was optimized by directly starting from R 2,2 , rather than the R 3,3 adjacency relation, as in the IA data structure.
This approach was extended to terrain datasets in [Fellegara et al., 2014], which also introduced a persistencebased simplification algorithm for noise removal. Figure 19 highlights the 73K triangles in the largest 2-manifold of the Morse complex for the MAUI terrain dataset (in red) and the blocks of the Stellar tree that were visited when extracting this region (blue-green squares).
While manifold extraction and persistence-based simplification operations were slightly more expensive than their IA counterparts for volumetric [Weiss et al., 2013] and terrain [Fellegara et al., 2014] datasets, the Stellar tree's vast memory savings and hierarchical encoding open the door to efficient parallel implementations on huge datasets, which we hope to explore in future work.

Concluding remarks
We have introduced the Stellar decomposition as a model for topological data structures over Canonical-Polytope (CP) complexes, a class of complexes that includes simplicial complexes and certain classes of cell complexes, like quadrilateral and hexahedral meshes. Stellar decompositions cluster the vertices of a complex into regions that contain sufficient information to locally reconstruct the star of their vertices. The model is agnostic about the domain of the complex (e.g., manifold, pure, non-manifold) and we have demonstrated the scalability of this model to large mixed-dimensional datasets in high dimension.
We introduced the Stellar tree as a concrete realization of the Stellar decomposition model over spatially embedded CP complexes. Stellar trees couple a spatial index H decomposing the complex's embedding space with a simple tuning parameter that limits the number of vertices indexed by a leaf block.
Stellar trees effectively exploit the spatial coherence of a CP complex Σ by using the clustering structure of H to reorder the arrays of top cells of Σ and to compress the resulting ranges of sequential indexes within the lists of vertices and top cells in the leaf blocks of H. We have demonstrated over a wide range of datasets that this process typically produces COMPRESSED Stellar trees that are only 1-10% larger than the original indexed base mesh for Σ while still retaining sufficient information to efficiently reconstruct all topological connectivity relations. The source code for our Stellar tree implementation will be released in the public domain.
In terms of storage size, Stellar trees compare quite favorably with state-of-the-art topological data structures. They are consistently half the size of their IA * data structure counterparts [Canino et al., 2011] and one to two orders of magnitude smaller than their Simplex tree counterparts [Boissonnat and Maria, 2014]. This is especially notable for high dimensional Vietoris-Rips complexes, a target application for the Simplex tree, for which Stellar trees have very low overhead. While Stellar trees support a much broader class of complexes, they have similar storage requirements as the dimensionspecific SOT data structure Rossignac, 2009, Gurung andRossignac, 2010], which supports only static pseudo-manifold triangle and tetrahedral complexes. In future work, it would be interesting to compare the Stellar tree against top-based extensions of the Simplex tree, such as the M ST and the SAL [Boissonnat et al., 2017], if public-domain implementations become available.
Despite the simplicity of their leaf block representation, Stellar trees provide a great deal of flexibility to customize the structure and layout of their expanded topological data structures to meet the needs of a given application. Such data structures are typically constructed by composing several local topological incidence and adjacency relations. We described efficient algorithms for reconstructing these relations within the subcomplex indexed by the leaf blocks of a Stellar tree and demonstrated the advantages of this approach compared to similar algorithms on the IA * and Simplex tree data structures. Stellar trees can also be used as an intermediary representation to generate topological data structures in a memory-constrained environment. For example, we used Stellar trees to generate IA * and Simplex tree representations for several of our larger complexes in Section 8 (as we discuss in the Appendix C). We also provided an overview of several mesh processing applications, ranging from mesh validation, to topology and shape preserving simplification and morphological analysis that have benefited from the Stellar trees representation.
One direction of future work would involve extending the Stellar tree representation to support a broader class of cell complexes. For example, it would not be difficult to extend support to indexed polyhedral cell complexes which define their cells in terms of their boundary polyhedral faces which are, in turn, defined by oriented lists of vertex indices [Muigg et al., 2011].
Another avenue for investigation is to extend our processing algorithms for parallel, distributed and/or out-of-core environments, which could be used for applications like multicore homology computation [Lewis and Zomorodian, 2014] on point cloud data. The Stellar tree's compact leaf block representation is already geared towards a parallel execution pattern since each block already has sufficient resources to query the connectivity of its local subcomplex. Preliminary results along this line look promising. A simple unoptimized OpenMP [OpenMP, 2015] adaptation of boundary and restricted vertex co-boundary queries yielded a 3-4x speedup compared to our serial approach on our 6 core machine.
Finally, while Stellar trees require their underlying complex to be spatially embedded, there is no such restriction on the Stellar decomposition model. Thus, we plan to investigate Stellar decompositions for abstract CP complexes, such as simplicial complexes representing social networks. Social network representation and processing poses new challenges in the social big data domain, such as the identification of key-players and communities in the dataset, as well as extracting topological properties of the network, like its homology or k-connectivity. Due to the irregularities of non-spatial datasets, one key challenge would be to define efficient decompositions (i.e., with a low average spanning number χ) using only the complex's connectivity information. A preliminary attempt for geolocalized social networks can be found in [Fellegara et al., 2016a], where the social network was represented in terms of its maximal cliques, i.e., sets of mutually related entities, corresponding to top simplices over the network's flag complex. The Stellar tree was built over the 2D embedding provided by the geospatial locations of the entities and simplified using homology-preserving edge-contractions [Fellegara et al., 2020b], enabling a study of the network's topological structure on a significantly reduced dataset.
As it can be noted, the generation process ignores how the partitioning on the vertices is obtained, since this step is defined by the data structure instantiating a Stellar decomposition. The reindexing of the vertices follows a traversal of the regions of ∆ in such a way that all vertices mapped to a region have a contiguous range of indices in the reindexed global vertex array Σ V (as detailed in A.1). Figure 20 illustrates a reindexing of the vertices of a triangle mesh in the plane into a decomposition defind by four rectangular regions (dashed lines).
We then insert each top CP k-cell σ, with index i σ in Σ T k , into all the regions of ∆ that index its vertices. This is done by iterating through the vertices of σ and inserting i σ into the r T array of each region r whose vertex map Φ V ERT (r) contains at least one of these vertices. As such, each top CP k-cell σ appears in at least one and at most |R k,0 (σ)| regions of ∆. Due to the vertex reindexing of step 1, this operation is extremely efficient. Determining if a vertex of a given cell lies in a block requires only a range comparison on its index i v .
Finally, we reindex the top CP cell arrays Σ T to better exploit the locality induced by the vertex-based partitioning and compress the local r T arrays using a sequential range encoding over this new index. The reindexing and the compression of the top CP cells is obtained following a traversal of the regions of ∆ in such a way that all top CP cells mapped from the same set of regions have a contiguous range of indices in the reindexed arrays Σ T . This last step is detailed in A.2 and in B. As we demonstrate in the paper, this compression yields significant storage savings in a wide range of mesh datasets.

A.1 Reindexing and compressing the vertices
After generating the partition ∆ and the vertex map Φ V ERT for the Stellar decomposition, we reindex the vertex array Σ V to better exploit the coherence induced by ∆. At the end of this process, each region of ∆ has a consecutive range of indices within the global vertex array Σ V , and thus it trivially compresses under This reindexing procedure is organized into three major steps, as outlined in Algorithm 4. The first step performs a traversal of the regions, which generates new indices for the vertices in Σ. For a region r, it generates a contiguous range of indices for the vertices in r. For example, in Figure 20, after executing step 1 on region b, we have v s = 4 and v e = 7.
The new indexes are then incorporated into mesh Σ by updating the vertex indices in R k,0 relations for all top k-cells in Σ T k (see step 2 of Algorithm 4) and then permuting the vertices (detailed in step 3 of Algorithm 9 in B). These updates take place in memory without requiring any extra storage.

A.2 Reindexing and compressing the top CP cells
After inserting the top CP cells of Σ T into ∆, we reorder the top CP cells array Σ T based on the partitioning and apply SRE compaction to the region arrays to generate our COMPRESSED encoding. This reindexing exploits the coherence of top CP cells that are indexed by the same set of regions, translating the proximity in ∆ into index proximity in Σ T . This procedure is organized into four main phases, as shown in Algorithm 5. A detailed descrip-  The EXTRACT TUPLES procedure (see Algorithm 6 in B), traverses the regions of ∆ to find the tuple of regions t b = (r 1 , . . . , r n ) in ∆ that index a top CP cell σ. Inverting this relation provides the list of top cells from Σ mapped to each such tuple of regions. As we iterate through regions, we ensure that each top CP cell in the complex is processed by only one region r, by skipping the top CP cells whose minimum vertex index i v is not in Φ V ERT (r). For example in Figure 21(a), triangle 5 is indexed by region a and b, and, thus, its tuple is t b = (a, b). The complete list of triangles in tuple (a, b) is {2, 5, 12}.
We use this inverted relation in EXTRACT CELL INDICES (see Algorithm 7 in B), to generate a new coherent order for the top CP cells of Σ T . Specifically, the prefix sum of the tuple cell counts provides the starting index for cells in that group. For example, when considered in lexicographic order, the first three region tuples, (a), (a, b) and (a, b, d) in Figure 21(b), with 1, 3 and 1 triangles, respectively, get starting indices 1, 2 and 5. We then assign increasing indices to the top CP cells of each group. Thus, e.g., the three triangles belonging to tuple (a, b) get indices {2, 3, 4} after this reindexing.
Finally, in COMPRESS CELLS and PERMUTE ARRAY (Algorithm 8 and 9 in B), we reorder and SRE-compact the r T region arrays and the global top CP cells array Σ T .
Experimental results In Table 7, we compare the time, and storage requirements to generate an IA or IA * data structure, depending on whether the complex has a manifold or non-manifold domain, using the Stellar tree or directly extracting it from the indexed representation. For each dataset, we compare the Stellar trees generated by using thresholds k S and k L and by using a local and a global algorithm against the direct approach on the original indexed representation of the complex. For the manifold (triangular, quadrilateral, tetrahedral and hexahedral) and pure (probabilistic) datasets, where all top cells have dimension d, we used Algorithm 10 to compute the adjacencies.
When comparing execution times, we find that the global Stellar tree approach is about 25% faster than the direct approach in most cases. However, due in part to the redundant lookups in the adjacency calculation, the local approach is slightly slower than the global approach, but still 10% faster than the direct approach in most cases. For example, it is almost twice as fast on F16, on par on Lucy and slower on the 5D probabilistic dataset. Considering the effects of the bucket threshold k V , we observe little discernible difference on the global Stellar tree approach. However, a larger bucketing threshold (k L ) yielded up to a 25% speedup in the local approach on our larger datasets, compared to its smaller (k S ) counterpart.
Lastly, we consider the storage requirements for generating the IA / IA * data structure. For both the local and global Stellar tree tree approaches, the auxiliary storage requirements are limited to the complexity of each leaf block, requiring only a few KBs of auxiliary storage for the manifold and non-manifold datasets, and a few MBs for the pure (probabilistic) datasets. In contrast, the direct approach requires hundreds of MBs for the medium sized datasets. We were not able to generate the IA * data structures using the direct approach on our largest datasets, which ran out of memory (OOM) on our workstation, despite its 64 GB of available RAM. Table 7: Generation times (seconds) and storage (number of references) for the IA * data structure from Stellar trees (k S and k L ) and the direct approach (dir.) on the original indexed representation of the complex. With the exception of V-RIPS complexes, the IA * is equivalent to the IA representation on these datasets. OOM --