Compact Representations of Simplicial Meshes in Two and Three Dimensions (cid:3)

.


Introduction
For many applications the space required to represent large unstructured meshes in memory can be the limiting factor in the size of a mesh.Standard representations of tetrahedral meshes, for example, can require 300-500 bytes per vertex.One option for using larger meshes is to maintain the mesh in external memory.To avoid thrashing, this requires designing algorithms for which the access to the mesh is carefully orchestrated.Although several such external memory algorithms have been designed [18,13,11,28,42,2,41,1], these algorithms can be much more complicated than their mainmemory counterparts, and can be significantly slower.
Another option for using larger meshes is to try to compress the representation within main memory.There has in fact been significant interest in compressing meshes [12,22,40,32,33,39,25,23,17].In three dimensions, for example, these methods can compress a tetrahedral mesh to less than a byte per tetrahedron [39]-about 6 bytes/vertex (not including vertex coordinates).These techniques, however, are designed for storing meshes on disk or for reducing transmission time, not for representing a mesh in main memory.They therefore do not support dynamic queries or updates to the mesh while in compressed form.
We are interested in compressed representations of meshes that permit dynamic queries and updates to the mesh.The goal is to solve larger problems while using standard random-access main-memory algorithms.In this paper we present data structures for representing two and three dimensional simplicial meshes.The data structures support standard operations on meshes including traversing among neighboring simplices, inserting and deleting simplices, and the ability to store data on simplices.For a class of well shaped meshes [30] with bounded degree these operations all take constant time.The precise definition of our interface is described in Section 4.Although our data structures are not as compact as those designed for disk storage, they still save a factor of between 5 and 10 over standard representations.
Our data structures are described in Section 5.They take advantage of the separator properties of well-shaped meshes [30] using recent results in graph compression [4,5].In particular our technique uses separators to relabel the vertices so that vertices that share a simplex are likely to have labels that are close in value.Pointers are then difference encoded using variable length codes.We use this technique to radially store the neighboring vertices around each vertex in 2D and around a subset of the edges in 3D.A query need only decode a single vertex in 2D or vertex and edge in 3D.For applications that need to generate new vertices, e.g., Delaunay refinement, we leave extra space in the label space and assign new labels based on the labels of the neighbors.
Section 6 describes an implementation of our data structure and Section 7 presents experimental results.The implementation uses about 5 bytes per triangle in 2D and about 7.5 bytes per tetrahedron in 3D when measured over a range of mesh sizes and point distributions.We present experiments based on using our representation as part of incremental Delaunay algorithms in both 2D and 3D.We use a variant of the standard Bowyer-Watson algorithm [9,43] and the exact arithmetic predicates of Shewchuk [38] for all geometric tests.We also present experiments based on a Delaunay refinement algorithm that removes triangles with small angles by adding new points at their circumcenters.All space is reported in terms of the total space including the space for the vertex coordinates and all other data structures required by the algorithm.The results for 1 Gbyte of memory are summarized as follows.
• We can generate a 2D Delaunay mesh with 110 million triangles (.47 Gbytes for the mesh, .44 Gbytes for the vertex coordinates, and about .1 Gbytes for auxiliary data used by the algorithm).
Compared to the Triangle code [37] (the most efficient we know of) our algorithm uses a factor of 3 less memory.It is about 10% slower than Triangle's divide-and-conquer algorithm and much faster than its incremental algorithm.
• We can generate a 3D Delaunay mesh with 100 Million tetrahedrons (.75 Gbytes for the mesh, .17Gbytes for the vertex coordinates, and .08Gbytes for auxiliary data).Compared to the Pyramid code [36], our algorithm uses a factor of 3.5 less memory, and is about 30% faster.
• We can generate a refined 2D Delaunay mesh with 80 million triangles with no angle less than 26%.This version dynamically generates new labels, and uses an extra level of indirection in our datastructure.
Our data structure can be used in conjunction with external memory algorithms.Also, although we only describe our implementation for 2D and 3D simplicial meshes, the ideas extend to higher dimensions.These topics are discussed, briefly, in Section 8.

Standard Mesh Data Structures
There have been numerous approaches for representing unstructured meshes in 2 and 3 dimensions.Some are specialized to simplicial meshes and others can be used for more general polytope meshes.For the purpose of comparing space usage, we review the most common of these data structures here.A more complete comparison for 2D structures can be found in a paper by Kettner [26].
In two dimensions most approaches are based on either triangles or edges.The simplest data structure is based on triangles.Each triangle has three pointers to the neighboring triangles, and three pointers to its vertices.Assuming no data needs to be stored on triangles or edges, this data structure uses 6 pointers per triangle.Storing data requires extra pointers.Shewchuk's Triangle code [37], and the CGAL 2D triangulation data structure [8] both use a triangle-based data structure.To distinguish the three neighbors/vertices of a triangle, a handle to a triangle typically needs to include an index from 1 to 3. The data structure used by Triangle, for example, includes such an index in the pointer to each neighbor (in the low 2 bits) so that a neighbor query not only returns the neighbor triangle, but returns in which of three orders it is held.
There are many closely related data structures based on edges, including the doubly connected edge list [31], winged-edge [3], half-edge [44], and quadedge [21] structures.In addition to triangulated meshes, these data structures can all be used for polygonal meshes.In these data structures each edge maintains pointers to its two neighboring vertices and to neighboring edges cyclically around the neighboring faces and vertices.Each edge might also maintain pointers to the neighboring faces and to edge data.The most space efficient of these data structures can maintain for each edge a pointer to the two neighboring vertices and to just two neighboring edges, one around each face and vertex.Assuming no data needs to be stored on a face or edge, this requires 4 pointers per edge, which for a manifold triangulation is equivalent to the 6 pointers per triangle used by the triangle structure (|E| = 3/2|T |).The half-edge data structure [44], used by CGAL [26], LEDA [29] and HGAM [19], maintains two structures per edge, one in each direction.These half-edges are cross referenced, requiring an extra two pointers per edge.The wingededge and quad-edge structures maintain pointers to all four neighboring edges, requiring 6 pointers per edge (9 per triangle).
In three dimensions there are analogous data structures based either on tetrahedrons or on faces and edges.Again the simplest data structure is to use a structure per tetrahedron.Each tetrahedron has 4 pointers to adjacent tetrahedrons, and 4 to its corner vertices.Assuming no data this requires 8 pointers per tetrahedron.This data structure is used by Pyramid [36] and CGAL [8].The face and edge data structures are often called boundary representations (b-reps).Such boundary representations are more general than the tetrahedron data structures, allowing the representation of polytope meshes, but tend to take significantly more space.Dobkin and Laszlo [14] suggest a data structure based on edge-face pairs, which in general requires 6 pointers per edge-face.For tetrahedral meshes this data structure can be optimized to 9 pointers per face (6 to the adjacent faces rotating around its 3 edges, and 3 to the corner vertices).This corresponds to 18 pointers per tetrahedron.Weiler's radial-edge representation [45], Brisson's cell-tuple representation [10], and Linehard's G-map representation [27] all take more space.
In summary, the most efficient standard data structures of simplicial meshes use 6 pointers per triangle in 2D and 8 pointers per tetrahedron in 3D.At least one extra pointer is required to store data on triangles in 2D or tetrahedrons in 3D.

Preliminaries
In this section we review some basic notions of combinatorial topology used in this paper.For a more detailed discussion the reader can refer to [15] and [34] among others.
An (abstract) simplicial complex K is a non-empty collection of finite sets which is closed under taking nonempty subsets.The elements of K are called simplices.The underlying set ∪K is called the vertex set and its elements are called vertices.The dimension of a simplex with d vertices is d − 1.The dimension of K is the maximum dimension among its simplices.A simplex τ is a face of a simplex γ iff τ ⊆ γ, and iff τ = γ we say that τ is a proper face of γ.We say that K is pure if every simplex is a face of a simplex of highest dimension.Let S be a subset of K. We call the collection of all simplices in S together with all their faces, Cl(S), the closure of S. The star of a simplex is the union of its superfaces, St(σ) = {γ : σ ⊆ γ}.The link of a simplex σ is the set of simplices in the closure of its star that do not intersect it, Let E be a mapping from the vertices of K to R m .We let |σ| denote the convex hull of the images of their vertices of σ under E, and let |K| = ∪ σ∈K |σ|.We say that |K| is an embedding of K iff for all simplices σ and τ it holds that |σ| ∩ |τ | = |γ| where γ is their maximum common face (which may be empty).We say that K is a d-manifold (with boundary) iff |K| is a d-manifold (with boundary).If K is a manifold of dimension d then the link of every (d − 2)-simplex is a cycle of edges and vertices (i.e., a 1-manifold).If K is a manifold with boundary, then the link of every (d − 2)-simplex is either a cycle or a path, i.e., a 1-manifold with or without boundary (see Figure 1 (a)).We will make use of this fact in our representation described in section 5.An ordering, s d , of a d-simplex, s d , is a total ordering of its vertices.An orientation, s d , of a simplex, s d , is a maximal set of orderings which are even permutations of each other1 .Every ordering s d on a simplex implies an orientation s d on the simplex, and every d-simplex, d > 0, has two possible orientations.The orientation s d of a simplex induces an orientation s d−1 on every d−1 subsimplex-i.e., for all s d−1 ∈ s d−1 there exists s d ∈ s d such that s d−1 is a prefix of s d .
For our purposes, a d-pseudomanifold is a pure dcomplex where every (d − 1)-simplex is contained in at most two d-simplices and where the dual graph is connected.The vertices of the dual graph are the dsimplices and the edges are the (d − 1)-simplices.A d-pseudomanifold is orientable if its d-simplices can be given orientations in such a way that when they meet at a (d−1)-simplex s, they induce opposite orientations on s.Every orientable d-pseudomanifold has two possible orientations, which can be specified by the orientation of one of its d-simplices.If K is a d-pseudomanifold then the link of every (d − 2)-simplex is a collection of disjoint cycles and/or paths (see Figure 1 (b)).
In this paper we will use the term simplicial mesh to refer to a pseudomanifold abstract simplicial complex with a given orientation.

Interface
In this section we present the interface for simplicial meshes that our data structure implements.It is a simplified version of an interface described in [6].The inter-face supports standard operations on meshes including, a mechanism to systematically traverse a mesh (e.g., reflect across a face, or rotate around a vertex), and for updating the mesh, including inserting and deleting simplices and associating data with simplices.The interface consists of four operations on ordered simplices (empty, up, down), and faces, and three operations on simplicial meshes (add,delete, and findUp).
Let s k = (v 1 , . . ., v k+1 ) be an ordered simplex.The empty operation creates an empty simplex empty() → ().The up operation adds a vertex v to s: Intuitively this means that it returns every (d − 1)-faces in the opposite orientation than the one given in s d .This can be easily implemented using the above operations.
Let M be a d-dimensional simplicial mesh.The add operation takes M and a highest dimension ordered simplex s d and returns a new mesh M that results from adding s d to M .We require that s d has consistent orientation with M .Note that when we add s d , we don't have to store all even permutations of s d -just storing s d is enough to determine the orientation of s d .The delete operation takes M and a highest dimension ordered simplex s d and returns the mesh M that results from removing s d from M .The findUp takes an ordered simplex s k , 0 ≤ k ≤ d, and M .If s k is not in M it returns null.Otherwise it returns an ordered simplex s d such that s d ∈ M and s k is a prefix of s d , or null if none exists.In the special case where k = d − 1, then there is at most one s d that can be returned.
In addition to the core interface, we also provide two operations to associate and retrieve data from simplices in a complex.The addData operation takes M , an ordered simplex s k , 0 ≤ k ≤ d and some user supplied data u.It associates u with s k in M .The operation findData takes M and an ordered simplex s k and returns the user data ud associated with s k in M .If there is no associated data then null is returned.
The interface as described can be used for most applications that traverse and update a simplicial mesh.Figure 2 gives an example of code that traverses a dsimplicial mesh with boundary and returns the (d − 1) boundary mesh.It recursively traverses the mesh in depth-first order storing flags on the d simplices when visited.Whenever a boundary d − 1 simplex is found (s d = null in the code), it is added to the output mesh.The code assumes the boundary is a simplicial mesh.

Data Structure
Here we describe our 2D and 3D data structures for simplicial meshes (simplicial orientable pseudo manifolds).We first describe uncompressed versions of the data structures and then describe how to compress them.Our data structures are based on storing the link for a set of (d − 2)-simplices.In 2D this is similar to the half-edge structure [44], and in 3D it is similar to the Dobkin and Laszlo [14] structure.We note, however, that all references are to vertex labels instead of pointers to other higher-dimensional simplex structures, allowing us to compress based on vertex labels.Our data structures have the property that if the degree of all vertices is bounded all queries take constant time.We first describe a version for manifold complexes.Our 2D data structure maps each vertex to its link, represented as a cycle of the labels of its neighboring vertices.The cycle is ordered radially around the vertex in the orientation of the complex, e.g., clockwise.A findUp query on the ordered edge (v 1 , v 2 ) can be answered by looking up the link for v 1 , finding v 2 in the link, and returning the next vertex in the link.A findUp on a vertex can be answered by selecting the first two vertices off of its link.
The link can be stored as a list of labels starting at an arbitrary point on the cycle.If the vertex has bounded degree, the lookup takes constant time.To analyze the space note that each edge appears in two cycles, and each appearance requires two pointers, one to the vertex label and one to the next element in the list.The total space is therefore 4 pointers/edge + 1 pointer/vertex.This is identical in space usage to the triangle-based structure, assuming that it also maintains a pointer from each vertex to one of its incident triangles.Our data structure is similar to the half-edge structure since there are effectively two structures per edge, one pointing in each direction.It differs, however, in that there are no direct cross pointers between the matching half edges.
In 3D the data structure maps a subset of all ordered edges to their link, represented as a cycle of vertex labels.The cycle is maintained in a consistent orientation, e.g., obeying a right-hand rule with respect to the order (direction) of the edge.The representative subset E is selected to include only the edges {v 1 , v 2 } for which either the labels of v 1 and v 2 are both odd, or they are both even.Furthermore an edge is only stored in one of its two orders, chosen using a fixed rule, e.g., lower labeled vertex first.Since for any triangle (2simplex) at least two labels have to be either odd or even, this sampling of the edges guarantees that every triangle has at least one representative ordered edge in E .The data structure also needs to supply a way to access the link given the vertex labels of any edge in E .This can be implemented using an adjacency list for each v ∈ V of all outgoing representative edges (v, v ) ∈ E .Each element of the list stores v and a pointer to the link of (v, v ).A findUp on an ordered triangle (v 1 , v 2 , v 3 ) works as follows.It first finds a representative ordered edge (v a , v b ) from the triangle.Let's call the third vertex on the triangle v c .It looks up the link of (v a , v b ) in the adjacency list for v a , and searches for v c in the link.If (v 1 , v 2 , v 3 ) and (v a , v b , v c ) have the same orientation (are an even permutation of each other) findUp returns the next vertex in the link, otherwise it returns the previous vertex in the link.A findUp on a vertex can be implemented by selecting any of its outgoing edges, and selecting the first two vertices of the edge's link.A vertex, however, might have no outgoing edges in E .For such a vertex v the data structure can store (v 1 , v 2 ) for any triangle (v, v 1 , v 2 ).The triangle can be used to find the tetrahedron.To support findUp on edges requires storing all edges (in one direction), but not necessarily their links.For edges not in E (i.e., odd-even edges), the data structure needs only store a single vertex in their link.
To analyze the space for this data structure we assume that the links of representative edges are stored as lists of vertices.Each list element has two pointers: one to the vertex and one to the next element in the list.For an edge e ∈ E there is a one-to-one correspondence between the triangles for which e is a face and list elements in the link of e.Since a triangle has 3 edges, and on average half the edges will appear in E2 , every triangle will contribute and average of 3 * .5 = 1.5 list elements to the overall data structure.Since there are twice as many triangles as tetrahedrons, each tetrahedron will contribute an average of 3 list elements, which corresponds to 6 pointers.We also need to store the vertex adjacency lists for out-edges in E .Each edge (v 1 , v 2 ) ∈ E will appear as an element in one list (v 1 ), and will require three pointers: one to v 2 , one to the link of (v 1 , v 2 ), and one to the next element in the list.Additionally a pointer from each vertex to its list is required.The total space to support findUp on triangles is therefore 6|T | + 3/2|E| + |V |.For a typical reasonably shaped mesh |E| ≈ 7/6|T | and |V | ≈ 1/6|T |, giving approximately 8 pointers/tetrahedron.This is the same as the data structure based on tetrahedrons.
The additional space to support findUp on vertices is trivial since most vertices already have an outgoing edge.To support findUp on edges, we need to separately store the excluded edges (v 1 , v 2 ) that are not in E in either direction.These can be stored off of v 1 using a linked list with 3 pointers per edge-one for v 2 , one for some v in the link of (v 1 , v 2 ), and one for the next pointer.This comes to about 3 * 1/2 * |E| = 7/6 * 3/2|T | = 7/4|T |.Many applications will not need findUp on edges, so in these cases this extra data need not be stored.
For manifolds with boundaries, the link might be a path of vertices instead of a cycle.We can simply keep the path starting at the first element.For pseudomanifolds the link of singular vertices (2D) or edges (3D) can consist of a set of cycles and/or paths.We call this set the link set and it can be represented as multiple lists.
A d-simplex s can be deleted by finding the representative (d − 2)-simplices that faces of s, and splitting a cycle or path of each of their links.For example, in Figure 1 when the triangle is deleted from (a) going to (b), the path for the link of vertex v is split into two paths.Similarly the cycles for the other two vertices on the triangle are each slit into a path.If splitting a link leaves the link set empty, then the (d − 2)-simplex is deleted.A d-simplex s can be added by finding the representative (d − 2)-simplices of s, and extending each of their link sets.This extension might add a new path to the set (if neither of the two new vertices are in the set), it might extend an existing path (if one vertex is in the set), it might join two existing paths (if the two vertices are in separate paths), or it might joint a path into a cycle (if the two vertices are the ends of the same path).If the data structure is restricted to manifolds with boundary, then the single path must either be extended by one (on either side), or jointed into a cycle.
Data can be added to the d simplices (or d − 1 simplices) by adding a data field to each element of the link.Since a d simplex will appear in multiple links, the data only needs to be stored on one of them (chosen in a fixed manner to make lookup easy).We make use of this in the compressed data structure.Compressed Data Structure: We first discuss how to compress the data structure in 2D.Compression in 3D is similar.We make use of difference coding, in which each element in a vertex's link is represented by its difference from the original vertex.If these differences are small, then a variable-length prefix code (such as the Gamma code of Elias [16]) can represent them efficiently.An additional sign bit can be added to allow for negative differences.To ensure that the differences are small, our algorithm relabels the vertices in a preprocessing phase which we will discuss later.
Once the vertices are relabeled, the link of a vertex can be represented by concatenating the code for its degree to the codes for the differences of its neighbors.(See Figure 3 for an example.)If a vertex has a link consisting of multiple cycles or paths (as can occur in a pseudomanifold), this link set can be represented by putting the cycles/paths one after the other with a count before each.If data is associated with some of the simplices, this can be interleaved with the codes for the neighbors.The resulting vertex encodings are stored in fixed-length blocks; if an encoding is larger than will fit in one block, multiple blocks may be formed into a linked list to hold the encoding.Our data structure makes use of a hashing technique to minimize the size of the pointers used in these linked lists.
When the data structure is queried, the code for the corresponding vertex is decompressed.When an update is made, the code for the corresponding vertices is decompressed, modified, and then compressed again.
Compression of a 3D data structure is similar except that the data structure stores the link around representative edges rather than around vertices.For each vertex the data structure stores a list of that vertex's representative out-edges, with pointers to the links of those out-edges.These pointers are compressed using the same hashing technique as above.
Generating Labels.If all the vertices are known before the algorithm begins, our algorithm can relabel them using a technique based on x-y cuts.Given a set of points, the technique first finds which of the x and y axes has the greatest diameter.It finds the approximate median in that coordinate and partitions the points on either side of that median.The points on one side are labeled first, then the points on the other side.This is done recursively to produce a labeling in which points that are near each other have similar labels.This is similar to a separator-based technique for graph compression through relabeling [4] except that it occurs before any edges have been added to the mesh.
If not all vertices are known before the algorithm begins, our algorithm can assign a sparse labeling to the initial vertices.When a new vertex is added, it is assigned a label that is close to the labels of its neighbors.It would be inefficient to allocate storage for every possible label; instead, our algorithm uses an extra level of indirection to map vertex labels to memory blocks.
6 Implementation 2D Triangulation.Our 2-dimensional compressed data structure is implemented as follows.
For difference encoding our structure uses the nibble code, a code of our own devising that stores integers using 4-bit "nibbles".Each nibble contains three bits of data and one "continue" bit.The continue bit is set to 0 if the nibble is the last one in the representation of a integer, and 1 otherwise.We find that this code is much faster than the gamma code while being almost as space-efficient.
It is sometimes necessary to store an extra bit b with a value v.This is accomplished with a shift operation: v ← 2v + b.In particular, if any value might be negative, our difference coder stores its absolute value plus a sign bit: v ← 2|v| + sign(v).
A vertex is represented with a nibble code for the degree of the vertex, followed by nibble codes for the differences to each of the vertex's neighbors.Our implementation stores two additional "special-case" bits with each neighbor to provide information about the triangle that precedes it in the link.One bit is set to indicate a gap in the link set: it indicates that there is no triangle preceding that neighbor in the mesh.The other bit is set when data is associated with the triangle preceding that neighbor.In this case, the code for that neighbor is followed with a nibble code representation of the data.
As an optimization, note that for many vertices none of the special-case bits will be set.Our implementation stores a bit with the degree of each vertex to indicate if none of its special-case bits are set; if this is so, those bits are omitted in the encoding of that vertex.
Our implementation stores the nibble codes for each vertex in an array containing one seven-byte block per vertex.If a block overflows (that is, if the storage needed is greater than seven bytes), additional space is allocated from a separate pool of seven-byte blocks.The last byte of the block stores a pointer to the next block in the sequence.Our implementation uses a hashing technique to ensure that the pointer never needs to be larger than one byte.This requires a hash function that maps (address, i) pairs to addresses in the spare memory pool.Our implementation tests values of i in the range 0 to 127 until the result of the hash is an unused block.It then uses that value of i as the pointer to the block.Under certain assumptions about the hash function, if the memory pool is at most 75% full, the probability that this technique will fail to find an i ≤ 127 is at most .75 12810 −16 .
If the vertices are labeled sparsely (so that new labels can be generated dynamically), our implementation also makes use of a hash mapping between labels and vertex data blocks.One byte of memory is allocated per label; if the label is in use, this byte contains a hash pointer to the first data block for that vertex.
One bit is stored with each block to indicate whether the current block is the last in the sequence.For the first block this bit is stored with the degree of the vertex; for subsequent blocks it is stored as the eighth bit of the one-byte pointer to that block.
There is a tradeoff in the sizes of the blocks used.Large blocks are inefficient since they contain unused space; small blocks are inefficient since they require space for pointers to other blocks.In addition, there is a cost associated with computing hash pointers by searching for unused blocks in the memory pool.Figure 4 shows the tradeoff between these factors for our Delaunay triangulation algorithm run on 2 20 uniformly distributed points in the unit square.We chose a block size of 7 since it gives the most efficient use of space.
To improve the efficiency of lookups our implementation use a caching system.When a query or update is made, the blocks associated with the appropriate vertex are decoded.The information is represented in uncompressed form as a list with one vertex in the link per element of the list.The lists are kept in a FIFO cache with a maximum capacity of 2000 nodes.Update operations may affect the lists while they are in the cache.Figure 4: The number of extra blocks needed for 2 20  vertices on a uniform distribution in 2D, and the total space required if we allocate 30% more blocks than are needed.
flushed from the cache.3D Triangulation.Our 3-dimensional structure is implemented as a slight generalization of our 2dimensional structure.Recall that our 3D data structure keeps a map from each vertex v to all of its representative out-edges.This is stored as a difference coded list of the corresponding neighbors.The code for each neighbor v is followed by a code for the number of nibbles in the encoding of the representative edge (v, v ), and a pointer to the first block containing the data for that edge.(The pointer is stored using the same hash trick as above to keep pointer sizes small.)Every representative edge has its own block allocated from the memory pool, with the capability to allocate additional blocks if needed.
When an edge is queried, our implementation loads only the list for one vertex and for the edge itself into the cache.It does not need to decompress the other edges adjoining that vertex.
Since the number of nibles needed per representative edge is quite variable, our data structure allocates from pools of 2, 4, 6, 8, or 10-byte blocks to reduce wasted space.The number of blocks in each pool was determined experimentally and is shown in Figure 5.The data structure ensures that each pool always has at least 10% free space; if a block cannot be allocated from a given pool, the data structure looks for a larger one.The initial block for each vertex comes from a separate array containing blocks of size 7.
Dynamic point generation.To support dynamic point generation we use an expanded label space.If a total of n vertices are to be generated, we allow for 2n possible labels.Each label receives a one-byte hash pointer which, if the label is in use, points to the initial data block for the corresponding vertex.The initial vertices are spread evenly across the label space.Figure 5: The number of blocks of each size that are allocated for an n-vertex 3D mesh, and the percentage of blocks that were used for n = 2 10 , 2 15 , and 2 20 .
Incremental Delaunay: We implemented a Delaunay triangulation algorithm in two and three dimensions using our compressed data structure.We employ the well-known Bowyer-Watson kernel [9,43] to incrementally generate the mesh.During the course of the algorithm a Delaunay triangulation of the current pointset is maintained.An incremental step inserts a new vertex into the mesh by determining the elements that violate the Delaunay condition.Those elements form the Delaunay cavity.The faces that bound the cavity are called the horizon.The mesh is modified by removing the elements in the cavity and connecting the new vertex to the horizon.
The cavity is connected, thus can be found by a local search on the current mesh.When a point p is inserted, the cavity is determined starting from any element that will get removed by the insertion of p.To archieve optimal runtime bounds we use the idea of Guibas et al. [20] and maintain an association of every point p not yet inserted into the mesh with the element t p that contains p.The search for the cavity of p will start at t p .In their algorithm the history of the mesh is kept and at the time p is inserted that history is used to locate the element t p in the current mesh.In contrast we do not keep the mesh history but maintain the association of noninserted points p to containing elements t p on the current mesh.
At each incremental step all points on cavity elements have to be reassociated with new elements using lineside tests in 2D and planeside tests in 3D, which accounts for the dominant cost of the algorithm.We have carefully implemented the bulldozing idea described in [6] and extended it to three dimensions.
Our implementation does not require extra memory for the lists of points since at any time a point is either a vertex in the mesh or in one such list.The memory that will be used to store the vertex in the mesh can first be used as a list node.
The algorithm maintains a work queue of elements whose interiors contain points.When no elements contain points (i.e., all have been added to the mesh), the algorithm terminates.
In this scenario all points are known at the beginning.We relabel the input points using cuts along coordinate directions as described above.The runtimes reported in the next section include this preprocessing step.
Delaunay Refinement: To test our implementation's performance for the case when new points are dynamically generated at runtime, we implemented a 2D Delaunay refinement code in the style of Ruppert [35].We augment a Delaunay triangulation by adding circumcenters of badly shaped triangles while maintaining the Delaunay property.When the initial triangulation is built we walk through the mesh once and check the quality of each element, queuing the ones not satisfying a preset minimum angle bound.The same work queue used in the triangulation phase of the algorithm is used to store the list of triangles to be split.
Whenever a new point p is generated the algorithm assigns as new label by considering the horizon vertices H of the cavity created by p and calculating the value v that minimizes the sum of the log norms to H.It then finds the closest label to v that is not yet used.
In the pure triangulation code, all vertices are known at the beginning, so we can store the point coordinates and the first level vertex arrays densely.In the refinement code we can only fill these arrays up to about 85% before the open address hashing takes prohibitively long.We also require extra memory for the additional map form the label space to the vertices.
2D Delaunay: We tested our 2D implementation on data drawn from several distributions to assess its memory needs for non-uniform data sets.We ran tests on the following distributions: Uniformly random, normal, kuzmin, and a line singularity.Details on these distributions can be found in [7].In Figure 6 we report the number of extra (overflow) 7-byte blocks used to store Delaunay meshes of various point distributions and the runtime of our implementation.It  Next, we compare runtime and memory usage of our implementation to Shewchuk's Triangle [37] code which is the most efficient code reported by Boissonnat et.al. [8].In Figure 7 we report the runtime of our (incremental) code vs. Triangle's divide-and-conquer and its incremental implementation.We report the total memory use of both codes in Figure 8 and break down our memory use for the simplicial mesh, point coordinates and the work queue in Figure 9.While using just about a third of the memory our code runs about 10% slower than Triangle's divide-and-conquer implementation and is about an order of magnitude faster than Triangle's incremental implementation.In our code 50% of the memory is used to represent the mesh, 40% to store the coordinates, and 10% for the work queue.
3D Delaunay: As in 2D we tested our 3D implementation on the same four point distributions.In our 3D structure we allocate memory blocks of different size.To compare the memory needs for various point distribution, we report the number of bytes used to store occupied blocks in Figure 10.As in 2D the runtimes differ, but the memory needed is nearly independent of the distribution.
We compare our 3D implementation with uniform random data to Shewchuk's Pyramid code [36]  Figure 13 breaks down the memory usage of our code. 3We note that the version of Pyramid we are using is a Beta release.In comparison our implementation runs slightly faster and uses only about one third of the memory.In 3D the representation of the mesh uses about 75% of the total memory; point coordinates and work queue account for 18% and 7%, respectively.
2D Delaunay refinement: We tested our 2D Delaunay refinement code and compare runtime and memory use to our pure 2D Delaunay code, see Figures 14  and 15.The Figures show problem size in terms of the final number of elements in the mesh.In the pure Delaunay code, all n points are known initially, whereas in the refinement code only n/2 points are known initially, the other n/2 are generated and labeled on the fly as described in Section 6.We refine the mesh up to a minimum angle of 26.85 • .
The runtimes for the two versions are almost identical.We need about 30% more memory in the refinement code.Additional memory is needed for the map from labels to vertices and for slack in the point coordinate array and the first level vertex array needed for our hashing technique.

Discussion
The representation we described can be used as an alternative to external memory (out-of-core) representations, when the mesh is within a factor of five or so of fitting in memory relative to a standard representation.Our representation has the advantage that it allows random access to the mesh without significant penalty, and can therefore be used as part of standard in-memory algorithms (or even code) by just exchanging the mesh interface.
In conjunction with external memory.For very large problems our representation can be used in conjunction with external-memory techniques.Since in our representation the ordering of the vertices is designed to be local (it is based on the quad/oct tree    decomposition), and the blocks of memory for vertices are laid out in this ordering, nearby vertices in the mesh will most likely appear on the same page.One problem is that if the data for a vertex overflows we now assign the overflow data to the extra blocks using a hash, which has no locality.To make sure that the overflow data has some spatial locality one could be more careful about assigning the extra blocks (e.g.preferentially within the same page as the original block).Based on this representation, algorithms that have a strong bias to accessing the mesh locally (e.g., see the recent work of Amenta, Choi and Rote [1]) will tend to have good spatial locality and work well with virtual memory when it does not fit into physical memory.
Generalizations to d-dimensions.The idea of storing the link of every d − 2 dimensional simplex generalizes to arbitrary dimension.The compression technique also generalizes to arbitrary dimension, but is likely to be ineffective for large dimensions.This is because the size of the difference codes depends on the separator sizes [4], which in turn depends on the dimension.Choosing an effective way to select the representative subset of the d − 2 dimensional simplices will depend on the dimension and would need to be considered to use our representation on dimensions greater than three.We have not done any experimentation to analyze the effectiveness of our techniques on dimensions greater than three, or to compare our representations to other representations.

Figure 1 :
Figure 1: Example of a 2D manifold complex with boundary (a) and a pseudomanifold complex (b) along with the link of a vertex v.The link is a single path in (a) and two paths in (b).

Figure 2 :
Figure 2: Pseudocode for computing the d−1 boundary M b of a d simplicial mesh M .

Figure 3 :
Figure 3: The neighborhood and corresponding difference code data for vertex 314.The first entry, 6, is the degree of the vertex.Other entries are the offsets of the neighbors.
3 .Figures 11 and 12 show the runtime and the memory usage.

Figure 15 :
Figure 15: Memory use in 2D, pure Delaunay vs. Delaunay refinement The lists are encoded back into blocks when they are can be seen that the runtime varies by about 40% while the number of extra blocks varies by about 10%.Furthermore the number of extra blocks used comes to only about 28% of the number of default blocks needed, which is one