# «POINTERLESS IMPLEMENTATION OF HIERARCHICAL SIMPLICIAL MESHES A N D EFFICIENT N E I G H B O R F I N D I N G IN A R B I T R A R Y DIMENSIONS* F. BETUL ...»

International Journal of Computational Geometry

\\J

^ World Scientific

&c Applications

www.worldscientific.com

Vol. 17, No. 6 (2007) 595-631

© World Scientific Publishing Company

## POINTERLESS IMPLEMENTATION OF

## HIERARCHICAL SIMPLICIAL MESHES A N D

## EFFICIENT N E I G H B O R F I N D I N G IN A R B I T R A R Y DIMENSIONS*

## F. BETUL ATALAY

Mathematics and Computer Science Department Saint Joseph's University, Philadelphia, PA 19131, USA fatalay@sju. edu## DAVID M. MOUNT

Department of Computer Science and Institute for Advanced Computer Studies University of Maryland, College Park, MD 20742, USA mountQcs. umd. edu Received 10 February 2005 Revised 24 May 2006 Communicated by JSB Mitchell## ABSTRACT

Keywords: Pointerless data structures; neighbor finding; hierarchical simplicial meshes.

1. I n t r o d u c t i o n Hierarchical simplicial meshes have been widely used in various application areas such as finite element computations, scientific visualization and geometric modeling. There has been considerable amount of work in simplicial mesh refinement, particularly in 2- and 3-dimensions, and a number of different refinement techniques have been proposed. 1 ' 2 ' 3 ' 4 ' 5 ' 6 ' 7 ' 8 Because of the need to handle data sets with temporal components, there is a growing interest in higher dimensional meshes. In this paper, we build on a bisection refinement method proposed by Maubach. 7 *This material is based upon work supported by the National Science Foundation under Grant No. 0098151. This work was done while the first author was a graduate student at the University of Maryland, College Park.

596 F. B. Atalay & D. M. Mount A hierarchical mesh issaid to be regular if the vertices of the mesh are regularly distributed and the process by which a cell is subdivided is identical for all cells.

Maubach developed a simple bisection algorithm based on a particular ordering of vertices and presented a mathematically rigorous analysis of the geometric structure of the hierarchical regular simplicial meshes in any dimension d? Each element of such a mesh is a d-simplex, that is, the convex hull of d + 1 affinely independent points. 9 The mesh is generated by a process of repeated bisection applied to a hypercube that has been initially subdivided into d\ congruent simplices. The subdivision pattern repeats itself on a smaller scale at every d levels. Whenever a simplex is bisected, some of its neighboring simplices may need to be bisected as well, in order to guarantee that the entire subdivision is compatible. Intuitively, a compatible subdivision is a subdivision in which pairs of neighboring cells meet along a single common face. A compatible simplicial subdivision is also referred to as a simplicial complex.10 (See Fig. 1 for a 2-dimensional example.) Compatibility is important, since otherwise, cracks occur along faces of the subdivision, which in turn present problems when using the mesh for interpolation.

** Fig. 1. Compatible simplicial mesh in the plane.**

In computer graphics, adaptively refined regular meshes in 2- or 3-dimensions have been of interest for their use in realistic surface and volume rendering. 11 ' 12 ' 13 In many such applications, efficiency of various operations such as traversal and neighbor finding on the mesh is most desired. Based on the 3-dimensional version of Maubach's method, Hebert 14 presented a more efficient symbolic implementation of regular tetrahedral meshes by introducing an addressing scheme that allows unique labeling of the tetrahedra in the mesh, and he showed how to compute face neighbors of a tetrahedron based on its label. Hebert's addressing scheme could be generalized to higher dimensions, however the neighbor finding algorithms are quite specific to 3-dimensions, and a generalization to higher dimensions is a definite challenge. In this paper, we present such an algorithm that works in arbitrary dimensions.

Our interest in higher dimensions is motivated by another computer graphics application that accelerates ray-tracing 15 through multi-dimensional interpolation.

In this application, rays in 3-space are modeled as points in a 4-dimensional parameter space, and each sample ray is traced through a scene to gather various geometric attributes that are required to compute an intensity value. (A ray is actually a 5-dimensional entity, but when rays are shot from a distant origin, it is more efficient to ignore the ray's origin and model it by its supporting line in space.) Pointerless Implementation of Hierarchical Simplicial Meshes 597 Since tracing a ray through a complex scene can be computationally intensive, our approach is to instead collect and store a relatively sparse set of sample rays in a fast data structure and associate a number of continuous geometric attributes with each sample. We then interpolate among these samples to reconstruct the value at intermediate rays. 16 ' 17 Because of variations in the field values, it is necessary to sample adaptively, with denser sampling in regions of high variation and sparser sampling in regions of low variation. An adaptively refined simplicial mesh is constructed over the 4-dimensional domain of interest, and the field values are sampled at the vertices of this subdivision. Given a query point, we determine which cell of the subdivision contains this point, and the interpolated value is an appropriate linear or multi-linear combination of the field values at the vertices of this cell.

For interpolation purposes, compatibly refined simplicial meshes are preferable over quadtree-like subdivisions, since they guarantee C° continuous interpolants.

The problem with decompositions based on hyp errect angles (such as quadtrees and kd-trees) is the problem of "cracks." These arise when neighboring leaf cells are refined to different refinement levels, and new sample vertices are inserted along the boundary of one cell but not the other. Although it is possible to eliminate cracks by further subdividing a quadtree subdivision to produce a simplicial complex, 18 ' 19 ' 20 these approaches do not scale well to higher dimensions due to the exponential increase in the number of vertices in each cell and combinatorial explosion of different cases that need to be considered. Another advantage with simplicial decompositions is that interpolation is simpler and more efficient because linear interpolations can be used, which are based on a minimal number of samples (d + 1 samples for ddimensional simplices) rather than multilinear interpolations (involving 2d samples) for hyperrectangles. To illustrate these advantages more concretely, consider the images generated from our ray-tracing application in Fig. 2. Images (a) and (c) show the result of an interpolation based on kd-trees, 17 and images (b) and (d) show the results of using the hierarchical simplicial decomposition described in this paper.

In addition to our own motivation, higher dimensional meshes are of interest for visualization of time-varying fields, and efficient algorithms for performing travers a l and neighbor finding is required.

Thus, our main objective in this paper is to present an efficient implementation of hierarchical regular simplicial meshes in any dimension d. Rather than representing the hierarchy explicitly as a tree using parent and child pointers, we use a pointerless representation in which nodes are accessed through an index called a location code. Location codes 21 ' 22 have arisen as a popular alternative to standard pointer-based representations, because they separate the hierarchy from its representation, and so allow the application of very efficient access methods, such as hashing. The space savings realized by not having to store pointers can be quite significant for large multidimensional meshes. Each node would nominally be associated with d + 4 pointers. These are its pointers to its parent and each of its two children in the hierarchy, and pointers to each of its d + 1 face-neighboring simplices.

598 F. B. Atalay & D. M. Mount

(c) (d) Fig. 2. Results of a ray-tracing application to produce an 800 x 800 image based on 4-dimensional interpolations using (a) a kd-tree based on 14,492 samples (96 CPU seconds) and (b) a simplex decomposition tree based on 6,072 samples (97 CPU seconds). Details of these images are shown in (c) and (d), respectively. Note the blocky artifacts in the kd-tree approach (c).

We store the mesh in a data structure called a simplex decomposition tree. We present a location code, called the LPT code, which can be used to access nodes of this tree. Our hierarchical decomposition is based on the same bisection method given by Maubach. 7 (Note that Maubach's representation is not pointerless.) In addition to efficient computation of neighbors, we show how to perform tree travers a l, point locations, and answer interpolation queries efficiently through the use of these codes.

The remainder of the paper is organized as follows. In the next section we present prior work and discuss pointerless representations for hierarchical structures. In Section 3 we present basic notation and definitions. In Section 4 we introduce the LPT code and in Sections 5 and 6 we explain how to use the LPT code to perform tree traversals and compute neighbors.

2. Pointerless Representations and Prior Work Regular subdivisions have the disadvantage of limiting the mesh's ability to adapt to the variational structure of the scalar field, but they provide a number of significant advantages from the perspectives of efficiency, practicality, and ease of use. The number of distinct element shapes is bounded (in our case by the dimension d), Pointerless Implementation of Hierarchical Simplicial Meshes 599 and hence it is easy to derive bounds on the geometric properties of the cells, such as aspect ratios and angle bounds. The regular structure relieves us from having to store topological information explicitly, since this information is encoded implicitly in the tree structure. Additionally, the hierarchical structure provides a straightforward method for performing point location, which is important for answering interpolation queries.

One very practical advantage of regularity involves performance issues arising from modern memory hierarchies. It is well know that modern memory systems are based on multiple levels, ranging from registers and caches to main memory and disk (including virtual memory). The storage capacity at each level increases, and so too does access latency. There are often many orders of magnitude of difference between the time needed to access local data (which may be stored in registers or cache) versus global data (which may reside on disk). 23 Large dynamic pointer-based data structures are particularly problematic from this perspective, because node storage is typically allocated and deallocated dynamically and, unless special care is taken, simple pointer-based traversals suffer from a nonlocal pattern of memory references.

This is one of the principal motivating factors behind I/O efficient algorithms 24 ' 25 and cache conscious and cache oblivious data structures and algorithms. 23,26 In contrast with pointer-based implementations, regular spatial subdivisions support pointerless implementations. Pointerless versions of quadtree and its variants have been known for many years. 27 ' 21 The idea is to associate each node of the tree with a unique index, called a location code. Because of the regularity of the subdivision, given any point in space, it is possible to compute the location code of the node of a particular depth in the tree that contains this point. This can be done entirely in local memory, without accessing the data structure in global memory. Once the location code is known, the actual node containing the point can be accessed through a small number of accesses to global memory (e.g., by hashing).

Prior work on pointerless regular simplicial meshes has principally been in 2and 3-dimensions. Lee and Samet presented a pointerless hierarchical triangulation based on a four-way decomposition of equilateral triangles. 22 Evans, Kirkpatrick and Townsend28 considered triangulations based on bisection of right triangles.

(This corresponds to the 2-dimensional case of the regular simplicial meshes we consider.) They developed a location code for this triangulation and provided an efficient neighbor finding method based on bit manipulation. Hebert presented a location code for bisection-based hierarchical tetrahedral meshes and a set of rules to compute neighbors efficiently in 3-space.14 Lee, De Floriani and Samet developed an alternative location code for this same tetrahedral mesh, and presented algorithms for efficient neighbor computation. 29 Both approaches are based on an analysis of specific cases that arise in the 3-dimensional setting, and so do not readily generalize to higher dimensions.