SHAPES AND IMPLEMENTATIONS IN THREE-DIMENSIONAL GEOMETRY Ernst Peter M"ucke Department of Computer Science University of Illinois at Urbana-Champaign September, 1993 Abstract Frequently, data in scientific computing is in its abstract form a finite point set in space, and it is often useful or required to compute what one might call the "shape" of the set. For that purpose, this thesis deals with the formal notion of the family of alpha shapes of a finite point set in threedimensional space. Each shape is a well-defined polytope, derived from the Delaunay triangulation of the point set, with a real parameter controlling the desired level of detail. Algorithms and data structures are presented that construct and store the entire family of shapes, with a quadratic time and space complexity, in the worst case. Implementations of the algorithms are discussed, with an emphasis on the robust construction of three-dimensional Delaunay triangulations. A general-purpose programming technique, called Simulation of Simplicity, is used to cope with degenerate input data. This method relieves the programmer from the task of providing a consistent treatment for every single special case that can occur. The programs that use the technique tend to be considerably smaller and more robust than others. iii Acknowledgments If it were not for Herbert Edelsbrunner, my advisor, mentor, and friend, I would never have pursued a PhD at the University of Illinois. I have to thank him for all his patience and confidence over the years. Many thanks also go to Don Hearn, C L Liu, Jean Ponce, and Parvin Vaidya, for taking time out of their busy schedules to serve on my committee. Thanks to past and present members of our research group, Nataraj Akkiraju, Ulrike Axen, Cecil Delfinado, Tamal Dey, Mike Facello, Pat Moran, Edgar Ramos, Harald Rosenberger, Nimish Shah, Tiow-Seng Tan, and last but not least, Roman Waupotitsch. We shared many constructive discussions, group meetings, and lunches. Special thanks to Riccardo Bettati, Ping Fu, and Nany Hasan, for their support and encouragement. I would also like to take the opportunity to thank all the people who, for whatever reason, granted me their love and friendship, and helped me emotionally through difficult times; this includes my parents, Ernst and Marianne, and my sister Andrea. Finally, thanks a mill to Raiza Kolia for her love: : : and for cracking the whip over the last year. iv Contents Chapter 1 Introduction : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 1 1.1 Points and Triangulations : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 3 1.2 Shapes of Point Sets : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 7 1.3 Robust Implementation of Geometric Algorithms : : : : : : : : : : : : : : : : : : : : : 9 2 Three-Dimensional Alpha Shapes : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 13 2.1 Alpha Shapes, Hulls, and Diagrams : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 16 2.2 Delaunay Triangulations and Alpha Complexes : : : : : : : : : : : : : : : : : : : : : : 18 2.3 Intervals and Face Classification : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 20 2.4 Combinatorial Analysis : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 23 2.5 Extensions : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 24 3 Data Structures : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 25 3.1 Triangle-Edge Data Structure : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 26 3.2 Alpha Shape File and Data Structure : : : : : : : : : : : : : : : : : : : : : : : : : : : 36 4 Algorithms : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 40 4.1 Constructing Delaunay Triangulations : : : : : : : : : : : : : : : : : : : : : : : : : : : 41 4.1.1 A Randomized Incremental Flip Algorithm : : : : : : : : : : : : : : : : : : : : 44 4.1.2 A Simple Point Location Scheme : : : : : : : : : : : : : : : : : : : : : : : : : : 49 v 4.2 Computing Alpha Intervals : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 50 4.3 Geometric Primitives : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 52 5 Simulation of Simplicity : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 56 5.1 Symbolic Perturbation : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 58 5.2 Implementing Predicates and Other Primitive Operations : : : : : : : : : : : : : : : : 62 5.2.1 Primitives for Delaunay Triangulations : : : : : : : : : : : : : : : : : : : : : : : 63 5.2.2 Primitives for Alpha Shapes : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 66 5.2.3 Evaluating Epsilon Polynomials : : : : : : : : : : : : : : : : : : : : : : : : : : : 69 6 Signatures : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 71 7 Performance : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 82 7.1 Randomized Incremental Flip with SoS : : : : : : : : : : : : : : : : : : : : : : : : : : 83 7.2 Constructing the Alpha Shape File and Data Structure : : : : : : : : : : : : : : : : : 93 8 Conclusion : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 95 References : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 97 Vita : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 102 vi Chapter 1 Introduction The aim of this work is to further reduce the gap between theory and practice in computational geometry. With its first appearance in the early seventies, the field is relatively young, but it has since blossomed enormously within the framework of algorithm design, data structures, and complexity theory. In the mid eighties, it achieved a somewhat mature status, as indicated by the publication of several books on this subject (for example, Mehlhorn [58], Preparata and Shamos [67], Toussaint [77], or Edelsbrunner [21]). Computational geometry seemed to provide "an exceedingly fruitful arena for the more theoretical computer scientist" (Forrest [34]); and now, in the nineties, it even found its way into the introductory texts of the computer science curriculum (like Cormen, Leiserson, and Rivest [13], or, more recently, Nievergelt and Hinrichs [63]), and it is time for the field to develop a stronger applied component. However, there is still a huge gap between the world of geometric "paper algorithms" and the "real world" of applications and implementations, and this in spite of wide areas of possible applications, such as, pattern recognition, image processing, geometric modeling, computer-aided design, robotics, computational graphics, and scientific visualization. This thesis studies shapes and triangulations. These are geometric concepts with a topological flavor which have connections to many application areas (which is not always true for the more "theoretical" problems and structures studied in computational geometry). Their construction and actual implementation is therefore important; this is the emphasis of this thesis. 1 Shapes. Scientific computing often deals with data which is, in its abstract form, a finite point set in three-dimensional space, and it is sometimes useful or necessary to compute what one might call the "shape" of such a set. For that purpose, this thesis implements the formal notion of the "family of alpha shapes" of a finite spatial point set. Each shape is a polytope, derived from the Delaunay triangulation of the given points, with a real parameter ff controlling the desired level of detail. The shapes are mathematically well-defined and can be computed and applied. They facilitate the computation of features, such as volume, surface area, connected components, voids, tunnels, and their sizes. All this makes alpha shapes amenable to applications in science and engineering. Implementation. The major obstacle towards implementing geometric textbook algorithms is the principal difficulty of combining numerical or geometric data with topological information. This results in possible inconsistencies due to numerical errors, leading to incorrect results, or even, program crashes. Until a few years ago, computational geometers, being mostly theoreticians, rarely addressed this problem. Even worse, simplifying assumptions are commonly made for the sake of a simpler argument and proof: "Without loss of generality," the input data is assumed to be "sufficiently well-behaved" such that the algorithm works correctly. On paper, such assumption are useful, in order not to distract the reader from the main concepts of the algorithm by discussing too many of the tedious implementation details. However, these "details" are crucial for the programmer and have to be dealt with at some point, or the algorithm will never leave the text processing facility of its designer (or publishers). This thesis uses a symbolic perturbation method called "Simulation of Simplicity" to achieve a robust and stable implementation of alpha shapes. This method simulates an infinitesimal perturbation of the points on the level of geometric predicates and relieves the programmer from the otherwise necessary case analysis (see section 5.1). To summarize, the aim of this work is two-fold: 1. Convey the concept of three-dimensional alpha shapes to the solid modeling, computational science, and computational engeneering communities. If it is true that "a picture says more than a thousand words," then it should also be true that "a running program, with threedimensional graphics, is more convincing than a theoretical paper." Note that everything described in this thesis is implemented and publicly available. 2. Advocate symbolic perturbation as a paradigm for implementing geometric algorithms. Threedimensional alpha shapes and Delaunay triangulations are both challenging programming problems. A successful implementation should be a convincing argument that the method is indeed practical. 2 1.1 Points and Triangulations This and the remaining subsections of the introduction list some of the basic terminology, definitions, and assumptions occurring throughout this thesis. They also give some background to the topics discussed in this thesis, and provide references to the related literature. Point sets. The discussion is about problems on a given finite point set S in d-dimensional Euclidean space IRd. The notation S 2 \Gamma IR d n \Delta is used for such a set S = fp1; p2; : : : ; png, where n = jSj is the size of the set. Usually, the points are specified by vectors, pi = (ssi;1; ssi;2; : : : ; ssi;d), of d coordinates in a (right-handed) Cartesian coordinate system. A point and its corresponding vector are often used interchangeably. When convenient, IRd is seen as a vector space in which vector addition and scaling have the usual meaning. As a matter of fact, this thesis deals mainly with point sets in IR3; however, this will be extended by attaching weights, or by lifting points in IR3 to IR4, whenever appropriate. A point x is an affine combination of a point set T = fpi1; pi2; : : : ; pikg if x = Pk`=1 *i`pi` for suitable *i` 2 IR with Pk`=1 *i` = 1. If none of the *i` is negative, then x is called a convex combination of T . The set of all convex combinations of T is the convex hull of T , denoted by conv(T ). Furthermore, T is said to be affinely independent if no point pi` 2 T is an affine combination of T \Gamma fpi`g. Orientation. The notion of the "orientation" of a sequence of d + 1 points in IRd will be very useful when dealing with triangulations. It is defined as follows. The orientation of (pi0; pi1; : : : ; pid) is positive if det 0BBB BB@ ssi0;1 ssi0;2 \Delta \Delta \Delta ssi0;d 1 ssi1;1 ssi1;2 \Delta \Delta \Delta ssi1;d 1 ... ... . .. ... ... ssid;1 ssid;2 \Delta \Delta \Delta ssid;d 1 1CCC CCA ? 0; and negative if the determinant is negative. If the points are not affinely independent, the determinant evaluates to zero, and their orientation is not defined. It is important to note that the orientation of a sequence depends only on the relative position of the points to each other and not on their absolute positions. However, recall from linear algebra that the determinant of a matrix changes its sign if two rows are exchanged. Thus, the orientation of a permutation of (pi0; pi1; : : : ; pid) is the same as the orientation of the sequence itself if the number of transpositions is even, otherwise, its orientation is the opposite. 3 As an alternative to the definition in terms of the above determinant, it is often helpful to keep the following intuition in mind, which explains the concept by examples. If the dimension d = 1, then the orientation of (pi0; pi1) is positive if pi0 ? pi1 and it is negative if pi0 ! pi1 (compare with figure 1.1, (a) and (b)). If d = 2, then (pi0; pi1; pi2) has positive orientation if the three points define a left-turn (or counterclockwise turn) in the plane, that is, pi2 lies to the left of the directed line that passes through pi0 and pi1 in this order. If (pi0; pi1; pi2) defines a right-turn, then its orientation is negative. Note that the orientation of (pi0; pi1; pi2) is the same as the orientation of (pi1; pi2) as "seen" from pi0 . Indeed, the line through pi1 and pi2 can be identified with IR1 as soon as we choose a direction of the line. This direction is provided by the location of pi0: it goes from left to right as "seen" from pi0 (see figure 1.1, (c) and (d)). If d ? 2, then the orientation of (pi0; pi1; : : : ; pid) is the same as the orientation of (pi1; : : : ; pid) as "seen" from pi0 . For example, (pi0; pi1; pi2; pi3) in IR3 has positive orientation if pi0 observes (pi1; pi2; pi3) making a left-turn. In most situations where the concept of orientation is used, the interest is in the position of one point, pi0 , relative to d other points, pi1; pi2; : : : ; pid. We thus say that pi0 lies on the positive side of (pi1; : : : ; pid) if (pi0; pi1; : : : ; pid) has positive orientation, and pi0 lies on the negative side of (pi1; : : : ; pid) if (pi0; pi1; : : : ; pid) has negative orientation. Simplices. If a subset T ` S of size jT j = k + 1, is affinely independent, its convex hull, conv(T ), is a k-simplex , denoted by oeT . The dimension of oeT is k = jT j \Gamma 1. For convenience we call a 0-simplex a vertex , a 1-simplex an edge, a 2-simplex a triangle, and a 3-simplex a tetrahedron. -qpi1 qpi0 (a) -qpi0 qpi1 (b) -qpi2 qpi1 qp i0\Theta \Theta \Theta \Theta (c) -qpi1 qpi2 qp i0B BB B (d) Figure 1.1: The orientation of d + 1 points in d dimensions. For d = 1, positive (a) and negative (b), and for d = 2, positive (c) and negative (d). 4 Simplicial complexes. A simplicial complex is a collection C of k-simplices, for 0 ^ k, that satisfies the following two properties. (i) If oeT 2 C then oeT 0 2 C for every T 0 ` T . In other words, with every simplex oeT , C contains all faces of oeT as well. (ii) If oeT ; oeT 0 2 C, then either oeT " oeT 0 = ;, or oeT " oeT 0 = oeT "T 0 = conv(T " T 0). Note that (i) implies that oeT "T 0 is also in C. In other words, the intersection of any two simplices in C is either empty or a face of both. The highest dimension of any simplex oeT 2 C is called the dimension of C. A subset C0 ` C is a subcomplex of C if it is also a simplicial complex. The underlying space of C, denoted by jjjCjjj, is the union of all simplices of C, or, in other words, the part of IR3 covered by C. The underlying space of a simplicial complex C, is also called the polytope of C. Delaunay Triangulations. A simplicial complex C is called a (geometric) triangulation of S 2 \Gamma IR 3 n \Delta , all vertices of C are points of S, and jjjCjjj = conv(S). It is often useful to restrict the simplices of a triangulation T by certain criteria. For example, a tetrahedron oeT is part of T if and only if its circumsphere does not enclose nor contain any points of S \Gamma T . It turns out that this criterion specifies a unique triangulation of S, provided the points are in general position, which we assume. It is called the (closest-point) Delaunay triangulation of S. The construction of Delaunay triangulations is a popular topic in the area of computational geometry (see for example [21, 67]), mainly because Delaunay triangulations, and their dual Voronoi diagrams (also known as Dirichlet tessellations), have numerous applications in several different disciplines of science. For relevant surveys and bibliographies consult [1, 16, 54]. Various approaches for their construction are described in the literature. Naturally, it was the planar case, which was solved first, but extensions to IR3 and higher dimensions followed. Incremental methods (for example [4, 39]) compete with divide-and-conquer algorithms (for example [8, 42]). Newer research studies randomization [41]. One approach uses a lifting map (see below) to transform the triangulation problem in IRd to the problem of constructing the convex hull in IRd+1. This idea goes back to [7]; details on the constructions of convex hulls in d dimensions can be found in [21]. Another approach is based on local transformations or "flips" [28, 46, 52]. A variant of this method will be discussed in section 4.1 of this thesis. Lifting map. Identify IR3 with the x1x2x3-space in IR4, that is, the subspace x4 = 0. The lifting map is a geometric transform that projects points p = (ss1; ss2; ss3) in IR3 along the x4-axis onto the paraboloid of revolution U : x4 = P3i=1 x2i in IR4. Let pU = (ss1; ss2; ss3; P3i=1 ss2k) be the image of 5 p, and define SU = fpU j p 2 Sg. The convex hull of SU , conv(SU ), is a convex polytope with n vertices in IR4. A facet belongs to the lower boundary of this polytope if the polytope lies on the side of the positive x4-axis of the hyperplane that contains the facet. Otherwise, it belongs to the upper boundary of conv(SU ). If we project all facets of the lower boundary of conv(SU ) parallel to the x4-axis into IR3, along with their subfaces, then we obtain the Delaunay triangulation of S (see [7, 21]). The upper boundary of conv(SU ) also corresponds to a triangulation of S in IR3: the furthest-point Delaunay triangulation consists of tetrahedra whose circumspheres contain all points of S. Note that this generally does not use all points of S as vertices, rather only the points of S that lie on the boundary of conv(S). Regular triangulations. The notion of regular triangulations [53] generalizes Delaunay triangulations. In general, a point set has a large number of regular triangulations, and a particular one can be chosen by specifying weights of the points. Since it is not a major topic of this thesis, but will be mentioned occasionally, regular triangulations are only be briefly defined in the following. Assume that each point p of S has a weight wp 2 IR. The regular triangulation of the weighted point set S consists of k-simplices oeT , with 0 ^ k ^ 3, such that 9z 2 IR3; wz 2 IR with d(z; p)2 \Gamma wp \Gamma wz = 0 for p 2 T and d(z; p)2 \Gamma wp \Gamma wz ? 0 for p 2 S \Gamma T ; where d(z; p) denotes the Euclidean distance between points p and z. Assuming general position of the points and weights, this criterion indeed specifies a unique triangulation. In general, it uses only a subset S0 ` S, with conv(S0) = conv(S), as vertices; the points in S \Gamma S0 are called redundant. Note that Delaunay triangulations are a special case of regular triangulations: just set all the weights to zero. Further details are omitted. General position. Throughout this thesis the underlying assumption is that the points of S are in general position. For the time being, this means that ffl no 4 points lie on a common plane, ffl no 5 points lie on a common sphere, and ffl for any fixed ff, the smallest sphere through any 2, 3, or 4 points of S has a radius different from ff. 6 The first condition implies that all subsets, T ` S with jT j = k + 1 and 1 ^ k ^ 3, are affinely independent, and thus form proper k-simplices. The second condition assures that the Delaunay triangulation of S is unique. The third condition is not very intuitive at this point, but will make the definitions of alpha shapes (chapter 2) a lot easier. Although the general position assumption is already rather restrictive, it will be extended when convenient (for example, in section 4.3). It simplifies forthcoming definitions, discussions, and algorithms tremendously, and is justified when using a symbolic perturbation scheme (chapter 5). 1.2 Shapes of Point Sets The geometric notion of "shape" has no associated formal meaning. This is in sharp contrast to other geometric notions, such as diameter, volume, convex hull, etc. In spite of the subjectivity, people still overlap enough in what they mean by the "shape of an object" such that the term is meaningful and commonly used. The goal of this work is to offer a concrete and formal definition of "shape" that can be computed and applied. It is not an attempt to cover the entire range of meanings the term "shape" has in contemporary language, even if it is restricted to geometric contexts. Nevertheless, it is sufficiently flexible to facilitate a wide range of applications including the study of molecular structures and the distribution of galaxies in our universe. Chapter 2 will discuss alpha shapes as a formal definition of "shape" for a finite point sets S in three-dimensional Euclidean space, IR3. Intuitively, think of a set S as a cloud of points and consider the shape of this cloud. A peculiar aspect of the common usage of the word "shape" is that its meaning varies with the degree of detail intended. This aspect will be reflected by defining a oneparametric family of shapes ranging from "fine" and "local" to "crude" and "global". Indeed, much of the research on pattern classification and structural analysis of point sets talks about the difference between macroscopic properties (features) and microscopic properties (texture). The notion of alpha shapes unifies the two seemingly separate worlds, offering a meaningful framework that cuts across the "small" and the "large." In IR2.In IR2.In IR2. A fair amount of related work has been done for point sets in IR2, and some for point sets in IR3. Jarvis [44] was one of the first to consider the problem of computing the shape as a generalization of the convex hull of a planar point set. A mathematically rigorous definition of shape was later introduced by Edelsbrunner, Kirkpatrick, and Seidel [26]. They proposed a natural generalization of convex hulls, called ff-hulls, and their combinatorial variant, ff-shapes (see figure 1.2). Their notion of ff-shapes is the two-dimensional analogue of the spatial notion described in chapter 2. Two7 dimensional ff-shapes are related to the dot patterns (Fairfield [32, 33]) and the circle diagrams used in bivariate cluster analysis (see for example Moss [61]). Different graph structures that serve similar purposes are the Gabriel graph [56], the relative neighborhood graph [76], or their parameterized version, the fi-skeleton [48]. In IR3.In IR3.In IR3. For spatial point sets, Boissonnat [3] suggested the use of Delaunay triangulations in connection with certain heuristics to "sculpture" a single connected shape of a "cloud" of points. The concept of three-dimensional alpha shapes of section 2.1 is more general and mathematically welldefined. More recently, Veltkamp [78] generalized neighborhood graphs and their variants to the so-called fl-neighborhood graph, emphasizing applications in IR2 and IR3. This two-parameter geFigure 1.2: The complete family of alpha shapes of a set of 7 points in IR2. The shape range from "crudest" (the convex hull) down through "finest" (the point set itself). Although a shape will be defined as a polytope, we display them with their natural decomposition into simplices inherited from the Delaunay triangulation. (This makes a difference in the first 5 examples.) 8 ometric graph structure forms a continuous spectrum that ranges from the empty to the complete graph, including the convex hull, the Delaunay triangulation, the Gabriel graph, and the fi-skeleton. Observe also the superficial similarity between alpha shapes and isosurfaces in IR3. The latter is a popular concept in volume visualization (see for example [18, 55]). In contrast to alpha shapes, isosurfaces are approximate by definition, with their exact form depending on the algorithm used for their construction. Alpha shapes, fi-skeletons, fl-neighborhood graphs, : : : , there is a whole variety of theoretical geometric concepts, trying to capture vague and ill-defined notions like "proximity," "structure," and "shape," and all share a surprisingly large number of potential application areas in science and engineering. Refer to Jaromczyk and Toussaint [43] for a survey and bibliography on relative neighborhood graphs and their variants, including a list of applications in computational morphology, computer vision, geographic analysis, and pattern recognition. Three-dimensional alpha shapes are particularly close to "real-life" applications. For example, in the framework of automatic mesh generation and geometric modeling, a well-defined and concrete notion of shape, which can be computed and applied, is often all that is needed. The problem of modeling and studying molecules in computational biology is more specifically related to shapes: (weighted) ff-shapes [23] are in a strict geometric sense dual to the so-called space filling diagrams (see [24], and also section 2.1). These diagrams can be defined as a union of balls (different atoms are represented by balls of different size) and ff-shapes can be be used to compute their topological and metric properties [25]. The computation of "molecular volume" and "surface" is recognized as especially useful in the highly-funded area of protein folding [37, 70]. Computational tools for such measurements exist for several years; refer, for example, to the works of Richards [69] and Connolly [9-12]. Through the use of ff-shapes these measurements can be improved in terms of computational speed and accuracy. It is thus hoped that ff-shapes prove to be a useful tool in future protein research. All this makes the design of robust algorithms and programs for their construction an important but challenging problem. 1.3 Robust Implementation of Geometric Algorithms In the mid eighties, Forrest [34] was among the first to recognize that correct implementations of even the simplest geometric operations are by no means trivial. In fact, he found it doubtful whether such implementations can ever exist. At that time, both theoretical and practical computational geometers were surprisingly willing to implement systems on rather shaky foundations, mainly by ignoring numerical issues at all, or at best, providing ad hoc solutions, like bounding the range of the input data for "save" computations (for an example see Beretta [2]). Over time, the situation has 9 changed, and a reasonable amount of research has been published addressing the practical issues of computational geometry. A seemingly straightforward geometric algorithm for the point-in-polygon problem is often used to illustrate what the problem is (Forrest [34]). The Parity Algorithm works as follows: let r be the horizontal half-line whose left endpoint is the query point q; count the number of intersections between r and the edges of the polygon; then, q lies inside the polygon if that number is odd; otherwise, q lies outside. It is not a trivial matter to implement this algorithm, even if we assume that q does not lie on the boundary of the polygon. There are only two "nondegenerate" or "simple" cases: either the intersection between r and an edge e of the polygon is empty, or r crosses e (see Figure 1.3, (a) and (b)). There are, however, four "degenerate" or "special" cases (as illustrated in Figure 1.3, (c) through (f)) that have to be taken into account. A correct answer is obtained if cases (c) and (e) are counted as one crossing and cases (d) and (f) are not counted at all. We realize that, even for this extremely simple problem, a substantial amount of the effort is required to cope with degenerated cases. Observe also that there are several seemingly plausible ways to treat degeneracies, and that some of them lead to incorrect implementations. The difficulties with degenerate cases are worsened if coupled with the ambiguities implied by round-off errors in common non-exact arithmetic. It is not clear what constitutes a "correct implementation" of a geometric algorithm. The literature commonly uses two terms, "robustness" and "stability." The former usually means that the program handles "all" input occurring in practice, whereas the latter often refers to how the implementation copes with numerical errors. Here, both notions are used synonymously, mostly following Fortune's definition [35]: A robust or stable implementation of a geometric algorithm is a program whose output is always "correct" in the sense that, although it might not produce the exact result for the given input S, it guarantees the correct result for a set S0 very close to the original. One school of researchers proposes the use of hardware-supported floating-point arithmetic including some sort of geometrical and numerical reasoning to achieve numerical stability. The work done in this area consists mainly of problem-specific solutions, but unfortunately, not always including implementations and empirical studies. Ottman, Thiemt, and Ullrich [64] present a robust algorithm for intersecting line-segments using "best possible evaluation" of arithmetic expressions and the q q q AA AA (a) q q q JJ JJ J (b) q -q q q @@\Gamma \Gamma \Gamma \Gamma (c) q q q q BB BB B \Delta \Delta \Delta \Delta (d) q q q q q BB BB B A AAA (e) q q q q q AA AA \Theta \Theta \Theta \Theta \Theta (f) Figure 1.3: The different cases in the Parity algorithm. 10 assumption of a precise implementation for the scalar product. Milenkovic [59, 60] achieves stable results for line arrangements by perturbing the data with an a priory constant as perturbation bound. A more general approach is due to Guibas, Salesin, and Stolfi [40]. They use interval arithmetic and backward error analysis, combined with a "good deal" of geometric reasoning. Their method results in a correct solution for a perturbed version of the input, plus a computed bound on the size of the implicit perturbation. Another approach solves numerical problems once and for all by giving up on fast floating-point arithmetic, and using exact integer arithmetic instead. Karasick, Lieber, and Nackman [47] use "adaptive-precision" arithmetic for computing the sign of determinants with integer or rational elements. They use these method to implement two-dimensional Delaunay triangulations employing the popular divide-and-conquer algorithm of Guibas and Stolfi [42]. They provide empirical tests and claim a factor of 4 for the computational overhead with respect to unstable floating-point implementations. A more rigorous approach is to not only tackle the numerical problems by using exact arithmetic but also cope with degeneracies as such by using a symbolic perturbation scheme. These methods provide the user with a black box making geometric data appear to be simple (that is non-degenerate) although in reality it is not. Two methods were developed at about the same time. The one proposed by Yap [81, 82] works for geometric primitives expressed by arbitrary polynomials. Its symbolic perturbation is implicit to the algorithm, and the consistency of the method has to be proven for each application. Our method, the Simulation of Simplicity, or SoS [27], on the other hand, achieves consistency by the proven existence of a physical equivalent to the symbolic perturbation. In contrast to Yap, we provide an efficient implementation: for example, three-dimensional Delaunay triangulations can be computed with a slowdown factor of 3 to 4, when compared with error-prone floating-point implementations (see section 7.1). This computational overhead is mostly due to the necessary use of exact arithmetic, and not to the negligible computational overhead resolving degeneracies.1 Symbolic perturbation schemes relieve both the designer and the implementer of geometric algorithms of problems with ffl numerical errors during the computation, and ffl degeneracies in the input data. Note that the occurrence of degneracies or special cases has nothing to do with round-off errors in numerical computations, rather, they are inherent to the structure of the input data. However, 1Emiris and Canny [29, 30] suggest a variant of SoS to reduce the computational steps necessary for resolving degeneracies. They do not, however, address the overhead due to exact arithmetic. 11 numerical errors are usually more serious for data configurations at or "near" special cases. A method like SoS can be thought of as a black box which breaks ties consistently, thus providing the user with the illusion of working with simple data, even if it is really in arbitrary, and therefore possibly degenerate position. SoS is convenient for theoreticians, who can now freely use the common simplifying assumptions and refer to our method. However, it is less clear that the overhead cost of the method is still practical. Especially, its usage of exact arithmetic is debatable. Several implementations will be presented in the form of a case study in order to convince the reader that the computational overhead of the simulation is feasible while giving the advantage of the reduced time spent for implementation and program maintenance. 12 Chapter 2 Three-Dimensional Alpha Shapes This section formally defines three-dimensional ff-shapes and discusses some of the related concepts and extensions. All this is supported by frequent illustrations. Conceptually, ff-shapes are a generalization of the convex hull of a point set. They constitute a whole family of shapes, ranging from the original point set up to the convex hull. They are parameterized by a real ff between 0 and 1. The ff-shape of a point set S 2 \Gamma IR 3 n \Delta is a polytope that is neither necessarily convex nor necessarily connected. For ff = 1, the ff-shape is identical to the convex hull of S. However, as ff decreases, the shape shrinks by gradually developing cavities. These cavities may join to form tunnels, and even closed holes (or voids) may appear (see figure 2.1). For ff = 0, the shape consists only of the points of S. Intuitively, a piece of the polytope disappears when ff becomes small enough so that a sphere with radius ff, or several such spheres, can occupy its space without enclosing any of the points of S. Think of IR3 filled with styrofoam and the points of S made of a more solid material, such as tiny rocks of granite. Now imagine a spherical eraser with radius ff. It is omnipresent in the sense that it carves out styrofoam at all positions where it does not enclose any of the sprinkled rocks, that is, points of S. The resulting object will be called the ff-hull. To make things more feasible we straighten the surface of the object by substituting straight edges for circular ones and triangles for spherical caps. The thus obtained object is the ff-shape of S (see figure 2.1 for some examples). It is a polytope in a fairly general sense: it can be concave or even disconnected, it can contain two-dimensional patches of triangles and one-dimensional strings of edges, and its components can be as small as single points. The parameter ff controls the maximum "curvature" of any cavity of the polytope. 13 Figure 2.1: Three-dimensional ff-shapes. (See table 2.1 for the legends to the above plates.) 14 1: Six different shapes of a set with points that were randomly generated on the surface of two linked tori. The first shape is the convex hull, for ff = 1, the last shape is the point set itself, for ff = 0. The ff-value used in the forth frame neatly separates the two tori. Further decreasing ff disassembles the shape. Singular triangles are shown in a darker color. 2: Six shapes of a set based on demo data for Silicon Graphics' Solidview program; of course, without any of the original connectivity information. The erasing sphere is shown to the right of the shape. Apart from a dense conglomerate of points representing part of the person's brain, the set is basically hollow with most points representing skin. 3, 4: These point sets are based on demo data for modeling and rendering programs. All connectivity information of the modeling data (that is, edges and triangles) was removed, leaving only the vertex coordinates. As described in section 2.2, each ffshape is triangulated by the tetrahedra of the corresponding ff-complex. This might be useful in the automatic generation of meshes for objects with nonconvex surfaces. 5: This data represents a time-averaged molecular dynamics structure of gramicidin A, a peptide that forms a channel for ion and water movement across lipid membranes. The major structural motif is a right-handed beta-bonded helix. The tunnel of the macro-structure can be detected using relatively large ff values. Smaller values of ff result in shapes with larger numbers of isolated triangles and edges. These ff-shapes reveal the helix of the micro-structure. 6: The "universe" data represents a simulation of the positions of galaxies within our universe. The theory is that galaxies first clustered into sheet-like structures, then progressed to filament-shaped structures at the intersection of multiple sheets. As filaments began to intersect, global clusters appeared. It is interesting to investigate the macro- and micro-structure of the galaxies, including the detection of large voids and local or global clusters. The full spectrum of ff-shapes promises to be useful in this study. Table 2.1: Legends to ff-shape plates. (For more information on the data sets displayed in figure 2.1 refer to chapter 7.) 15 2.1 Alpha Shapes, Hulls, and Diagrams For 0 ! ff ! 1, let an ff-ball be an open ball with radius ff. For completeness, a 0-ball is a point, and an 1-ball is an open half-space. For S 2 \Gamma IR 3 n \Delta , an ff-ball b is called empty if b " S = ;. Alpha shapes. For 0 ^ k ^ 2, a k-simplex oeT , with T ` S, is said to be ff-exposed if there is an empty ff-ball b with T = @b " S, where @b is the sphere or plane bounding b. A fixed ff thus defines sets Fk;ff of ff-exposed k-simplices for 0 ^ k ^ 2. The ff-shape of S, denoted by Sff, is the polytope whose boundary consists of the triangles in F2;ff, the edges in F1;ff, and the vertices in F0;ff (see figure 2.1). The k-simplices in Fk;ff are also called the k-faces of Sff. But which connected components of IR3 \Gamma @Sff are interior and which are exterior to Sff? Fix the value of ff and notice that for each ff-exposed triangle oeT there are two (not necessarily empty) ff-balls, b1 6= b2, so that T ` @b1 and T ` @b2. If both ff-balls are empty then oeT does not belong to the boundary of the interior of Sff. Otherwise, assume that b1 is empty and that b2 is not. In this case, oeT bounds the interior of Sff. More specifically, the interior of Sff and the center of b1 lie on different sides of oeT . The definition of interior and exterior of Sff is possibly more natural in the context of Delaunay triangulations and ff-complexes as described in section 2.2. There are quite a few natural geometric concepts that are closely related to ff-shapes. Some of them are discussed in the remainder of this section. In each case, the emphasis is on how the concept is related to ff-shapes and how this relation can enrich our understanding of ff-shapes. Section 2.2 briefly reviews Delaunay triangulations and their dual incarnations known as Voronoi diagrams. The relevance of the Delaunay triangulation of a point set is that each ff-shape of the set is the underlying space of a subcomplex of the triangulation. These subcomplexes are termed ff-complexes and defined in section 2.2. Extensions of these notions are mentioned in section 2.5. Alpha hulls. For 0 ^ ff ^ 1, the ff-hull of S, denoted by Hff, is defined as the complement of the union of all empty ff-balls. This is the formal counterpart of the styrofoam object described earlier. Sample members of the continuous family of ff-hulls are the convex hull of S, for ff = 1, and S itself, for ff sufficiently small. Observe that Hff1 ` Hff2 if ff1 ^ ff2. Alpha diagrams. Another interesting concept defined by ff-balls is what we call the ff-diagram of S, denoted by Uff. For 0 ! ff ! 1, Uff is the union of all ff-balls whose centers are points in S. In chemistry and biology, ff-diagrams are known as space filling diagrams. However, there they are usually not restricted to equally large balls. This restriction can be removed when merging the concept of weighted ff-shapes [23] with ff-diagrams (see section 2.5). 16 Observe that a point x 2 IR3 belongs to Uff if the ff-ball centered at x is not empty. Denote this ff-ball by bx. This implies the following close relationship between Hff and Uff. x 2 Uff () bx " Hff 6= ;; and x 2 Hff () bx ` Uff: Consider the boundary of Uff. It consists of spherical caps, circular arcs, and vertices which we call corners. These are the 2-, 1-, and 0-faces of Uff. These caps, arcs, and corners are in close correspondence with the vertices, edges, and triangles of Sff. Some definitions are needed to describe this correspondence. Let ff be fixed, with 0 ! ff ! 1, and let T be a subset of S of size jT j = k + 1, with 0 ^ k ^ 2. Define KT = Tp2T @bp, where bp is the ff-ball centered at p, as before. Besides general position assume KT 6= ;. For jT j = 1, KT is a sphere, for jT j = 2, KT is a circle, and for jT j = 3, KT is a pair of points. It follows from the definitions that T is ff-exposed if and only if KT contains at least one face of the boundary of Uff. If jT j = 1 this face is a cap, if jT j = 2 it is an arc, and if jT j = 3 it is a corner. This fact can be expressed as follows. oeT is a k-face of Sff, with 0 ^ k ^ 2 () KT contains at least one (2 \Gamma k)-face of Uff. Moreover, the number and structure of (2 \Gamma k)-faces contained in KT are reflected by oeT and the way it is embedded in Sff. For example, if jT j = 3, then KT contains no, one, or two corners of Uff. First, it contains no corner if and only if oeT is not a triangle of Sff. This is mentioned above. Second, KT contains one corner if and only if oeT bounds the interior of Sff. Third, both points of KT are corners of Uff if and only if both sides of oeT face the outside of Sff, that is, oeT is a triangle of Sff that does not bound its interior. Such a triangle will be called singular in section 2.3. Similarly, if jT j = 2 and oeT is an edge of Sff, then every angular interval between two incident triangles of Sff that faces the outside corresponds to an arc of Uff contained in KT . If there is no incident triangle (in this case, oeT is an singular edge) then the entire circle KT belongs to the boundary of Uff. Analogous statements can be made for jT j = 1 where spherical caps on KT correspond to solid angles around the vertex oeT that face the outside of Sff. All this amounts to a one-to-one correspondence between the (2\Gamma k)-faces of Uff and the k-faces of Sff, provided the latter are interpreted with multiplicities reflecting the number of exposed sides, angular intervals, or solid angles. The one-to-one correspondence also preserves (or reverses) incidences. This suggests that Uff, or more specifically, the boundary of Uff, be represented by Sff, or more specifically, the faces of Sff and their incidences. 17 2.2 Delaunay Triangulations and Alpha Complexes A finite point set S 2 \Gamma IR 3 n \Delta defines a special triangulation known as the Delaunay triangulation of S (see for example [21, 67]). Assuming general position of the points, this triangulation is unique and decomposes the convex hull of S into tetrahedra. The triangulation is named after the Russian geometer Boris Delaunay (also Delone) who introduced it in his seminal paper [14] in 1934. As explained below, the Delaunay triangulation of S is dual to another complex defined by S known as the Voronoi diagram [79, 80]. Both complexes are related in an interesting way to the family of all ff-shapes of S. The relationship between ff-shapes and Delaunay triangulations is of particular importance here. Delaunay triangulations. For 0 ^ k ^ 3, let Fk be the set of k-simplices oeT = conv(T ), T ` S and jT j = k + 1, for which there are empty open balls b with @b " S = T . Notice that F0 = S. The Delaunay triangulation of S, denoted by D, is the simplicial complex defined by the tetrahedra in F3, the triangles in F2, the edges in F1, and the vertices in F0. By definition, for each simplex oeT 2 D, there exist values of ff * 0 so that oeT is ff-exposed. Conversely, every face of Sff is a simplex of D. This implies the following relationship between the Delaunay triangulation and the boundary of Sff. Fk = [ 0^ff^1 Fk;ff; for 0 ^ k ^ 2. We take advantage of this relationship by representing the family of ff-shapes of S implicitly by the Delaunay triangulation of S. This will be described in detail in chapter 4. Voronoi diagrams. For a point p 2 S, define V (p), the Voronoi cell of p, as the set of points x 2 IR3 so that the Euclidean distance between x and p is less than or equal to the distance between x and any other point of S. Each Voronoi cell is a convex polyhedron, and the collection of all Voronoi cells, one for each point of S, defines the Voronoi diagram of S, denoted by V. We call a Voronoi cell also a 3-cell of V. Each 2-cell of V is the intersection of two Voronoi cells, each 1-cell or edge is the intersection of three 3-cells, and each 0-cell or vertex is the intersection of four 3-cells. There is a natural one-to-one correspondence between the k-simplices of D and the (3 \Gamma k)-cells of V. Let T be a subset of S of size jT j = k + 1, with 0 ^ k ^ 3, and define VT = Tp2T V (p). oeT is a k-simplex of D () VT is a (3 \Gamma k)-cell of V; for 0 ^ k ^ 3. This correspondence preserves (or reverses) incidences which implies that D and V are indeed dual to each other. 18 Observe that VT is the set of all points x for which there exists an empty open ball bx centered at x with T ` @bx " S. Equality holds if and only if x belongs to the relative interior of VT . It follows that oeT is ff-exposed if and only if there is a point x in the relative interior of VT whose distance from the points in T is ff. Since VT is convex there is a single interval so that oeT is ff-exposed if and only if ff belongs to this interval. We will exploit this fact later when we discuss how to decide when a simplex of D is ff-exposed. Since all faces of Sff are simplices of D, it follows that the interior of Sff is naturally triangulated by the tetrahedra of D. This idea leads to the concept of ff-complexes as defined below. Alpha complexes. Each k-simplex oeT of D defines an open ball bT bounded by the smallest sphere @bT that contains all points of T . Let %T be the radius of bT . For k = 3, @bT is the circumsphere of oeT , for k = 2, the circumcircle of oeT is a great circle of @bT , and for k = 1, the two points in T are antipodal on @bT . Call @bT the smallest circumsphere and %T the radius of oeT . For 1 ^ k ^ 3 and 0 ^ ff ^ 1, define Gk;ff as the set of k-simplices oeT 2 D for which bT is empty and %T ! ff. Furthermore, define G0;ff = S, for all ff. The sets Gk;ff do not necessarily define a simplicial complex because it can happen that G3;ff contains a tetrahedron but not all triangles of this tetrahedron belong to G2;ff. Similarly for triangles and edges. With this in mind, we define the ff-complex of S, denoted by Cff, as the simplicial complex whose k-simplices are either in Gk;ff or they bound (k + 1)-simplices of Cff. By definition, Cff1 is a subcomplex of Cff2 if ff1 ^ ff2. As defined earlier, the underlying space of Cff, denoted by jjjCffjjj, is the union of all simplices of Cff. It is a polytope in the sense specified in section 1.1, and identical to Sff, a property which is rigorously proven for the more general weighted case in appendix A of [23]: For all 0 ^ ff ^ 1, Sff = jjjCffjjj: This can be considered an alternative definition of ff-shapes. It makes it easy to specify the intervals of the simplices of D alluded to in the above discussion of Voronoi diagrams. For example, let oeT be a k-simplex of D. If bT is empty then oeT belongs to Cff if and only if ff 2 (%T ; 1].1 In this case, set *T = %T . Otherwise, bT is not empty and oeT belongs to Cff if and only if ff 2 (*T ; 1], where *T is the minimum *R over all (k + 1)-simplices oeR of D with T ` R. 1Because of the general position assumption we have ff 6= %T . It is therefore irrelevant whether the interval is open or closed at its left endpoint. 19 2.3 Intervals and Face Classification For each simplex oeT 2 D there is a single interval so that oeT is a face of the ff-shape Sff if and only if ff is contained in this interval. This was mentioned in section 2.2. It will be convenient to study these intervals for the ff-complex Cff rather than the ff-shape. Also, we break each interval into three (possibly empty) parts that correspond to values of ff for which the simplex is an interior, regular, or singular simplex of Cff. A simplex oeT 2 Cff is said to be8 ??!??: interior if oeT 62 @Sff, regular if oeT 2 @Sff and it bounds some higher-dimensional simplex in Cff, and singular if oeT 2 @Sff and it does not bound any higher-dimensional simplex in Cff. Notice that there are Delaunay edges and triangles that can never be singular because their smallest circumsphere encloses other points of S.2 Therefore, we call a simplex oeT 2 D( attached if jT j = 2; 3 and bT " S 6= ;, and unattached otherwise: Recall that %T is the radius of the smallest circumsphere of oeT . In order to break up the interval for which oeT belongs to Cff, we introduce values _T and _T for which oeT changes from singular to regular and from regular to interior, respectively. Before that, however, let up(oeT ) be the set of all simplices in D that contain a simplex oeT 2 D, with jT j ^ 3, as a proper face, that is, up(oeT ) = foeT 0 2 D j T ae T 0g: If oeT is a tetrahedron, define _T = _T = %T . Otherwise, _T = minf%T 0 j oeT 0 2 up(oeT ); unattachedg and _T = maxf%T 0 j oeT 0 2 up(oeT )g: 2It is convenient to extend the general position assumption so that no smallest circumsphere of two or three points of S contains another point of S. A slightly more general assumption is the following. If a sphere of radius ff contains three points of S then no two of them are antipodal, and if it contains four points then no three lie on a great-circle of the sphere. 20 It is sufficient to consider only the set up1(oeT ) = foeT 0 2 up(oeT ) j jT 0j = jT j + 1g; that is, all faces incident to oeT whose dimension is one higher than that of oeT , in order to derive the values _T and _T : _T = min f%T 0 j oeT 0 2 up1(oeT ); unattachedg [ f_T 0 j oeT 0 2 up1(oeT ); attachedg ! ; and _T = maxf_T 0 j oeT 0 2 up1(oeT )g: Specifying intervals. The intervals of ff values in which oeT is an interior, regular, or singular simplex of Cff are shown in the table 2.2. Because of the general position assumption, ff is different from all % values and therefore also from all _ and _ values. We can thus define all intervals as open, except at endpoints 0 and 1. It is necessary to distinguish simplices that bound the convex hull of S from the others. The next paragraph briefly explains the entries of table 2.2 for the case of triangles that do not bound the convex hull of S. The arguments for tetrahedra, triangles bounding the convex hull, edges, and vertices are similar. Consider a triangle oeT 2 D, T = fpi; pj; pkg, that does not bound the convex hull of S; we denote this by oeT 62 @conv(S). Let oeT 0 and oeT 00 be the two incident tetrahedra in D, and assume T 0 = T [ fpug, and T 00 = T [ fpvg. Furthermore, be 0 ! %T 0 ! %T 00 ! 1; in other words, _T = %T 0 and _T = %T 00. Now, fix a value for ff. If %T 00 ! ff ^ 1 then the triangle oeT is not ff-exposed. It will, however, be oeT is : : : singular regular interior tetrahedron (%T ; 1] edge or triangle, 62 @conv(S), unattached (%T ; _T ) (_T ; _T ) (_T ; 1] 62 @conv(S), attached (_T ; _T ) (_T ; 1] 2 @conv(S), unattached (%T ; _T ) (_T ; 1] 2 @conv(S), attached (_T ; 1] vertex, 62 @conv(S) [0; _T ) (_T ; _T ) (_T ; 1] 2 @conv(S) [0; _T ) (_T ; 1] Table 2.2: Intervals of ff values for which oeT 2 D belongs to the ff-complex Cff. 21 part of the interior of Sff, because both incident tetrahedra are in Cff. If %T 0 ! ff ! %T 00 then the triangle is ff-exposed and oeT 0 is in Cff but oeT 00 is not. This means that oeT 0 is a regular triangle of Cff. For ff ! %T 0, neither oeT 0 nor oeT 00 are tetrahedra of Sff, but oeT can still be a singular triangle, that is, if and only if %T ! ff and neither pu nor pv are inside bT . If one of the two points is inside bT , then oeT is attached, and oeT can never be a singular triangle of Cff, no matter what ff value is selected. The ff-complex consists of all interior, regular, and singular simplices for a given ff value. The interior of the ff-shape is triangulated by the interior simplices. The boundary of the interior is formed by the set of regular triangles and their edges and vertices. We refer to the endpoints of the intervals in table 2.2 as ff-thresholds. This does not include 0 and 1. Since all _T and _T values are % values of other simplices, each ff-threshold is the radius of a simplex in D. More specifically, the set of ff-thresholds is exactly the set of radii of all unattached k-simplices for 1 ^ k ^ 3. Define the ff-spectrum as the sorted sequence of ff-thresholds. This concept will appear again in section 3.2. Computing intervals. Assume that each simplex oeT 2 D is marked as either "2 @conv(S)" or "62 @conv(S)" after the construction of D. With this, the above intervals can be computed by classifying oeT as attached or unattached, and by computing %T , _T , and _T , whenever applicable. We said that oeT is attached if and only if one of the simplices that contain oeT has a vertex in bT , the open ball bounded by the smallest circumsphere of oeT . This implies that oeT can be classified in time proportional to jup1(oeT )j. The time it takes to classify all simplices is proportional to the number of simplices in D, because each simplex has only a constant number of faces. In other words, assuming that constant time suffices to decide whether or not a point belongs to bT (see section 4.3), a simplex can be classified in constant amortized time. Furthermore, assume that, given T with oeT 2 D, %T can be computed in constant time (again, see section 4.3). By processing tetrahedra before triangles before edges before vertices, we can get _T and _T simply as the minimum and maximum of the values %T 0, _T 0, and _T 0, for oeT 0 2 up1(oeT ). This also takes only constant amortized time per simplex. 22 2.4 Combinatorial Analysis In contrast to ff-hulls and ff-diagrams, the ff-shapes of a finite point set form a discrete family, even though they are defined for all real numbers ff, with 0 ^ ff ^ 1. Indeed, Sff1 6= Sff2 if and only if S3k=0 Gk;ff1 6= S3k=0 Gk;ff2. Thus, Sff1 6= Sff2 if and only if there is an empty open ball bounded by a smallest circumsphere of an edge, triangle, or tetrahedron of D whose radius lies between ff1 and ff2. Recall that such a radius is referred to as an ff-threshold because it separates two ff-shapes. The number of ff-shapes exceeds the number of ff-thresholds by one. It follows that one plus the total number of k-simplices of D, for 1 ^ k ^ 3, is an upper bound on the number of different ff-shapes. An upper bound on the number of simplices of D can be obtained using a relationship between Delaunay triangulations in IR3 and certain convex polytopes in IR4. We briefly describe this relationship below and refer to [21, 74] for details. The lifting map of section 1.1 connects Delaunay triangulations in IR3 with convex polytopes in IR4 by projecting the points of S onto the paraboloid of revolution U : x4 = P3i=1 x2i . According to the upper bound theorem for convex polytopes, the maximum numbers of 1-, 2-, and 3-faces of a convex polytope with n * 5 vertices in IR4 are at most 12 (n2 \Gamma n), n2 \Gamma 3n, and 12(n2 \Gamma 3n), respectively (see for example [6, 57]). The lifting map implies the same upper bounds for jFkj, with 1 ^ k ^ 3. As a matter of fact, the upper bound for jF3j is one less than for the number of 3-faces of conv(SU ) because at least one 3-face belongs to the upper boundary of conv(SU ). By a result of [74], these bounds are tight even though the vertices of conv(SU ) are constrained to lie on a second degree surface in IR4. We summarize these results for n = jSj. jF0j = n; jF1j ^ 12 (n2 \Gamma n); jF2j ^ n2 \Gamma 3n; and jF3j ^ 12(n2 \Gamma 3n \Gamma 2): By adding one to the sum of the bounds for jF1j, jF2j, and jF3j, we obtain the following result. Lemma 2.4.1. S 2 \Gamma IR 3 n \Delta has at most 2n2 \Gamma 5n different ff-shapes. This bound is too pessimistic for two reasons. First, although the upper bounds on the number of simplices of D are tight, there are many fewer simplices for most point sets. For example, if the n points are uniformly distributed in the unit ball then the expected number of simplices is only O(n) (see Dwyer [19]). Second, not all edges and triangles of D have a smallest circumsphere that bounds an empty open ball. However, since the circumsphere of every tetrahedron of D bounds an empty open ball by definition, the number of different ff-shapes is always at least proportional to the number of simplices of D. 23 2.5 Extensions The definitions presented in sections 2.1 and 2.2 above can be extended in various ways. So far we refrained from mentioning these extensions in order to avoid unnecessary complications and to be faithful to the currently available implementations. For completeness, negative values of ff, weighted points, and higher dimensions are now briefly discussed. Negative alpha values. This extension has been described in [26] for the two-dimensional case.3 For negative ff, ff-complexes are most naturally defined as subcomplexes of the so-called furthestpoint Delaunay triangulation of S (see section 1.1 or [21, 67]). For T ` S and jT j = 4, the tetrahedron oeT belongs to this triangulation if and only if bT contains all points of S \Gamma T . The ff-complex is defined by all tetrahedra whose circumspheres have radius larger than \Gamma ff. The ff-shape is, again, the underlying space of the ff-complex. These shapes exhibit far less interesting geometric and topological properties than the ones for positive ff, and are thus less interesting for applications. We omit further details. Weighted points. Recall the relationship between ff-shapes and ff-diagrams described in section 2.1. It is interesting to consider diagrams for different ball sizes, and this is indeed done in chemistry and biology where space filling diagrams are usually defined as unions of balls with arbitrary (and possibly different) radii. In order to represent such diagrams by polytopes similar to ff-shapes it is necessary to introduce weighted ff-shapes. These can be defined using subcomplexes of so-called regular triangulations (see section 1.1 or [22, 53]). Given a finite set of points, each with a real weight, the regular triangulation is a unique simplicial complex whose underlying space is the convex hull of the point set. If all weights are the same then it equals the Delaunay triangulation of the points. Details can be found in [23]. Higher dimensions. It is fairly straightforward to generalize all important concepts of this section (like ff-shapes, ff-hulls, ff-diagrams, ff-complexes, Delaunay triangulations, and Voronoi diagrams) to finite point sets S in IRd, for arbitrary dimension d. This generalization, combined with an extension to weighted points is developed in [23]. Note however, that the implementation details become progressively hairier with increasing dimension, and the worst-case complexity of the problem grows exponentially. For example, if S is a set of n points in IRd then the Delaunay triangulation of S can consist of up to \Theta (nb d+1 2 c) faces (see [74]). Although the running time of the programs constructing ff-shapes will get substantially worse as d increases, there might be applications that warrant implementations in low dimensions higher than three. 3The definitions in [26] are slightly different from the ones in this thesis. In particular, their ff is the same as minus one over ff in this text. Thus, our negative values of ff correspond to their positive values. 24 Chapter 3 Data Structures Two different data structures are used to store the whole family of ff-shapes of a given data set S. One represents the connectivity and order among the simplices of the three-dimensional Delaunay triangulation D of S. The other is used to implement the face classification of section 2.3. A trianglebased data structure is used to store D. This is described in section 3.1. The collection of intervals assigned to the simplices of D is stored in a few linear arrays. This is discussed in section 3.2. 25 3.1 Triangle-Edge Data Structure The data structure presented in this section is designed to represent the topology and order of a three-dimensional triangulation. The oldest and most widely used representation for the topological structure of a cell complex is the incidence graph (see for example [21]). This graph contains a node for every cell (of any dimension) in the complex, with arcs connecting nodes if the corresponding cells are incident to each other and their dimensions differ by one. This approach stores the cells explicitly but does not provide any direct access to ordering information. A fairly different approach is to represent the cells implicitly by a set of basic elements and operators. This idea was introduced with the quad-edge data structure by Guibas and Stolfi [42], proposed for modeling two-dimensional manifolds and their duals. Faces and vertices are represented implicitly by rings of directed and oriented edges. A set of traversal functions defines an abstract edge algebra. The edge-facet structure due to Dobkin and Laszlo [17] extends this idea to three-dimensions. Brisson [5] generalized this idea to subdivided manifolds of arbitrary dimensions. The triangle-edge data structure presented here is a simplification of Dobkin and Laszlo's edge-facet structure. By reducing the scope to exactly what is needed, namely, simplicial complexes without duals, the compactness as well as the response time of the structure can be improved dramatically. Data structure. The triangle-edge data structure is triangle based. Assume a unique index f for each triangle in a triangulation T of a set of points S = fp1; p2; : : : ; pngIR3. The data structure's atomic unit is the so-called triangle-edge pair a = hf; vi, with 0 ^ v ^ 5, which identifies one of six versions of triangle f = oefpi;pj;pkg. Each version corresponds to one of the six directed edges incident to f (see figure 3.1). Sometimes it will be convenient to associate edges or triangles with triangle-edge pairs. We write a = hf; vi j ei;j j fi;j;k, where ea = ei;j is the directed edge from pi to pj, and fi;j;k is the oriented triangle determined by the sequence of vertices (pi; pj; pk). Let fi;j;k correspond to a = hf; vi. The positive halfspace of a triangle-edge a, denoted by a+, is the set of all points x 2 IR3 which see the endpoints of fi;j;k make a counterclockwise turn; that is, (x; pi; pj; pk) has positive direction. Each triangle f defines two edge rings, one traversing the edges of f in counterclockwise direction, and one in clockwise direction, as seen from hf; 0i+. Similarly, each edge e defines two triangle rings traversing the triangles incident to e in the two opposite directions. Each triangle-edge a belongs to exactly one edge ring and exactly one triangle ring (see figures 3.1 and 3.2). 26 The following basic triangle-edge functions are defined: ffl a:org = i if pi is the origin of the directed edge ea. ffl a:sym = b if ea = ei;j and eb = ej;i, that is, ea and eb are the two directed versions of the same undirected edge. ffl a:enext = b if ea and eb are both in the same edge ring and eb is the successor of ea. ffl a:fnext = b if ea and eb are in the same triangle ring and eb is the successor of ea. q q q CC CC C'' '' 'pi pj pk hf; 0i j ei;j enext\Gamma ! q q q CC CC C'' '' 'pi pj pk CC CC CO hf; 2i j ej;k enext\Gamma ! q q q CC CC C'' '' 'pi pj pk' '' ''/ hf; 4i j ek;i enext\Gamma ! hf; 0i l sym l sym l sym l sym q q q CC CC C'' '' 'pi pj pk oe hf; 1i j ej;i enext\Gamma q q q CC CC C'' '' 'pi pj pkC CC CCW hf; 3i j ek;j enext\Gamma q q q CC CC C'' '' 'pi pj pk '' '' '7 hf; 5i j ei;k enext\Gamma hf; 1i Figure 3.1: The six versions of triangle fi;j;k, connected in two edge rings. qp i (a:org = i) qpj qpk0 qpk1 qpk 2 qpk` \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi * \Theta \Theta \Theta \Theta \Theta \Theta \Theta \Theta \Theta \Theta \Theta \Theta \Theta \Theta aaa aaa aaaa aa A AA AA AA AAA iiii ii CC CC CC C a j fi;j;k0 a:enext j fj;k0;i a:fnext j fi;j;k1 \Delta \Delta \Delta Figure 3.2: The basic triangle-edge functions. The eyes of the reader are in half-space a+. The triangle ring of a is hhfi;j;k0; fi;j;k1; fi;j;k2; : : :; fi;j;k`ii. The triangle ring of a:sym j fj;i;k0 is in opposite direction: hhfj;i;k0; fj;i;k`; : : :; fj;i;k2; fj;i;k1ii. 27 Figure 3.2 illustrates these functions. The postfix form "a:func" and its equivalent functional form "func(a)" are used as synonyms throughout this text. Exponents are used to express functional concatenation; for example, "a:enext2" is equivalent with "a:enext:enext." Also, "a:dest" is used as an abbreviation for "a:sym:org." Consistency. Table 3.1 lists some of the axioms that hold for the edge-triangle functions. These relations can be used to test the data structure for consistency. Of course, axioms involving enext traversals are trivially true since, as we will see below, our implementation has the edge rings "hardwired." However, all axioms concerning fnext must hold if the facet rings are properly linked. These relations, together with, for example, Euler formulas and their derivatives, are useful to check whether the current state of the data structure is indeed a valid representation of a triangulation. For example, checking axioms (T13) and (T14) makes sure that the facets of all tetrahedra are properly linked. Implementation. The internal representation of the triangle-edge functions takes advantage of the fact that the structure of all edge rings is identical. Therefore, it can be implemented by a global lookup table ve[0::5] with six entries ve[v], 0 ^ v ^ 5, such that, a:enext = ht; ve[v]i, for any triangle-edge pair a = ht; vi. However, actual pointers are needed for the successor relations in the (T0) a 2 T () a:org 2 S and a:sym; a:enext; a:fnext 2 T (T1) a:sym2 = a (T2) a:enext3 = a (T3) a:enexti 6= a for i = 1; 2 (T4) a:sym:enexti 6= a for i = 1; 2; 3 (T5) a:sym:enext:sym:enext = a (T6) a:fnexti = a for some i ! 1 (T7) a:sym:fnexti 6= a for all i (T8) a:sym:fnext:sym:fnext = a (T9) a:org 6= a:sym:org (T10) a:enext:org = a:sym:org (T11) a:enext2:org 6= a:org (T12) a:fnext:org = a:org (T13) a+ " jjjT jjj 6= ; =) a = a:turn3; with a:turn = a:enext:sym:fnext:enext (T14) a+ " jjjT jjj 6= ; =) pa:fnext:enext2:org 2 a+ (T15) x 2 a+ =) sequence (x; pa:org; pa:enext:org; pa:enext2:org) has positive orientation Table 3.1: Some axioms for the triangle-edge functions. 28 triangle rings. Note that there are six versions of a triangle but only three vertices (that is, possible origins). To save memory, a second lookup table vo[0::5] is used to map the indices corresponding to the origins of the six different versions of a triangle f to its base triangle hf; 0i. See table 3.2 for the internal declaration of the data structure, and the pseudocode of the basic triangle-edge functions. In total, there are nine fields per record for a triangle f : three origin fields for the indices, i, j, and k, where hf; 0i j fi;j;k, and six next.f fields for the pointers to the next triangle in the triangle rings of hf; vi, for 0 ^ v ^ 5. For convenience, assume a (dynamically allocated) linear array tr[0::\Lambda ] to type Triangle = record origin: array [0::2] of Integer; next.f: array [0::5] of Integer; end; const vo: array [0::5] of Integer f0; 1; 1; 2; 2; 0g; ve: array [0::5] of Integer f2; 5; 4; 1; 0; 3g; var tr: dynamic array [0::\Lambda ] of Triangle; function org(a) is assert a = hf; vi; return tr[f ]:origin[vo[v]]; end; function sym(a) is return a + (if a is odd then \Gamma 1 else 1); end; function enext(a) is assert a = hf; vi; return hf; ve[v]i; end; function fnext(a) is assert a = hf; vi; return tr[f ]:next.f[v]; end; Table 3.2: The basic triangle-edge functions. 29 store all the triangle records. Moreover, we map triangle-edge pairs to single integers, which allows us to handle all references to triangle pairs, and not just triangles, as integers: hf; vi j 8f + v. Note how this conversion formula enables the usage of fast shift operations instead of multiplication. This would not be possible for the otherwise more natural conversion formula 6f + v. Since lemma 2.4.1 states that triangulations of S have O(n2) triangles in the worst case, it makes sense to assign only half as many bits for each origin field than for each next.f field.1 This gives a total of 30 bytes for a Triangle record, assuming the usual 32 bits integers for edge-triangle pairs. Using 32-bit origin fields increases the record size to 36 bytes. This lifts the restriction on the maximum number of vertices for all practical purposes, provided that the complexity of the triangulation does not show quadratic worst-case behavior. The sign bits of the triangle-edge fields can be used to encode markers. For the time being, assume only the sign bits of the origin fields to be used for that purpose. This gives the user the possibility to mark each triangle with up to eight different "colors," where three of them are predefined: BLANK for unmarked triangles, and HULL or SYMHULL for triangles f bounding the convex hull of S, depending on whether hf; 0i+ " conv(S) or hf; 1i+ " conv(S) is ;. Creating and deleting triangles is easy when the records of the linear array tr[0::\Lambda ] are managed in a freelist. To delete a triangle, mark its record as "not active" by simply clearing its origin fields, and add it to the freelist. To create a new triangle fi;j;k, retrieve an unused index f from the list, set its origin fields to i, j, and k, and initialize all its triangle rings to be trivial , that is, set the next.f fields such that hf; vi:fnext = hf; vi; for 0 ^ v ^ 5. Ring manipulation. A most useful operation is fsplice, which modifies triangle rings. It takes two triangle-edge pairs a and b as arguments and returns no value. However, assuming that hha, a:fnext, a:fnext2, : : : , a:fnext`\Gamma 1ii denotes the triangle ring of a, with a:fnext` = a, the call of operation fsplice(a; b) has the following side effects on the rings of a and b: hh: : : : : :-- -z "" i ; a; iiz ""-- -: : : : : : : : :; b; : : :: : :-- -z "" iii ii fsplice(a;b)() hha; : : :: : :-- -z "" iii ; : : : : : :-- -z "" i ii 6j hhb; iiz ""-- -: : : : : : : : :ii (a) If the two rings are distinct, it combines them into one ring. (b) If the rings are identical, it breaks them into two distinct rings. 1Let M be the maximum value for Integer numbers. To be save also in the worst case, n has to be restricted so that 8(n2 \Gamma 3n) + 5 ^ M. For M = 231 \Gamma 1 (integers), this gives n ^ 16385. For M = 232 \Gamma 1 (unsigned integers), this is n ^ 23171. Both bounds are smaller than 215 \Gamma 1. 30 It is worthwhile to mention, that fsplice is the inverse to itself. The implementation of fsplice can be found in table 3.3. Also, observe how fmerge(a; b) = fsplice(a; b:sym:fnext:sym) combines aessoe,rr aessoe,r r aessoe,rr aessoe,r r aessoe,r r aessoe,rr aessoe,r r aessoe,rr6 ? -oe ? 6 oe 6 ? -oe ? 6 oe a b fmerge(a; b)=) aessoe,rr aessoe,r r aessoe,rr aessoe,r r aessoe,r r aessoe,rr aessoe,r r aessoe,rr6 ? -oe oe -oe ? 6 oe -oe oe a b aessoe,rr aessoe,r r aessoe,r r aessoe,rr aessoe,rr 6 ? -oe ? 6 oe a c fmerge(c; c)(= fmerge(a; c)=) aessoe,rr aessoe,r r aessoe,r r aessoe,rr aessoe,rr 6 ? -oe oe oe ?oe a c Figure 3.3: The effects of fmerge(a; b) = fsplice(a; b:sym:fnext:sym). procedure fsplice(a; b) is swap pointers a; b; swap pointers a:fnext:sym; b:fnext:sym; end; refinement swap pointers a; b is assert a = hfa; vai and b = hfb; vbi; tr[fa].next.f[va] $ tr[fb]:next.f[vb]; end; procedure delete(a) is assert a = hf; vi and b = a:enext and c = a:enext2; fmerge(a; a); fmerge(b; b); fmerge(c; c); for v = 0; 1; 2 do tr[f ]:origin[v] 0 end for; add t to freelist; end; Table 3.3: More pseudocode for triangle-edge functions. 31 the two rings of a and b if they are different (see figure 3.3). When fmerge is applied to two identical arguments, like in fmerge(a; a), where a denotes some triangle-edge in a nontrivial triangle ring, it unlinks a from this ring. This fact is used in the code of delete, which can be found in table 3.3. On the other hand, if a = hf; vi denotes a trivial triangle ring, and b is some other ring, then fmerge(a; b) adds the triangle f to the ring of b. procedure collect(ph; a; a0; L) is assert a = hf; vi j fi;j;k and (ph; pi; pj; pk) has positive orientation; if f is marked with HULL or SYMHULL then unmark f ; L L [ fag; collect dfs(ph; a:fnext:sym; a0; L); collect dfs(ph; a:enext:fnext:sym; a0; L); collect dfs(ph; a:enext2:fnext:sym; a0; L); else L fa; a:fnext:sym; a:enext:fnext:sym; a:enext2:fnext:symg; end if; end procedure collect dfs(ph; a; a0; L) is assert a = hf; vi j fi;j;k; if f is marked with HULL or SYMHULL then if (ph; pi; pj; pk) has positive orientation then unmark f ; L L [ fag; collect dfs(ph; a:enext:fnext:sym; a0; L); collect dfs(ph; a:enext2:fnext:sym; a0; L); else a0 a; end if; end if; end Table 3.4: Collecting triangles visible from point ph. 32 Usage. This section closes with a couple of sample programs using the triangle-edge functions of the data structure. These procedures will be used for the construction of Delaunay triangulations in section 4.1. The idea is to add a point ph, which is not yet a vertex of the current triangulation T , to the triangulation by either attaching it to the outside of the convex hull or by inserting it inside the appropriate tetrahedron. We assume that the current data structure is a valid representation of T , and that all triangles bounding (the convex hull of) T are properly marked with either HULL or SYMHULL; call such a triangle a hull triangle. Furthermore, a facet a = oeT is said to be visible from pi, if oeT [fpig does not intersect any triangle of T , and pi 2 a+, the positive side of a. Procedure collect(ph; a; a0; L) is to collect all triangles that are visible from ph, and to store them in a set L. Table 3.4 lists the code. We assume that an oracle provided an initial triangle-edge a which corresponds to one such triangle. In other words, a j fi;j;k, the sequence (ph; pi; pj; pk) has positive orientation, and one of the following two cases applies. 1. ph 62 jjjT jjj. Then a lies on the boundary of the convex hull, and collect will use depth-first search to add all visible hull triangles to L. As a side effect, their hull bits will be deleted. This will tear a hole into the boundary of jjjT jjj, which has to be fixed afterwards. For that purpose, we maintain a hull triangle a0 which is not visible from ph, but lies on the boundary to the region covered by the triangles of L (see below). 2. ph 2 jjjT jjj. In this case, ph lies inside the tetrahedron oeT = ta+ , that is fpi; pj; pkg ae T and oeT " a+ is nonempty. In this case, the output parameter a0 will be undefined, and L will simply consist of the triangles of oeT . In order to insert the new vertex ph into the current triangulation, collect(ph; a; a0; L) is used in conjunction with a procedure connect(ph; a0; L) which adds the tetrahedron th;i;j;k to T for every triangle fi;j;k 2 L. The full pseudocode of this procedure is omitted; however, some crucial details are provided in the following. For the time being, assume that ph 62 jjjT jjj. Furthermore, let the triangles of L be processed in the same order as they were gathered by collect (that is, in depth-first order). With this, fi;j;h 62 T conveniently implies that fj;k;h 62 T and fk;i;h 62 T . Consulting figure 3.4, five case can be distinguished: 33 A. fi;j;h 62 T (which implies fj;k;h 62 T and fk;i;h 62 T ). All three missing triangles have to be added to T . Also, we need to perform appropriate modifications to the corresponding triangle rings. B. fi;j;h 2 T and fj;k;h 62 T and fk;i;h 62 T . Two missing triangles have to be added, and their triangle rings modified. C. fi;j;h 2 T and fj;k;h 2 T and fk;i;h 62 T . One new triangle has to be created, with the appropriated triangle rings. D. fi;j;h 2 T and fj;k;h 62 T and fk;i;h 2 T . This case is the mirror image of case C. E. fi;j;h 2 T and fj;k;h 2 T and fk;i;h 2 T . In this case, all adjacent triangles have already been created at an earlier stage. However, it is possible that the triangle rings are not yet properly closed. In some sense, only after this is done, the tetrahedron th;i;j;k is really added to T . Obviously, it is not enough to just create the new triangles, they also have to be linked into the data structure by appropriate modifications to the corresponding triangle rings. Since the details of this procedure are tedious and repetitive, only case B will be considered as an example. Table 3.5 lists the appropriate code segment. Two triangles, fj;k;h and fk;i;h, have to be created, and five triangle rings, corresponding the edges ej;k, ek;i, ei;h, ej;h, and ek;h, have to be altered. Note that the triangle rings of the edges of the new triangles are initially trivial. Thus, x:sym:fnext:sym = x, and calls to fmerge in table 3.5 can be replaced by calls to fsplice with the same parameters. pi pj pk ph (a) pi pj pk ph (b) pi pj pk ph (c) pi pj pk ph (d) pi pj pk ph (e) Figure 3.4: The five cases when connecting ph to a j fi;j;k. 34 To reiterate, connect(ph; a0; L) scans through L, the list of triangles visible from ph, and processes them in the same order as they were generated by collect(ph; a; a0; L), starting with some visible hull triangle a. For every visited triangle, it then distinguishes the five cases of figure 3.4 and executes a series of fmerge or fsplice operations, following the pattern of table 3.5. After this, the facets of every tetrahedron th;i;j;k, for fi;j;k 2 L, are properly added to the triangulation T . Since we were adding the point ph outside jjjT jjj we still have to fix the hull bits of the new triangles on the convex hull. A simple depth-first traversal starting at the hull triangle a0, which was provided by collect as not visible from ph, but initially adjacent to some triangle in L, can do the job. Finally, observe that the above procedure works independently of whether ph 62 jjjT jjj or ph 2 jjjT jjj. In the latter case, ph lies inside some tetrahedron, and L, as generated by procedure connect, consist of its four triangles. They are processed by connect, which adds the new triangles after encountering the cases A, B, C, and E. Overall, the time complexity of collect and connect is proportional to the number of triangles visible from ph. This concludes the discussion of the edge-triangle data structure. This structure and its traversal functions will be used in section 4.1, where we construct three-dimensional Delaunay triangulations. connect refinement case B(ph; a) is assert a j fi;j;k 2 L and fi;j;h 2 T and fj;k;h 62 T and fk;i;h 62 T ; b fj;k;h; c fk;i;h; fmerge(a:enext; b); f ej;k g fmerge(a:enext2; c); f ek;i g fmerge(a:fnext:sym:enext; c:enext); f ei;h g fmerge(a:fnext:enext:sym; b:enext2); f ej;h g fmerge(b:enext; c:enext2:sym); f ek;h g end Table 3.5: Connecting ph to a j fi;j;k; case B. 35 3.2 Alpha Shape File and Data Structure According to chapter 2, it is possible to represent the whole family of ff-shapes of S by the faces or simplices of the Delaunay triangulation D, augmented by a set of intervals. How to store these intervals is the topic of this section. Figure 3.5 shows the overall layout of the data structure, which is discussed in more detail below. Since it is often convenient to have the interior of a shape triangulated, we represent ff-complexes, rather than ff-shapes; see section 2.2. As already discussed, there is a single interval for each simplex oeT 2 D so that oeT belongs to the ff-complex Cff if and only if ff is contained in this interval. Each face is classified as either interior, regular, or singular, which breaks down the interval for each face into three (possibly empty) subintervals; see table 2.2, section 2.3. Face indexing. The data structure of section 3.1, which is used to store D, is triangle based; vertices, edges, and tetrahedra are only represented implicitly. However, some sort of indexing for faces of all dimensions is needed in order to conveniently store additional information for each face (that is, the above mentioned intervals) in linear arrays. For vertices and triangles this is easy. SoS already provides a unique index for each point or vertex. Triangles can be indexed as in the triangle-edge data structure. Hash tables are used to index edges and tetrahedra. For each edge and tetrahedron, the hash function provides a unique index between 1 and me ? ne = jF1j and between 1 and mt ? nt = jF3j, respectively.2 Note that there will be unused indices, corresponding to empty records in the hash tables. The number of unused indices depends on the given load factor (for example, the load factor may be set to nm = 0:8). Let E(oeT ) be the set of all edge-triangle pairs that correspond to a simplex oeT 2 D. For a tetrahedron oeT j ta+ these are the corresponding four edge-triangle pairs; and for an edge oeT j ea, we set E(oeT ) = fa:fnexti; a:sym:fnexti j i = 1; 2; 3; : : :g; that is, the edge-triangle pairs of the facets in the two triangle rings (in opposite direction) corresponding to ea. We define the minimum edge-triangle pair representing an edge oeT j ea, or a tetrahedron oeT j ta+ , as amin = minfhf; vi j hf; vi 2 E(oeT )g.3 2The actual ranges of the hash indices for edges and tetrahedra are 0 : : : me \Gamma 1 and 0 : : : mt \Gamma 1. For convenience, we add 1 to hash indices to get face indices starting at 1. 3It is easy to order the edge-triangle pairs. After all, we are already representing them as integers: hf1; v1i ! hf2; v2i if and only if 8f1 + v1 ! 8f2 + v2; compare this with section 3.1. 36 1 mt %\Gamma rank Tetrahedron rank table 1 nf %\Gamma rank _\Gamma rank _\Gamma rank Triangle rank table 1 me %\Gamma rank _\Gamma rank _\Gamma rank Edge rank table 1 n _\Gamma rank _\Gamma rank Vertex rank table 1 n Vertex coordinates 1 nf Triangle-edge data structure 1 n SoS parameter matrix 0 mt \Gamma 1 amin Tetrahedron hash table 0 me \Gamma 1 amin Edge hash table 1 r Master pointer ff-threshold Spectrum 1 e Master list : : : : : : z ""-- - Face index Face type Rank type EoSL flag Master entry -- -z "" Sublist for rank r -- -z "": : : r + 1 Rank r 6 ffrr 6 ffr+1r 6 : : : : : : r r ? r ? r ? ? roe roe oe Figure 3.5: The ff-shape file and data structure. 37 Note that amin can be computed in constant time for each tetrahedron t+a , and in constant amortized time for each edge ea in D. Since amin is unique for each face, it can be used as a key for the hash functions. Moreover, it also allows the implementation of inverse indexing: given a face index i, we can obtain an edge-triangle pair representing the face; either by the amin reference in the (i \Gamma 1)-st entry of the corresponding hash table, or, in the case of triangles, by directly referring to the i-th triangle record. Thus, we achieve a unique (and reversible) indexing scheme for all faces of D. This will be used in the following to build the linear arrays storing the intervals for the face classification of section 2.3. Alpha ranks. Recall that the set of ff-thresholds consists of all end points of the intervals in table 2.2. The sorted sequence of ff-thresholds, ff1; ff2; : : : ; ffr, with ffr ! ffr+1 for 1 ^ r ^ r, forms the ff-spectrum. For convenience, we always assume ff1 = 0 and ffr = 1; if 0 and 1 do not occur in the set of ff-thresholds, these values are added to the spectrum. We say the threshold ffr is of rank r. Each rank stands for an end point of a certain value. This is used to store the intervals of the faces of D in four so-called rank tables, one for each dimension; see figure 3.5. For edges and triangles we need to store the ranks of the %, _, and _ values, for vertices the ranks of _ and _, and for tetrahedra only the rank of %; compare this with the intervals in table 2.2. Note that certain intervals are empty. A "void" rank 0 is used to encode these; again, compare this with table 2.2. ffl Attached edges and triangles have %-rank = 0. ffl Faces bounding conv(S) have _-rank 6= 0 and _-rank = 0. ffl Furthermore, redundant (or dumped) vertices and empty entries in the edge and tetrahedra rank tables are marked with _-rank = _-rank = 0.4 Master list. In general, end points of more than one ff-interval will be mapped to the same ff-threshold in the spectrum. For many application, such as computing "Betti numbers" of ffcomplexes [15], or computing the volume, surface, and other metric features of ff-diagrams or space filling diagrams [25], it is desirable to be able to retrieve, given a rank r, all faces whose intervals have ffr as an end point. In other words, for each ff-threshold ffr in the spectrum, we need to store the list of faces whose intervals have this threshold value as an end point. For convenience, we merge these lists in order of increasing rank to a list of sublists, the so-called master list; see figure 3.5. Each entry stores the corresponding face index, the face type (that is, the dimension of the face), the rank type (the entry may have type %, _, or _), and a flag indicating the end of a sublist. The faces in each sublist are sorted in order of non-decreasing dimension. 4The concept of "redundant" points or vertices occurs in the context of weighted ff-shapes and regular triangulations; see section 4.1.1. Even in the case of Delaunay triangulations, our actual implementation chooses to dump duplicate points. Neither redundant nor dumped vertices occur in D. 38 The master list can be implemented as a linear array, provided that the spectrum stores an index for each threshold ffr, pointing to the corresponding sublist within the array. Note that, if we scan the spectrum and master list in order of decreasing ranks, we get the faces in the same order as they are removed from the ff-complex when "shrink wrapping" the ff-shape from the convex hull conv(S), for ffr = 1, to the point set S itself, for ff1 = 0. Also notice that the sublists of ff1 = 0 and ffr = 1 are empty in the nondegenerate case. Furthermore, observe that we do not explicitly store the % = 0 values of the vertices. Thus, for ff = 0, the set of vertices is initialized with S, by default.5 Section 4.3 will list the formulas needed to compute the ff-thresholds, and section 5.2 will discuss their robust implementation using the SoS perturbation scheme. Here, we only mention that it is neither necessary nor desirable to maintain the perturbation for the values of the ff-spectrum: simple floating-point values will do. In fact, SoS uses exact arithmetic, and the resulting ff-thresholds are rationals formed as ratios of long-integer numbers. These exact threshold values are used to sort the spectrum consistently. After this, however, they are discarded, and mapped to floating-point values. This is done so that ffl different actual values (that is, rationals) for ff-thresholds with the same floating-point value are mapped to different ranks, while ffl ff-thresholds, which are only different with respect to the perturbation, are mapped to the same rank. To conclude the discussion on the data structures used to represent the family of ff-shapes for a point set S 2 IR3, it should be noted that interval trees (see for example [67]) can be used as an alternative to (or in addition to) the master list. In particular, rendering shapes of arbitrary ff will be more efficient with interval trees. However, the master list is more compact, and thus better suited for storage purposes; besides, it is always possible to reconstruct the interval trees from the master list, if desired. 5This is no longer valid for weighted ff-shapes. There, the %-values of vertices are not necessarily all zero and their actual values have to be stored in the vertex rank table as well. The sublist for ff1 = 0 is no longer necessarily empty in the non-degenerate case. 39 Chapter 4 Algorithms The major challenge in computing three-dimensional ff-shapes is the construction of Delaunay triangulations of spatial point sets. This problem is discussed in section 4.1. Section 4.2 then shows how to compute the additional information which is needed for the face classification in section 2.3. Finally, a list of all the required geometric primitive operations, including a revised general position assumption, will be given in section 4.3. We assume general position of the point sets throughout this section. 40 4.1 Constructing Delaunay Triangulations The algorithm discussed in this section, is based on local transformations or flips (as defined below) and goes back to Lawson [51, 52]. He introduced this method (also known as Lawson's flip method) for constructing two-dimensional Delaunay triangulations in 1972. Given a finite point set S 2 IR2, the method initially constructs an arbitrary triangulation T of S. This triangulation is then altered step by step through a sequence of edge flips until all edges are Delaunay edges, and thus the Delaunay triangulation D is obtained. Flips in IR2.Flips in IR2.Flips in IR2. Consider an edge ei;j = oefpi;pjg in T . If ei;j belongs to the boundary of conv(S) then it also belongs to the Delaunay triangulation of S. Otherwise, it is incident to two triangles, fi;j;k and fi;j;u, in T . We call ei;j locally Delaunay if it belongs to the Delaunay triangulation of the points pi, pj, pk, and pu. If ei;j is not locally Delaunay then it cannot belong to the Delaunay triangulation of S either. In this case, we can replace ei;j by ek;u = oefpk;pug since the union of fi;j;k and fi;j;u is convex; see figure 4.1. This operation is called an edge flip in IR2. It is easy to see that ek;u is locally Delaunay. Again, it helps to use the lifting map of section 1.1 to illustrate this procedure; of course, this time, it is a projection of points in IR2 onto U , the paraboloid of revolution in IR3. The triangulation T is the Delaunay triangulation D if and only if its projection TU is identical to the lower boundary of conv(SU ); call it DU . If an edge ei;j is not locally Delaunay, then its projection is not part of DU . Therefore, it is possible to attach a tetrahedron underneath the edge and its two incident triangles. In other words, flipping edges in IR2 corresponds to attaching tetrahedra in IR3. Clearly, this process has to end sometimes. A classical lemma by Delaunay [14] implies that this happens only when TU = DU , that is, when T is a Delaunay triangulation. Also note, that the upper bound theorem implies that no more that O(n2) tetrahedra can be attached in the worst case, which makes Lawson's flip method a quadratic algorithm. pi pi pj pj pk pk pu pu Figure 4.1: Edge flipping in IR2. 41 Flips in IR3.Flips in IR3.Flips in IR3. In order to generalize the flip algorithm to three dimensions, we have to first specify what we mean by a "flip" in IR3. Consider a set S0 of five points in IR3 in general position. If there exists a point p 2 S0 such that conv(S0) = conv(S0 \Gamma fpg), then S0 has only one possible triangulation with four tetrahedra, all incident to p. This triangulation is trivially the Delaunay triangulation of S0. Otherwise, all points of S0 lie on the boundary of conv(S0), and S0 has only two possible triangulation, one with two tetrahedra, and one with three tetrahedra: one of them has to be the Delaunay triangulation; see figure 4.2. Similar as for the two-dimensional case, the lifting map is used for illustration; this time again as a map from IR3 to IR4 like in section 1.1. The two triangulations of S0 correspond to the lower and upper boundary of conv(S0U ), which is a four dimensional simplex. The triangulation corresponding to the lower boundary is the (closest-point) Delaunay triangulation. The other one is furthest-point Delaunay triangulation. A flip is the operation that replaces the triangulation that corresponds to the upper boundary by the one that corresponds to the lower boundary. Observe that a flip is defined with a direction, so that it tends toward the Delaunay triangulation of the set S. This idea can be extended to arbitrary dimensions [28]. pi pi pj pj pk pk pu pu pv pv Figure 4.2: The two triangulations of set S = fpi; pj; pk; pu; pvg. The triangulation to the left consists of two tetrahedra: ti;j;k;u and ti;j;k;v, while the triangulation to the right consists three tetrahedra: tu;v;i;j, tu;v;j;k, and tu;v;k;i. 42 Local Delaunayhood. A triangle oeT = fi;j;k, so T = fpi; pj; pkg, is called locally Delaunay if either (a) fi;j;k 2 conv(S), or (b) fi;j;k is the shared facet of two tetrahedra oeT 0 and oeT 00 in T , with T 0 = T [ fpug and T 00 = T [ fpvg, and pv lies outside @bi;j;k;u, the circumsphere of oeT 0. Note that, pv 62 bi;j;k;u if and only if pu 62 bi;j;k;v. As illustrated in figure 4.2 we distinguish two different ways to flip faces that are not locally Delaunay: 1. Let fi;j;k be a triangle in T which is not locally Delaunay, and let ti;j;k;u and ti;j;k;v be the two incident tetrahedra. If the union of the two tetrahedra is convex, then they can be replaced by three tetrahedra: tu;v;i;j , tu;v;j;k, and tu;v;k;i, and all of them are locally Delaunay. This is called a triangle(-to-edge) flip. 2. Let eu;v be an edge in T which is not locally Delaunay and has three incident tetrahedra tu;v;i;j , tu;v;j;k, and tu;v;k;i. They can be replaced by two tetrahedra ti;j;k;u and ti;j;k;v which are locally Delaunay.1 This is called a edge(-to-triangle) flip. Incremental flip. Unfortunately, the straightforward generalization of Lawson's flip method to IR3, using the above local transformations, does not work. Starting with an arbitrary initial triangulation, it is indeed possible to get stuck with faces that are not Delaunay but are still not transformable. An example of such a non-transformable face is a triangle whose two incident tetrahedra form a nonconvex union. This was first observed by Joe [45]; however, he also proves the following lemma [46]. Lemma 4.1.1. If a single point p is added to the Delaunay triangulation of a point set S in IR3, then there exists a sequence of flips that results in the Delaunay triangulation of S [ fpg. This result is the basis for the incremental flip algorithm in three dimensions. It was extended to arbitrary dimensions in a paper by Rajan [68]. However, his arguments require a priority queue to sort the flips, which takes logarithmic time per flip. Edelsbrunner and Shah [28] eliminate this requirement by proving that incremental flipping works even for an arbitrary sequence of flips. As a matter of fact, they show that it even works for the more general problem of constructing regular triangulations in IRd. 1Note that the non-Delaunayhood of eu;v implies that ti;j;k;u [ ti;j;k;v is convex. 43 Randomization. A crucial step in all incremental flip algorithms is the point location problem, which occurs when a point p is added to the triangulation of S. Joe [46] avoids this problem by adding the points in a sequence sorted along a fixed direction, for example, the x1-axis, and shows that the resulting algorithm has a time complexity of O(n2), if S is a set of n points in IR3. This is optimal in the worst case. Guibas, Knuth, and Sharir [41] study the complexity of a similar algorithm in IR2 where points are added in a random sequence. They show that, while \Theta (n2) flips are required in the worst case, the expected number of flips is only O(n). Edelsbrunner and Shah [28] extend this result to a randomized incremental flip algorithm in IRd and prove a time complexity of O(nd(d+1)=2e) in the worst case, and O(n log n + ndd=2e) in the expected case, for the worst point distribution. The second term is of the same order of magnitude as the maximum number of possible simplices in the Delaunay triangulation. For IR3, this gives again time O(n2). However, it is possible to give a probabilistic time analysis that depends on the number of simplices in the final Delaunay triangulation. In [28] it is shown that the expected number of flips is about the same as the expected number of simplices in the final Delaunay triangulation. Empirical results with an existing implementation justify this; see [31]. Edelsbrunner and Shah's method also adds points in a random sequence, and uses the history of performed flips to guide the point location. Sections 4.1.1 and 4.1.2 present a variant of this algorithm. The main difference between our method and the original algorithm in [28] is that it takes no additional storage for solving the point location problem. In fact, point location is implemented in a straightforward, yet efficient, way. Empirically, this saves around 80% of memory over the method in [28], and is still 10-20% faster, when compared to the implementation of [31] (at least for a certain range of n; see section 7.1). An extension to regular triangulations is straightforward. 4.1.1 A Randomized Incremental Flip Algorithm The basic idea of the incremental flip algorithm is the following. S be a set of n points in IR3. Let 4 ! i ^ n and assume that the Delaunay triangulation of the first i \Gamma 1 points in S 2 \Gamma IR 3 n \Delta is already constructed; call it Di\Gamma 1. Add the i-th point pi 2 S to the triangulation, and restore Delaunayhood by flipping; this results in Di. Repeat this process until i = n. If implemented correctly, this strategy always leads to the Delaunay triangulation of S; see [28, 46]. As mentioned earlier, there is the problem of locating point pi in triangulation Di\Gamma 1. Joe [46] avoids this by adding the points in a sequence sorted along the x1-axis. Thus, pi lies always outside jjjDi\Gamma 1jjj, and it is easy to find the faces visible from pi. This approach does no longer work, when the algorithm is to be randomized. Rather, one has to be prepared to add points pi inside and outside jjjDi\Gamma 1jjj. Edelsbrunner and Shah [28] achieve this as follows. First, they assume an artificial "infinite" tetrahedron enclosing all points of S, and make it the initial triangulation. Therefore, all points of 44 S can be added inside the existing triangulation. Second, they use a directed acyclic graph (called the history dag) to store the history of all performed flips. This generalizes the two-dimensional technique presented in [41]. This history dag serves as a data structure for point location. The resulting algorithm is elegant and efficient in the expected case, but it has certain drawbacks. The infinite tetrahedron complicates the robust implementation of the geometric primitive operations considerably. Moreover, the resulting triangulation includes artificial triangles connecting each face of conv(S) with one of the points at infinity. These triangles will have to be removed in a post processing step. The history dag has the additional disadvantage in that it can grow quite large. Empirical studies with an existing implementation suggest that one has to expect 4 to 5 times as many triangles in the history dag than in the final Delaunay triangulation; see section 7.1 and [31]. The randomized algorithm presented below does neither assume an initial tetrahedron at infinity nor does it require any additional data structure for point location. Given a triangulation, it adds points that can lie lie inside simplices of the triangulation or outside all simplices. This is only slightly more complicated than adding points inside only, while simplifying the implementation of the geometric primitives. Section 4.1.2 will discuss how to perform point location without any additional storage. This method is simple and turns out to be fairly efficient; see section 7.1. Global algorithm. Table 4.1 shows the global structure of the algorithm in pseudocode notation. Let (p1; p2; : : : ; pn) be a random sequence of the points of S, and let T denote the evolving triangulation, as it is stored in the edge-triangle data structure of section 3.1. Initially, T consists of t1;2;3;4 and its faces. The main loop incrementally adds points pi to T . procedure delaunay(S; n) is let (p1; p2; : : :; pn) be a random sequence of the points in S; initially, T consist of t1;2;3;4 and its faces; for i 5; : : : ; n do assert T = Di\Gamma 1; a search(pi); collect(pi; a; a0; L); connect(pi; a0; L); flip(L; pi); end for; assert T = D; end Table 4.1: The global structure of the randomized incremental flip algorithm. 45 Again, we need some more terminology for further discussion. The terms triangles, facets, and edge-triangles (in the context of the data structure) are used synonymously. However, facets or edge-triangles are always meant to be oriented. As in section 3.1, a facet a = oeT is said to be visible from pi, if oeT [fpig does not intersect any triangle of T , and pi 2 a+, the positive side of a. Furthermore, call a = oeT a link facet of pi if oeT [fpig is a tetrahedron in T ; denote this by a 2 L(pi). Function search(pi) implements the point location step returning an edge-triangle a 2 T visible from pi; see section 4.1.2 for more details. If pi lies outside jjjT jjj, then a will lie on the boundary of jjjT jjj, oriented such that a+ " jjjT jjj = ;. Otherwise, pi 2 jjjT jjj, and pi lies inside ta+ , the tetrahedron incident to a on the positive side of a. In both cases, procedure collect(pi; a; a0; L) collects all facets in T visible from pi. It uses depth-first search, starting with a; see table 3.4, section 3.1. Assume these facets are stored in the list L. With this, pi can be added to T by connecting the facets of L to pi, and L then becomes the set of initial link facets of pi. This is done by calling connect(pi; a0; L) as discussed in section 3.1. Flipping link facets. In general, after adding pi, there will be triangles in T which are not locally Delaunay. However, these triangles stay in some sense local around pi. They are either incident to pi, or part of the link L(pi). Table 4.2 lists the pseudocode of the procedure that restores the Delaunayhood of T by flipping the edges or triangles that are not locally Delaunay, until none such exist, and T = Di. The results of [28] prove the following: 1. The while-loop in table 4.2 always terminates with T = Di. 2. The correctness is independent of the order in which the facets are flipped. 3. It is sufficient to consider link facets only. A geometric test for local Delaunayhood (sphere test, see section 4.3) has to be performed, whenever a facet is selected from L. By the mechanism that administers link facets it is possible to run into facets that are no longer part of the triangulation. So, we first have to test if the facet still exists in T . During the flip procedure, triangles are continually created and destroyed. Deleted triangles are marked as such, and discarded when they later show up in L. Assume now a facet a 2 L was found that is not locally Delaunay. The algorithm tries to get rid of it by either flipping one of its edges e*, * = 1; 2; 3, or by flipping the triangle a itself. As we will see shortly (lemma 4.1.2), it is possible to test whether an edge is transformable without any geometric tests. Since geometric tests are computational expensive, we check transformability for edges first; this reduces the overal number of tests. If none of the e* is transformable, a series of geometric tests (orientation tests, see section 4.3) decides whether ta+ [ ta\Gamma is convex and triangle a can be flipped. 46 If so, a is flipped; if not, a is (temporarily) discarded as non-transformable. In any case, whenever a flip is performed, the link L(pi) changes and new link facets have to be added to L. As it turns out, there cannot be more than three new link facets per flip; see figure 4.3. Therefore, the total time required by the while-loop is proportional to the number of flips performed. procedure flip(L; pi) is assert L = L(pi) and T is current triangulation; while 9 a 2 L such that a 2 T and a is not locally Delaunay do remove a from L; let e* denote the *-th edge of a, * = 1; 2; 3; if 9 e* such that e* is transformable then flip edge e* in T ; else if a is transformable then flip triangle a in T ; end if; add new link facets to L, if any; end while end Table 4.2: Restoring Delaunayhood of triangulation T after adding pi. pu pi pv pj pw pk pi pu pj pv (a) (b) Figure 4.3: Face flips in IR3. (a) The triangle flip from fu;v;w 2 L(pi) to ei;j corresponds to deleting fu;v;w, and adding fi;j;u, fi;j;v, fi;j;w. The triangles fu;v;j, fv;w;j, fw;u;j are added to L(pi). (b) The edge flip from eu;v to fi;j;k corresponds to deleting three triangles fi;u;v, fj;u;v, and fk;u;v, (the last two are part of L(pi)), and adding fi;j;k. The triangles fj;k;u and fj;k;v need to be added to L(pi). 47 Call the number of triangles that are incident to an edge e the degree of e. An edge e = oeT is transformable if it is not locally Delaunay and has degree 3. If this is the case, an edge flip is performed, replacing three tetrahedra by two. With respect to the data structure storing the triangulation, this means that the three incident triangles, call them oeTj , with j = 1; 2; 3, are deleted, and triangle oeT 0, T 0 = (S Tj) \Gamma T , is created and inserted. This can be done with a sequence of fsplice operations, as discussed in section 3.1. Since e was an edge of a former link facet, pi 2 T 0 and there are two new link facets to be added to L; see figure 4.3 (b). Lemma 4.1.2. Consider a degree-3 edge in T which does not lie on the boundary of jjjT jjj. Either all three incident triangles are locally Delaunay or all three triangles are not locally Delaunay. Let fi;u;v, fj;u;v, and fk;u;v be the three triangles incident to eu;v. Lemma 4.1.2 follows directly from the definition of local Delaunayhood (again, it helps to look at figure 4.2 or 4.3): fi;u;v is locally Delaunay () pk 62 bi;u;v;j () pj 62 bi;u;v;k = bk;u;v;i fj;u;v is locally Delaunay () pi 62 bj;u;v;k () pk 62 bj;u;v;i = bi;u;v;j fk;u;v is locally Delaunay () pj 62 bk;u;v;i () pi 62 bk;u;v;j = bj;u;v;k Consider now an edge e* of a link facet a that is not locally Delaunay. If e* has degree 3 and does not lie on the boundary of jjjT jjj, lemma 4.1.2 implies that all three incident triangles are not locally Delaunay. Thus, e* is transformable, and no additional geometric test (in this case, orientation tests) is necessary to determine this. An edge flip is performed. If none of the edges of a non-Delaunay link facet are transformable, we have to try to flip the triangle itself. A triangle a = oeT is transformable if it is not locally Delaunay, and the union of the two incident tetrahedra oeT 0 and oeT 00 is convex. In this case a triangle flip can be performed, replacing two tetrahedra by three. There are three new link facets to be added to L; see figure 4.3. One last remark with respect to the list L that stores the link facets. It can be implemented as a stack, or queue, or whatever. There is no need to delete a specific edge-triangle from L. If we delete triangles during the flipping process, we rather mark them as "deleted" and later check the marks when selecting facets from L. However, in order to avoid duplicates in L, we also mark edge-triangles as "on the list" when they are added to L. Regular triangulations. It is relatively easy to modify the code in table 4.1 such that it works for regular triangulations as well. All we need to do is to identify redundant points and remove them from the triangulations. Observe that a non-locally Delaunay triangle and one of its edges cannot 48 be transformable at the same time. This is because the edge has degree 3 and thus the union of the two tetrahedra incident to the triangle is non-convex. Moreover, it is true that, in the context of Delaunay triangulations, there can be at most one transformable edge per triangle. In the context of regular triangulations, however, there can one or two such edges. In the latter case, they are incident to a redundant point. Further details are omitted. 4.1.2 A Simple Point Location Scheme Procedure delaunay calls a function search to locate a point ph in the current (Delaunay) triangulation T , assuming ph is not already a vertex of T . The algorithm whose implementation is studied in this thesis employs the following simple method: (i) Choose a random facet a 2 T such that ph 2 a+. (ii) If jT j " a+ = ;, then ph lies outside the current triangulation, and a is visible from ph. Otherwise, there is a unique tetrahedron ta+ 2 T . If ph 2 ta+ , then ph is located inside (the tetrahedron ta+ of) T . In both cases we are done. If ph 62 ta+ , then there exists some facet b of the tetrahedron ta+ , with the same orientation as a, and ph 2 b+. Repeat step (ii) with a = b. function search t(ph; a) is assert a j fi;j;k and ta+ = tu;i;j;k and (ph; pi; pj; pk) has positive orientation; if f is marked HULL then return a; else if (ph; pi; pj; pu) has positive orientation then return search t(ph; a:fnext); else if (ph; pj; pk; pu) has positive orientation then return search t(ph; a:enext:fnext); else if (ph; pk; pi; pu) has positive orientation then return search t(ph; a:enext2:fnext); else return a; end if; end if; end Table 4.3: Tail-recursive point location in a triangulation. 49 Function search t(ph; a), as listed in table 4.3, uses tail recursion to implement step (ii). Let (a0; a1; : : : ; a`) be the sequence of edge-triangles a with ph 2 a+ for which search t(ph; a) is called. Since T is a Delaunay triangulation, the volume of T`*=0 a+* is strictly decreasing, for increasing `. Thus, the function always terminates successfully. This simple idea of walking through the triangulation is not new, it was already mentioned in [4], with a reference to [39], where it was used for computing Delaunay triangulations in the plane. It is difficult to give a good estimate for the time complexity of search. In the worst case, it is quadratic in n, the number of points in T . A more realistic estimate is that the number of steps required to locate ph is proportional to n 1 3 (imagine T to be a cubic grid; see also [4, 39]). Empirically, it performs even better. This is especially so, if we choose the initial triangle from a random sample of size k. For example, k = c \Delta log(n), for some constant c. Then, we choose a "good" starting triangle according to some fast heuristics.2 With this, the time complexity of search t seems to be logarithmic in n, and roughly two times faster than the history dag method (at least for the observed range of n; see section 7.1). 4.2 Computing Alpha Intervals There is really not much to say about the algorithmic issues involved in computing the intervals of table 2.2, which are then used for the face classification for ff-shapes in section 2.3. It is pretty much straigthforward, and this section only outlines the involved steps, giving relevant cross-references to other sections in this thesis. The goal is to build the interval structure of section 3.2, which augments the triangle-edge data structure of section 3.1 storing the Delaunay triangulation D. Together, these two data structures provide a complete and compact representation of the whole family of ff-shape for a given point set S. The most difficult step in constructing the intervals is the computation of the radii of the faces in D, and the test whether or not a given face is attached. The relevant formulas are explained in section 4.3. Another aspect that has to be considered is the robust implementation of these formulas; this is discussed in section 5.2. Given all that, it is quite straightforward to construct the ff-shape file and data structure, using O(m log m) time and O(m) space, m denoting the number of simplices in D. 2Here, "fast" means faster than the actual orientation tests performed in the body of search t. Our actual implementation uses a floating-point code of the Euclidean distance function. In comparison with the exact long-integer arithmetic of SoS, this is computational inexpensive. 50 1. The hash tables implementing reversable index functions for edges and tetrahedra have to be computed. We use universal hashing (see for example [13]) to achieve this. An appropriate load factor of, for example, 80%, has to be specified. Since disk space is usually thight, we reconstruct the hash tables whenever the ff-shape files are read by an application program, rather than keeping them on file. 2. The %-values for each unattached simplex oeT 2 D are computed, using the formulas of section 4.3. For robustness, this is done employing SoS. However, the perturbation is only used for the attachment test. 3. The sorted sequence of %-values, that is, the spectrum, is computed. Although the spectrum stores the threshold values as simple floating-point numbers, we have to keep them temporarily as (exact) long-integer rationals, in order to achieve a consistent sorting. This operation can be quite memory intensive; see section 7.2. After the sorting, the long-integer ratios are mapped to floatingpoint numbers with the same total ordering, and the additional memory can be freed. All further computations are based on the ff-ranks of the spectrum; see section 3.2. 5. The _- and _-values are computed, based on the %-ranks of the individual faces. This is done for each Delaunay simplex, in order of decreasing dimension. This completes the construction of the rank tables that encode the intervals of section 2.3, and can be done in constant time for each simplex. 6. In a final step, we sort the temporary list of %-, _-, and _-values, keeping the references to the corresponding faces. The result is the master list. In addition, the entries of the sublists of each rank are sorted with respect to the dimension of the corresponding faces. As it turns out, the memory management in the above construction process is crucial. To cut down the notable space requirements, it is useful to compress the involved temporary data structures as well as all the structures that are finally put on file. The major source of savings are the type fields and flags in the master list; see section 3.2. One possibility is to break up the record structures and use bit fields provided by the C programming language. This results in fairly ugly code, but also saves roughly 18% of main memory, and 33% in file size.3 This completes the short overview of how to build the ff-shape file and data structure of section 3.2. As already mentioned, the crucial part are the formulas for the ff-thresholds. They are discussed in the next section. 3In total, we can save 40 \Gamma 50% in file size if we further compress the file structure using standard compression tools, like Unix's compress. 51 4.3 Geometric Primitives The set of primitive operations needed to compute ff-shapes in IR3 can be separated into two groups: one for constructing Delaunay triangulations, and one for computing the intervals for face classifications. Constructing Delaunay triangulations requires essentially two geometric tests. One test decides on which side of a plane, spanned by three points, a fourth point lies. The other test decides on which side of a sphere, spanned by four points, a fifth point lies. By the "side" of a sphere, we mean the inside or the outside. In order to generate the intervals and face classifications of section 2.3 (table 2.2), we need to compute the radius of the smallest circumsphere of a tetrahedron, triangle, or edge, and test whether a point lies inside or outside this sphere. While the two tests for Delaunay triangulations are fairly common in geometric algorithms and computer graphics, the operations involving smallest circumspheres of triangles and edges are rather specialized. All operations share the problem of degenerate cases, which we can ignore because of the general position assumption. This section gives a formula for each of the primitive operations mentioned above. Their actually implementation using SoS is discussed in chapter 5.2. Assume that the points of S are labeled as p1; p2; : : : ; pn, and that each point pi is given by the vector (ssi;1; ssi;2; ssi;3) of its three coordinates. To simplify the notation in the remainder of this section, we use minors, which are determinants of submatrices of a given matrix. For convenience, define ssi;0 = 1 for all i and use the following notation for minors: Mi1;i2;:::;ikj1;j2;:::;jk = det 0BBB BB@ ssi1;j1 ssi1;j2 : : : ssi1;jk ssi2;j1 ssi2;j2 : : : ssi2;jk ... ... . . . ... ssik;j1 ssik;j2 : : : ssik;jk 1CCC CCA : Orientation test. Let T = fpi; pj; pkg and define hT as the unique plane that contains all three points of T . This plane can be oriented if we replace the set T by the sequence T ; for example, T = (pi; pj; pk). A point pu lies on the positive side of hT if the sequence (pu; pi; pj; pk) has positive orientation; otherwise, it is said to lie on the negative side of hT . We get the following. T = (pi; pj; pk): pu lies on the positive side of hT () Mu;i;j;k1;2;3;0 ? 0; (4.3-1) and pu lies on the negative side if the determinant is negative. Notice that the determinant equals zero if and only if the points are in degenerate position, that is, they lie on a common plane. 52 Sphere test. Given a set T = fpu; pi; pj; pkg, we need to decide whether another point pv lies inside or outside the sphere @bT . We can assume that the degererate case where pv 2 @bT does not occur. A possible implementation of this test is discussed in [27] using an extension of the lifting map (as mentioned in section 2.4) to three-dimensional spheres. Consider the paraboloid of revolution U : x4 = P3`=1 x2` in IR4. A sphere @b with center c = (fl1; fl2; fl3) and radius ae is mapped to the hyperplane @bU : x4 = P3`=1(2fl`x` \Gamma fl2` ) + ae2. This hyperplane has the property that U " @bU projected along the x4-axis into the x1x2x3-space yields @b. Moreover, a point p lies inside (outside) @b if and only if pU lies vertically below (above) @bU. The resulting formula assumes that each point pi 2 S has a fourth coordinate ssi;4 = P3j=1 ss2i;j. T = fpu; pi; pj; pkg: pv lies inside @bT () Mu;i;j;k1;2;3;0 \Delta Mu;i;j;k;v1;2;3;4;0 ? 0 (4.3-2) The first minor, Mu;i;j;k1;2;3;0, is a corrective term that is necessary because the sphere @bT does not change if the points of T are permuted. The second minor, Mu;i;j;k;v1;2;3;4;0, expresses the fact that the lifting map transforms a sphere test in IR3 to a hyperplane test in IR4. Radius of a smallest circumsphere. Next, consider computing the radius %T of @bT , the smallest circumsphere of oeT , for all k-simplices oeT 2 D, with 1 ^ k ^ 3. The formulas for the square of %T are given in (4.3-3) through (4.3-5). Note that computing %2T will be sufficient for our purposes since %T can never be negative. We distinguish the cases when k = jT j \Gamma 1 is 1, 2, or 3. T = fpi; pjg: %2T = P 3`=1 iMi;j`;0j2 4 (4.3-3) T = fpi; pj; pkg: %2T = (P 3`=1 iMi;j`;0j2) \Delta (P3`=1 iMj;k`;0j2) \Delta (P3`=1 iMk;i`;0j2) 4 \Delta (iMi;j;k2;3;0j 2 + iMi;j;k 1;3;0j 2 + iMi;j;k 1;2;0j 2) (4.3-4) T = fpi; pj; pk; pug: %2T = iM i;j;k;u 2;3;4;0j 2 + iMi;j;k;u 1;3;4;0j 2 + iMi;j;k;u 1;2;4;0j 2 + 4 \Delta Mi;j;k;u 1;2;3;0 \Delta M i;j;k;u 1;2;3;4 4 \Delta iMi;j;k;u1;2;3;0j 2 (4.3-5) In order to explain these three formulas we introduce some notation. Let a; b; c; d be points in IR3. We write jabj for the length of the edge conv(fa; bg), and jabcj for the area of the triangle conv(fa; b; cg). In the case of two points, %T is the same as half the distance between pi and pj. Equation (4.3-3) follows because jpipjj2 = iMi;j1;0j 2 + iMi;j 2;0j 2 + iMi;j 3;0j 2: 53 To handle the case of a triangle spanned by T = fpi; pj; pkg, we use the formulas %T = jpipjj \Delta jpjpkj \Delta jpkpij4 \Delta jp ipjpkj and jpipjpkj2 = 14 \Delta (iMi;j;k2;3;0j 2 + iMi;j;k 1;3;0j 2 + iMi;j;k 1;2;0j 2) which can be found in any good mathematical handbook. Finally, we obtain (4.3-5) using the extension of the lifting map mentioned above. If @b is the sphere through points pi; pj; pk; pu then @bU is the hyperplane through the four points pi;U ; pj;U ; pk;U; pu;U . The equation for the hyperplane can be computed directly from the coordinates of the lifted points. From this equation it is easy to compute the center and the radius of @b. Attached and unattached edges and triangles. We still have to consider the problem of deciding whether an edge or triangle oeT 2 D is attached or not. By definition, oeT is attached if there is a oeR 2 up1(oeT ) so that the point in R \Gamma T belongs to bT . If oeT is an edge, say, T = fpi; pjg and R \Gamma T = fpkg, then this can be done by comparing %T with the distance between pk and pi+pj2 . Straightforward algebraic manipulations lead to the following equation. T = fpi; pjg: pk 2 bT () 3X `=1 i Mi;j`;0j 2 \Gamma 3X `=1 i Mi;k`;0 + Mj;k`;0j 2 ? 0 (4.3-6) Now let oeT be a triangle, for example, T = fpi; pj; pkg and R \Gamma T = fpug. To see whether or not the point pu belongs to bT , we compute the center c of the circumsphere @bR of the tetrahedron oeR. Observe that pu 2 bT if and only if c and pu do not lie on the same side of the plane through pi; pj; pk. In other words, we need to test whether or not the sequences (pu; pi; pj; pk) and (c; pi; pj; pk) have different orientation. Some rather tedious algebraic manipulations are needed to derive the following equation which expresses the test in terms of minors. T = fpi; pj; pkg: pu 2 bT () 8??!??: Mi;j;k;u2;3;4;0 \Delta Mi;j;k2;3;0 + Mi;j;k;u1;3;4;0 \Delta Mi;j;k1;3;0 + Mi;j;k;u1;2;4;0 \Delta Mi;j;k1;2;0 \Gamma 2 \Delta Mi;j;k;u1;2;3;0 \Delta Mi;j;k1;2;3 ? 0 (4.3-7) 54 General position revisited. The general position assumption used in this paper assures that no geometric test is ambiguous. We summarize and revise the necessary assumptions below and include pointers to the formulas for which the assumptions are relevant. ffl No 4 points lie on a common plane; compare with (4.3-1). ffl No 5 points lie on a common sphere; compare with (4.3-2). ffl No smallest circumsphere of 2, 3, or 4 points has a radius equal to any given ff; compare with (4.3-3), (4.3-4), and (4.3-5). ffl No point lies on the smallest circumsphere of 2 or 3 other points; compare with (4.3-6) and (4.3-7). These assumptions are indeed very restrictive and rarely true for real-life data. We will deal with this apparent shortcoming in section 5.1. 55 Chapter 5 Simulation of Simplicity Simulation of Simplicity (SoS) is a paradigm that allows implementations of geometric algorithms to work on objects in arbitrary position as if they were in simple position. A set of objects is called to be in simple or nondegenerate position if it does not contain any degeneracies; otherwise, it is said to be in arbitrary position. The terms "simplicity" and "degeneracy" can only be defined with respect to the primitives used by the algorithm in question. For example, to compute the convex hull of a two-dimensional point set, an algorithm might use the so-called ccw test as a primitive operation: applied to an ordered triplet of points, ccw determines whether or not they define a left- or rightturn. A degeneracy occurs when points are collinear. Topological primitives. A ternary topological primitive is an operation that classifies some given configuration of input objects as one out of three possible cases, represented by +1, 0, and \Gamma 1. Here, 0 indicates a degeneracy, and +1 and \Gamma 1 distinguish between two complementary nondegenerate cases. Topological operations with more cases can be obtained by combining several ternary tests. In theory, the scope of SoS is restricted to ternary topological primitives whose result can be computed by evaluating the sign of some polynomial expression in the parameters representing the objects involved. For practical purposes (read, efficiency), we further restrict the scope to primitives expressible in (polynomial expressions of) determinants. As it turns out, this covers many, if not most, of the primitive operations that occur in computational geometry. Although, this is sufficient for the algorithms discussed in this thesis, it certainly does not include all possible geometric operations, especially those that are not even topological; for example, the computation of distances or intersections. Nevertheless, SoS remains a powerful tool since a large variety of geometric algorithms are either already using only primitives expressible in determinants, or can be rewritten to do so. 56 Intuition. The basic idea of SoS is to perturb the given objects ever so slightly such that all degeneracies disappear. However, it is important that the perturbation is small enough in order not to change the nondegenerate position of objects relative to one another. Since it is usually rather difficult and costly to actually come up with such a perturbation, SoS performs it symbolically by replacing the parameters specifying each object by polynomials in ". These polynomials will be chosen in such a way that it is sufficient to assume " positive but arbitrarily small to guarantee nondegeneracy, while enabling the usage of " as an indeterminant for a symbolic evaluation of primitive operations on the perturbed set. This means that the perturbation is never actually performed. In other words, SoS can be viewed as a careful way of breaking ties: a black box that consistently assigns a simple result (namely, +1 or \Gamma 1) to a primitive operation, even if it really evaluates to 0, the degenerate case.1 Section 5.1 outlines the perturbation method for point sets in IRd, the points given in Cartesian coordinates. This restriction can easily be lifted, as mentioned in the original paper [27]. 1The "Systematic treatment Of Special cases" is yet another possible extension of the acronym SoS, and so is "Save Our Software." It is amusing to quote from dictionaries. For example, A S Hornby's Oxford Student's Dictionary of American English: SOS (letters of international code signal of distress) n [C] (fig) urgent appeal for help. Or the Dictionary of American Slang by H Wentworth and S B Flexner: SOS S.O.S. n. 1 Lit. and fig., a distress signal; fig., a request for aid, as in finishing a task successfully or in extracting oneself from an embarrassing situation. From the Morse Code telegraph and radio distress signal. Usually thought to be the initials from "save our souls" or "save our ship," the letters actually represent no words, but were chosen because they are easy to transmit and understand without confusion. 2 [taboo] The usual routine tasks, procedures, or glib task; the usual unfairness, deception, or bull; the usual shit. Some W.W.II Army use. From the initials for "same old shit," reinforced jocularly by the Morse Code distress signal. 57 5.1 Symbolic Perturbation Let S be a finite set of points in IRd, with a unique index i assigned to each point pi, which itself is given by d coordinates: pi = (ssi;1; ssi;2; : : : ; ssi;d) The parameter matrix P is the matrix formed by all coordinates ssi;j, assume that the indices are assigned contiguously from 1 through n. It is not necessary, though convenient, that the indices be contiguous. In fact, not even the size of S has to be fixed. The only prerequisite is that each point pi gets its unique serial number i. Conceptually, the symbolic perturbation replaces each coordinate ssi;j, for 1 ^ i ^ n and 1 ^ j ^ d, by a polynomial in ": ssi;j(") = ssi;j + "(i; j) (5.1-1) Call ssi;j(") the perturbed or "-expanded coordinate. The "-expanded parameter matrix, the set of perturbed points, and the perturbed points themselves, will be denoted as P("), S("), and pi("), respectively. The polynomial "(i; j) is called the "-expansion for the original coordinate ssi;j. The product of zero or more "-expansions is called an "-product; if the number of factors is zero, the "-product is defined to be equal to 1. Requirements and axioms. The choice of the polynomials "(i; j) will be guided by the three requirements SoS has to meet, for sufficiently small " ? 0: (R1) S(") must be simple. (R2) S(") must retain all nondegenerate properties of the original set S. (R3) The computational overhead caused by simulating S(") must be feasible. In order to meet these requirements, there is no need to further specify the "(i; j) expansion, if the "-expansions satisfy the following axioms: (E1) "(i; j) ? 0 but arbitrarily small. (E2) "(i; j) AE "(u; v) if and only if (i; j) OE (u; v). 58 Here, "OE" denotes some linear order on the index pairs, and "AE" means "arbitrarily much larger" or "orders of magnitude larger." In other words, coordinates with smaller index pairs are perturbed much more than those with larger index pairs. A more rigorous way to express the fact that "- expansions differ by "orders of magnitude" for different index pairs is to state that, the following condition always holds for fixed (u; v): Y (i;j)OE(u;v) "(i; j) ? "(u; v) (5.1-2) Correctness. In order to show that this axiomatic definition of "(i; j) is sufficient to meet the requirements of SoS, some more definitions and notations are needed. First, define a polynomial expression by a finite number of applications of the following two inductive rules: (i) All original coordinates ssi;j are (elementary) polynomial expressions. (ii) If A and B are polynomial expressions and c is some constant, then c \Delta A, A + B, and AB form polynomial expressions. In a similar fashion, an "-expanded polynomial expression, or "-polynomial for short, can be defined by a finite number of applications of the inductive rule (ii), but with "-expanded coordinates ssi;j +"(i; j) as elementary "-poynomials. Note that, for our purposes, it makes sense to consider only polynomial expressions not identical to zero. Let I(e) denote the set of index pairs associated with an "-product e, where I(1) is the empty set, and "\Gamma " denotes the usual set difference. For I(e1) 6= I(e2) call I(e1) smaller than I(e2) if (i) I(e1) \Gamma I(e2) is the empty set, or (ii) (i; j) OE (u; v), for (i; j) the largest index pair in I(e1) \Gamma I(e2), and (u; v) the largest index pair in I(e2) \Gamma I(e1). With this, the following lemma holds. Lemma 5.1.1. Let c1 and c2 be two positive constants, and let e1 6= e2 be two "-products. Then, for sufficiently small " ? 0, c1 \Delta e1 ? c2 \Delta e2 if and only if I(e1) is smaller than I(e2). Lemma 5.1.1 is an immediate consequence of condition (5.1-2) and axiom (E1) since the "(i; j) factors in e1 and e2 can be assumed to be small enough to compensate for the influence of any constants c1 59 and c2. Lemma 5.1.1 is crucial for showing that the perturbed point set S(") is indeed simple and consistent with the original set S. Lemma 5.1.2. S(") is nondegenerate with respect to any topological primitive that can be expressed by evaluating the sign of some polynomial expression. Moreover, S(") retains all nondegenerate properties of S. Let C denote the polynomial expression associated with the primitive at hand, and let C(") be its "-expansion. It can be assumed that C is not identical to zero.2 The "-polynomial C(") can be rewritten as C(") = mX `=0 C` \Delta e` where the coefficients C` are polynomial expressions and the e` are "-products with I(e`) smaller than I(ek) for ` ! k. Axioms (E1) and (E2) imply that the absolute value of the first term with nonzero coefficient is larger than the sum of all other terms. Thus, the sign of the first nonzero coefficient is the sign of the entire "-polynomial. Note that e0 = 1 and C0 = C, the polynomial expression of the unperturbed primitive. In other words, S(") retains all nondegenerate properties of S. This also means, that there is no computational overhead for a priori simple cases. All coefficients C` are now polynomial expressions in the coordinates ssi;j. In order to proof the simplicity of P ("), it remains to be shown that there is always at least one term with nonzero coefficient in C("). This follows by induction over the length of "-polynomials, since Lemma 5.1.1 assures that terms with different "-products can never cancel each other. Notice, that the coefficients C` of C(") can be interpreted as "subexpressions" of the original expression C. Compare this with the partial derivatives of Yap's method in [82] and [81], or with the subdeterminants in the original SoS paper [27]. So, the exact definition of the "-expansion is irrelevant, and, indeed, not used in the implementation of SoS: it is sufficient to presume Axioms (E1) and (E2). However, the linear order on the index 2Obviously, primitives with a corresponding polynomial expression identical to zero can hardly be useful. Note that this usually implies that the indices of the involved objects are pairwise different; compare this with the primitive operations of section 5.2. 60 pairs has to be specified. For historic reasons it is defied similar, but not equal, to the lexicographical order: (i; j) OE (u; v) if ( i ! u, ori = u and j ? v. (5.1-3) It is still important to show that there is at least one suitable "-expansion. The existence of such an expansion implies the physical existence of an appropriately perturbed point set, which is the guarantee for the consistency of the method. Lemma 5.1.3. There exists an "-expansion "(i; j) which satisfies Axioms (E1) and (E2), assuming a linear ordering on index pairs as in (5.1-3). Let ffi be an integer larger than the maximum number of "-expansions with identical index pairs within any of the "-products that occur during the simulation. For example, ffi = 2 if the used "-polynomials contain only "-products with pairwise different "(i; j) factors, or set ffi = 3 when the "-products contain squares of "(i; j) but no higher power. With this, choose "(i; j) = "ffi d\Delta i\Gamma j for all (i; j) and d, the number of parameters per object. Obviously, Axiom (E1) holds if " ? 0 and sufficiently small. The appropriate choice of ffi, with ffi * 2, makes "(i; j) satisfy condition (5.1-2), which implies Axiom (E2). This proves lemma 5.1.3. 61 5.2 Implementing Predicates and Other Primitive Operations Since P (") is simple, an implementation employing SoS can neglect any degeneracies of the input data and assume that topological primitive operations always return a non-zero result. Thus, such primitives can be seen as (boolean) predicates, returning only TRUE or FALSE by mapping their "degenerate result" consistently onto the two complementary non-degenerate cases. However, the objects are now internally represented by symbolic expressions, and even for the simplest tasks (for example, comparing two coordinates) SoS has to provide special procedures to handle the "- expansion. The aim of this section is to illustrate how SoS predicates can be implemented, especially the ones needed for the construction of three-dimensional Delaunay triangulations and ff-shapes, as they are listed in section 4.3. We assume axioms (E1) and (E2) for the "-expansion, a globally available parameter matrix P, and the index pairs linearly ordered as in (5.1-3); see section 5.1. Furthermore, we extend the notation of minors, as introduced in section 4.3, to subdeterminants of the "-expanded parameter matrix. Assume a column 0 which is not perturbed, namely, ssi;0(") = ssi;0 = 1, and define: Mi1;i2;:::;ikj1;j2;:::;jk (") = det 0BBB BB@ ssi1;j1(") ssi1;j2 (") : : : ssi1;jk (") ssi2;j1(") ssi2;j2 (") : : : ssi2;jk (") ... ... . .. ... ssik;j1 (") ssik;j2(") : : : ssik;jk (") 1CCC CCA : Sign of determinants. The original paper on SoS [27] presents a whole list of predicates involving points in IRd as well as planes and hyperplanes, all expressed in terms of the sign of determinants. In fact, the original version of SoS considered only two types of determinants: det \Lambda d(") = Mi1;i2;:::;id1;:::;d\Gamma 1;0(") and det \Delta d(") = Mi1;i2;:::;id1;:::;d\Gamma 1;d("). The main reason for restricting the scope to determinants is pragmatic. On the one hand, there are many algorithms in computational geometry that are based solely on primitive operations expressed by such determinants. On the other hand, determinants can be evaluated more efficiently than general polynomials of the same complexity, and there is a particularly clean way to generate the coefficients (that is, the subdeterminants) of det \Lambda d(") and det \Delta d(") in the right order; see the next.v procedure in [27]. It should be noted that an extension to arbitrary polynomials is possible, in principle. Software packages for symbolic mathematics can be used to generate the coefficients of arbitrary "-polynomials, followed by some semi-automatic procedure to transform the sequence of coefficients into program code. Alternatively, refer to Rosenberger [71] for an algorithm that generates the coefficients of arbitrary "-polynomials directly (although he uses a slightly different for of "-perturbation). 62 5.2.1 Primitives for Delaunay Triangulations As stated in section 4.3, the incremental flip algorithm requires two geometric primitives: the orientation test and the sphere test. Their implementation is discussed below. Orientation test. Providing a predicate testing a sequence of points for their positive or negative orientation is essential for constructing as well as storing a triangulation. The predicate will also serve as an illustration of how SoS implements topological primitive operations. However, to keep things simple, we start with the orientation test in IR2; an extension to IR3 will follow. Predicate ccw(i; j; k) takes three pairwise different indices, i, j, and k, and returns TRUE if and only if the sequence of the three points (pi("); pj("); pk(")) in IR2 has positive orientation. In other words, it returns TRUE if and only if pi("), pj("), and pk(") make a counterclockwise turn. Without loss of generality, assume i ! j ! k. With this, ccw(i; j; k) can be implemented by evaluating the sign of Mi;j;k1;2;0("). The corresponding "-polynomial looks as follows. Mi;j;k1;2;0(") = det 0BB@ ssi;1 ssi;2 1 ssj;1 ssj;2 1 ssk;1 ssk;2 1 1CCA \Gamma det ssj;1 1ss k;1 1 ! \Delta "(i; 2) + det ssj;2 1ss k;2 1 ! \Delta "(i; 1) + det ssi;1 1ss k;1 1 ! \Delta "(j; 2) + 1 \Delta "(j; 2)"(i; 1) \Sigma \Delta \Delta \Delta Note that terms not explicitly listed are irrelevant because the fifth coefficient is the constant 1 and therefore always different from 0. The sequence of the five relevant terms of Mi;j;k1;2;0("), as well as of any higher-dimensional "-determinant, can be generated using the next.v algorithm described in [27]. This is a procedure that returns the sorted sequence of subdeterminants, while eliminating irrelevant ones, that is, terms with repeated coefficients and terms smaller than the largest term with guaranteed nonzero coefficient. The sequence of relevant terms can then be used to generate the code 63 for computing the sign of the "-polynomial in a more or less straightforward manner. Table 5.2 shows the resulting code for signMi;j;k1;2;0("), which then can be used to implement the ccw predicate as in table 5.1. Extending the orientation test to IRd, d ? 2, is straightforward. Here, we only consider the case d = 3, because this is what is needed for the construction of three-dimensional Delaunay triangulations. Predicate positive(u; i; j; k) assumes four pairwise different indices u, i, j, and k and returns TRUE if and only if the sequence of points (pu("); pi("); pj("); pk(")) in IR3 has positive orientation. In other words, it returns TRUE if and only if pu(") sees (pi("); pj("); pk(")) making a counterclockwise turn. function ccw(i; j; k) is sort indices so that i ! j ! k, using s swaps; r "signMi;j;k1;2;0(")"(i; j; k) if s is even then return r ? 0 else return r ! 0 end Table 5.1: The orientation test in IR2. function "signMi;j;k1;2;0(")"(i; j; k): \Sigma 1 is assert i ! j ! k; r signMi;j;k1;2;0; if r 6= 0 then return r; r \Gamma signMj;k1;0; if r 6= 0 then return r; r signMj;k2;0; if r 6= 0 then return r; r signMi;k1;0; if r 6= 0 then return r; return +1 end Table 5.2: Computing signMi;j;k1;2;0("): 64 The pseudocode for positive can be found in table 5.3. It implements the formula (4.3-1) in section 4.3, and is a straightforward extension of the code for ccw, provided there is an implementation for Mu;i;j;k1;2;3;0("). This list of the coefficients of this "-polynomial can either be computed using the above mentioned next.v algorithm or it can be found in form of a table in [27]. There are plenty of algorithms for point set problems which are based on computing the orientation of a sequence of points. Prime examples are algorithms for the construction of convex hulls; see for example [66, 72, 73] or [21, 54, 67]. References to other applications can be found in [27], which also extends the concept of orientation to points with homogeneous coordinates in arbitrary dimensions. Sphere test. The crucial primitive operation for constructing Delaunay triangulations is to check if a triangle is locally Delaunay or not. This is equivalent to testing whether a given point lies inside or outside the circumsphere sphere of four other points; see section 4.1. The two-dimensional version of the sphere test is a common operation in computational geometry (see for example [42, 50]). For IR3, we state the predicate as follows. Predicate in sphere(u; i; j; k; v) assumes five pairwise different indices u, i, j, k, and v and returns TRUE if and only if the perturbed image of the fifth point, pv, lies inside the sphere through the perturbed images of the first four points, pu, pi, pj , and pk. Table 5.4 implements the in sphere predicate. As mentioned in section 4.3, this assumes that the point set S is lifted to SU 2 IR4, by attaching a fourth coordinate ssi;4 = P3j=1 ss2i;j. The implementation uses two orientation tests, one in IR3 and one in IR4, and deduces the result according to formula (4.3-2). function positive(u; i; j; k) is sort indices so that u ! i ! j ! k, using s swaps; r "signMu;i;j;k1;2;3;0(")"(u; i; j; k) if s is even then return r ? 0 else return r ! 0 end Table 5.3: The orientation test in IR3. 65 Remark. In our implementation, SoS perturbs the sum of squares, P3j=1 ss2i;j, as if it was a forth coordinate, ssi;4. This "cheap trick" (Emiris [30]) is used for reasons of compactness and efficiency. However, since primitives in IR3 and IR4 are used together, inconsistencies might occur, and some care is in order. Recall that the lifted points pi 2 SU lie on the paraboloid of revolution U : x4 = P3`=1 x2` in IR4. This is no longer true for SU ("). In particular, SU (") is no longer necessarily a convex point set. Unfortunately, this happens to be the underlying assumption for the correctness of the algorithm constructing the Delaunay triangulation of S. Fortunately, there is only one (special) case for which SU (") is not convex, namely, if S includes coinciding points pi = pj, i 6= j; and it is easy to eliminate such points in a preprocessing step, without changing the overall time complexity to the worse. If, on the other hand, pi 6= pj, for all i 6= j, then SU (") is convex, because the unperturbed SU is convex, and the " of the perturbation is assumed to be arbitrarily small. Thus, the "cheap trick" works, and the in sphere predicate can be implemented using SoS as shown in table 5.4. Refer to [30] or [50] for alternative (though, not necessarily more efficient) approaches. 5.2.2 Primitives for Alpha Shapes Using SoS to implement the primitive operations necessary for the construction of alpha shapes turns out to be considerably more difficult than using SoS for Delaunay triangulations. Some of the details are quite tedious. Radius of a smallest circumsphere. For the face classification of section 2.3 we need to compute the radius %T of the smallest circumsphere of a k-simplex oeT 2 D, with 1 ^ k ^ 3; see formulas (4.3-3), (4.3-4), and (4.3-5) in section 4.3. We face the following problems. function in sphere(u; i; j; k; v) is sort indices so that u ! i ! j ! k; r1 "signMu;i;j;k1;2;3;0(")"(u; i; j; k); sort indices so that u ! i ! j ! k ! v, using s swaps; r2 "signMu;i;j;k;v1;2;3;4;0(")"(u; i; j; k; v); if s is even then return r1r2 ? 0 else return r1r2 ! 0 end Table 5.4: The sphere test in IR3. 66 First, computing the radius of circumspheres is not a predicate operation. What we really need is to compare the radii of two simplices, possibly of different dimension. We approach this by providing primitives computing the radius of a simplex, and implementing a separate function that compares two arbitrary radii. Second, under SoS, the radius of a circumsphere is a fraction of two "-polynomials. Fortunately, there is no need to keep the perturbation when representing ff-shapes, once the Delaunay triangulation is built and all topological decisions are made; compare this with sections 3.2 and 4.2. A simple mapping of fractions of "-polynomials to normal rational numbers is sufficient. For an "-polynomial A(") let cA and "A be its first non-zero coefficient and the corresponding "-product. Then, map A(")B(") to cA cB for "A = "B, 0 for I(eB) is smaller than I(eA), and 1 otherwise. Note that, although the ff-spectrum of section 3.2 stores only the floating-point representation of the radii %T , we temporarily need to maintain them as ratios of long-integer numbers to keep the spectrum consistent. After sorting, however, the exact rationals can be mapped into floating-point numbers with the same total ordering. Third, the primitive operations are now based on polynomial expressions of determinants. This means that the original framework for generating the "-coefficients in [27], which was based on single determinants, needs to be extended accordingly. For example, generate the "-polynomials of single determinants using the next.v algorithm, and implement the elementary arithmetic operations on "-polynomials, that is, A(")+B("), A(")B("), and c\Delta A("), symbolically, while maintaining the proper ordering of the "-terms. Finally, notice that the determinants in formulas (4.3-3), (4.3-4), and (4.3-5) no longer fit in the model of det \Lambda d(") = Mi1;i2;:::;id1;:::;d\Gamma 1;0(") and det \Delta d(") = Mi1;i2;:::;id1;:::;d\Gamma 1;d("). The new primitives require "-expanded minors with arbitrary column indices, although pairwise different and sorted (except a possible zero index for the last column). However, since only the relative order of the involved indices is of any relevance for the comparison of "-products, an appropriate extension to the next.v algorithm of [27] is trivial. For completeness, let us now state the definition of the primitive operation computing the radii. The listings of the corresponding pseudo codes are omitted, for they provide little or no additional insight. 67 Primitives rho1(i; j), rho2(i; j; k), and rho3(i; j; k; u) return the ratio %2 = ab identical to the square of the radius of the smallest circumsphere of a given 1-, 2-, or 3-simplex. An edge oeT is given by T = fi; jg, a triangle oeT is given by T = fi; j; kg, and a tetrahedron oeT is given by T = fi; j; k; ug. In all cases, the indices are assumed to be pairwise different. The resulting ratio is either a normal rational ab , with a 6= 0 and b ? 0, 0, or 1.e Attached and unattached edges and triangles. Two more operations remain to be implemented. To test whether an edge or a triangle is attached or unattached, we need predicates that decide whether a given point is inside the smallest cicumsphere @bT of a given edge or triangle oeT . Formulas (4.3-6) and (4.3-7) express this in terms of minors. Predicate in ball3(i; j; k; u) assumes pairwise different indices i, j, k, and u and returns TRUE if and only if pu(") 2 bT with T = fpi("); pj("); pk(")g. Predicate in ball2(i; j; k) assumes pairwise different indices i, j, and k, and returns TRUE if and only if pk(") 2 bT with T = fpi("); pj(")g. SoS versions of both predicates can be implemented in the same way as all the operations above: generate the sequence of terms for the "-expanded minors, and compute the sequence for the entire expression using polynomial arithmetic. However, in the case of the in ball predicates there is one bothersome detail. It is now no longer possible to permute the indices without totally changing the meaning of the predicate. In other words, the indices can no longer be sorted; and this was the crucial property for the approach to generate the terms of the "-polynomials without knowing the actual values of the involved indices. One solution would be to sort the terms on the fly, at run time, using the actual index values, but this is too inefficient. The other extreme would be generate the sequence of points for each possible perturbation of the indices. Fortunately, ther is a feasible alternative. Observe that the indices of the points defining bT actually can be sorted without changing the meaning of the predicate. For the point not in T , we need to distinguish all possible cases. For example, assuming i ! j, there are 3 cases for in ball2(i; j; k): i ! j ! k, i ! k ! j, and k ! i ! j. Analogously, we have 4 different cases for in ball3(i; j; k; u). For each case, we generate (and provide code for) the corresponding sequence of relevant terms. 68 5.2.3 Evaluating Epsilon Polynomials Some remarks on low-level details of the SoS implementation are in order. How do we generate the sequence of significant terms? How is its evaluation implemented? Generating significant terms. The current approach to generate the sequence of relevant terms for a particular "-polynomial can be summarized as follows. 1. Use the next.v algorithm of [27] to generate the significant terms of "-expanded determinants, that is, minors of the parameter matrix. 2. If the primitive operation is a polynomial expression of such determinants, use polynomial arithmetic to generate the significant terms for the entire expression, while maintaining the proper ordering. Use commonly available software for symbolic mathematics (for example: Mathematica) to simplify the expressions for each coefficient. 3. Translate the sequence of relevant terms into C code that evaluates the "-polynomial term by term until a nonzero coefficient is reached. This process can be automated. 4. Bind the resulting code for the new predicate into the SoS Library and use it as a black box in the application code. This process resembles a compilation of symbolic "-expressions. Another approach would be to build an "interpreter" for "-polynomials by providing procedures for an on-line generation of the significant terms of "-expanded minors, polynomial expressions thereof, or, in the most general case, arbitrary "-polynomials. Executing such a procedure would certainly cost more time then just scanning compiled "-polynomials. However, this idea still seems to be worthwhile, since the compiled code for "-expressions can get fairly large,3 and the vast majority of the code is practically never executed. This is because degenerate cases can be expected to occur with very low probability (see for example [71]), and the probability that a term of depth t needs to be evaluated drops dramatically with increasing t (see section 7.1). Exact Arithmetic. Obviously, SoS requires the coefficients of "-polynomials to be evaluated exactly; this disables the use of hardware supported and fast floating-point arithmetic. The current implementation uses long-integer arithmetic, with a cost of O(`) elementary operations for adding and subtracting two numbers with ` digits each; multiplication requires O(`2) operations. An alternative would be modular arithmetic based on the Chinese remainder theorem. This would speed up multiplication to only O(`) operations. However, there are at least two reasons why modular arithmetic is not the method of choice for our purposes. First, when evaluating determinants, mod3For example: the "-polynomial in the numerator of formula (4.3-5) has 213 significant terms. The corresponding C code produces an object file of 127k bytes. 69 ular arithmetic does not seem to give an advantage over long-integer arithmetic; see [36, 71]. The situation is even worse for sign computations (that is, comparisons with 0) for they require the expensive modular inverse. Second, in contrast to long-integer arithemtic, modular arithmetic requires always maximum length for object representation. This is bad for the ff-shape application: recall that we have to (temporarily) store the radii of all faces in the Delaunay triangulation. Therefore, we would need an enormous amount of memory. For further details on long-integer as well as modular arithmetic see [49]. SoS employs a "home-made" implementation for long-integer arithmetic: the Lia library,4 which compares favorably with packages that are commonly available.5 Lia consists of carefully tuned C code, providing functions for long-integer addition, subtraction, and multiplication. It uses linear arrays a[0::`] of 30-bit integer digits to store the long-integer objects in the usual radix representation: a = P`i=1 a[i]ri\Gamma 1. All digit-by-digit operations are done with normal 32-bit integer arithmetic; the multiplication routine temporarily splits the Lia digits into half-digits to avoid overflow. Note that using a power of two as radix (here: r = 230) enables shift-operations to accelerate digit normalization. A recent paper by Fortune and Van Wyk [36] presents several ideas which would speed up Lia even further; for example: delayed normalization, or doing the digit-by-digit operations with hardware supported floating-point arithmetic. It should be stated that some sort hardware support for exact arithmetic would be highly desirable: either for the special case of determinants, or, even better, for long-integer arithmetic in general. Hardware "long-integer accelerators" would not only speed up SoS considerably, they would also be very useful in the field of symbolic computation, which employs exact arithmetic heavily. For a first step in this direction, see Rosenberger's thesis [71], in which he presents the design of a chip computing the sign of 3-by-3 long-integer determinants. Also, refer to Yap [83] for a discussion on the general benefits of exact computation, and how to make it a viable alternative, even without hardware support. Evaluating determinants. We use expansion by minors (also known as cofactor expansion) for the exact evaluation of determinants; however, it is implemented using dynamic programming. Assume, a k-by-k determinant has to be evaluated. Evaluate all \Gamma k2\Delta 2-by-2 minors of the first two rows, then all \Gamma k3\Delta 3-by3 minors of the first three rows, and so on and so forth, until the entire determinant is evaluated. (This is cofactor expansion taking advantage of common subexpressions.) The time complexity of this method is exponential in k, but although there are asymptotically faster algorithms, expansion by minors with dynamic programming is the method of choice for small k; that is, for k ^ 10, according to [71]. 4An earlier version of Lia is described in [62]. 5For example, the mp library which comes with the Unix operating system, or the BigNum package, as distributed by INRIA [75]. 70 Chapter 6 Signatures Although individual ff-shapes are quite interesting geometric objects themselves, it is even more useful to have also efficient algorithms that can compute geometric and topological properties of individual (as well as entire families of) shapes. One of the major advantages of the sound concept of ff-shapes is that such efficient computations are indeed possible. This short chapter presents first steps in this direction, gives relevant references, and discusses how geometric and topological properties can be used to find the "right shape." 71 Features. We refer to a geometric or topological property of a shape as a feature of the shape. The computation of elementary metric features, like the volume or the surface area of a given shape Sff, are easy. For example, the volume of Sff is simply the sum of the volumes of all the tetrahedra in the corresponding ff-complex Cff. Computing elementary combinatorial features, such as the number of surface triangles or the number of singular points, is similarly straightforward. The computation of topological properties, such as the connectivity of the ff-shape as expressed by its first three homology groups (see for example [38]) is more difficult (but also more useful), especially when studying the change in a feature for varying ff. More challenging and very useful, is the computation of metric properties connected to the dual ff-diagrams (see section 2.1). Such ball diagrams, also known as space filling diagrams, are used to model molecules in computational biology, and features of interest include the volume or surface area of the union of balls, or the volume of the holes or "voids" inside the structure. All these, and more, can conveniently be computed within the framework of (weighted) ff-shapes. For more details, see [23-25]. Signatures. The history of a certain feature over all ff, 0 ^ ff ^ 1, is called a signature. Consider some specific feature, such as, the number of connected components in the ff-shape. The corresponding signature is an integer function c: [0; 1] ! IR, so that c(ff) is the number of components of Sff. This function represents the number of components as ff changes continuously from 0 to 1. Given the ff-spectrum and the master list, as discussed in section 3.2, it is fairly easy to compute c. Start at ff = 0, and maintain a union-find data structure (see for example [13]) storing the connected components of Sff as the threshold values are processed in increasing order. The number of components corresponds to the 0-homology group of Sff. Signatures for higher-order homology groups (0, 1, and 2, for simplicial complexes in three dimensions) are also of interest, and a more challenging task to compute. According to common terminology in algebraic topology, their corresponding signature values are called the Betti numbers. A fair amount of algebraic terminology is necessary to properly define the Betti numbers of a simplicial complex. In three dimensions, the following intuitive interpretation of the Betti numbers b0, b1, b2 apply: b0 is the number of connected components of the complex or shape; b1 is the number of independent tunnels, where "independent" roughly means that one tunnel cannot be obtained as a sum of other tunnels (addition respects orientation and oppositely oriented simplices erase each other); finally, b2 is the number of voids of the complex or shape.1 Further details are omitted; refer to Delfinado and Edelsbrunner [15] who give an efficient algorithm computing b0, b1, and b2 for simplicial cell complexes. An implementation for ff-shapes, based on the master list, is available. 1Here, the term "void" is used for a bounded connected component of the complement of the shape. 72 The author believes that the possibility of feature and signature computations is the major advantage of the ff-shapes over less rigorously defined concepts such as, for example, iso-surfaces. Signature graphs are particularly handy for the selection of an appropriate ff-value, which typically depends on the application and the user's momentary focus of attention. Moreover, such graphs provide us with a new form of quantification for the external and internal structure of point sets with possible applications to shape and feature recognition. Computing signatures. In the following, we discuss how the ff-spectrum and master list (see figure 3.5 in section 3.2) can be used to efficiently compute elementary metric signatures. Here, "efficiently" means, in constant time for each face in the Delaunay triangulation. Take, for example, the volume, V : [0; 1] ! IR, and assume that V (ff), the volume of Sff, is to be stored in an linear s 0; for r 1 to r do let e be the first entry of the sublist for spectrum[r]; while e exists do if face type of e = ``Triangle" then let fijk denote the triangle corresponding to entry e; if rank type of e = ``%" then f fijk does not belong to Sffr\Gamma 1 and is singular in Sffr g s s + 2 \Delta area(fijk); else if rank type of e = ``_" then if fijk is attached then f the attached fijk does not belong to Sffr\Gamma 1 and is regular in Sffr g s s + area(fijk); else f fijk is singular in Sffr\Gamma 1 and regular in Sffr g s s \Gamma area(fijk); end if else if rank type of e = ``_" then f fijk is regular in Sffr\Gamma 1 and interior in Sffr g s s \Gamma area(fijk); end if end if; let e be the next entry of the sublist for rank r; end while; A[r] s; end for; assert A[r] = A(ffr), for 1 ^ r ^ r; Table 6.1: Computing a signature: A(ff), the area of the triangles bounding Sff. 73 array V [1::r] with one entry for each ff-rank. The basic procedure is: initially, set V [0] = 0, then, scan through the ff-spectrum and master list in direction of increasing rank r, and add the volume of each tetrahedron to V [r], as it becomes part of Sffr . Table 6.1 gives the details of the analogous procedure of computing A: [0; 1] ! IR, where A(ff) is the surface area of Sff. Assuming a linear array A[1 : : : r], allocated to store the signature values, we scan through the entries of the master list, ignoring all faces but the triangles. Note that we have to distinguish between attached and not attached triangles, and when a triangle becomes singular, regular, or interior (according to table 2.2 in section 2.3). We set A(ff) = Arg(ff) + 2Asg(ff); where Arg(ff) and Asg(ff) denote the portions of the surface area contributed by the regular and singular triangles. Examples. This chapter closes with figures 6.1 through 6.6, showing selected signature graphs and shapes of some of the data sets used in the empirical experiments of chapter 7; the legends to these figures are in table 6.2. We show the graphs of the elementary metric signatures functions: volume, surface area (of the total, regular, and singular boundary), and (singular) edge length. Often, weighted combinations of these turn out to be useful as well. For example, we can consider w(ff) = maxf0; 4A(ff) 1 2 \Gamma 8V (ff) 1 3 \Gamma 2L(ff)g; where L(ff) denotes the total length of all singular edges in Sff. Such functions are mostly heuristic, and depend on the application at hand. Here, our aim was to take demo data for modeling and rendering programs, strip it from all connectivity information (that is, remove all edges and triangles), and re-compute the model shapes using the ff-shape software on the vertex coordinates. Maximizing the above signature function seems to be a promising heuristic to find "good shapes." In general, point sets with "similar" features share "similar" charateristics in their signature functions, but this subject is still wide open for future research. 74 Figure 6.1: "tori" -- The points were randomly generated on the surface of two linked tori. Four shapes are shown: the point set itself, and three other shapes maximizing one signature function each: the singular edge length, the surface of the singular boundary, and the heuristic signature w. The last shape seems to be the best. It neatly separates the two tori. Note that, in this, and the following figures, singular triangles are displayed in a darker color than the regular triangles. Figure 6.2: "spiral" -- This data set consists of points generated on a spiral, and on a separate inner ring. The four displayed shapes maximize the singular edge length, singular boundary, total boundary, and regular boundary area, respectively. Note that, although, we are able to achieve ff-shapes with similar and recognizable features as the original model, we cannot totally separate the spiral from the ring. Also observe how maximizing the singular edge length reveals how the original model was generated. Figure 6.3: "liss-5-8" -- This data set consists of points generated on the surface of a self-intersecting "worm" around a sphere. As in the "spiral" data, the original model consists of "moving rings." This is seen best in the shape maximizing the singular edge length. The three other shapes maximize the singular boundary, the heuristic map w, and the regular boundary. Note the steep drop in volume, and the corresponding hike in the surface area, when we decrease ff around ff56000. This is exactly when we hit the radius of the sphere enclosed by the "worm" and a void appears inside the shape. Figure 6.4 illustrates this by displaying two shapes which look quite the same from the outside, but have a huge difference in volume. The difference of the two, reveals the inner void.2 Figure 6.5: "teapot" -- This point set consists of the vertices of a model of the famous "Utha teapot." The first two shapes are based on the singular edge length: one at the global maximum, one at a local minimum. The third shape maximizes w, and seems to capture the "teapot" best. The fourth shape lies in the ff-spectrum close to the rank where the volume and area signatures show a dramatic change and spike around ff13500. The difference shape in figure 6.6 reveals a rather bizarre void growing inside the "teapot" as we decrease ff. Table 6.2: Legend to following figures displaying ff-shapes and signatures. (Refer to chapter 7 for more information on the data sets of figures 6.1 through 6.6.) 75 S Sff727 Sff1866 Sff3389 727 1866 3389 7811 Length Rank Singular edges 727 1866 3389 7811 Area Rank Boundary triangles Regular triangles Singular triangles 727 1866 3389 7811 Volume Rank 727 1866 3389 7811 w Rank w Figure 6.1: Shapes and signatures of the "tori" data set. 76 Sff232 Sff1987 Sff2785 Sff5853 232 1987 2785 5853 7939 Length Singular edges 232 1987 2785 5853 7939 Area Boundary Regular Singular 232 1987 2785 5853 7939 Volume 232 1987 2785 5853 7939 w w Figure 6.2: Shapes and signatures of the "spiral" data set. 77 Sff4222 Sff30710 Sff33859 Sff45858 4222 30710 33859 45858 56000 Length Rank Singular edges 4222 30710 33859 45858 56000 Area Rank Boundary triangles Regular triangles Singular triangles 4222 30710 33859 45858 56000 Volume Rank 4222 30710 33859 45858 56000 w Rank w Figure 6.3: Shapes and signatures of the "liss-5-8" data set. 78 Sff56003 Sff55860 Sff56003 \Gamma Sff55860 Figure 6.4: The difference of two "liss-5-8" shapes. 79 Sff4221 Sff6195 Sff8652 Sff13500 4221 6195 8652 13500 15434 Length Rank Singular edges 4221 6195 8652 13500 15434 Area Rank Boundary triangles Regular triangles Singular triangles 4221 6195 8652 13500 15434 Volume Rank 4221 6195 8652 13500 15434 w Rank w Figure 6.5: Shapes and signatures of the "teapot" data set. 80 Sff137003 Sff13500 Sff13700 \Gamma Sff13500 Sff13700 \Gamma Sff13500 (bottom view) Figure 6.6: The difference of two "teapot" shapes. 81 Chapter 7 Performance This chapter discusses empirical results based on the actual implementations of the algorithms and methods of this thesis. Hardware. The experiments were run on Silicon Graphics Indigo and Crimson workstations, featuring 50 Mhz R4000 RISC processors and 48 Mb and 192 Mb of main memory, respectively. The machines were chosen such that their main memory was always large enough for the data and page swapping is not a factor. Software. It should be noted that all software written for this thesis is publically available, including the data sets used in the experiments.1 1The Silicon Graphics binaries can be downloaded from ftp.ncsa.uiuc.edu:/SGI/Alpha-shape/alpha-1.0 or, alternatively, from ftp.cs.uiuc.edu:/pub/edels/geometry. The latter ftp server also contains the C sources of the incremental flip code for three-dimensional Delaunay triangulations. 82 7.1 Randomized Incremental Flip with SoS This section takes a look at the performance of detri 1.2. This program constructs the Delaunay triangulation of a given point set in IR3, employing the randomized incremental flip algorithm with its "walk through" point location scheme as discussed in sections 4.1.1 and 4.1.2. In the following, we will address the following questions about the program. What is the actual performance of the method? How efficient is the "walk through" point location scheme? How large is the speedup due to randomization? How large is the penalty for the exact arithmetic in SoS? The data. Table 7.1 lists all sample data sets used in the experiment. The data consists of points in IR3, given by their x1-, x2-, and x3-coordinates. For this discussion, we restrict ourselves to 32- bits integer coordinates.2 As it turns out, the performance of algorithms that construct Delaunay triangulations depend on the distribution of the points. For this reason, we categorize the data sets with the goal in mind to collect sets with similar distribution characteristics in each category. The categories are as follows: C : : : Chaotic (or strange) attractors. They are often used in computer graphics to model natural phenomena, such as, smoke.3 Chaotic attractors are essentially "very complex" surfaces with the points rendered in a random fashion (see for example [65]). G : : : Point sets from geometric models. Such data sets consist of point coordinates only; all the original connectivity information (that is, edges and triangles) was removed. The geometric models were found in public domain data bases for solid modelers. M : : : Molecular data. These point sets consists of the coordinates of atom centers. The geometric structure of many proteins is available in public domain data bases. O : : : Other, miscellaneous data. For example, "grid1" with its points on a perfect grid, or the "universe" data, obtained from a simulation of the positions of galaxies within our universe.4 R : : : Range data. It consists of random samples from images obtained with a laser range finder (in this case, the images were of composers' busts).5 The data consists of surfaces of the form x3 = f(x1; x2). S : : : Points on one or more surfaces. These include points generated randomly on certain surfaces. V : : : Volumetric data. These are (random) threshold samples of "voxel" data used in volume visualization, where each grid point (x1; x2; x3) has some value f(x1; x2; x3).6 r : : : Random coordinates. 2This restriction was necessary to compare the performance of detri 1.2 with other programs; detri1.2 as such works also for coordinates given in fix-point format. 3The chaotic attractor data used here was obtained from (or generated with code of) Randy Hudson, University of Chicago, and Paul Haeberli, Silicon Graphics. 4The "universe" data is due to Marc Dyksterhouse who took first steps to investigate the micro- and macro-structure of galaxies using ff-shapes [20]. 5Range data source: Fridtjof Stein, USC Institute for Robotics and Intelligent Systems. 6The volumetric data is due to Pat Moran, NCSA. 83 D SoS depth memory CPU S jSj jF0;1j jF2j jF3j mean Mb seconds molecule M 318 34 4000 1984 0.000052 0.391 4.34 (\Sigma 1:59%) tori S 800 193 12197 6003 0.000000 1.781 14.23 (\Sigma 1:27%) grid1 O 1000 272 12192 5961 0.993537 * 1.133 16.19 (\Sigma 1:31%) spiral G 1248 191 19301 9556 0.017412 * 2.781 21.87 (\Sigma 0:92%) molecule2 M 1366 70 17974 8953 0.000000 1.609 23.01 (\Sigma 1:16%) universe O 1717 77 21321 10623 0.000000 1.945 26.55 (\Sigma 1:18%) fract S 1849 147 23925 11890 0.057878 * 2.133 31.59 (\Sigma 1:20%) pdb9pap M 1908 54 25318 12633 0.000000 2.219 30.92 (\Sigma 0:78%) ma2 S 2206 871 46427 22779 0.016603 * 5.312 40.88 (\Sigma 1:11%) body O 2630 248 35184 17469 0.014703 * 3.055 44.93 (\Sigma 5:53%) fract2 S 2704 141 35203 17532 0.047149 * 3.102 46.84 (\Sigma 1:22%) f15 G 3059 227 42941 21358 0.129221 * 3.625 60.20 (\Sigma 8:37%) br M 3762 89 51559 25736 0.000000 4.437 75.67 (\Sigma 9:05%) liss 5 8 G 4200 534 58486 28977 0.000068 4.961 87.39 (\Sigma 8:77%) teapot G 4668 1247 55529 27142 0.017828 * 5.156 91.73 (\Sigma 8:44%) ku2 V 5292 462 74836 37188 0.362767 * 6.312 112.45 (\Sigma 7:40%) phone 1 G 6070 105 80131 40014 0.011200 * 7.031 123.46 (\Sigma 0:58%) mozart 07k R 7000 259 86775 43259 0.381461 * 7.844 128.40 (\Sigma 1:47%) r 075 r 7500 111 99409 49650 0.000000 8.664 142.36 (\Sigma 0:59%) bach 08k R 8000 192 101944 50877 0.388049 * 9.062 154.76 (\Sigma 0:93%) rat T2a V 9600 263 125933 62836 0.043155 * 11.016 168.19 (\Sigma 1:35%) ra C 10000 566 130172 64804 0.000038 11.430 197.13 (\Sigma 1:08%) r 105 r 10500 131 139951 69911 0.000000 12.148 199.52 (\Sigma 1:25%) r 115 r 11500 149 153193 76523 0.000000 13.281 225.04 (\Sigma 1:09%) smoke 12 C 12300 616 157688 78537 0.000000 13.984 249.51 (\Sigma 1:04%) r 125 r 12500 136 166258 83062 0.000000 14.422 236.84 (\Sigma 1:19%) r 135 r 13500 137 179805 89835 0.000000 15.578 261.74 (\Sigma 2:03%) mozart 14k R 14000 310 173918 86805 0.418923 * 15.648 266.79 (\Sigma 0:86%) smoke 14 C 14300 653 183767 91558 0.000000 16.250 291.39 (\Sigma 1:26%) r 145 r 14500 117 193601 96743 0.000000 16.758 280.81 (\Sigma 2:30%) smoke 15 C 15000 672 192898 96114 0.000000 17.023 315.41 (\Sigma 1:27%) r 155 r 15500 148 206806 103330 0.000000 17.891 300.78 (\Sigma 2:20%) bach 16k R 16000 220 203784 101783 0.440613 * 18.070 314.84 (\Sigma 1:34%) smoke 16 C 16300 707 209931 104613 0.000000 18.492 334.69 (\Sigma 1:36%) r 165 r 16500 147 220087 109971 0.000001 19.039 319.32 (\Sigma 1:50%) r 175 r 17500 143 233323 116591 0.000000 20.180 336.36 (\Sigma 1:30%) smoke 18 C 18300 768 235782 117508 0.000000 20.805 391.08 (\Sigma 2:12%) r 185 r 18500 137 247015 123440 0.000000 21.352 362.63 (\Sigma 1:30%) r 195 r 19500 155 260007 129927 0.000000 22.477 379.73 (\Sigma 1:60%) smoke 20 C 20300 817 261141 130163 0.000000 23.008 437.05 (\Sigma 1:79%) r 215 r 21500 161 287217 143529 0.000000 24.805 417.25 (\Sigma 1:49%) smoke 22 C 22300 854 286810 142979 0.000000 25.273 471.38 (\Sigma 1:49%) smoke 24 C 24300 910 312882 155987 0.000000 27.539 530.65 (\Sigma 1:85%) rat T2b V 26304 424 347042 173310 0.024887 * 30.125 485.98 (\Sigma 1:85%) Table 7.1: The sample data, its combinatorial complexity, and the performance of detri 1.2. C : : : Chaotic (or strange) attractors. G : : : Points from geometric models. M : : : Molecular data. O : : : Other, miscellaneous data. R : : : Range data. S : : : Points on one or more surfaces. V : : : Volumetric data. r : : : Random coordinates. * : : : Highly degenerated data. 84 The experiment. The code of detri 1.2 was run several times on each data set (sample size N = 13). The tables and figures in this section show the sample means Y of the corresponding measurements. Their confidence intervals are described in relative terms as a percentage of Y . The confidence level used was 90%. We say a sample mean has a confidence interval of \Sigma H% if the interval [Y \Gamma H%; Y + H%] contains the the real mean value _ with a probability of 90%. Let us now take a closer look at table 7.1. For each data set S, it lists its size, jSj, and the number of points on the boundary of the convex hull, F0;1. The complexity of the Delaunay triangulation is given by the number of triangles, jF2j, and tetrahedra, jF3j. The number of edges can then be determined by Euler's formula: jSj \Gamma jF1j + jF2j \Gamma jF3j = 1. It is interesting to observe how degenerate a data set is, with respect to the detri 1.2 code, of course.7 This is, in some sense, measured by the mean "depth" of SoS evaluations. Recall that SoS transforms the polynomial expression C corresponding to a certain geometric primitive into a sequence of subexpressions C`, with ` = 0; 1; 2; : : :, and C0 = C. If ` = t is the first index for which C` 6= 0, we say t is the "depth" of this particular evaluation; see section 5.1. For each run we compute the mean value of the t values with respect to all evaluation of the primitive. Table 7.1 then takes the sample mean of the N runs, and lists it in the "SoS depth, mean" column; the corresponding confidence intervals are neglible. Data sets with a particular high mean depth are marked as "highly degenerate." Note the exceptional high spike for the "grid1" data; these are 1000 points on a perfect 10\Theta 10\Theta 10 grid. Most degeneracies in our experiment are due to coplanarities; for example, the "ku2" data contains volumetric data that comes in "slices." Data sets with a mean depth of 0:0 are in general position. Table 7.1 also lists the corresponding memory requirements and running times of detri 1.2. Whereas the memory requirement for different runs of a particular data set were more or less the same (\Sigma 0:01%), there was slight variance in the running times (the table lists the corresponding confidence intervals in parentheses). This is due to the randomization of the algorithm, and the differences in the system load. 7Note that degeneracies are always "defined" with respect to the algorithm at hand, and the set of geometric primitives used. 85 CPU cycle distribution. The next question to ask is: "Where do all the cycles go?" Table 7.2 lists the major procedures of the incremental flip code of tables 4.1 and 4.2 (section 4.1.1) and their relative share of the total number of CPU cycles. The percentages were obtained by run-time profiles generated with pixie, a tool that comes with the Unix operating system on Silicon Graphics workstations. It partitions an executable program into "basic blocks" and writes an equivalent program containing additional code fpr counting the execution of each basic block.8 The new program, including the code for the basic block count, is then used to generate the run-time profiles. Although, the counts will change slightly from data set to data set (and from randomized run to the next), the variance is only marginal. Table 7.2 lists the CPU cycle distribution when executing detri 1.2 for the "spiral" data. The quoted percentages are rounded to the next larger integer value and can be considered as typical for all the sample data sets. It is not surprising that the majority of the time is spent in procedure flip, see table 4.2. Within this procedure, the tests in the head of the while loop are most time-consuming. This is where the sphere tests are performed which check incoming triangles for local Delaunayhood. Recall the code of the in sphere predicate in table 5.4 (section 5.2). While it is listed evaluating two "-deteriminants, 8A basic block is a region of the program that can be entered only at the beginning and exited only at the end. For more information on the pixie tool refer to the Unix man pages. Procedure Share of Share of SoS total CPU cycles primitives operations delaunay 98% search 10% 29% of positive collect 1% 3% of positive connect 1% flip 86% get next facet in L 60% 100% of in sphere triangle-to-edge flip 23% 68% of positive edge-to-triangle flip 3% Table 7.2: CPU cycle distribution while executing procedure delaunay inside detri 1.2. Procedure Share of total CPU cyles in sphere 57% positive 28%P 85% Procedure Share of total CPU cyles mul 39% psub 13% padd 4%P 56% Table 7.3: CPU cycle distribution for SoS and long-integer arithmetic inside detri 1.2. 86 one 5-by-5 and one 4-by-4 (see also section 4.3), the actual implementation eliminates the evaluation of the 4-by-4 determinant. It was only a corrective term to check for the orientation of the involved points. However, it is easy to implement the incremental flip algorithm such that the orientation of the points, for which the in sphere predicate is called, is always positive. Another substantial amount of time is spent on performing orientation tests; they are needed for the point location and the triangle flips.9 As seen in table 5.3 (section 5.2) the implementation of the positive predicate consists of evaluating a 4-by-4 determinant. One has to keep in mind that both predicates are implemented with SoS, and that subdeterminants are sometimes evaluated in order to resolve degeneracies. However, the mean "depth" column in table 7.1 shows us that degenerate cases occur with low probability. In fact, the evaluation of k-by-k determinants, with k ! 4, typically account for less than 0:5% of the total CPU time. Overall, SoS is responsible for roughly 85% of the CPU time; see table 7.3. The long-integer arithmetic done within SoS accounts for 56% of the CPU time. Exact arithemtic. We see, that the largest chunk of CPU time is spent for the robust implementation of the geometric primitives; especially, in doing long-integer arithmetic. Plot (a) in figure 7.1 shows the number of executed long-integer operations, distinguishing between addition/subtraction and multiplication. It is curious to note, that between 3 to 5 times as many pure subtractions (psub) occur than pure additions (padd); see figure 7.1 (b). We conclude that long-integer arithmetic is the real overhead of SoS. Although our implementation already uses fairly efficient code, there is certainly room for improvement. Empirical comparisons of SoS code versus fast but nonrobust floating-point implementations of the incremental-flip algorithm suggest a slow-down factor between 2.5 and 4.0, assuming input data with coordinates with no more than 32 bits. The following experiment with the detri 1.2 code was carried out. - Run detri and write the parameters for all the determinant evaluations to a file. - Read this file with a different program evaluating all these determinants using both long-integer as well as floating-point arithemtic. - Time the difference. Let t(detri), t(Lia), and t(FP) denote the CPU time needed for the overall execution of detri 1.2, and for the above long-integer and floating-point evaluations, respectively. With this, the formula t(detri*) = t(detri) - t(Lia) + t(FP) gives a rough estimate for the running time of a hypo9Recall that edge flips do not need any geometric tests; see section 4.1.1. 87 thetical floating-point version detri*, which (with some "magic") would produce the correct result. See figures 7.1 (c) and (d) for the corresponding slow-down factors, observed for the data sets in our sample with jSj ! 16k.10 The results show that long-integer arithemtic slows down the determinant evaluations by a factor of 11, which leads to an overall slow-down factor for detri 1.2 of roughly 3, when compared to the hypothetical floating-point implementation. We believe that this is an adequate price to pay for a compact and robust implementation of a rather involved geometric algorithm. A price, that can be reduced when more efficient implementations for long-integer arithmetic are available; see section 5.2.3 for a short discussion on this subject. 10Due to disk space limitations, the runs of the experiment for figures 7.1 (c) and (d) could not be performed for data with more than 16k points. 0 10000000 20000000 30000000 40000000 50000000 60000000 0 5000 10000 15000 20000 25000 Data Size muladd/sub (a) Long-integer operations. (\Sigma 1:1%) 0 10000000 20000000 30000000 40000000 50000000 0 5000 10000 15000 20000 25000 Data Size psubpadd (b) Subtractions versus additions. (\Sigma 1:1%) 0 2 4 6 8 10 12 14 0 5000 10000 15000 20000 25000 t(Lia) / t(FP) Data Size (c) Observed slow-down factors of long-integer versus floating-point arithmetic with respect to determinant evaluations. Sample mean: 10:79 (\Sigma 2:9%) 0 2 4 6 8 10 12 14 0 5000 10000 15000 20000 25000 t(detri) / t(detri*) Data Size (d) Hypothetical slow-down factors of long-integer versus floating-point arithemtic with respect to the overall program. Sample mean: 2:98 (\Sigma 3:1%) Figure 7.1: Long-integer arithmetic in detri 1.2. 88 Performance Evaluation. In the following, we use a series of plots to illustrate the actual performance of detri 1.2. The data constitutes the sample means Y of the measurements in question (sample size N = 13). The corresponding confidence intervals are usually too small to be drawn in the plot. Thus, we simply state the worst confidence interval, with respect to the whole sample, in the captions of the figures. Figure 7.2 (a) plots the actual run-time measurements of detri 1.2 for the sample data sets. Observe a certain variance between the data of different categories. In particular, the "r xxx" point sets, which were created with a random number generator, performed slightly better than, for example, the "smoke xx" samples of a chaotic attractor; compare this with table 7.1. In any case, a positive correlation between the required CPU seconds and the number of necessary primitive operations (that is, the required evaluations of "-determinants) is obvious; see figure 7.2 (b): this plot shows the total count of 5-by-5 and 4-by-4 minor evaluations. Note that all 5-by-5 determinants occur in in sphere tests, but 4-by-4 minors are due to positive tests and to degenerate in sphere tests. For example, observe the particular high number of evaluations of 4-by-4 determinants for a certain 0 100 200 300 400 500 600 0 5000 10000 15000 20000 25000 CPU Seconds Data Size detri 1.2 (a) Run-time measurements. (\Sigma 9:1%) 0 500000 1000000 1500000 2000000 0 5000 10000 15000 20000 25000 Data Size 5-by-5 (in_sphere)4-by-4 (positive) (b) Evaluations of "-determinants. (\Sigma 1:0%) 0 50000 100000 150000 200000 250000 300000 0 5000 10000 15000 20000 25000 Data Size triangle flips edge flips skips (c) Flips and skips. (\Sigma 2:1%) 0 0.5 1 1.5 2 2.5 0 5000 10000 15000 20000 25000 (Flips + Skips) / Triangles Data Size detri 1.2 (d) Steps per triangle. (\Sigma 2:1%) Figure 7.2: Performance of detri 1.2. 89 subgroup of the data sets. These spikes belong to the "ku2," "mozart xxk," and "bach xxk" data sets.11 They all share a high degree of degeneracy (see table 7.1). Plot (c) in figure 7.2 draws the number of triangle(-to-edge) and edge(-to-triangle) flips, including the unsuccessful flips (or "skips"). Adding these numbers up results in a run-time independent cost model for the incremental flip method. It measures the number of times the while loop of the flip procedure in table 4.2 is executed. Figure 7.2 (d) divides this number of "steps" of the algorithm by nf = jF2j, the number of triangles in the final Delaunay triangulation. We see that the performance of the algorithm is sensitive to the size of the output. Within the window of our sample data, we can say that the algorithm never needs more that 1:9nf steps to construct the Delaunay triangulation, and nf is the output size. Comparison with other flip algorithms. How does section 4.1.1's variant of the incremental flip algorithm compare with the other methods? In the following, detri 1.2 is compared with two other implementations: rdetri, a code implementing the incremental flip as described in Edelsbrunner and Shah [28],12 and detri 1.1, an earlier, non-randomized version of the algorithm, which essentially follows the method of Joe [46].13 Figure 7.3 shows four plots comparing the three programs, which, by the way, all use SoS to cope with degeneracies. The major difference between our variant of the incremental flip and the one of [28] is that they use a history dag to implement the point location; see section 4.1.1. In contrast to point location with the history dag, the simple "walk through" method of section 4.1.2 does not use any extra memory. This explains the difference in memory requirements between the two implementations; see figure 7.3 (a). It is also twice as fast in locating points in the triangulation, provided we employ the "random sample" heuristic of section 4.1.2 to find a starting triangle relatively close to the query point. Figure 7.3 (b) documents this by plotting the "average search length" counting the average number of positive tests required to locate a point with respect to one run. For the date in our experiment, we see that the history dag (of rdetri) always performs worse than the "walk through" method (of detri 1.2). However, the latter seems to have a higher variance between the different categories of data. In the range below n = 10k, we have poorer performance for data with many points on the boundary of the convex hull (for example, "ma2" and "teapot") and for data with a high count of degeneracies (for example, "ku2"). Beyond n = 10k, three groups 11Actually, "grid1" has a spike in 4-by-4 determinants as well. It is just not very visible in the plot. 12The rdetri code is due to Facello [31]. 13The detri 1.1 code is a predecessor of detri 1.2, both written by the author. 90 with different levels of average search length emerge. The random coordinates of "r xxx" perform best, followed by the chaotic attractors of the "smoke xx" data; the highly degenerated range data sets of Mozart's and Bach's busts perform the worst. A further investigation of the "walk through" strategy and its performance would be interesting. The same is true for a better tuning of the size of the random sample to determine a "good" starting point for the search. The current implementation simply chooses the closest triangle out of K = 15 log(ni) random triangles, with ni = jDij, the number of points in the current triangulation. This procedure yields reasonable results. The detection of the "good" initial triangle usually accounts for no more that 15% of the time spent in function search. Of course, search itself is not responsible for more that 10% of the total running time of the program; see table 7.2. Therefore, the overall performance of detri 1.2 is only slightly better than the one of rdetri: on average by 14%. Figure 7.3 (c) plots the speedup factor detri 1.2 versus rdetri; the average factor in our sample is 1:14. 0 20 40 60 80 100 120 140 160 0 5000 10000 15000 20000 25000 Megabytes Data Size detri 1.2rdetri (a) Memory requirements. (\Sigma 0:2%) 0 10 20 30 40 50 60 70 80 90 0 5000 10000 15000 20000 25000 Average Search Length Data Size detri 1.2rdetri (b) Average search length, point location.(\Sigma 1:2%) 0 2 4 6 8 10 0 5000 10000 15000 20000 25000 Data Size detri 1.2 vs. rdetri (c) Speedup: detri1.2 versus rdetri. 0 2 4 6 8 10 0 5000 10000 15000 20000 25000 Data Size detri 1.2 vs. detri 1.1 (d) Speedup: detri 1.2 versus detri 1.1. Figure 7.3: Comparing detri 1.2 with rdetri and detri 1.1. 91 When comparing detri 1.2 with the non-randomized version detri 1.1, we achieve a speedup of roughly 2; see figure 7.3 (d). Note even more dramatic speedup factors ( ? 6) for volumetric data ("ku2" and "rat T2[ab]"). The 1.1 version of detri essentially implements Joe's original algorithm [46]. However, there is one difference. Rather than inserting the points sorted along the x1-, or some other axis, detri 1.1 inserts them sorted along some arbitrary direction. We observed that Joe's algorithm performs particulary badly if the data points come in planes orthogonal to the chosen direction in which these points are to be sorted. By choosing this direction in an random fashion, this problem can be eased.14 Artifacts of SoS. Note that SoS does not allow the user to "selectively" remove certain types of degeneracies, while retaining others. Rather, it fully implements what theoreticians call "general position." Moreover, observe that algorithms employing SoS produce results for the perturbed data, S("), and not for the original one, S. Some postprocessing can be used to remove or repair parts of the output that collapses to a degenerate state because of the degeneracy of the input data. This is mentioned in [27]. In the case of Delaunay triangulations consider, for example, coplanar points on the boundary of the convex hull. The perturbation will move some of these points off the boundary, while some will remain there. In the Delaunay triangulation of S("), the points, that moved from the boundary of conv(S) towards the inside, will be covered by "flat" tetrahedra. These tetrahedra have infinitesimal thickness in the perturbed setting, and are completely flat in the original setting. It is easy to identify and remove these tetrahedra in a postprocessing step. However, it cannot be denied, that artifacts like these contribute to the overhead of SoS. It certainly costs some time to (a) construct them, and then, (b) remove them afterwards, if necessary. For example, 26:6% of the tetrahedra in the Delaunay triangulation of the highly degenerate "grid1" data set are flat; for "mozart 07k" the percentage is 24:8%. These kind of data sets consist of a particularly large portion of coplanar points, many of which lie on the boundary of the convex hull (in the unperturbed setting). Usually, this is not the case, and the number of totally flat tetrahedra, due to SoS, is not a factor. 14It should be noted that, in order to implement Joe's algorithm, we need one more primitive not yet mentioned in this text: the smaller predicate, which compares coordinates. An SoS implementation of this predicate was stated in [27]. It is quite easy to modify this primitive operation such that it compares coordinates along a given direction (that is, it compares the distance from two points to a plane). Details are omitted. 92 7.2 Constructing the Alpha Shape File and Data Structure In section 4.2 we briefly discussed how to construct the data structure representing the whole family of ff-shapes of a given point set. This essentially means to augment the data structure storing the Delaunay triangulation with the face-classifying intervals of section 2.3. In addition, the ffspectrum, and the master list need to be constructed; see section 3.2. Let us now take a short look at the empirical behavior of mkalf 1.0, our actual implementation of this procedure. Plots (a) and (b) in figure 7.4 show the run-time behavior of the program. Its running time is roughly proportional to the size of the master list, which is of the same order as the number of faces in the Delaunay triangulations. The real bottleneck, however, is mkalf's memory requirements. While it is true that we need no more than 60 bytes for each entry in the master list, it definitely sums up in the end; see plots (c) and (d) in figure 7.4. For data sets with n * 25k, we need at least 70 Mb 0 200 400 600 800 1000 0 5000 10000 15000 20000 25000 CPU Seconds Data Size mkalf 1.0 (a) Run-time measurements. 0 0.2 0.4 0.6 0.8 1 0 5000 10000 15000 20000 25000 Milliseconds / |Master List| Data Size mkalf 1.0 (b) CPU time per entry in master list. Sample mean: 0:53 (\Sigma 4:9%). 0 10 20 30 40 50 60 70 80 0 5000 10000 15000 20000 25000 Megabytes Data Size total main memorylong-integer ratios alpha-shape (file) (c) Memory requirements. 0 10 20 30 40 50 60 70 80 0 5000 10000 15000 20000 25000 Bytes / |Master List| Data Size mkalf 1.0 (d) Bytes per entry in master list. Sample mean: 52:72 (\Sigma 1:2%). Figure 7.4: Performance of mkalf 1.0. 93 of main memory, and 30 Mb of file space (this includes the file storing the Delaunay triangulation). The total main memory requirements of mkalf include the hash-tables for face indexing, temporary lists to construct and sort the master list, and the exact long-integer ratios of the threshold values in the ff-spectrum, none of which is stored on file at the end. At this stage of the implementation, considerable care has already been taken to keep most record types in the data structure as compact as possible, but further work is needed to cut down mkalf's appetite for memory. 94 Chapter 8 Conclusion Three-dimensional alpha shapes are an approach to formalize the intuitive notion of "shape" for spatial point sets data which occurs frequently in the computational sciences. The alpha shape is a concrete geometric concept that is uniquely defined. It thus stands in sharp contrast to many common concepts in computer graphics, such as for example, the isosurface, which is approximate by definition and whose exact form depends on the algorithm used to construct it. Alpha shapes are a generalization of convex hulls. Given a finite point set S and a real parameter ff, the ff-shape of S is a polytope which is not necessarily convex or connected. The set of all real values ff leads to a discrete family of shapes capturing the intuitive notion of "crude" versus "fine" shapes of a point set. The parameter ff controls the maximum "curvature" of any cavity of the polytope. Three-dimensional alpha shapes have immediate benefits to scientists in many areas. In biology, for instance, it can be used to construct three-dimensional models of proteins, and to identify holes and tunnels inside the protein; this is useful for determining the possible location and number of water molecules. In astronomy, the concept can help in analyzing the distribution and shapes of galaxies in the universe. In chemistry and physics, it can, for example, be employed to find a meaningful definition for the volume of nano-crystals, which are small clusters of atoms attached to each other by weak forces. Other possible application areas include solid modeling, mesh generation, and feature extraction. 95 A major contribution of this thesis is the actual implementation of the discussed theoretical concepts. This includes a program for constructing three-dimensional Delaunay triangulations, which is based on a variant of the incremental-flip algorithm. Robustness is achieved via Simulation of Simplicity (SoS), a symbolic perturbation scheme that copes with degenerate input data. SoS assumes exact arithmetic, and an efficient implementation of thus necessary long-integer arithmetic is provided as well. The source code for the Delaunay triangulation, SoS, and long-integer arithmetic is freely available via anonymous ftp from cs.uiuc.edu, directory: /pub/edels/geometry. The Silicon Graphics executables of "Alvis -- A 3D Alpha Shape Visualizer" are also available from ftp.ncsa.uiuc.edu, directory SGI/Alpha-shape. The distribution on both ftp servers includes most of the data sets cited in this thesis. The Alvis tool allows the user to interactively select ff-values and display the corresponding shape. A small collection of signature graphs provides combinatorial, metric, and topological information about individual shapes, and also the entire family of shapes. The signatures turn out to be useful in quickly selecting interesting shapes in the typically huge family. We see the program as an extension to the theoretical work, as an "animated figure," so to speak. One of its purposes is to effectively convey the concept of ff-shapes to the engeneering and scientific computing community, where we believe it will find many applications. 96 References [1] F Aurenhammer. Voronoi diagrams -- a survey of a fundamental data structure. ACM Computing Surveys, 23(3):345-405, 1991. [2] G B Beretta. An Implementation of a Plane-Sweep Algorithm on a Personal Computer. PhD thesis, Eidgen"ossische Technische Hochschule, Z"urich, 1984. [3] J D Boissonnat. Geometric structures for three-dimensional shape representation. ACM Transactions on Graphics, 3(4):266-286, 1984. [4] A Bowyer. Computing Dirichlet tessellations. The Computer Journal, 24(2):162-166, 1981. [5] E Brisson. Representing geometric structures in d dimensions: Topology and order. Discrete and Computational Geometry, 9(4):387-426, 1993. [6] A Bro/nsted. An Introduction to Convex Polytopes. Graduate Texts in Mathematics. Springer-Verlag, New York, 1983. [7] K Q Brown. Voronoi diagrams from convex hulls. Information Processing Letters, 9(5):223-228, 1979. [8] P Cignoni, C Montani, and R Scopigno. A merge-first divide-and-conquer algorithm for Ed Delaunay triangulations. Internal Report C92/19, Institute of the National Research Council (CNUCE), Pisa, 1992. [9] M L Connolly. Analytical molecular surface calculation. Journal of Applied Crystallography, 16:548-558, 1983. [10] M L Connolly. Computation of molecular volume. Journal of the American Chemical Society, 107(5):1118-1124, 1985. [11] M L Connolly. Molecular surface triangulation. Journal of Applied Crystallography, 18:499-505, 1985. [12] M L Connolly. Molecular interstitial skeleton. Computers & Chemistry, 15(1):37-45, 1991. [13] T H Cormen, C E Leiserson, and R L Rivest. Introduction to Algorithms. MIT Press, Cambridge, Massachusetts, 1990. [14] B Delaunay. Sur la sph`ere vide. Izvestiya Akademii Nauk SSSR, Otdelenie Matematicheskii i Estestvennyka Nauk, 7:793-800, 1934. [15] C Delfinado and H Edelsbrunner. An incremental algorithm for Betti numbers of simplicial complexes. In Proceedings of the Nineth Annual Symposium on Computational Geometry, pages 232-239, 1993. 97 [16] D P Dobkin. Computational geometry and computer graphics. Proceedings of the IEEE, 80(9):1400-1411, 1992. [17] D P Dobkin and M J Laszlo. Primitives for the manipulation of three-dimensional subdivisions. Algorithmica, 4(1):3-32, 1989. [18] R A Drebin, L Carpenter, and P Hanrahan. Volume rendering. Computer Graphics, 22(4):65-74, 1988. ACM Siggraph '88. Conference Proceedings. [19] R A Dwyer. Higher-dimensional Voronoi diagrams in linear expected time. Discrete and Computational Geometry, 6(4):343-367, 1991. [20] M D Dyksterhouse. An alpha-shape view of our universe. Master's thesis, Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, Illinois, 1992. [21] H Edelsbrunner. Algorithms in Combinatorial Geometry. Springer-Verlag, Berlin, 1987. [22] H Edelsbrunner. Geometric algorithms. To appear in P Gruber and J Wills, editors, Handbook of Convex Geometry, North-Holland, 1992. [23] H Edelsbrunner. Weighted alpha shapes. Technical Report UIUCDCS-R-92-1760, Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, Illinois, 1992. [24] H Edelsbrunner. The union of balls and its dual shape. In Proceedings of the Nineth Annual Symposium on Computational Geometry, pages 218-231, 1993. [25] H Edelsbrunner and P Fu. Measuring space filling diagrams. NCSA Technical Report UIUCNCSA-TR93-010, National Center for Super Computing Applications, Urbana, Illinois, 1993. [26] H Edelsbrunner, D G Kirkpatrick, and R Seidel. On the shape of a set of points in the plane. IEEE Transactions on Information Theory, IT-29(4):551-559, 1983. [27] H Edelsbrunner and E P M"ucke. Simulation of Simplicity: A technique to cope with degenerate cases in geometric algorithms. ACM Transactions on Graphics, 9(1):66-104, 1990. [28] H Edelsbrunner and N R Shah. Incremental topological flipping works for regular triangulations. In Proceedings of the Eighth Annual Symposium on Computational Geometry, pages 43-52, 1992. [29] I Emiris and J Canny. A general approach to removing degeneracies. In Proceedings of the 32-nd Annual IEEE Symposium on Foundation of Computer Science, pages 405-413, 1991. [30] I Emiris and J Canny. An efficient approach to removing geometric degeneracies. In Proceedings of the Eighth Annual Symposium on Computational Geometry, pages 74-82, 1992. [31] M A Facello. Constructing Delaunay and regular triangulations in three dimensions. Technical Report UILU-ENG-93-1712, Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, Illinois, 1993. [32] J Fairfield. Contoured shape generation: Forms that people see in dot patterns. In Proceedings of the IEEE Conference on Cybernetics and Society, pages 60-64, 1979. [33] J Fairfield. Segmenting dot patterns by Voronoi diagram concavity. IEEE Transactions on Pattern Analysis and Machine Intelligence, 5(1):104-110, 1983. [34] A R Forrest. Computational geometry in practice. In R A Earnshaw, editor, Fundamental Algorithms for Computer Graphics, pages 707-724. Springer-Verlag, Berlin, 1985. 98 [35] S Fortune. Stable maintenance of point-set triangulations in two dimensions. Unpublished manuscript, 1989. [36] S Fortune and C J Van Wyk. Efficient exact arithmetic for computational geometry. In Proceedings of the Nineth Annual Symposium on Computational Geometry, pages 163-172, 1993. [37] C Gh'elis and J Yon. Protein Folding. Academic Press, 1982. [38] P J Giblin. Graphs, Surfaces, and Homology. Second edition. Chapman and Hall, London, 1981. [39] P J Green and R Sibson. Computing Dirichlet tessellations in the plane. The Computer Journal, 21(2):168-173, 1978. [40] L Guibas, D Salesin, and J Stolfi. Epsilon geometry: Building robust algorithms from inprecise computations. In Proceedings of the 5th Annual Symposium on Computational Geometry, pages 208-217, 1989. [41] L J Guibas, D E Knuth, and M Sharir. Randomized incremental construction of Delaunay and Voronoi diagrams. In G Goos and J Hartmanis, editors, Lecture Notes in Computer Science 443: ICALP 90 Proceedings, pages 414-431. Springer-Verlag, New York, 1990. [42] L J Guibas and J Stolfi. Primitives for manipulation of general subdivisions and the computation of Voronoi diagrams. ACM Transactions on Graphics, 4(2):74-123, 1985. [43] J W Jaromczyk and G T Toussaint. Relative neighborhood graphs and their relatives. Proceedings of the IEEE, 80(9):1502-1571, 1992. [44] R A Jarvis. Computing the shape hull of points in the plane. In Proceedings of the IEEE Computing Society Conference on Pattern Recognition and Image Processing, pages 231-241, 1977. [45] B Joe. Three-dimensional triangulations from local transformations. SIAM Journal on Scientific and Statistical Computing, 10(4):718-741, 1989. [46] B Joe. Construction of three-dimensional Delaunay triangulations using local transformations. Computer Aided Geometric Design, 8(2):123-142, 1991. [47] M Karasick, D Lieber, and L R Nackmann. Efficient Delaunay triangulation using rational arithmetic. ACM Transactions on Graphics, 10(1):71-91, 1991. [48] D G Kirkpatrick and J D Radke. A framework for computational morphology. In G T Toussaint, editor, Computational Geometry, pages 234-244. Elsevier North Holland, New York, 1985. [49] D E Knuth. The Art of Computer Programming. Volume 2: Seminumerical Algorithms. Addison-Wesley, Reading, Massachusetts, 1969. [50] D E Knuth. Axioms and Hulls. Lecture Notes in Computer Science 606. Springer-Verlag, Berlin, 1992. [51] C L Lawson. Generation of a triangular grid with applications to contour plotting. Memo 299, Jet Propulsion Laboratory, Pasadena, California, 1972. [52] C L Lawson. Software for C1 surface interpolation. In J R Rice, editor, Mathematical Software III, pages 161-194. Academic Press, New York, 1977. [53] C Lee. Regular triangulations of convex polytopes. In P Gritzmann and B Sturmfels, editors, Applied Geometry and Discrete Mathematics: The Victor Klee Festschrift, pages 443-456. American Mathematical Society, Providence, RI, 1991. 99 [54] D T Lee and F P Preparata. Computational geometry -- a survey. IEEE Transactions on Computers, 33(12):1072-1101, 1984. [55] W E Lorensen and H E Cline. Marching cubes: A high resolution 3D surface construction algorithm. Computer Graphics, 21(4):163-169, 1987. ACM Siggraph '87. Conference Proceedings. [56] D W Matula and R R Sokal. Properties of Gabriel graphs relevant to geographic variation research and the clustering of points in the plane. Geographical Analysis, 12(3):205-222, 1980. [57] P McMullen and G C Shepard. Convex Polytopes and the Upper Bound Conjecture. Cambridge Unversity Press, London, 1971. [58] K Mehlhorn. Data Structures and Algorithms 3: Multidimensional Sorting and Computational Geometry. Springer-Verlag, Berlin, 1984. [59] V J Milenkovic. Verifiable implementations of geometric algorithms using finite precision arithmetic. Artificial Intelligence, 37:377-401, 1988. [60] V J Milenkovic. Calculating approximate curve arrangements using rounded arithmetic. In Proceedings of the 5th Annual Symposium on Computational Geometry, pages 197-207, 1989. [61] W W Moss. Some new analytic and graphic approaches to numerical taxonomy, with an example from the Dermanyssidae (Acari). Systematic Zoology, 16:177-207, 1967. [62] E P M"ucke. Sos -- a first implementation. Master's thesis, Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, Illinois, 1988. [63] J Nievergelt and K Hinrichs. Algorithms and Data Structures with Applications to Graphics and Geometry. Prentice-Hall, 1993. [64] T Ottman, G Thiemt, and C Ullrich. Numerical stability of geometric algorithms. In Proceedings of the 3rd Annual Symposium on Computational Geometry, pages 119-125, 1987. [65] C Pickover. A note on rendering 3-D strange attractors. Computers and Graphics, 12(2):263-267, 1988. [66] F P Preparata and S J Hong. Convex hulls of finite sets of points in two and three dimensions. Communications of the ACM, 20(2):87-93, 1977. [67] F P Preparata and M I Shamos. Computational Geometry -- An Introduction. Springer-Verlag, New York, 1985. [68] V T Rajan. Optimality of the Delaunay triangulation on IRd. In Proceedings of the 7th Annual Symposium on Computational Geometry, pages 357-363, 1991. [69] F M Richards. Calculation of molecular volumes and areas for structures of known geometry. Methods in Enzymology, 115:440-464, 1985. [70] F M Richards. The protein folding problem. Scientific American, 264(1):54-63, 1991. [71] H Rosenberger. Degeneracy Control in Geometric Programs. PhD thesis, Department of Computer Science, University of Illinois at Urbana-Champaign, Ubana, Illinois, 1990. [72] R Seidel. A convex hull algorithm optimal for point sets in even dimensions. Technical Report 81-14, Department of Computer Science, University of British Columbia, Vancouver, British Columbia, 1981. [73] R Seidel. Constructing higher-dimensional convex hulls in logarithmic cost per face. In Proceedings of the 18th Annual ACM Symposium on Theory of Computing, pages 484-507, 1986. 100 [74] R Seidel. Exact upper bounds for the number of faces in d-dimensional Voronoi diagrams. In P Gritzmann and B Sturmfels, editors, Applied Geometry and Discrete Mathematics: The Victor Klee Festschrift, pages 517-530. American Mathematical Society, Providence, RI, 1991. [75] B Serpette, J Vuillemin, and J C Herv'e. BigNum: A Portable and Efficient Package for ArbitraryPrecision Arithmetic. INRIA. [76] G T Toussaint. The relative neighborhood graph of a finite planar set. Pattern Recognition, 12(4):261- 268, 1980. [77] G T Toussaint (Editor). Computational Geometry. Elsevier North Holland, New York, 1985. [78] R C Veltkamp. Closed Object Boundaries from Scattered Points. PhD thesis, Erasmus Universiteit, Rotterdam, 1992. [79] G Voronoi. Nouvelles applications des param`etres continus `a la th'eorie des formes quadratiques. Premier M'emoire: Sur quelques propri'et'es des formes quadratiques positives parfaites. Journal f"ur die Reine und Angewandte Mathematik, 133:97-178, 1907. [80] G Voronoi. Nouvelles applications des param`etres continus `a la th'eorie des formes quadratiques. Deuxi`eme M'emoire: Recherches sur les parall'ello`edres primitifs. Journal f"ur die Reine und Angewandte Mathematik, 134:198-287, 1908. [81] C K Yap. Symbolic treatment of geometric degeneracies. In Proceedings of the 13-th IFIP Conference on System Modeling and Optimization, Chuo University, Tokyo, 1987. [82] C K Yap. A geometric consistency theorem for a symbolic perturbation scheme. In Proceedings of the 4th Annual Symposium on Computational Geometry, pages 134-142, 1988. [83] C K Yap. Towards exact geometric computation. In Proceedings of the Fifth Canadian Conference on Computational Geometry, pages 405-419, 1993. 101 Vita Ernst M"ucke was born in Graz, Austria, on August 13, 1962. He received the Dipl.-Ing. degree (graduated engineer) in information and data processing from the Technical University of Graz in 1985, and the M.S. degree in computer science from the University of Illinois at Urbana-Champaign in 1988. He completed a Ph.D. degree in computer science at the same university in 1993. His research interests include computational geometry and scientific visualization, data structures and algorithms, software design and computer graphics. 102