Minkowski Sum Construction and other Applications of Arrangements of Geodesic Arcs on the Sphere
aa r X i v : . [ c s . C G ] J un TEL AVIV UNIVERSITY @ ביבא-לת תטיסרבינוא
Raymond and Beverly SacklerFaculty of Exact SciencesThe Blavatnik School of Computer Science
Minkowski Sum Construction andother Applications of Arrangementsof Geodesic Arcs on the Sphere
Thesis submitted for the degree of “Doctor of Philosophy”by
Efraim Fogel
This work has been carried out under the supervision ofProf. Dan HalperinSubmitted to the Senate of Tel-Aviv UniversityOctober 2008
Acknowledgements
This research project would not have been possible without the support of many people.I wish to express my gratitude to my advisor, Prof. Dan Halperin who was abundantlyhelpful and offered invaluable assistance, support, insights, and guidance. This sentense isfar short of describing the magnitude of positive influence Danny had on this research andthis researcher.I would like to thank Ron Wein and Ophir Setter from Tel-Aviv University and Eric Berberichfrom Max-Planck-Insitut f¨ur Informatik, for ongoing and fruitful collaboration. I would alsolike to thank all other members of the applied computational geometry group at the computerscience school in Tel-Aviv University who made the study duration festive.I would like to thanks the members of the
Cgal
Editorial Board in particular and allmembers of the
Cgal developers’ community in general for sharing their wisdom throughmany useful advises and rich discussions.During the years I participated in the EU-funded projects
Ecg (Effective ComputationalGeometry for Curves and Surfaces; contract No. IST-2000-26473),
Movie (Motion Planningin Virtual Environments, contract No. IST-2001-39250) and
Acs (Algorithms for ComplexShapes; contract No. IST-006413) projects. I wish to thank all the other participants inthese projects, with whom I enjoyed working.Finally, I would like to express my love and gratitude to my beloved families; for theirunderstanding and endless love, through the duration of my studies.iii
Abstract
We present two exact implementations of efficient output-sensitive algorithms that compute
Minkowski sums of two convex polyhedra in R . We do not assume general position. Namely,we handle degenerate input, and produce exact results. We provide a tight bound on theexact maximum complexity of Minkowski sums of polytopes in R in terms of the number offacets of the summand polytopes. The complexity of Minkowski sum structures is directlyrelated to the time consumption of our Minkowski sum constructions, as they are outputsensitive. We demonstrate the effectiveness of our Minkowski-sum constructions throughsimple applications that exploit these operations to detect collision between, and answerproximity queries about, two convex polyhedra in R .The algorithms employ variants of a data structure that represents arrangements embed-ded on two-dimensional parametric surfaces in R , and they make use of many operationsapplied to arrangements in these representations. We have developed software componentsthat support the arrangement data-structure variants and the operations applied to them.These software components are generic, as they can be instantiated with any number type.However, our algorithms require only (exact) rational arithmetic. These software componentstogether with exact rational-arithmetic enable a robust, efficient, and elegant implementationof the Minkowski-sum constructions and the related applications. These software compo-nents are provided through a package of the Computational Geometry Algorithm Library( Cgal ) [5] called
Arrangement on surface 2 [WFZH07a]. The code of
Cgal in general,the
Arrangement on surface 2 package in particular, and all the rest of the code developedas part of this thesis adhere to the
Generic Programming paradigm and follow the
ExactGeometric Computation paradigm.We also present exact implementations of other applications that exploit arrangementsof arcs of great circles, also known as geodesic arcs, embedded on the sphere. For example,we implemented robust polyhedra central-projection and Boolean set-operations applied topoint sets embedded on the sphere bounded by geodesic arcs. We use them as basic blocks inan exact implementation of an efficient algorithm that partitions an assembly of polyhedrain R with two hands using infinite translations. This application makes extensive use ofMinkowski-sum constructions and other operations on arrangements of geodesic arcs embed-ded on the sphere. It distinctly shows the importance of exact computation, as imprecisecomputation might result with dismissal of valid partitioning-motions.We have produced three movies that explain some of the concepts portrayed in thisthesis [20]. The first movie [FFHL02] explains what Minkowski sums are and demonstrateshow they are used in various applications. The second movie [FH05] demonstrates the firstmethod we have developed to construct Minkowski-sums of convex polyhedra. The thirdmovie [FSH08b] illustrates exact construction and maintenance of arrangements induced bygeodesic arcs and applications that exploit such arrangements.Additional information is available at the following web sites: http://acg.cs.tau.ac.il/projects/internal-projects/gaussian-map-cubical Throughout the thesis a number in brackets (e.g., [20]) refers to the link list starting on page 118, andan alphanumeric string in brackets (e.g., [FFHL02]) is a standard bibliographic reference. v http://acg.cs.tau.ac.il/projects/internal-projects/arrangements-on-surfaces , http://acg.cs.tau.ac.il/projects/internal-projects/exact-complexity-of-minkowski-sums , and http://acg.cs.tau.ac.il/projects/internal-projects/arr-geodesic-sphere Auxiliary programs, source code, data sets, and documentation can be downloaded from . ontents CONTENTS k = 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 674.2 The Lower Bound for k = 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 694.2.1 Constructing P . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.2.2 Constructing P i , i ≥ P . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.3 Maximum Complexity of Minkowski Sums of Many Polytopes . . . . . . . . 744.3.1 The Lower Bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.3.2 The Upper Bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 ONTENTS vii5.3.6 Pairwise Minkowski Sum Projection . . . . . . . . . . . . . . . . . . . 865.3.7 Motion-Space Construction . . . . . . . . . . . . . . . . . . . . . . . 885.3.8 Motion-Space Processing . . . . . . . . . . . . . . . . . . . . . . . . . 885.4 Additional Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 895.5 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
A Software Components, Libraries, and Packages 101
A.1 Visual Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101A.2 Software Availability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102iii
CONTENTS ist of Figures G and G rotated about the Y axis . . . . . . . . . . . . . . 704.3 Different views of P . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 714.4 Views of P and P . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.5 The Minkowski sum of M , and M , rotated about the Y axis . . . . . . 734.6 The overlay of the Gaussian maps of three rotated tetrahedra . . . . . . . . 745.1 The Split Star assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77ix LIST OF FIGURES ist of Tables
Cgm planar maps . . . . . . . . 553.3 Complexities of primal and dual representations . . . . . . . . . . . . . . . . 563.4 Complexities of primal and dual Minkowski-sum representations . . . . . . . 613.5 Time consumption of the Minkowski-sum computation . . . . . . . . . . . . 615.1 Split Star partitioning directions and corresponding subassemblies . . . . . . 885.2 Time consumption of assembly partitioning of the Split Star . . . . . . . . . 90xiii
LIST OF TABLES he Guide is definitive. Re-ality is frequently inaccurate.
Douglas Adams Introduction
Let P and Q be two closed convex polyhedra in R d . The Minkowski sum of P and Q is theconvex polyhedron M = P ⊕ Q = { p + q | p ∈ P, q ∈ Q } . A polyhedron P translated by avector t is denoted by P t . Collision Detection is a procedure that determines whether P and Q overlap. The Separation Distance π ( P, Q ) and the Penetration Depth δ ( P, Q ) defined as π ( P, Q ) = min {k t k | P t ∩ Q = ∅ , t ∈ R d } ,δ ( P, Q ) = inf {k t k | P t ∩ Q = ∅ , t ∈ R d } are the minimum distances by which P has to be translated so that P and Q intersect orbecome interior disjoint, respectively. The problems of finding the distances above can alsobe posed given a normalized vector r that represents a direction, in which case the minimumdistance sought is in direction r . The Directional Penetration-Depth , for example, is definedas δ r ( P, Q ) = inf { α | P α~r ∩ Q = ∅ , α ∈ R } . We present two exact and robust implementations of efficient output-sensitive algorithmsto compute the Minkowski sum [FFHL02] of two convex polyhedra, polytopes for short, in R [FH05, FH07, FSH08a, BFH + R .1 Chapter 1. Introduction
We compare our Minkowski-sum constructions with three other methods that produceexact results we are aware of. One is a simple method that computes the convex hull of thepairwise sums of vertices of two polytopes. The second is based on Nef polyhedra embeddedon the sphere, and the third is based on linear programming. We conducted experimentswith a broad family of polytopes and compared the performance of our methods with theperformance of the other. The results reported in Table 3.5 clearly shows that both ourmethods are significantly faster.Each method we have developed uses a different variant of Gaussian maps, also knownas normal diagrams or slope diagrams, to maintain dual representations of polytopes. Eachmethod employs a different variant of a generic data-structure that represents arrangement embedded on two-dimensional parametric surfaces in R to maintain the dual representations.(Arrangements embedded on two-dimensional parametric surfaces are subdivisions of thesurface, as induced by curves embedded on the surface. They play a central role in thisthesis; see Chapter 2 for a formal definition and Figure 2.1 for an illustration.) Each methodmakes use of many operations applied to arrangements in the corresponding representations.We present a tight bound on the exact maximum complexity of Minkowski sums ofpolytopes in R [FHW07]. In particular, we prove that the maximum number of facetsof the Minkowski sum of k polytopes with m , m , . . . , m k facets respectively is boundedfrom above by P ≤ i 5) + P ≤ i ≤ k m i + (cid:0) k (cid:1) . Given k positive integers m , m , . . . , m k , we describe how to construct k polytopes with corresponding number offacets, such that the number of facets of their Minkowski sum is exactly P ≤ i 5) + P ≤ i ≤ k m i + (cid:0) k (cid:1) . When k = 2, for example, the expression above reduces to4 m m − m − m + 26.We also present an exact implementation of an efficient algorithm that partitions an as-sembly of polyhedra in R with two hands using infinite translations [FH08]. This applicationmakes extensive use of Minkowski-sum constructions and other operations on arrangementsof arcs of great circles, also known as geodesic arcs, embedded on the sphere, such as poly-hedra central-projection and Boolean set-operations of point sets embedded on the spherebounded by geodesic arcs. The assembly partitioning demonstrates the importance of exactcomputation, as imprecise computation might result with dismissal of valid partitioning-motions.Both methods that construct Minkowski sums and the related applications are imple-mented on top of the Computational Geometry Algorithms Library ( Cgal ) [FT07], andare mainly based on the arrangement package of the library [FWH04, WFZH07b, FHK + + 07, BFH + .2. Background: Minkowski Sums geometry traits class that handles thespecific family of curves they are interested in, which mainly involves algebraic computation.The separation is enabled by a modular design and conveniently implemented within the generic-programming paradigm. The separation is a key aspect of the package, as well as ofother central Cgal components, such as the various triangulation packages [BDTY00] andconvex-hull algorithms (see [cga07] for more details), has been forced since its early stages,and heightened by our new design.The package comes with many geometric traits classes that handle all kinds of curvesorganized in a structured hierarchy. In particular, we mention the geometric traits class thathandles continuous piecewise linear curves, referred to as polylines and the geometric traitsclass that handles geodesic arcs embedded on the sphere [FSH08b, FSH08a, BFH + − , ∗ , and / in unlimited precision over the rationals, such as therational number type CGAL::Gmpq based on GMP — Gnu’s Multi Precision library [12]. Minkowski sums are closely related to proximity queries [LM04]. For example, the mini-mum separation distance between two polytopes P and Q is equal to the minimum distancebetween the origin and the boundary of the Minkowski sum of P and the reflection of Q through the origin [CC86]. Computing Minkowski sums, collision detection and proximitycalculation constitute fundamental tasks in computational geometry [HKL04, LM04, Sha04].These operations are ubiquitous in robotics, solid modeling, design automation, manufac-turing, assembly planning, virtual prototyping, and many more domains; see, e.g., [BdLT97,EOR92, KR91, Lat91, VKSM05].The wide spectrum of ideas expressed in the massive amount of literature published aboutthe subject during the last three decades has inspired the development of quite a few usefulsolutions. For an extensive overview about the subject and a comprehensive list of packagessee [LM04]. However, only recent advances in the implementation of computational-geometryalgorithms and data structures made our exact, robust, and efficient implementation possi- Commonly referred to as a field number type. Chapter 1. Introduction ble. The exact geometric-computation paradigm [Yap04] designed for implementing compu-tational geometry algorithms particularly prevails in Cgal . While in general the underlyingarithmetic is highly time consuming compared to machine floating-point arithmetic, majorefficiency is gained by computing predicates only to sufficient precision to evaluate themcorrectly; see Section 1.3.2 and e.g., [PF06].Various methods to compute the Minkowski sum of two polyhedra in R have beenproposed. The goal is typically to compute the boundary of the sum provided in somerepresentation. The combinatorial complexity of the Minkowski sum of two polyhedra of m and n features respectively can be as high as Θ( m n ). One common approach to computeit is to decompose each polyhedron into convex pieces, compute pairwise Minkowski sums ofpieces of the two, and finally the union of the pairwise sums. Computing the exact Minkowskisum of non-convex polyhedra is naturally expensive. Therefore, researchers have focused oncomputing an approximation that satisfies some criteria, such as the algorithm presented byVaradhan and Manocha [VM06]. They guarantee a two-sides Hausdorff distance bound onthe approximation, and ensure that it has the same number of connected components as theexact Minkowski sum. Computing the Minkowski sum of two convex polyhedra remains akey operation, and this is what we focus on. The combinatorial complexity of the sum can beas high as Θ( mn ) when both polyhedra are convex. For the complexity of the intermediatecase, where only one polyhedron is convex, cf. [AS97, Sha04].Convex decomposition is not always possible, as in the presence of non-convex curvedobjects. In these cases other techniques must be applied, such as approximations usingpolynomial/rational curves in 2D [kLsKE98]. Seong at al. [SKS02] proposed an algorithmto compute Minkowski sums of a subclass of objects; that is, surfaces generated by slope-monotone closed curves. Flato and Halperin [AFH02] presented algorithms based on Cgal for robust construction of planar Minkowski sums. While the citations in this paragraphrefer to computations of Minkowski sums of non-convex polyhedra, and we concentrate onthe convex cases, the latter is of particular interest, as our method makes heavy use ofthe same software components, in particular the Cgal arrangement package, which wentthrough a few phases of improvements [FWH04, WFZH07b, BFH + 07, BFH + kinetic framework in two dimensions introduced byGuibas et al. [GRS83] was the definition of the convolution operation in two dimensions, asuperset of the Minkowski sum operation, and its exploitation in a variety of algorithmicproblems. Basch et al. extended its predecessor concepts and presented an algorithm tocompute the convolution in three dimensions [BGRR96]. An output-sensitive algorithmfor computing Minkowski sums of polytopes was introduced in [GS87]. Gritzmann andSturmfels [GS93] obtained a polynomial time algorithm in the input and output sizes forcomputing Minkowski sums of k polytopes in R d for a fixed dimension d , and Fukuda [Fuk04]provided an output-sensitive polynomial algorithm for variable d and k . Ghosh [Gho93]presented a unified algorithm for computing 2D and 3D Minkowski sums of both convex andnon-convex polyhedra based on a slope diagram representation. Computing the Minkowskisum amounts to computing the slope diagrams of the two objects, merging them (see detailsin Section 3.2.2,) and extracting the boundary of the Minkowski sum from the mergeddiagram. Wu et al. [WSD03] introduced an improved version of Ghosh’ algorithm for convex .3. Background: Programming Several definitions of the term generic programming have been proposed since it was firstcoined around the early sixties, along with the introduction of the Lisp programming lan-guage. Here we confine ourself to the classic notion first described by Musser et al. [MS88],who considered generic programming as a discipline that consists of the gradual lifting ofconcrete algorithms abstracting over details, while retaining the algorithm semantics andefficiency. Within this context, several approaches have been put into trial through the in-troduction of new features in existing computer languages, or even new computer languagesall together. The software described in this thesis is written in C ++ , a programming languagethat is well equipped for writing software according to the generic-programming paradigmthrough the extensive use of class templates and function templates. Concepts and Models One crucial abstraction supported by all contemporary computer languages is the subroutine(also known as procedure or function, depending on the programming language). Anotherabstraction supported by C ++ is that of abstract data typing, where a new data type isdefined together with its basic operations. C ++ also supports object-oriented programming,which emphasizes packaging data and functionality together into units within a runningprogram, and is manifested in hierarchies of polymorphic data-types related by inheritance.It allows referring to a value and manipulating it without needing to specify its exact type. Asa consequence, one can write a single function that operates on a number of types within aninheritance hierarchy. Generic programming identifies a more powerful abstraction (perhapsless tangible than other abstractions). It is a formal hierarchy of polymorphic abstractrequirements on data types referred to as concepts , and a set of classes that conform preciselyto the specified requirements, referred to as models . Models that describe behaviors arereferred to as traits classes [Mye98]. Traits classes typically add a level of indirection intemplate instantiation to avoid accreting parameters to templates.A generic algorithm has two parts: The actual instructions that describe the steps of thealgorithm, and a set of requirements that specify which properties its argument types mustsatisfy. The following swap() function is an example of the first part of a generic algorithm. Chapter 1. Introduction template Assignable [Aus99]. The int data type, for example, is a model of thisconcept, so it can be used to instantiate the function template [Aus99] [25].A concept is a set of requirements divided into four categories, namely, associated types,valid expressions, invariants, and complexity guarantees. When a type meets all requirementsof a concept, the type is considered a model of the concept. When a concept extends therequirements of another concept, the former is said to be a refinement of the latter. Associated Types are auxiliary types. For example, a type that represents a two-dimensionalpoint, namely Point 2 , is required by the arrangement geometry traits concept; seeSection 2.5. Valid Expressions are C ++ expressions that must compile successfully. For example, p =q , where p and q are both objects of type Point 2 . Valid expressions identify the setof operations a model of the concept must be able to perform. Invariants are run-time characteristics such as time and space complexity bounds. In ourcontext invariants typically take the form of preconditions and postconditions, whichmust always be satisfied. For example, a condition that requires that an input point p lies on an input curve c on invocation to a predicate that accepts both p and c asparameters. Having preconditions typically minimizes the concept, as the operationsprovided by a model must operate only on restricted arguments. Formally, removingpreconditions from, and introducing postconditions to, a requirement set results witha refined concept. Complexity Guarantees are maximum limits on the computing resources consumed bythe various expressions. Traits Classes The name “traits class” comes from a standard C ++ design pattern [Mye98], which providesa way of associating information with a compile-time entity (typically a type). For example,the standard class-template std::iterator traits Iso rectangle 2 (and implicitly also a point typesay Point 2 ). Moreover, it should supply two three-valued predicates that compare twopoints by their x -coordinates and by their y -coordinates, respectively. It should also supportthe construction of an axis-parallel iso-rectangle from four points that specify its extremal x and y -coordinates. Note, however, that the actual representation of points and rectan-gles (the coordinate system, the number-type used to represent the coordinates, etc.) andthe implementation of the traits-class operations is entirely decoupled from the function bounding rectangle() we have introduced.Consider an imaginary generic implementation of a data structure that handles ge-ometric arrangements embedded on two-dimensional parametric surfaces in space called Arrangement on surface 2 . Its prototype is listed below. This template class must be This notation means that pts begin points to the first point in the range, while pts end points afterthe range ends (it is therefore called a past-the-end iterator , and need not point to any valid point object). The predicate return value is SMALLER , EQUAL , or LARGER . Chapter 1. Introduction instantiated with a traits class that in turn defines a type that represents a certain family ofcurves, and some functions (or function objects; see [10] for an exact definition) that operateon curves of this family. template Cgal are parameterized by a model of the Kernel concept. Thechoice of a particular model determines, among the other, the type of arithmetic used, asexplained in the following sections. One can easily switch between different models of the Kernel concept, but here lies a trap, as Section 1.3.2 reveals. A kernel model that supportsexact arithmetic must be used to ensure robustness, although inexact arithmetic could beused at a certain risk. Libraries Alexander Stepanov began exploring the potential of compile-time polymorphism for rev-olutionizing software development in 1979. With the help of several other researchers hiswork evolved into a prime generic-programming library — the Standard Template Library( STL ). This library had became part of the C ++ standard library in 1994, approximatelyone year before early development of Cgal started; see Section 1.3.3 for details about theevolution of Cgal .Through the years a few other generic-programming libraries emerged. One notable li-brary in the context of computational geometry is Leda (Library of Efficient Data Types andAlgorithms), a library of combinatorial and geometric data types and algorithms [MN00] [17].Early development of Leda started in 1988, ten years before the first public release of Cgal became available. In some sense Leda is a predecessor of Cgal , although the two librariesare headed in different directions. While Leda is mostly a large collection of fundamentalgraph related and general purpose data-structures and algorithms, Cgal is a collection oflarge and complex data-structures and algorithms focusing on geometry.A noticeable influence on generic programming is conducted by the Boost online commu-nity, which encourages the development of free C ++ software gathered in the Boost librarycollection [3]. It is a large set of portable and high quality C ++ libraries that work wellwith, and are in the same spirit as, the C ++ Standard Template Library. The Boost GraphLibrary ( BGL ) [SLL02], which consists of generic graph-algorithms, serves a particularlyimportant role in our context. An arrangement instance, for example, can be adapted as a .3. Background: Programming BGL graph, and passed as input to generic algorithms already implemented in the BGL ,such as the Dijkstra shortest path algorithm. We use the BGL to compute the strongly con-nected components of a directed graph — a phase in the assembly partitioning operation;see Section 5.3.8 for more details about the application of this operation. Implementing geometric algorithms and data structures is notoriously difficult, much harderthan may seem when just considering the algorithm as described in a paper or a book.In the traditional computational-geometry literature two assumptions are usually made tosimplify the design and analysis of geometric algorithms. First, inputs are in “general po-sition”. That is, degenerate input (e.g., three curves intersecting at a common point) isprecluded. Secondly, operations on real numbers yield accurate results (the “real Ram ”model [PS85], which also assumes that each basic operation takes constant time). Unfortu-nately, these assumptions do not hold in practice, as degenerate input is commonplace inpractical applications and numerical errors are inevitable. Thus, an algorithm implementedwithout keeping this in mind may yield incorrect results, get into an infinite loop, or justcrash, while running on a degenerate, or nearly degenerate, input (see [KMP + 08, Sch00] forexamples). These pitfalls have become well known, and have been the subject of intensiveresearch [Sch00, Yap04].Indeed, the last decade has seen significant progress in the development of software forcomputational geometry. The mission of such a task, which Kettner and N¨aher [KN04]call geometric programming , is to develop software that is correct, efficient, flexible (namelyadaptable and extensible ), and easy to use. Separation of Topology and Geometry The use of the generic-programming paradigm enables a convenient separation between thetopology and the geometry of data structures. This is a key aspect in the design of ge-ometric software, and is put into practice, for example, in the design of Cgal polyhedra, Cgal triangulations, and our Cgal arrangements. This separation allows the convenientabstraction of algorithms and data structures in combinatorial and topological terms, re-gardless of the specific geometry of the objects at hand and the algebra used to representthem. This abstraction is realized through class and function templates that represent spe-cific data-structures and algorithmic frameworks, respectively. Consider again our imaginary Arrangement on surface 2 class template from the previous section; its improved prototypeis listed below. It is instantiated with two classes. The first, referred to as a geometric traits class, defines the set of geometric-object types and the operations on objects of these types.The second, defines the topological-object types and the operations required to maintain the Adaptability refers to the ability to incorporate existing user code, and extendibility refers to the abilityto enhance the software with more code in the same style. In this context, we sometimes say combinatorics instead of topology, and say algebra or numerics insteadof geometry. We always mean the same thing: The separation between the abstract, graph-like structure(the topology) from the actual embedding on the surface (the geometry). Chapter 1. Introduction incident relations among objects of these types. template Exact Geometric Computation The need for robust software implementations of computational-geometry algorithms hasdriven many researchers over the last decade to develop variants of the classic algorithms thatare less susceptible to degenerate inputs. The approaches taken to overcome the difficultiesin robustly implementing geometric algorithms can be roughly divided into two categories:(i) Exact computing, and (ii) fixed-precision approximation. In the latter approach thealgorithms are modified so that they can consistently cope with the limited precision ofcomputer arithmetic. In the former, which is the approach taken by Cgal in general and thearrangement package in particular, ideal computer arithmetic is emulated for the specific typeof objects being manipulated, and the code is prepared for successfully handling degenerateinput.Advances in computer algebra enabled the development of efficient software librariesthat offer exact arithmetic manipulations on unbounded integers, rational numbers ( GMP — Gnu’s multi-precision library [12]), and algebraic numbers (the Core library [KLPY99] [6]and the numerical facilities of Leda [MN00, Chapter 4] [17]). These exact-number typesserve as fundamental building-blocks in the robust implementation of many geometric ap-plications in general (see [Yap04] for a review) and of those that employ arrangements inparticular. .3. Background: Programming Cgal , such as numerical filtering.Here, computation is carried out using a number type that supports only inexact arithmetic(e.g., double-precision floating-point arithmetic), while applying a filter that checks whetherthe computation has reached a stage of uncertainty, an event referred to as a filter failure in the hacker’s jargon. When a filter failure occurs, the computation is re-done using exactarithmetic.Switching between number types and exact computation techniques, and choosing theappropriate components that best suit the application needs, is conveniently enabled throughthe generic-programming paradigm, as it typically requires only a minor code change reflectedin the instantiating of just a few data types. Cgal Chronicles Several research groups in Europe started to develop small geometry libraries on their ownin the early 1990s. A consortium of several sites in Europe and Israel was founded in1995 to cultivate the labor of these groups and gather their produce in a common librarycalled Cgal — the Computational Geometry Algorithms Library. The goal was to promotethe research in computational geometry and translate the results into useful, reliable, andefficient programs for industrial and academic applications [Ove96, cga07, KN04, FGK + Cgal development efforts to date.An INRIA startup, Geometry Factory [11] was founded in January 2003. The com-pany sells Cgal commercial licenses, support for Cgal , and customized developments basedon Cgal .In November 2003, when Version 3.0 was released, Cgal became an Open SourceProject [5], allowing new contributions from various resources. Most parts of Cgal arenow distributed under the GNU Lesser General Public License (GNU LGPL). Cgal has evolved through the years and is now representing the state-of-art in imple-menting computational geometry software in many areas. The implementations of the Cgal software modules described in this thesis are complete and robust, as they handle all de-generate cases. They rigorously adhere to the generic-programming paradigm to overcomeproblems encountered when effective computational geometry software is implemented. Aglimpse at the structure of Cgal is given in the following subsection. Cgal Content Cgal is written in C ++ according to the generic-programming paradigm described above.It has a common programming style, which is very similar to that of the STL . Its ap-2 Chapter 1. Introduction plication programming-interface (API) is homogeneous, and allows for a convenient andconsistent interfacing with other software packages and applications. The library consistsof about 900,000 lines of code divided among approximately 4,000 files. Cgal also comeswith numerous examples and demos. The manual comprises about 3,500 pages. There areapproximately 65 chapters arranged in 14 parts. The Arrangement on surface 2 package,for example, consists of about 140,000 lines of code divided among approximately 300 files,described in about 300 pages of a didactic manual.One distinguished piece consists of the geometric kernels [FGK + Cgal withvarious visualization tools (i.e., input and output streams). An important contribution ofthis piece is the number-type support. This piece also contains extensive debugging utilitiesthat handle warnings and errors that may result from unfulfilled conditions.The rest of the library offers a collection of geometric data structures and algorithms suchas convex hull, polygons and polyhedra and operations on them (Boolean operations, polygonoffsetting), 2D and 3D triangulations, Voronoi diagrams, surface meshing and surface subdi-vision, search structures, geometric optimization, interpolation, and kinetic data-structures.The 2D arrangements and its related data-structures naturally fit in. These data structuresand algorithms are parameterized by traits classes that define the interface between themand the primitives they use. In many cases, a kernel can be used as a traits class, or at leastthe subtypes of a kernel can be used as components of traits classes for these data structuresand algorithms. Cgal Arrangement Package History Cgal contains an elaborate and efficient implementation of a generic data-structure thatrepresents an arrangement induced by general types of curves embedded on a two-dimensionalparametric surface in R , but it has not been like this from the beginning. The first versionof the Cgal arrangement package supported only line segments, circular arcs, and restrictedtypes of parabolas embedded in the plane.While the first version supported only limited types of curves, it was originally designedwith the vision of supporting general curves [FHH + 00, Han00, HH00]. This vision wasreflected, among the other, through the separation between the topological and the geometricaspects enabled by the generic-programming paradigm (see Section 1.3.2). Most of theprinciples related to the topology, e.g., the use of a doubly-connected edge list ( Dcel ) tomaintain the incident relations between the arrangement cells (i.e., vertices, halfedges, andfaces) were conceived from the start. However, the types of curves that induce arrangements(see Section 2.5) gradually expanded. A couple of years after the introduction of Cgal arrangement Wein extended its implementation to support arcs of conics [Wei02]. Thearsenal of geometric traits continues to grow to date. .3. Background: Programming , together with Weinwe improved and refined the software design of the Cgal arrangement package [FWH04].This new design formed a common platform for a preliminary comparison between differentapproaches to handle arcs of conics [FHW + Cgal Version 3.2 released in 2006 included an arrangement package that supported onlybounded curves in the plane. This forced users to clip unbounded curves before insertingthem into the arrangement; it was the user responsibility to clip without loss of information.However, this solution is generally inconvenient and outright insufficient for some applica-tions. For example, representing the minimization diagram defined by the lower envelope ofunbounded surfaces in R [Mey06] generally requires more than one unbounded face, whereasan arrangement of bounded clipped curves contains a single unbounded face. Cgal Version 3.3 released a year later in 2007 already included an arrangement pack-age that handled unbounded planer curves. As a matter of fact it included much more.Together with Berberich, Melhorn, and Wein we observed the possibility to maximize codereuse by generalizing the various algorithms applied to arrangements and their implementa-tions so that they could be employed on a large class of surfaces and curves embedded onthem [BFH + 07, BFH + Cgal , expected to be released in 2010, is planned to include an ar-rangement package that constructs, maintains, and operates on arrangements embedded oncertain two-dimensional orientable parametric surfaces. The package already exists as a pro-totypical Cgal package under the new name Arrangement on surface 2 to better reflectits capabilities. For example, it includes (i) a geometric traits that handles geodesic arcs em-bedded on the sphere [FSH08b, BFH + + 07, BFH + + Cgal , such as the Boolean set operations 2 package. Only little effortis now required to support Boolean set-operations on point sets bounded by general curvesembedded on two-dimensional parametric surfaces. ECG is a Shared-Cost RTD (FET Open) Project of the European Union devoted to Effective Compu-tational Geometry for curves and surfaces [7]. Chapter 1. Introduction The rest of this thesis is organized as follows. In Chapter 2 we give an overview of the Cgal arrangement package, which provides the common infrastructure for all software solutionsdescribed in this thesis. This chapter contains selected sections from several papers andfrom a book chapter we have co-authored, and it provides new material that has not beenpublished yet. In particular, the chapter contains a description of the architecture of the ar-rangement generic data-structure, part of which also appeared in the chapter Arrangement of the book Effective Computational Geometry for Curves and Surfaces [FHK + Computational Geometry — Theory and Application (specialissue on Cgal ) [WFZH07b]. Preliminary results were first introduced in (i) the proceed-ings of the 12 th Annual European Symposium on Algorithms (ESA) [FWH04] and (ii) the1 st Workshop of Library-Centric Software Design [WFZH05]. The chapter also presents ar-rangements induced by general curves and embedded on parametric surfaces, and it offersa detailed description of a particular type of arrangement, namely arrangements of geodesicarcs on the sphere. The presentation of general arrangements embedded on parametric sur-faces is based on a joint work with E. Berberich, K. Melhorn, and R. Wein. Preliminaryresults of this work were first published at the 23 rd European Workshop on ComputationalGeometry (EWCG) [BFHW07]. Improved results appeared in the proceeding of the 15 th Annual European Symposium on Algorithms (ESA) [BFH + + th European Workshop on Computational Geometry (EWCG) [FSH08a].A movie rendering the results of this work was presented at the 24 th Annual Symposium onComputational Geometry (SoCG) . The proceeding of this symposium contains the relatedextended abstract [FSH08b]. Mature results were presented in a manuscript [BFH + R . We describehow the input polytopes in polyhedral-mesh representation both methods accept are con-verted into the corresponding internal representations unique to each method. We providethe theoretical concepts both methods rely on, and detailed descriptions specific to eachmethod. The first method was also published in the journal Computer Aided Design [FH07].A preliminary version of this paper appeared in the proceedings of the 8 th Workshop onAlgorithm Engineering and Experimentation (Alenex’06) [FH06]. The second method ap-peared in the publications [FSH08a, FSH08b, BFH + k polytopes in R in terms of the number of facets of the polytope summands. It is basedon collaborative work with C. Weibel. The results of this work were introduced in the pro-ceedings of the 23 rd annual Symposium on Computational Geometry (SoCG) [FHW07], andwere accepted for publication in the journal Discrete and Computational Geometry [FHW]. .4. Thesis Outline and Related Publications R — a solution to a problem in the domainof assembly planning. This application uses several types of operations on arrangements ofgeodesic arcs embedded on the sphere as basic blocks. In this context the chapter introducesexact implementations of additional applications that exploit geodesic arcs embedded on thesphere, such as polyhedra central-projection and Boolean set-operations applied to point setsembedded on the sphere bounded by geodesic arcs. Great parts of this chapter are extractsfrom a paper recently published in the 8 th International Workshop on Algorithmic Founda-tions of Robotics (WAFR) [FH08]. Specific background of assembly planning is provided atthe beginning of the chapter.We refer the reader to some ongoing research and future prospects and conclude inChapter 6.The software access-information along with some further design details are provided inthe Appendix.6 Chapter 1. Introduction common mistake that peo-ple make when trying todesign something completelyfoolproof is to underestimatethe ingenuity of completefools. Douglas Adams Arrangements on Surfaces Given a finite collection C of geometric objects (such as lines, planes, or spheres) the arrange-ment A ( C ) is the subdivision of the space where these objects reside into cells as inducedby the objects in C . In this thesis we deal only with arrangements embedded on certaintwo-dimensional orientable parametric surfaces in R , i.e., planes, cylinders, spheres, tori,and surfaces homeomorphic to them. In this case the objects in C embedded on the surface S are curves that divide S into a finite number of cells of dimension 0 ( vertices ), 1 ( edges )and 2 ( faces ). Figure 2.1 shows various types of arrangements embedded on two-dimensionalparametric surfaces.The Arrangement on surface 2 package of Cgal [WFZH07a] is a generic implemen- As a convention, Cgal prescribes the suffix for all data structures of planar objects and the suffix for all data structures of 3D objects. In the case of arrangements on surfaces the suffix indicates thedimension of the parameter space of the embedding surface. (a) (b) (c)Figure 2.1: (a) An arrangement of circles in the plane, (b) an arrangement of lines in the plane, and(c) an arrangement of geodesic arcs on the sphere. Chapter 2. Arrangements on Surfaces tation of a complete software package that constructs and maintains arrangements em-bedded on two-dimensional parametric surfaces [FHK + 07, FWH04, WFZH07b, BFH + + Arrangement on surface 2 package.As mentioned in Section 1.3.2 Cgal in general and the Arrangement on surface 2 pack-age of Cgal in particular follow the Exact Geometric Computation paradigm. The devel-oped code covers all cases to successfully handle degenerate input, while ideal computerarithmetic is emulated. While the Cgal package supports arrangements induced by generalalgebraic curves of arbitrary degree, the constructions and uses of arrangements describedin this thesis require only rational arithmetic. This is a key property that enables efficientimplementations of all the algorithms presented in the thesis. Cgal arrangements and their related components are the results of ongoing collaborativeresearch we were, and still are, deeply involved with; see Section 1.3.3 for the evolution ofthe arrangement package. We, in particular, had a significant impact on specific topicswithin this research area as follows: We came up with great parts of the geometry-traitsconcept hierarchy, and the corresponding sets of minimal requirements (see Section 2.5); wecontributed the geometry-traits module for polylines (see Section 2.5.4) and the geometryand topology traits modules for geodesic arcs embedded on the sphere (see Section 2.6); weled the design of the Boolean set-operation architecture (See Section 2.7.1), and we providedsolutions to numerous issues encountered during the development of all these components.This chapter provides a comprehensive overview of these components, with a slight emphasison asspects related to our work. It explains how arrangements induced by any type of curvescan be constructed, maintained, and used by other applications. Both closely related Mapc [KCMK00] [19] and Esolid [CKF + 04] [8] libraries consist ofan arrangement-construction module for algebraic curves. However, these implementationsmake some general-position assumptions on the input curves. The Leda library [MN00]includes geometric facilities that allow the robust construction and maintenance of planarmaps of line segments that may contain degeneracies. However, the resulting planar mapsare represented as simple graphs that cannot fully describe the topological structure of thearrangement. For example, it is impossible to encode the containment relation betweendisconnected components of the graph (i.e., to keep track of the holes contained in a face;see Section 2.3.2). Leda -based implementation of arrangements of conic curves and of cubiccurves were developed under the Exacus project [BEH + 05] [9]. Cgal ’s arrangement package was the first complete software implementation, designedfor constructing arrangements of arbitrary planar curves and supporting operations andqueries on such arrangements. The package was employed by many users to develop avariety of applications in various domains. For example, it was used to solve geometric op-timization problems [Rog03, COLYKT03], to construct Minkowski sums efficiently [AFH02, .2. Parametric Surfaces + 05] by considering the planar arrangements of theirprojected intersection curves. A better approach to compute such arrangements [BFH + + Arrangement on surface 2 package includes a genericimplementation of an elaborate version of this mechanism exploited by several higher-leveloperations supported by the package.The famous sweep-line algorithm of Bentley and Ottmann [BO79] was originally formu-lated for sets of non-vertical line segments, with the “general position” assumption that nothree segments intersect at a common point and no two segments overlap. Many generaliza-tions have been introduced ever since, such as the ability to handle more general curves [SH89]and to deal with degeneracies(see [dBvKOS00, Section 2.1] and [MN00, Section 10.7] for adiscussion about degeneracies induced by line segments).Effective algorithms for manipulating arrangements of curves have been a topic of con-siderable interest in recent years with an emphasis on exactness and efficiency of implemen-tation [FHK + Leda external package (LEP) SphereGeometry [18] handles geodesic arcs on thesphere using an implicit representation, which enables the use of exact rational arithmeticto handle objects of this type. The package contains implementations of basic algorithmsrelated to geodesic arcs on a sphere, such as computing the spherical convex hull, the unionof two spherical polygons, and the width of a three-dimensional set of points. It does not,however, support arrangements. A parametric surface S is defined by a continuous function f S : P → R , where the domain P = U × V is a rectangular two-dimensional parameter space with bottom, top, left, and rightboundaries, and the range S = f S ( P ). U and V are open, half-open, or closed intervals withendpoints in R ∪ {−∞ , + ∞} .We use u min , u max , v min , and v max to denote the endpoints of U and V , respectively. For example, the standard parameterization of the plane is f S ( u, v ) =( u, v, U = V = ( −∞ , + ∞ ), and the unit sphere is commonly parameterized as f S ( u, v ) = (cos u cos v, sin u cos v, sin v ), where P = [ − π, π ] × [ − π , π ].A contraction point p ∈ S is a singular point, which is the mapping of a whole boundary of0 Chapter 2. Arrangements on Surfaces the domain P . For example, if the top boundary is contracted, we have ∀ u ∈ U, f S ( u, v max ) = p ′ for some fixed point p ′ ∈ R . An identification curve C ⊂ S is a continuous curve,which is the mapping of opposite closed boundaries of the domain P . If the left and rightboundaries are identified, we have ∀ v ∈ V, f S ( u min , v ) = f S ( u max , v ), (and similarly for thebottom and top boundaries). For example, consider the sphere as parameterized above. Itscontraction points are (0 , , ± f S ( u, − π ) = (0 , , − 1) and f S ( u, π ) = (0 , , 1) for all u .Its identification curve is { f S ( π, v ) | − π ≤ v ≤ π } , as f S ( − π, v ) = f S (+ π, v ) for all v .A parameterizable curve γ is a continuous function γ : I → P where I is an open, half-open, or closed interval with endpoints 0 and 1, and γ is injective, except for at a finitenumber of points. If 0 I , lim t → γ ( t ) exists (in the closure of P ) and lies in an openside of the boundary. Similarly, if 1 I , lim t → − γ ( t ) exists and lies in an open side of theboundary. A curve C in S is the image of a curve γ in the domain.A curve is closed in the domain if γ (0) = γ (1); in particular, 0 ∈ I and 1 ∈ I . A curveis closed in the surface S (or simply closed) if f S ( γ (0)) = f S ( γ (1)). A curve γ has two ends ,the 0-end h γ, i and the 1-end h γ, i . If d ∈ I , the d -end has a geometric interpretation. Itis a point in P . If d I , the d -end has no geometric interpretation. You may think of it as apoint on an open side of the domain or an initial or terminal segment of γ . If d I , we saythat the d -end of the curve is open. Consider for example the equator curve on the sphereas parameterized above. This curve is given by γ ( t ) = ( π (2 t − , t ∈ [0 , γ is the point ( − π, 0) in P and a point on the equator of the sphere. It is closed on thesphere, but non-closed in P .A u -monotone curve is the image of a curve γ , such that if t < t , then u ( γ ( t ))
Arrangement on surface 2 package handles inducing curves that are decompos-able into parameterizable weakly u -monotone curves as defined above. Any two weakly u -monotone curves must intersect only a finite number of times or overlap only in a finitenumber of sections, if at all. The curves must be embedded on parameterizable surfacesas defined above. A curve can be unbounded or reach the boundaries of the embeddingsurface. A boundary may define a contraction point or an identification curve. We allownon-injectivity on the boundary, denoted ∂ P , and require bijectivity only in P \ ∂ P (theinterior of P ). More precisely, we require that f S ( u , v ) = f S ( u , v ) and ( u , v ) = ( u , v )imply ( u , v ) ∈ ∂ P and ( u , v ) ∈ ∂ P . Informally, we require that all geometric operationsdefined in Section 2.5 be applicable on our curves.Code reuse is maximized by generalizing the prevalent algorithms and their implementa-tions originally designed to operate on arrangements embedded in the plane. The generalizedcode handles features embedded in the modified surface e S : f e S = f S ( u, v ) | ( u, v ) ∈ P \ ∂ P defined over the interior of the parameter space, where identification curves, contractionpoints, and points at infinity are removed. Specific code that handles unbounded features u -monotone curves refer to weakly u -monotone curves hereafter. We do not support surfaces that contain a contracted identification curve. .3. The Arrangement Package Architecture The main class of the package, namely Arrangement on surface 2 , constructs and maintainsthe embedding of a set of continuous weakly u -monotone curves that are pairwise disjointin their interiors on a two-dimensional parametric surface in R . It provides the necessarycapabilities for maintaining the embedded graph, while associating geometric data with thevertices, edges, and faces of the graph. The embedded graph is represented using a doubly-connected edge list ( Dcel ) [dBvKOS00, Section 2.2], which maintains the incidence relationson its features [WFZH07b]. Each edge of the subdivision is represented by two halfedgeswith opposite orientation, and each halfedge is associated with the face to its left. It isbased on an implementation of a halfedge data-structure ( Hds ) [Ket07b] also used by thepolyhedral-surfaces package [Ket99, Ket07a].An important guideline in the design is to decouple the arrangement representation fromthe various algorithms that operate on it. Thus, the Arrangement on surface 2 class pro-vides only a restricted set of methods for locally modifying the arrangement; see Section 2.3.2.Non-trivial algorithms that involve geometric operations are implemented as free (global)functions that use the interface of the arrangement class; see Section 2.4. For example, thepackage offers free functions for incremental or aggregated insertion of curves that may notnecessarily be u -monotone, and the insertion location of which are unknown a priori . Eachinput curve is subdivides into several u -monotone subcurves before inserted using one of themember methods listed in Section 2.3.2. The Arrangement on surface 2 Dcel to the embedding modi-fied surface e S . It determines whether the embedded surface is bounded, or otherwisewhether a boundary defines a contraction point or an identification curve. If the induc-ing curves are confined to the modified parameter space, the tasks of the topology-traitsclass are minimal. However, in other cases it maintains the features that escape themodified parameter space e P .The underlying Dcel in turn associates a point with each vertex and a u -monotone curvewith each halfedge pair, where the geometric types of the point and the u -monotone curve2 Chapter 2. Arrangements on Surfaces are defined by the geometry-traits class. Users may extend the default Dcel data-structure,in order to attach additional data to the Dcel records, as explained in Section 2.3.3.The Arrangement on surface with history 2 class-template represents an arrangementof general curves embedded on a two-dimensional parametric surface, and maintains the con-struction history of the arrangement. Input curves that induce the arrangement are split into u -monotone subcurves that are pairwise disjoint in their interior, and these subcurves arethe embeddings of the arrangement halfedges. While using the Arrangement on surface 2 class we lose track of the connection between input curves and their final embeddings, inthe Arrangement on surface with history 2 data-structure each embedded u -monotonecurve is extended to store a pointer to the input curve associated with it, or a container ofcurve pointers in case the embedded u -monotone curve is an overlapping section of severalinput curves.The Arrangement on surface with history 2 class is a simple decorator for Arrangement on surface 2 . It inherits from an Arrangement on surface 2 class-templateinstantiated with a geometry-traits class that extends the u -monotone curve type. It alsostores the set of input curves, and maintains a data structure that enables efficient traversalof all halfedges induced by given input curves. The cross-pointers between input curvesand arrangement halfedges are maintained using an observer (see Section 2.4.4) that keepstrack of each change that involves an arrangement halfedge and updates its underlying curveaccordingly.Users can traverse the original curves of each arrangement halfedge, or iterate over allhalfedges induced by a given input curve. Tracing back the curve (or curves) that inducedan arrangement edge is essential in a variety of applications that use arrangements, such asrobot motion planning; see, e.g., [HH03]. It is possible, for example, to efficiently remove acurve from the arrangement by deleting all edges it induces.Arrangements embedded in the plane are very common and, at least as far as the ar-rangement package of Cgal is concerned, have a longer history than their generalizationfor two-dimensional surfaces in R . The Arrangement 2 class-template represents a planarsubdivision. It maintains the embedding of continuous weakly u -monotone curves in the xy plane, parameterized the natural way. That is, the two parameters u and v are directlymapped to x and y , respectively. Thus, u -monotonicity implies x -monotonicity and viceversa. The Arrangement 2 class-template is parameterized with a geometry-traits class andwith a Dcel data-structure. It inherits from an Arrangement on surface 2 class-templateinstantiated with the geometry traits template parameter and with a specific topology-traitsclass suitable for the plane. The dedicated topology traits is instantiated with the Dcel template parameter. Similarly the Arrangement with history 2 class-template representsa planar subdivision, and maintains the construction history of the arrangement.The package offers various operations on arrangements stored in these representations,such as point location, insertion of curves, removal of curves, and overlay computation. A decorator attaches additional responsibilities to an object dynamically [GHJV95]. .3. The Arrangement Package Architecture The interface of Arrangement on surface 2 consists of various methods that enable thetraversal of arrangement features. The class supplies iterators over its vertices, halfedges, orfaces. The classes Vertex , Halfedge , and Face , nested in the Arrangement on surface 2 class, supply in turn methods for local traversals. For example, it is possible to visit allhalfedges incident to a specific vertex. Halfedges stored in doubly-connected lists formchains. The chains define the inner and outer connected components of the boundary (CCB)of each face. A bounded face in the Arrangement 2 data structure has a single outer CCBrepresenting the outer boundary of the face, and may have several inner CCBs representingits holes. However, a face in the general Arrangement on surface 2 data structure mayhave several inner and outer CCBs; see Section 2.6. Naturally, it is possible to traverse allthe halfedges along the inner and outer boundaries of a given face. (a) (b)(c) (d)Figure 2.2: The arrangementimmediate insertion methods.The newly inserted curve isdrawn in bright blue. Verticescreated as a result of the inser-tion are drawn in bright red. Alongside with the traversal methods, the arrangementclass also supports several methods that modify the ar-rangement. The functions insert in face interior(C,f) (Figure 2.2 (a)), insert from left vertex(C,v) (Figure 2.2(b)), insert from right vertex(C,v) ) (Figure 2.2 (c)), and insert at vertices(C,v1,v2) (Figure 2.2 (d)) create an edgethat corresponds to a u -monotone curve C , whose interior isdisjoint from existing edges and vertices. The choice of whichone to use depends on whether the curve endpoints are asso-ciated with existing non-isolated arrangement vertices: (i) Ifboth curve endpoints do not correspond to any existing ver-tex, insert in face interior() is used to generate a newhole inside an existing face. (ii) If exactly one endpoint cor-responds to an existing Dcel vertex, one of the functions insert from left vertex() or insert from right vertex() is called, depending on which endpoint is associated with anexisting vertex. It forms an “antenna” emanating from an ex-isting connected component. (iii) Otherwise, both endpointscorrespond to existing vertices, and insert at vertices() iscalled to connect these vertices using a pair of twin halfedges. These functions hardly involveany geometric operations, if at all. They accept topologically related parameters, and usethem to operate directly on the Dcel records, thus saving algebraic operations, which areespecially expensive when high-degree curves are involved.Other modification methods included in the arrangement class enable users to split anedge into two, to merge two edges incident to a common vertex, and to remove an edge fromthe arrangement. It is also possible to insert a point in the interior of a given face, creatingan isolated vertex that corresponds to this point, or to remove an isolated vertex from thearrangement. Unless we force checking preconditions. In this case the precondition evaluation involves geometriccomputation. Chapter 2. Arrangements on Surfaces As mentioned above the Arrangement on surface 2 is parameterized by a topological traits,which in turn is parameterized by a Dcel class. Users may extend the default Dcel data-structure, in order to attach additional data to the Dcel records. The default Dcel modelsimply associates a point with each Dcel vertex and a u -monotone curve with each halfedgepair. Although it is possible to store auxiliary data with the curves or points by extendingtheir respective types (see Section 2.5.3), it is sometimes necessary to extend the vertex,halfedge, or face topological features of the Dcel . Many times it is desired to associateextra data with the arrangement faces only. For example, when an arrangement representsthe subdivision of a country into regions associated with their population density. In thiscase there is no alternative other than to extend the Dcel face, as there is no geometry-traitsclass entity that corresponds to an arrangement face. A similar mechanism for extendingtopological features with auxiliary attributes can be found in other components of Cgal ,such as the triangulation packages [PY07] and the polyhedral-surfaces package [Ket07a]. The Arrangement on surface 2 package offers a generic implementation of the sweep-linealgorithm [dBvKOS00, Section 2.1] in form of a class template called Sweep line 2 . Ithandles any set of arbitrary u -monotone curves, and serves as the foundation of a family ofconcrete operations, such as computing all intersection points induced by a set of curves,constructing an arrangement of curves, aggregately inserting a set of curves into an existingarrangement, and computing the overlay of two arrangements. A concrete algorithm isrealized through a sweep-line visitor, a template parameter of Sweep line 2 , which followsthe visitor design-pattern [GHJV95], and models the concept SweepLineVisitor 2 . In thiscase, a visitor defines an operation based on the sweep-line algorithm to be performed on anarrangement without the need to change the arrangement structure. Users may introducetheir own sweep based algorithms by implementing an appropriate visitor class. Another parameter of the Sweep line 2 class-template is the geometry-traits class, whichmust be instantiated with a model of the ArrangementXMonotoneTraits 2 concept; see Sec-tion 2.5 for the precise definition of this concept. It defines the minimal set of geometricprimitives, among the other, required to perform the sweep-line algorithm briefly describednext.An imaginary vertical curve is swept over the surface from left to right, transforming thestatic two-dimensional problem into a dynamic one-dimensional one. At each time duringthe sweep a subset of the input u -monotone curves intersect this vertical line in a certainorder. The subset of curves and their order along the sweep line may change as the linemoves along the u -axis, implying a change in the topology of the arrangement, only at a The Boost Graph Library, for example, uses visitors [SLL02, Section 12.3] to support user-definedextensions to its fundamental graph algorithms. .4. The Arrangement Facilities event points , namely intersection points of two curves and left endpoints orright endpoints of arcs of curves. The event points, namely endpoints and all the intersec-tion points that have already been discovered, are stored in a uv -lexicographic order in adynamic event queue, named the U -structure . The ordered sequence of segments intersect-ing the imaginary vertical line is stored in a dynamic structure called the V -structure . Bothstructures are maintained as balanced binary trees that enable their efficient maintenanceusing an advanced implementation of red-black trees [Wei05].During the sweep-line process the event objects in the U -structure are sorted lexico-graphically, and the subcurve objects are stored in the V -structure in the same order as thelexicographic order of their intersection with the imaginary sweep-line. The Sweep line 2 class performs only the operations required to maintain the U -structure and the V -structure,while the visitor class is responsible for producing the actual output of the algorithm. When-ever the sweep-line class handles an event point p , it sends a notification to its visitor. Usingthis information, the visitor can access the subcurves incident to p and the neighboringsubcurves from above and below. The map overlay of two planar subdivisions S and S is a planar subdivision S , such thatthere is a face f in S if and only if there are faces f and f in S and S respectively, suchthat f is a maximal connected subset of f ∩ f [dBvKOS00, Section 2.3]. The overlay oftwo two-dimensional subdivisions embedded on a surface is defined similarly.The overlay of two given arrangements, conveniently referred to as the “blue” and the“red” arrangements, is straightforwardly implemented using a sweep-line visitor. A consol-idated set of the “blue” and “red” curves is processed, while the imaginary vertical line isswept over the surface. The u -monotone curve type is extended with a color attribute (whosevalue is either BLUE or RED ); see Section 2.5.3. Using the extended type we filter out un-necessary computations. For example, we ignore monochromatic intersections, and computeonly red–blue intersection points (or overlaps). This way the arrangement of a consolidatedset of “blue” and “red” curves is computed efficiently.The overlay visitor needs to construct a Dcel that properly represents the overlay oftwo input arrangements, the Dcel ’s of which are potentially independently extended (seeSection 2.3.3). A face in the overlay arrangement corresponds to overlapping regions of theblue and red faces. An edge in the overlay arrangement is due to a blue edge, a red edge,or an overlap of two differently colored edges. An overlay vertex is due to a blue vertex,a red vertex, a coincidence of two differently colored vertices, or an intersection of a blueand a red curve. In each case, the data associated with the overlay Dcel feature shouldbe computed from the red and blue Dcel features that induce it. To this end, the overlayvisitor is parameterized by an overlay-traits module, which defines the merging operationsbetween various Dcel features, achieving maximum genericity and flexibility for the users.The instantiated overlay traits models the OverlayTraits concept. The concept requires theprovision of ten functions that handle all possible cases as follows:1. A new vertex v is induced by coinciding vertices v r and v b .6 Chapter 2. Arrangements on Surfaces 2. A new vertex v is induced by a vertex v r that lies on an edge e b .3. An analogous case of a vertex v b that lies on an edge e r .4. A new vertex v is induced by a vertex v r that is contained in a face f b .5. An analogous case of a vertex v b contained in a face f r .6. A new vertex v is induced by the intersection of two edges e r and e b .7. A new edge e is induced by the overlap of two edges e r and e b .8. A new edge e is induced by the an edge e r that is contained in a face f b .9. An analogous case of an edge e b contained in a face f r .10. A new face f is induced by the overlap of two faces f r and f b .We apply the overlay operations in four different ways in this thesis; see Sections 3.2.2, 3.3.2,5.3.6, and 5.3.7 for the different applications. Each application requires the provision of adifferent set of the ten functions above. The zone [Hal04] of a u -monotone curve C in an arrangement A is the set of cells of A ( C )intersected by the curve C .The Arrangement on surface 2 package includes the Arrangement zone 2 class-template,which computes the zone of an arrangement. Similar to the Sweep line 2 template, the Arrangement zone 2 template is parameterized with a zone visitor, a model of the concept ZoneVisitor 2 , and it serves as the foundation of a family of concrete operations, such as in-serting a single curve into an arrangement and determining whether a query curve intersectswith the curves of an arrangement.The zone of a curve C is computed by locating the left endpoint of C in the arrangement,and then “walking” along the curve towards the right endpoint, keeping track of the vertices,edges, and faces crossed on the way (see, for example, [dBvKOS00, Section 8.3] for thecomputation of the zone of a line in an arrangement of lines).It is sometimes necessary to compute the zone of a curve in an arrangement withoutactually inserting the curve. In other situations, the entire zone is not required, as in the caseof a process that only checks whether a query curve passes through an existing arrangementvertex; if the answer is positive, the process can terminate as soon as the vertex is located.While the sweep-line algorithm operates on a set of input u -monotone curves and its visitorscan just use the notifications they receive to construct their output structures, the zone-computation algorithm operates on an arrangement object and its visitors may modify thesame arrangement object as the computation progresses. This makes the interaction of themain class with its visitors slightly more intricate. Some arrangement-based algorithms and applications should be bound to a specific arrange-ment instance and receive notifications on various topological changes this arrangementundergoes. This is not just a convenience, but crucial to the usability of the package, as itmight be the only way for providing an algorithm with a certain input, such as data that .4. The Arrangement Facilities Arrangement on surface 2 package supports a notification mechanism, which fol-lows the observer design-pattern [GHJV95]. In this case of one-to-many dependency a set ofobserves depend on a single arrangement, so that when the arrangement changes state, allits dependents are notified and updated automatically. Using this mechanism it is possibleto attach any number of observer instances to a specific arrangement, such that all attachedobservers get notified on local and global changes the arrangement undergoes.The Arr observer Arr observer class-templateserves as a base class for other observer classes and defines a set of virtual notification func-tions, giving them all a default empty implementation. The interface of the base class isdesigned to capture all possible changes that arrangements can undergo, with a minimal setof topological events.The set of functions can be subdivided into three categories as follows:1. Notifiers of changes that affect the entire topological structure. Such changes occurwhen the arrangement is cleared or when it is assigned with the contents of anotherarrangement.2. Notifiers of a local change to the topological structure, such as the creation of a newvertex or an edge, the splitting of an edge or a face, the formation of a new hole insidea face, the removal of an edge, etc.3. Notifiers of a global change initiated by a free (global) function, and called by thefree function (e.g., incremental or aggregate insert; see Section 2.3). This categoryconsists of a single pair of notifiers, neither of them is called by methods of the Arrangement on surface 2 class-template itself. Issuing point-location queries (orany other queries for that matter) between the calls to the “before” and “after” func-tions of this pair is forbidden. See [WF05] for a detailed specification of the arrangement observer class sketched above.Each arrangement object stores a list of pointers to Arr observer objects, and wheneverone of the structural changes listed in the first two categories above is about to take place,the arrangement object invokes the appropriate function of each of its observers. It also doesso immediately after the change has taken place. In addition, a free function may choose totrigger a similar notification, which falls under the third category above. This constraint improves the efficiency of the maintenance of auxiliary data structures for the relevantpoint-location strategies, which have to update their data structures according to the changes the arrange-ment undergoes (see Section 2.4.5 for more details). Since no point-location queries are issued betweenthe invocation of before global change() and after global change() , it is not necessary to perform anupdate each time a local topological change occurs, and it is possible to postpone the updates until after theglobal operation is completed. Chapter 2. Arrangements on Surfaces In case the new observer is attached to a non-empty arrangement, its constructor mayextract the relevant data from the non-empty arrangement using various traversal methodsoffered by the public interface of the Arrangement on surface 2 class, and update anyinternal data stored in the observer. This is necessary, for example, in case of the non-stateless point-location strategies, as shown in the next section. Point location is defined as follows: Given a point, find the arrangement cell that contains it.The Arrangement on surface 2 package provides the means to answer this query. Typically,the result of the point-location query is one of the arrangement faces, but in degeneratesituations the query point can lie on an edge, or it may coincide with a vertex. Sincethe arrangement representation is decoupled from the algorithms that operate on it, the Arrangement on surface 2 class does not support point-location queries directly. Instead,the package provides a set of classes that are capable of answering such queries, all are modelsof the concept ArrangementPointLocation 2 . Each model employs a different algorithm or strategy for answering queries. A model of this concept must define the locate() function,which accepts an input query point and returns an object representing the arrangementcell that contains this point (a polymorphic CGAL::Object instance that can either be a Face handle , a Halfedge handle , or a Vertex handle ).The following models of the concept ArrangementPointLocation 2 are included in the Arrangement on surface 2 package. • Arr naive point location locates the query point naively, by exhaustively scanningall arrangement cells. It is the only strategy with unlimited support; see Section 2.6. • Arr walk along a line point location simulates a reverse traversal along an imag-inary vertical ray emanating from the query point toward infinity. It starts from theunbounded face of the arrangement and moves downward toward the query point untilit locates the arrangement cell containing it. • Arr landmarks point location ArrangementLandmarksTraits 2 concept, which adds two requirements to the basic Arrange-mentBasicTraits 2 concept: (i) Approximating the coordinates of a given point p using the double-precisionarithmetic, and (ii) constructing a u -monotone curve that connects two given points p and q , where p rep-resents a landmark point and q is the query point. Most traits classes included in the arrangement packageare models of this refined concept. .5. Geometry-Traits Concepts Dcel that represents the arrangement. For a complete survey see [HH08].Each of the “landmarks” point-location class and the trapezoidal point-location class usesan observer to receive notifications whenever the arrangement is modified. For example, thedefault generator employed by the “landmarks” strategy uses the arrangement vertices aslandmarks, so whenever a new vertex is created (by the insertion of a new edge, by thesplitting of an existing edge, or by the insertion of an isolated point), it should be insertedinto the nearest-neighbor search structure maintained by the respective landmark class. Theusage of the notification mechanism makes it possible to associate several point-locationobjects with the same arrangement simultaneously.The “landmarks” and the trapezoidal point-location strategies are both characterized byvery efficient query time at the cost of time-consuming preprocessing. Naturally, these strate-gies exhibit better overall performance when the number of arrangement updates is relativelysmall compared to the number of issued queries. For a report on extensive experiments withthe various point-location strategies see [HH08]. ArrangementBasicTraits 2ArrangementXMonotoneTraits 2ArrangementLandmarksBasicTraits 2ArrangementTraits 2ArrangementLandmarksXMonotoneTraits 2ArrangementLandmarksTraits 2 Figure 2.3: Refinement hierarchy of geometry-traits concepts. The implementations of thevarious algorithms that con-struct and manipulate ar-rangements are generic, asthey are independent of thetype of curves they handle.All steps of the algorithmsare enabled by the minimalset of geometric primitivesgathered in the geometry-traits class, a model of ageometry-traits concept. Thegeometry-traits concept isfactored into a hierarchy ofrefined concepts. The refine-ment hierarchy is defined according to the identified minimal requirements imposed by dif-ferent algorithms that operate on arrangements, thus alleviating the production of traitsclasses that handle complicated curves, and increasing the usability of the algorithms. Therequirements listed by the geometry-traits concepts include only the utterly essential typesand operations, and fully specify all the preconditions that the input must satisfy, as thesemay simplify the implementation of models of this concept even further.0 Chapter 2. Arrangements on Surfaces The following sections are dedicated to a detailed description of the hierarchy. We listthe minimal requirements of each layer in the hierarchy, and provide formal definitions forthe required operations. The letters x and y are used in the code to refer to the two surfaceparameters, as arrangements embedded in the xy -plane are more common and familiar.The names of required nested types (e.g., X monotone curve 2 ) and valid expressions (e.g., compare x ) are faithful to the original source code. However, we use the letters u and v in the formal definitions below to refer to the two surface parameters, as these definitionsapply to the general case of arrangements embedded on surfaces. Let cmp u () and cmp v ()denote two predicates that accept two points and compare them by their u -coordinates andby their v -coordinates respectively. We use the following notation. For a point p , ( u p , v p ),denotes a pre-image, and for a curve C , γ denotes a pre-image, that is, p = f S ( u p , v p ) and C ( t ) = f S ( γ ( t )) for all t ∈ I .The basic concept ArrangementBasicTraits 2 requires the definition of the types Point 2 and X monotone curve 2 . The latter represents a u -monotone curve, and the former is thetype of the endpoints of the curves, representing a point on the surface. This conceptlists the minimal set of predicates on objects of these two types sufficient to enable theoperations provided by the Arrangement on surface 2 class-template itself, namely theinsertion of bounded u -monotone curves that are interior disjoint from any vertex and edgein the arrangement. All points and curves in the set below are required to have an inversepre-image in P \ ∂ P . In particular all curves are u -monotone. Compare x 2 : Compare two points by their u -coordinates. Compare xy 2 : Compare two points lexicographically by their u and then by their v -coordi-nates. Construct min 2 : Return the lexicographically smaller (left) endpoint of a given curve. Construct max 2 : Return the lexicographically larger (right) endpoint of a given curve. Is vertical 2 : Determine whether a weakly u -monotone curve is vertical. Compare y at x : Given a point p and a curve C , such that the u p lies in the u -range of C ,determine whether p is above, below, or lies on C . More precisely, if C is vertical,determine whether p lies on C , or above or below C . Otherwise, since u ( γ (0)) ≤ u p ≤ u ( γ (1)) must hold and C is u -monotone, there must be a unique 0 ≤ t ′ ≤ 1, thatsatisfies u ( γ ( t ′ )) = u p . Return cmp v ( p, γ ( t ′ )). Compare y at x right : Given two curves C and C that share a common left endpoint p ,determine the relative position of the two curves immediately to the right of p . Moreprecisely, return cmp v ( γ ( ǫ ) , γ ( ǫ )), where ǫ , ǫ > Compare y at x left : Given two curves C and C that share a common right endpoint p ,determine the relative position of the two curves immediately to the left of p . Moreprecisely, return cmp v ( γ (1 − ǫ ) , γ (1 − ǫ )), where ǫ , ǫ > optional requirement with ramifications in case it is not fulfilled; see Sec-tion 2.5.1.The set of predicates listed above is also sufficient for answering point-location queries by thevarious point-location strategies, with a small exception of the “landmarks” strategy, whichrequires a traits class that models the refined concept ArrangementLandmarksTraits 2 . This .5. Geometry-Traits Concepts u -monotone curves that may intersect in theirinterior, requires an arrangement instantiated with a traits class that models the concept ArrangementXMonotoneTraits 2 . This concept refines the basic arrangement-traits conceptdescribed above, as it requires an additional method for computing intersections between u -monotone curves, among the other. An intersection point between two curves is also rep-resented by the Point 2 type. The refined traits concept also requires methods for splittingcurves at these intersection points to obtain pairs of interior disjoint subcurves and mergingpairs of subcurves. In summary, a model of the refined concept must provide the additionaloperations bellow. All points and curves in the set below are required to have an inversepre-image in P \ ∂ P . In particular all curves are u -monotone. Intersection 2 : Compute the intersections between two given curves C and C . Split 2 : Split a given curve C at a given point p , which lies in the interior of C , into twointerior disjoint subcurves. Merge 2 : Merge two mergeable curves C and C into a single curve C . Is mergeable 2 : Determine whether two curves C and C that share a common endpointcan be merged into a single continuous curve representable by the traits class.The further refined concept ArrangementTraits 2 enables the construction of arrange-ments induced by general curves. A model of the refined concept must define a third typethat represents a general (not necessarily u -monotone) curve, named Curve 2 . It also hasto supply a method that subdivides a given curve into simple u -monotone subcurves, andpossibly isolated points. We refer to the entire hierarchy of refinements defined above asa single “abstract” concept called NoBoundaryTraits , as it represents concepts the modelsof which handle curves that must have inverse pre-images in P \ ∂ P . We use this abstractconcept to simplify the description of the hierarchy defined below. NoBoundaryTraitsHasBoundaryTraitsBoundedBoundaryTraitsUnboundedBoundaryTraitsAllBoundaryTraits Figure 2.4: Abstract refinement hierarchy of geometry-traits concepts for arrangement on surfaces. The package introduces addi-tional concepts, models of which areable to handle unbounded curves orcurves that reach the boundaries, theendpoints of which coincide with con-traction points or lie on identifica-tion curves; see Figure 2.4. The“abstract” HasBoundaryTraits sub-hierarchy lists additional predicatesrequired to handle both curves thatreach or approach the boundaries ofthe parameter space. It has no mod-els. The refined BoundedBoundaryTraits and UnboundedBoundaryTraits sub-hierarchies listadditional predicates required to handle bounded and unbounded curves, respectively. The For example, the curve ( x + y )( x + y − 1) = 0 is comprised of two u -monotone circular arcs, whichtogether form the unit circle, and a singular isolated point at the origin. Chapter 2. Arrangements on Surfaces geometry-traits class that handles arcs of great circles models the BoundedBoundaryTraits concept, as the parameter space is bounded in all four directions. Finally, the AllBoundary-Traits sub-hierarchy refines all the above. A model of this concept can handle unboundedcurves in some directions and bounded curves in others.In the rest of this section all curves are required to be u -monotone. The HasBoundary-Traits concept requires the following additional operations: Parameter space in x 2 : Given a curve C and an index d ∈ { , } that identifies one of itsends, determine the location of its pre-image in the domain P along the u dimension.More precisely, determine whether u ( γ ( d )) is equal to u min , u max , or falls in between.In case of an unbounded curve, determine whether lim t → d u ( γ ( t )) is equal to −∞ or+ ∞ . Parameter space in y 2 : Given a curve C and an index d ∈ { , } that identifies one of itsends, determine the location of its pre-image in the domain P along the v dimension.More precisely, determine whether v ( γ ( d )) is equal to v min , v max , or falls in between.In case of an unbounded curve, determine whether lim t → d v ( γ ( t )) is equal to −∞ or+ ∞ . Compare x near boundary 2 : There are two predicates:1. Given a point p , the inverse of which is in P \ ∂ P , a curve C , and an index d ∈ { , } that identifies an end of C , compare the u coordinates of p and a pointalong C near its given end. More precisely, return cmp u ( p, γ ( | d − ǫ | )), where ǫ > C and C and two corresponding indices d , d ∈ { , } thatidentify two ends of C and C respectively, compare the u coordinates of twopoints along C and C respectively near their given ends. More precisely, returncmp u ( γ ( | d − ǫ | ) , γ ( | d − ǫ | )), where ǫ , ǫ > Compare y near boundary 2 : Given two curves C and C , and a single index d ∈ { , } thatidentifies two ends of C and C , compare the v coordinates of two points along C and C respectively near the given ends. More precisely, return cmp v ( γ ( | d − ǫ | ) , γ ( | d − ǫ | )), where ǫ , ǫ > UnboundedBoundaryTraits concept requires the following additional operations: Is bounded 2 : Given a curve C and an index d ∈ { , } that identifies an end of C , determinewhether the curve end is bounded.The BoundedBoundaryTraits concept requires the following additional operations: Is on x identification 2 : This predicate applies only to a parameterization that has avertical identification curve. Given a point p (respectively a curve C ), determinewhether p (respectively C ) lies on the vertical and identified sides of the boundary.More precisely, determine whether u p ∈ { u min , u max } . (Respectively, determine whether .5. Geometry-Traits Concepts u ( γ ( t )) ∈ { u min , u max } , ∀ t ∈ [0 , Is on y identification 2 : This predicate applies only to a parameterization that has ahorizontal identification curve. Given a point p (respectively a curve C ), determinewhether p (respectively C ) lies on the horizontal and identified sides of the boundary.More precisely, determine whether v p ∈ { v min , v max } for all pre-images of p . (Respec-tively, determine whether v ( γ ( t )) ∈ { v min , v max } , ∀ t ∈ [0 , Is on x contraction 2 : This predicate applies only to a parameterization that has a con-tracted vertical boundary. determine whether p coincides with a contraction point.More precisely, determine whether u p is equal to u min or u max . Is on y contraction 2 : This predicate applies only to a parameterization that has a con-tracted horizontal boundary. determine whether p coincides with a contraction point.More precisely, determine whether v p is equal to v min or v max . Compare x on identification 2 : This predicate applies only to a parameterization thathas a horizontal identified sides of the boundary. Given two points p and p that lieon the horizontal identification arc, compare their u -coordinates. Compare y on identification 2 : This predicate applies only to a parameterization thathas a vertical identified sides of the boundary.. Given two points p and p that lie onthe vertical identification arc, compare their v -coordinates.All traits-class operations are implemented as function objects ( functors ) according to Cgal ’s guidelines. This allows extending the geometric types above, without the need toredefine the methods that operate on them; see [HHK + 07] for details on the extensible kernel.For a detailed specification of the various concept requirements see [WF05]. The geometry-traits adaptor class-template implements geometric operations that are notprovided by a model of the geometry-traits concept itself, using the operations supplied bya model of the geometry-traits concept as basic blocks. It decreases the effort required todevelop geometry-traits models, and at the same time increases the usability of the geometry-traits models, adapting them for extended uses. A geometry-traits type is injected as atemplate parameter into the adaptor class, which inherits from it, centralizing all geometricoperations. In cases where the efficiency of methods is crucial, a developer has a way tooverride these methods with optimized ones.For example, in order to determine whether a point p is in the u -range of a u -monotonecurve C , the adaptor simply compares p to the endpoints of C . It checks whether p lies tothe right of the left endpoint and to the left of the right endpoint.In some cases, the geometry-traits adaptor class uses a tag-dispatching mechanism toselect the appropriate implementation of a geometry-traits class operation. Tag dispatchingis a technique that uses function overloading to dispatch a function at compile time , basedon properties of the types of the arguments the function accepts [4]. This mechanism enablesusers to implement their traits class with a reduced or alternative set of operations. Theadaptor respects the tags listed below every geometry-traits class must define.4 Chapter 2. Arrangements on Surfaces Has left category : A Boolean tag that indicates whether the traits class provides thepredicate compare y at x left , which compares two u -monotone curves to the left of a common right endpoint. This predicate is required only by some point-locationstrategies and by the zone-computation algorithm. While in some cases it is fairly easyfor the traits-class implementer to provide it, in other cases it can be rather difficult,or even quite impossible. When this tag is false, the traits-class adaptor resorts to asomewhat less efficient algorithm that uses (other) existing traits-class predicates. Has merge category A Boolean tag that indicates whether a model of the ArrangementX-MonotoneTraits 2 supports the merge of u -monotone curves. If the tag is true, thetraits class must provide the two operations merge 2 and is mergeable 2 . The mergeroperation is used to eliminate redundant features in the arrangement. For example,if we have a T-shaped structure formed by two line segments, and the vertical seg-ment forming the “leg” is removed, then it is possible to merge the two horizontalsub-segments. When the has-merge tag is false, the adaptor simply declares any pairof curves as non-mergeable. The only effect on the arrangement is that we cannotremove redundant vertices (of degree two) following the deletion of edges. Boundary category : A quadruple tag that categorizes the traits class according to the hi-erarchy described in Figure 2.4. The adaptor provides empty implementations of theoperations that are never invoked, yet required for smooth compilation. Table 2.1: Geometry-traits models Curve Family Degree Surface Boundness Arithmetic linear segment 1 plane bounded rationallinear segments, rays, lines 1 plane unbounded rationalpiecewise linear curves ∞ plane bounded rationalcircular arcs, linear segments ≤ ≤ ≤ ≤ ≤ n plane unbounded algebraicplanar B´ezier curves ≤ n plane unbounded algebraicunivariate polynomials ≤ n plane unbounded algebraicgeodesic arcs on sphere ≤ ≤ ≤ n Dupin cyclides bounded algebraicThe large number of geometry-traits models already implemented enables the construc-tion and maintenance of arrangements induced by many different types of curves. Thepackage itself contains several models of the geometry-traits concept. A few other modelshave been developed by other groups of researchers. Models are distinguished not only bythe different families of curve they handle, but also by their suitability for constructing and .5. Geometry-Traits Concepts + + + 02, EKP + + + + 07, BFH + + Traits-class decorators are used to extend the geometric entities defined by the traits classwith additional, possibly non-geometric, data. An alternative way to achieve this is to extendthe geometric types of the kernel, as the kernel is fully adaptable and extensible [HHK + u -monotonecurve types with distinct types of data, and there is a convenient one, where the data attachedto the u -monotone curve type is a set of objects, the type of which is attached to the (general)curve type. This set usually contains a single data object, unless the u -monotone curvecorresponds to an overlapping section of two curves or more. When a curve with a data field d is split into u -monotone subcurves, each subcurve is associated with a singleton set { d } .When two u -monotone curves overlap, the decorator takes the union of their data sets, andassociates it with the resulting overlapping subcurve. Polylines are of particular interest, as they can be used to approximate more complex curves.At the same time handling them is easier than handling higher-degree algebraic curves, asrational arithmetic is sufficient to carry out exact computations on polylines.The geometry-traits model that handles polylines is a class-template called Arr polyline traits 2 . It must be instantiated with a geometry-traits class that modelsthe concept ArrangementLinearTraits . This concept refines the ArrangementTraits concept,as it adds a variant — it must handle line segments. This variant cannot be enforced by the6 Chapter 2. Arrangements on Surfaces compiler, but rather be verified at run time. A polyline curve is represented as a vector of SegmentTraits::X monotone curve 2 objects (namely segments). The polyline-traits classdoes not perform any geometric operations directly. Instead, it solely relies on the function-ality of the instantiated segment-traits class. For example, when we need to determine theposition of a point with respect to a u -monotone polyline, we use binary search to locatethe relevant segment that contains the point in its u -range, then we compute the position ofthe point with respect to this segment. Thus, operations on u -monotone polylines of size m typically take O (log m ) time.Users are free to choose the underlying segment-traits class based on the number ofexpected intersection points (see discussion above in Section 2.5.2). Moreover, it is possibleto instantiate the polyline-traits class-template with a traits class that handles segmentswith some additional data attached to them (see Section 2.5.3). This makes it possible toassociate different data objects with the different segments that compose a polyline. In this section we concentrate on the particular category of arrangements embedded onsurfaces, where the embedding space is the sphere, and the inducing objects are geodesic arcs.There is an analogy between this class of arrangements and the class of planar arrangementsinduced by linear curves (i.e., segments, rays, and lines), as properties of linear curves in theplane can be often, but not always, adapted to geodesic arcs on the sphere.An arrangement of geodesic arcs embedded on the sphere is defined as an instance ofthe Arrangement on surface 2 class-template instantiated with appropriate geometry- andtopology-traits classes, namely Arr geodesic arc on sphere traits 2 and Arr spherical topology traits 2 , respectively. The geometry-traits class is tailored tohandle geodesic arcs as efficiently as possible, and defines the parameterization used: P =[ − π + α, π + α ] × [ − π , π ], f S ( u, v ) = (cos u cos v, sin u cos v, sin v ), where α is a variablethat must be set at compile time, and is by default 0. This parameterization induces twocontraction points p s = (0 , , − 1) and p n = (0 , , Dcel to have a single face, theembedding of which, is the entire sphere. It is designed to retain the variant that this facealways contains the north pole during modifications the arrangement may undergo. Thetopology-traits class is required, for example, to inform its users that the top and bottomboundaries of the parameter space are contracted and the left and right boundaries are iden-tified. It maintains a search structure of vertices that coincide with the contraction pointsor lie on the identification arc. .6. Arrangements of Geodesic Arcs on the Sphere The geometry-traits class for geodesic arcs on the sphere is parameterized with a geometrickernel [HHK + 07] that encapsulates the number type used to represent coordinates of geo-metric objects and to carry out algebraic operations on those objects. The implementationhandles all degeneracies, and is exact, as long as the underlying number type is rational,even though the embedding surface is a sphere. We are able to use high-performance ker-nel models instantiated with exact rational number-types for the implementation of thisgeometry-traits class, as exact rational arithmetic suffices to carry out all necessary alge-braic operations. The ability to robustly construct arrangements of geodesic arcs on thesphere, and robustly apply operations on them using only (exact) rational arithmetic is akey property that enables an efficient implementation.A point in our arrangement is defined to be an unnormalized vector in R , representingthe place where the ray emanating from the origin in the relevant direction pierces the sphere.An arc of a great circle is represented by its two endpoints, and by the plane that containsthe endpoint vectord and goes through the origin. The orientation of the plane and thesource and target points determine which one of the two great arcs is considered.The point type is extended with an enumeration that indicates whether the vector (i)pierces the south pole, (ii) pierces the north pole, (iii) intersects the identification arc, or (iv)is in any other direction. An arc of a great circle is extended with three Boolean flags thatindicate whether any one of the x, y, z coordinates of the normal of the plane that defines thearc vanishes. These flags are used to minimize the number of invocations of the geometry-traits operations, which has a drastic effect on the performance of arrangement operationsat the account of a slight increase in space consumption. This representation enables anexact yet efficient implementation of all geometric operations required by the geometry-traitsconcept using exact rational arithmetic, as normalizing vectors (that represent directions andplane normals) is completely avoided.We describe in details four predicates, namely Compare x 2 , Compare xy 2 , Compare x near boundary 2 , and Compare y near boundary 2 ; see Section 2.5 for the com-plete set of the concept requirements. The former compares two points p and p by their u -coordinates. The concept admits the assumption that the input points do not coincidewith the contraction points and do not lie on the identification arc. Recall that points arein fact unnormalized vectors that represent directions in R . We project p and p ontothe xy -plane to obtain two-dimensional unnormalized vectors ˆ p and ˆ p , respectively. Wecompute the intersection between the identification arc and the xy -plane to obtain a thirdtwo-dimensional unnormalized vector ˆ d . Finally, we test whether ˆ d is reached strictly before8 Chapter 2. Arrangements on Surfaces xy ˆ d ˆ p ˆ p ˆ p is reached, while rotating counterclockwise starting at ˆ p . Thisgeometric operation is supported by every geometric kernel of Cgal .In the figure on the right ˆ d is reached strictly before ˆ p is reached.Therefore, the u -coordinate of p is larger than the u -coordinate of p . The predicate Compare xy 2 compares two points p and p lexico-graphically. It first applies Compare x 2 to compare the u -coordinatesof the two points. If the u -coordinates are equal, it applies a pred-icate that compares the v -coordinates of two points with identical u -coordinates, referred to as Compare y 2 . This predicate first com-pares the signs of the z -coordinates of the two unnormalized inputvectors. If they differ, it concludes that the point with the positive z -coordinate has a v -coordinate that is larger than the v -coordinateof the point with the negative z -coordinate. If the signs are identical,it compares the squares of their normalized z -coordinates, essentiallyavoiding the square-root operation. If the sign of the (unnormalized) z -coordinates of bothpoints is positive (resp. negative), the point with the larger (resp. smaller) square of nor-malized z -coordinate has a larger v -coordinate.The predicates above accept points, the pre-images of which, liein the interior of the parameter space. However, there is also aneed to lexicographically compare the ends of arcs, the pre-imagesof which reach the boundary of the parameter space. The predicate Compare x near boundary 2 accepts either (i) a point, the pre-imageof which lies in the interior of the parameter space, and an arc end,or (ii) two arc ends. Such an arc end is provided by a vertical arcand an index that identifies one of the two ends of the arc, and mustcoincide with one of the contraction points. The first variant compares the u -coordinatesof the input point and a point along the input arc near its given end, whereas the secondvariant compares the u -coordinates of two points along the input arcs near their respectivegiven ends. Recall, that the u -coordinates of all points along a vertical arc are the same ( C and C in the figure above). Thus, we can compare the u -coordinate of an arbitrary pointon a vertical arc that lie inside the parameter space. For example, for the second case, wecompare the two vectors perpendicular to the normals to the planes that define the verticalarcs, respectively, e.g., the u -coordinate of a point on the arc C near its top end is smallerthan the u -coordinate of a point on the arc C near its top end, and in particular it is smallerthan the u -coordinate of the bottom end of C . The Compare y near boundary 2 predicatecompares the v -coordinate of two arcs ends, the pre-images of which lie on the same (left orright) identified side of the boundary of the parameter space. We use the aforementioned Compare y 2 predicate to compare the end points. If the points are equal, we compare thenormals to the plane that define the arcs. In our example, the left end of C is smaller thanthe left end of C , which is smaller than the left end of C .All the required geometric types listed in the traits concept are maintained using onlyrational numbers. All required geometric operations are implemented using only rational .7. Applications Degeneracies, such as overlapping arcs that occur during intersection compu-tation, are properly handled. The end result is a robust yet efficient implementation. Arrangement on surfaces have many applications this thesis falls short to list. However,we do list a few samples we were involved (or remotely involved) with the implementationof which, i.e., Regularized Boolean Set-Operations, Envelopes of Surfaces, and Voronoi di-agrams. Minkowski sum construction is covered in details in the following chapter. TheBoolean set-operation results, minimization diagrams, maximization diagrams and Voronoidiagrams, (see Section 2.7.2 for definitions), and Minkowski-sums are all represented as ar-rangements, and as such can be passed as input to consecutive operations on arrangementssupported by the Arrangement on surface 2 package and its derivatives. Together with R. Wein and B. Zuckerman we have developed a package that supports Booleanset-operations on point sets bounded by u -monotone curves embedded on two-dimensionalparametric surfaces in R [FWZH07]. In particular, it contains the implementation of reg-ularized Boolean set-operations, intersection predicates, and point containment predicates.A regularized Boolean set-operation op ∗ can be obtained by first taking the interior of theresulting point-set of an ordinary Boolean set-operation ( P op Q ) and then by taking theclosure [Hof04]. That is, P op ∗ Q = closure(interior( P op Q )). Regularized Boolean set-operations appear in constructive solid geometry (CSG), because regular sets are closed underregularized Boolean set-operations, and because regularization eliminates lower dimensionalfeatures, namely isolated vertices and “antennas” (namely, dangling edges), thus simplifyingand restricting the representation to physically meaningful solids. Ordinary Boolean set-operations, which distinguish between the interior and the boundary of a polygon, are notimplemented within this package. However, we implemented a specialized ordinary unionoperation as part of an assembly partitioning application; see Chapter 5. ArrangementBasicTraits 2ArrangementDirectionalXMonotoneTraits 2GeneralPolygonSetTraits 2 Figure 2.5: Refinement hierarchy of geometrytraits concepts for Boolean set-operations. The operands and results of the regularizedoperations are general polygons that may haveholes. The boundaries of a general polygonand of holes, if present, are general u -monotonecurves. The Arrangement on surface 2 class isemployed to represent a point set embedded ona two-dimensional parametric surface as an ar-rangement. A point set is typically constructedfrom a single general polygon or a collection ofinterior disjoint general polygons. The underly-ing arrangement must be instantiated with a geometry traits that models the concept Gener- Points are represented as unnormalized vectors; The coordinates of such points are converted into ma-chine floating-point only for rendering purposes. Chapter 2. Arrangements on SurfacesalPolygonSetTraits 2 . This concept refines the concept ArrangementDirectionalXMonotone-Traits 2 , which in turns refines the concept ArrangementBasicTraits 2 (see Section 2.5).The ArrangementDirectionalXMonotoneTraits 2 concept treats its u -monotone curvesas objects directed from one endpoint appointed to be the source to the other endpointappointed to be the target. Thus, it requires few additional operations on u -monotonecurves: Compare endpoints xy 2 : Given a u -monotone curve C , compare the source and targetpoints of C lexicographically. Construct opposite 2 : Given a u -monotone curve C , construct the opposite curve of C (namely, swap the source and target endpoints of C ). Intersection 2 : Compute the intersections between two given curves C and C . Merge 2 : Merge two mergeable curves C and C into a single curve C . Is mergeable 2 : Determine whether two curves C and C that share a common endpointcan be merged into a single continuous curve representable by the traits class.Most traits classes bundled in the Arrangement on surface 2 package and distributed with Cgal , are models of the concept ArrangementDirectionalXMonotoneTraits 2 . The GeneralPolygonSetTraits 2 concept requires its models to define a type that rep-resents a general polygon and another one that represents general polygon with holes inaddition to the Point 2 and X monotone curve 2 types that must be defined by all modelsof the generalized concept. It also requires the provision of several operations that operateon these two types listed below. Construct polygon 2 : Given a sequence C of u -monotone curves, construct a general poly-gon that has C as its outer boundary. Construct curves 2 : Given a general polygon P , obtain the sequence of u -monotone curvesthat comprise the boundary of P . Construct general polygon with holes 2 : Given a general polygon P and a (possiblyempty) set of holes H , construct a general polygon with holes that has P as its outerboundary and H as its holes. Construct outer boundary : Given a general polygon-with-holes P , obtain the general poly-gon that is its outer boundary. Construct holes : Given a general polygon-with-holes P , obtain the holes of P if any. Is unbounded : Given a general polygon-with-holes P , determine whether it has an outerboundary. The Arr polyline traits 2 traits class is not a model of the ArrangementDirectionalXMonotone-Traits 2 concept, as the u -monotone curve it defines is always directed from left to right. Thus, an oppositecurve cannot be constructed. .7. Applications Lower envelopes of functions on parametric surfaces are defined in a way similar to thestandard definition of lower envelopes of bivariate functions in space [Hal04]. Let S be atwo-dimensional parametric surface is R . Given a set of bivariate functions F = { f , . . . , f n } ,where f i : S → R , their lower envelope Ψ( u, v ) is defined to be their point-wise minimumΨ( u, v ) = min ≤ i ≤ n f i ( u, v ). The minimization diagram M ( F ) of the set F is the two-dimensional map obtained by the projection of the lower envelope onto S . Upper envelopesand maximization diagrams are defined analogously for the point-wise maximum of thefunctions. (a) (b) (c)Figure 2.6: Lower envelopes of various types of surfaces (a) The lower envelope of triangles. (b) Thelower envelope of spheres. (c) The lower envelope of planes. (The minimization diagrams is drawnabove the planes for clarity.) The Envelope 3 package of Cgal [Mey06, MWZ07] computes the lower (or the upper)envelope of a set of surfaces in R . It is based on the Arrangement on surface 2 package,and like the base package, it handles degenerate input, and produces exact results. An ar-rangement data-structure is used to represent the resulting minimization diagram [Hal04].The envelope computation is enabled by a traits class — a model of the concept Envelope-Traits 3 , which refines the ArrangementTraits 2 concept. The Envelope 3 package currentlycontains three models of the EnvelopeTraits 3 concept that can be used to compute the en-velope of triangles, spheres, and planes, respectively; see Figure 2.6 for an illustration. Othermodels of the EnvelopeTraits 3 concept have been developed, for example, a traits class thatenables the computation of the envelope of a set of quadric surfaces [BM07]. Voronoi diagrams were thoroughly investigated, and were used to solve many geometricproblems, since introduced by Shamos and Hoey to the field of computer science [SH75](although their origin dates back centuries ago; see [OBSC00]). The concept of computingcells of points that are closer to a certain object than to any other object, among a finitenumber of objects, was extended to various kinds of geometric sites, ambient spaces, anddistance functions, e.g., power diagrams of circles in the plane, multiplicatively weighted2 Chapter 2. Arrangements on Surfaces Voronoi diagrams, additively weighted Voronoi diagrams [AK00]. This space decompositionis strongly connected to arrangements [ES86], a property that yields a very general approachfor computing Voronoi diagrams.Given a set of n points P = { p , . . . , p n } , p i ∈ S , we define R ( P, p i ) = { x ∈ S | ρ ( x, p i ) <ρ ( x, p j ) , j = i } , where ρ ( p i , p j ) is some given distance function. R ( P, p i ) is the regionof all points that are closer to p i then to any other point in P . The Voronoi diagram of P over S is defined to be the regions R ( P, p ) , R ( P, p ) , . . . , R ( P, p n ) and their boundaries.Edelsbrunner and Seidel [ES86] observed the connection between Voronoi diagrams in R d and lower envelopes in R d +1 of the corresponding distance functions to the sites. From theabove definitions it is clear that if f i : S → R is set to be f i ( x ) = ρ ( x, p i ), for i = 1 , . . . , n ,then the minimization diagram of { f , . . . , f n } over S is exactly the Voronoi diagram of P over S . (a) (b)Figure 2.7: Voronoi diagrams on the sphere.(a) The Voronoi diagram of 14 random points.(b) The power diagram of 10 random circles. K. E. Hoff et al. developed a technique forcomputing generalized 2D and 3D Voronoi dia-grams using interpolation-based polygon raster-ization hardware [HKL + et al. developeda new framework to compute different types ofVoronoi diagrams embedded on certain paramet-ric surfaces in an exact manner [SSH08]. Theframework is based on the exact computation ofthe lower envelope of the site-distance functionsover the surface [Mey06]. It provides a reducedand convenient interface between the construction of the diagrams and the construction ofenvelopes, which in turn are computed using the Envelope 3 package [MWZ07]. Obtaininga new type of Voronoi diagrams only amounts to the provision of a geometry-traits classthat handles the type of bisector curves of the new diagram type. Essentially, every type ofVoronoi diagram, the bisectors of which can be handled by a geometry traits class, can beimplemented using this framework. In particular the geometry-traits class for geodesic arcsembedded on the sphere enables Voronoi diagrams of points on the sphere and their gener-alization, power diagrams, also known as Laguerre Voronoi diagrams, on the sphere, as thebisector curves between point sites on the sphere are great circles [NLC02, OBSC00], and soare the bisectors between circle sites on the sphere under the Laguerre distance [Sug02]; seeFigure 2.7. Power diagrams on the sphere have several applications similar to the applica-tions of power diagrams in the plane. For example, determining whether a point is includedin the union of circles on the sphere, and finding the boundary of the union of circles on thesphere [IIM85, Sug02]. In certain cases, the distance to a site may depend on various parameters associated with the site, e.g.,in the cases of M¨obius diagrams or anisotropic diagrams. .7. Applications v0 pi/2−piu0pi We implicitly construct envelopes of distance functions definedover the sphere to compute Voronoi diagrams. The image to theright illustrates the distance function from (0 , ∈ [ − π, π ] × [ − π , π ]on the sphere in the parameter space. The great circle bisector oftwo point sites on the sphere is the intersection of the sphere andthe bisector plane of the points in R (imposed by the Euclideanmetric).If a point on the sphere is given as a general vector in R , it must be normalized. Ifthe normalization results with a point with irrational coordinates, then it must be approx-imated to a point that lies exactly on the sphere [CDR92]. Once approximated, all furthercomputations are carried out using exact rational arithmetic.(a) (b)Figure 2.8: Arrangements on the sphere. Figure 2.8 (a) on the left shows an arrange-ment on the sphere induced by (i) the conti-nents and some of the islands on earth, and (ii)the institutions that hosted SoCG during thismillennium, which appear as isolated vertices.The sphere is oriented such that College Park,MD, USA is at the center. The arrangementconsists of 1054 vertices, 1081 edges, and 117faces. The data was taken from gnuplot [13] andgoogle maps [14]. Figure 2.8 (b) shows an ar-rangement that represents the Voronoi diagramof the nine cities, the institutions above are lo-cated at, namely College Park, Gyeongju, Sedona, Pisa, New York, San Diego, Barcelona,Medford, and Hong Kong.As mention above Voronoi diagrams, among the other, are repre-sented as arrangements and can be passed as input to consecutive oper-ations on arrangements supported by the Arrangement on surface 2 package and its derivatives. The figure on the right shows the overlayof the two arrangements shown in Figure 2.8.4 Chapter 2. Arrangements on Surfaces fficiency is just intelligentlaziness. Anonymous Minkowski Sum Construction We present two exact and robust implementations of efficient output-sensitive algorithms tocompute the Minkowski sum of two polytopes in R . We demonstrate the effectiveness of ourMinkowski-sum computations through simple applications that exploit these operations todetect collision, and compute the Euclidean separation distance between, and the directionalpenetration depth of, two polytopes in R . In Chapter 5 we show a more involved applicationof these operations.Each method we have developed uses a different variant of Gaussian maps to main-tain dual representations of polytopes. Each method employs a different variant of two-dimensional arrangements to maintain the dual representations, and it makes use of manyoperations applied to arrangements in the corresponding representations. The first methoduses the traditional (spherical) Gaussian map. The map is represented as an arrangementof geodesic arcs embedded on the unit sphere. The second method uses a data structurecalled Cubical Gaussian Map . It consists of six arrangements induced by linear segmentsembedded in the plane. The six arrangements correspond to the six faces of the unit cube— the parallel-axis cube circumscribing the unit sphere.A simple method to compute the Minkowski sum of two polytopes is to compute the con-vex hull of the pairwise sum of the vertices of the two polytopes. Although there are manyimplementations of various algorithms to compute Minkowski sums and answer proximityqueries, we are unaware of the existence of complete implementations of methods to computeexact Minkowski sums other than (i) the naive method above, (ii) a method based on Nefpolyhedra embedded on the sphere [HKM07], and (iii) an implementation by Weibel [28] ofFukuda’s algorithm [Fuk04]. Both our methods exhibit much better performance than theother methods in all cases, as demonstrated by the experiments reported in Table 3.5. Ourmethods well handle degenerate cases that require special treatment when alternative rep-resentations are used. For example, the case of two parallel facets facing the same direction,one from each polytope, does not bear any burden on our methods, and neither does the456 Chapter 3. Minkowski Sum Construction extreme case of two polytopes with identical sets of facet normals.The results of experimentation with a broad family of convex polyhedra are reported.The relevant programs, source code, data sets, and documentation are available at . A short movie [FH05] that describes some of the conceptsof the cubical Gaussian-map method can be downloaded from http://acg.cs.tau.ac.il/projects/internal-projects/gaussian-map-cubical/Mink3d.avi . Another shortmovie [FSH08b] that describes some of the concepts of the (spherical) Gaussian map methodamong the other can be downloaded from http://acg.cs.tau.ac.il/projects/internal-projects/arr-geodesic-sphere/movie/aos-xvid.avi .Both methods are implemented on top of Cgal , and are mainly based on the arrangementpackage of the library (see Chapter 2 and [FWH04, WFZH07b]), although other parts, suchas the polyhedral-surface package developed by L. Kettner [Ket99], are used as well. Insome cases it is sufficient to build only portions of the boundary of the Minkowski sumof two given polytopes to answer collision and proximity queries efficiently. This requireslocating the corresponding features that contribute to the sought portion of the boundary.As both methods we have developed employ two-dimensional arrangements implementedon top of the Cgal arrangement package, we harness the ability to answer point-locationqueries efficiently that comes along, to locate corresponding features of two given polytopes.The rest of this chapter is organized as follows. The Gaussian maps dual representationsof polytopes in R are described in Section 3.1 along with some of their properties. InSections 3.2 and 3.3 we show how 3D Minkowski sums can be computed efficiently, afterthe summands are converted to (spherical) Gaussian maps and cubical Gaussian maps,respectively. Section 3.4 presents an exact implementation of an efficient collision-detectionalgorithm under translation based on either of the dual representations. In the last section,(a) (b)(c) (d)Figure 3.1: (a) A tetrahedron, (b) the Gaus-sian map of the tetrahedron, (c) a cube, and (d)the Gaussian map of the cube. dedicated to experimental results, we highlightthe performance of our methods on variousbenchmarks. Suggestions for future directionsare provided in the conclusion chapter in Sec-tion 6.4. The software access-information alongwith some further design details are provided inthe Appendix. The Gaussian map G = G ( P ) of a com-pact convex polyhedron P in Euclidean three-dimensional space R is a set-valued functionfrom P to the unit sphere S , which assigns toeach point p on the boundary of P the set of out-ward unit normals to support planes to P at p .Thus, the whole of a facet f of P is mapped un-der G to a single point, representing the outwardunit normal to f . An edge e of P is mapped to .2. The (Spherical) Gaussian-Map Method 47a (geodesic) segment G ( e ) on S , whose length is easily seen to be the exterior dihedral angleat e . A vertex v of P is mapped by G to a spherical polygon G ( v ), whose sides are the imageunder G of edges incident to v , and whose angles are the angles supplementary to the planarangles of the facets incident to v ; that is, G ( e ) and G ( e ) meet at angle π − α whenever e and e meet at angle α . In other words, G ( v ) is exactly the “spherical polar” of the link of v in P . (The link of a vertex is the intersection of an infinitesimal sphere centered at v with P , rescaled, so that the radius is 1.) The above implies that G ( P ) is combinatorially dualto P and an arrangement embedded on the unit sphere [HRS92]. Extending the mappingabove, by marking each face f = G ( v ) of the arrangement with its dual vertex v , enables aunique inverse Gaussian mapping, denoted by G − , which maps an extended arrangementembedded on the unit sphere back to a polytope boundary. x d = 1 u ˆ u d o An alternative and practical definition follows. A direction in R can be represented by a point u ∈ S . Let P be a polytopein R , and let V denote the set of its boundary vertices. Fora direction u , we define the extremal point in direction u to be λ V ( u ) = arg max p ∈ V h u, p i , where h· , ·i denotes the inner product.The decomposition of S into maximal connected regions, so thatthe extremal point is the same for all directions within any region forms the Gaussian mapof P . For some u ∈ S the intersection point of the ray ~ou emanating from the origin withone of the planes listed below is a central projection of u denoted as ˆ u d , and illustrated onthe right. The relevant planes are x d = 1 , d = 1 , , 3, if u lies in the positive respectivehemisphere, and x d = − , d = 1 , , Cubical Gaussian Map ( Cgm ) C = C ( P ) of a polytope P in R is a set-valued function from P to the six faces of the unit cube whose edges areparallel to the major axes and are of length two. A facet f of P is mapped under C to acentral projection of the outward unit normal to f onto one of the cube faces. Observe that,a single edge e of P is mapped to a chain of at most four connected segments that lie infour adjacent cube-faces respectively, and a vertex v of P is mapped to at most five abuttingconvex dual faces that lie in five adjacent cube-faces, respectively. The decomposition of theunit-cube faces into maximal connected regions, so that the extremal point is the same forall directions within any region forms the Cgm of P . Likewise, the inverse Cgm , denotedby C − , maps the six extended arrangement embedded on the six faces of the unit cube tothe polytope boundary. Figure 3.2 shows the Cgm of a tetrahedron. Armed with the geometry-traits class for geodesic arcs on the sphere (see Section 2.6), wecompute Minkowski sums of convex polyhedra, by overlaying their respective Gaussian mapsrepresented as arrangements of geodesics on the sphere. Each face f of the Gaussian mapis extended with the coordinates of its dual vertex v = G − ( f ), resulting with a uniquerepresentation.8 Chapter 3. Minkowski Sum Construction (a) (b) (c)Figure 3.2: (a) A tetrahedron, (b) the Cgm of the tetrahedron, and (c) the Cgm unfolded. Thicklines indicate real edges. An input model of a polytope is provided as a polyhedral mesh in R . A polyhedral meshrepresentation consists of an array of boundary vertices and the set of boundary facets,where each facet is described by an array of indices into the vertex array. Constructing theGaussian map of a model given in this representation is done indirectly. First, the Cgal Polyhedron 3 [Ket99] data-structure that represents the model is constructed. This datastructure provides quick access to the incidence relations on the polytope features. Then,the Gaussian map is constructed from the accessible information stored in the Polyhedron 3 data-structure. Once the construction of the Gaussian map is complete, the Polyhedron 3 intermediate representation is discarded.Tetrahedron Octahedron Icosahedron Dioctagonal PyramidPentagonal Truncated Geodesic EllipsoidHexecontahedron Icosidodecahedron Sphere level 4 like polyhedronFigure 3.3: Gaussian maps of various polytopes. The Polyhedron 3 data-structure, like the arrangement Dcel , is based on the imple-mentation of an Hds ; see Section 2.3. It consists of extendible vertices, halfedges, and facets .2. The (Spherical) Gaussian-Map Method Polyhedron 3 data-structure is extended with a Boolean flagthat indicates whether the vertex or the halfedge respectively have already been processedduring the construction of the arrangement that represents the Gaussian map. Each facetis extended with the handle of an arrangement vertex. The handle extension of a facet f that has already been processed points to the dual vertex v = G ( f ). The procedure thatconverts an (extended) Polyhedron 3 data-structure P representing a polytope to an (ex-tended) Arrangement on surface 2 data-structure G ( P ) consists of two steps. First all fieldextensions of all vertices, halfedges, and facets of P are cleared. Then, a recursive functionprovided with an arbitrary vertex v of P as a single parameter is invoked. This function tra-verses the halfedges incident to v . When it encounters an unprocessed halfedge e , it obtainsthe normal n and the vertex-handle extension h ( v ) of the facet f adjacent to the left of e and the normal n and the vertex-handle extension h ( v ) of the facet f adjacent to the leftof the next halfedge in the cyclic chain of halfedges incident to v . Finally, the geodesic shortarc between n and n is constructed and inserted into the arrangement G ( P ), as explainedbelow, using one of the efficient insertion member-functions of Arrangement on surface 2 ;see Section 2.3.2. Once the insertion is complete, the halfedge e is marked as processed. v ismarked as processed once all its incident halfedges are processed. The function recursivelyinvokes itself providing an unprocessed vertex adjacent to v , and terminates when no suchvertex is found.Let C indicate the new geodesic arc to be inserted into the arrangement representing theGaussian map. Assume that C is u -monotone with respect to the parameterization definedby the geometry-traits class; see Section 2.6. That is, it does not intersect the identificationarc. Let v , v , f , and f be the two vertices and two facets as described above. There arefour cases to handle as follows.1. If v and v are both null, it implies that this is the first attempt to insert a geodesicarc into the arrangement. In this case we call insert in face interior( C , f ) , where f is the single face the Dcel was initialized with; see Section 2.6. The handle of thetwo new vertices associated with the endpoints of the newly created geodesic arc C arestored in the records of the corresponding facets f and f of P respectively for lateruse.2. If v is null but v is not, we call either insert from left vertex( C , v ) or insert from right vertex( C , v ) depending on whether the existing vertex v is tothe right or to the left of C , and update the vertex-handle field of the correspondingfacet f with the new vertex v .3. We handle the analogous case where v is null but v is not similarly.4. If both v and v are not null, we call insert at vertices( C , v , v ) . In this case novertex-handle field needs to be updated.If C intersects the identification arc, it is first split at the intersection point into two u -monotone arcs C and C . Then, C and C are inserted according to four cases similar tothe above. The handling is a bit more intricate. For example, consider the case where v is null but v already exists. Assume that C reaches the right boundary of the parameter0 Chapter 3. Minkowski Sum Construction space and C reaches the left boundary (they both meet at an identified point). If v isassociated with the left endpoint of C , we first call insert from left vertex( C , v ) , andthen insert from left vertex( C , v ′ ) , where v ′ is the new vertex associated with the rightendpoint of C introduced while C is inserted. Otherwise, v must be associated with theright endpoint of C . In this case we first call insert from right vertex( C , v ) , and then insert from right vertex( C , v ′ ) , where v ′ is the new vertex associated with the rightendpoint of C introduced while C is inserted. The other three cases are handled similarly.We have created a large database of models of polytopes. Table 3.3 lists, for a smallsubset of our polytope collection, the number of features in the arrangement of geodesic arcsembedded on the sphere that represents the Gaussian map of each polytope. Recall that thenumber of faces ( F ) of the Gaussian map is always equal to the number of vertices of thepolytope. However, the number of vertices ( V ) of the Gaussian map is either equall to, orgreater than, the number of facets of the primal representation due to intersections betweenGaussian-map edges and the identification arc. A similar argument holds for the edges.That is, the number of halfedges ( HE ) of the Gaussian map is either equall to, or greaterthan, twice the number of edges of the primal representation. An edge of the Gaussianmap that intersects the identification arc must be split at the intersection point into two u -monotone geodesic arcs. The table also lists the time in seconds ( t ) it takes to constructthe arrangement once the intermediate polyhedron is in place, on a Pentium PC clocked at1.7 GHz. (a) (b)Figure 3.4: (a) The Minkowski sum of a tetra-hedron and a cube and (b) the Gaussian map ofthe Minkowski sum. The overlay (see Section 2.4.2 for the exact def-inition) of the Gaussian maps of two polytopes P and Q respectively identifies all pairs of fea-tures of P and Q that have parallel support-ing planes, as they occupy the same space onthe unit sphere, thus, identifying all the pair-wise features that contribute to the boundaryof the Minkowski sum of P and Q . A facet ofthe Minkowski sum is either a facet f of Q trans-lated by a vertex of P supported by a plane par-allel to f , or vice versa, or it is a facet parallelto two parallel planes supporting an edge of P and an edge of Q , respectively. A vertex of the Minkowski sum is the sum of two vertices of P and Q respectively supported by parallel planes.When the overlay operation progresses, new vertices, edges, and faces of the resultingarrangement are created based on features of the two operands. When a new feature iscreated its attributes are updated. There are ten cases that arise and must be handled; seeSection 2.4.2 for the precise enumeration of the various cases. For example, a new face f isinduced by the overlap of two faces f and f of the two summands, respectively. The primalvertex associated with f is set to be the sum of the primal vertices associated with f and f , respectively. .3. The Cubical Gaussian-Map Method ⊕ Cube DP ⊕ ODP PH ⊕ TI El16 ⊕ OEl16Figure 3.5: Gaussian maps of Minkowski sums. DP — Dioctagonal Pyramid, PH — Pentagonal Hexe-contahedron, TI — Truncated Icosidodecahedron, GS4 — Geodesic Sphere level 4, El16 — Ellipsoid-likepolyhedron made of 16 latitudes and 32 longitudes. Table 3.4 lists the number of features ( V , HE , F ) in the arrangement that represents theGaussian map of the respective Minkowski sums. Table 3.5 shows the time in seconds ( t ) ittakes to construct the arrangement once the Gaussian maps of the summands are in place. While using the Cgm increases the overhead of some operations sixfold, and introducesdegeneracies that are not present in the case of alternative representations, it simplifies theconstruction and manipulation of the representation, as the partition of each cube face is aplanar map of segments, a well known concept that has been intensively experimented withduring recent years. Indeed, all the basic software components that the Cgm layer dependson are available in Cgal version 3.3 and higher, while many of the software componentsrequired by the (spherical) Gaussian map method are expected to appear only in a futurerelease of Cgal . The Cgm method, being more mature, exhibits better performance thanthe (spherical) Gaussian map method. One of the reasons for the performance gap is the lackof optimized primitives that operate on unnormalized vectors in R in case of the (spherical)Gaussian map method. Evidently, most of the methods of the geometry-traits class thathandle geodesic arcs embedded on the sphere project their input u -monotone curves andpoints onto one of the axis-aligned planes every time they are invoked, while all geometricobjects in case of the Cgm method are projected onto the plane a priori . This creates anopportunity for optimization; see Section 6.1.4. In addition, the Cgm data structure, beingbased on components confined to the plane, has a broader recognition, as it can be used inrestricted environments, e.g., 3D hardware accelerators; see Section 6.5. We use the Cgal Arrangement 2 data-structure to maintain the planar maps. The construc-tion of the six planar maps from the polytope features and their incident relations is similarto the construction of the arrangement of geodesic arcs that represents the Gaussian map;see Section 3.2.1. As in the case of the (spherical) Gaussian map, the Polyhedron 3 inter-mediate representation is discarded once the construction of the Cgm is complete. However,2 Chapter 3. Minkowski Sum Construction while a single edge of P is mapped to at most two u -monotone geodesic arcs in case of the(spherical) Gaussian map, it can be mapped to a chain of at most four connected segmentsthat lie in four adjacent cube-faces, respectively. In any case the constructions of the Gaus-sian maps of both methods respectively amount to the insertion of curves that are pairwisedisjoint in their interior into the arrangement, an operation that is carried out efficiently,especially when one or both endpoints are known. The construction of the Minkowski sum,described in Section 3.3.2, amounts to the computation of the six overlays of six pairs ofplanar maps, respectively, an operation well supported by the data structure as well.A related dual representation had been considered and discarded before the Cgm rep-resentation was chosen. It uses only two planar maps that partition two parallel planesrespectively instead of six, but each planar map partitions the entire plane. In this 2-maprepresentation facets that are near orthogonal to the parallel planes are mapped to pointsthat are far away from the origin. The exact representation of such points requires coordi-nates with large bit-lengths, which increases significantly the time it takes to perform exactarithmetic operations on them. Moreover, facets exactly orthogonal to the parallel planesare mapped to points at infinity, and require special handling all together.Features that are not in general position, such as two parallel facets facing the samedirection, one from each polytope, or worse yet, two identical polytopes, typically requirespecial treatment. Still, the handling of many of these problematic cases falls under the “gen-eral” case, and becomes transparent to the Gaussian-map layer (either cubical or spherical).Consider for example the case of two neighboring facets in one polytope that have parallelneighboring facets in the other. This translates to overlapping segments in case of the Cgm method and overlapping geodesic arcs in case of the (spherical) Gaussian-map method, onefrom each Gaussian map of the two polytopes, that appear during the Minkowski sum com-putation. The algorithm that computes it is oblivious to this condition, as the underlyingarrangement data structure, either embedded in the plane or on the sphere, is perfectly ca-pable of handling overlapping curves. However, as mentioned above, other degeneracies doemerge, and are handled successfully. One example, in case of the Cgm , is a facet f mappedto a point that lies on an edge of the unit cube, or even worse, coincides with one of the eightcorners of the cube. Figure 3.7(a,b,c) depicts an extreme degenerate case of an octahedronoriented in such a way that its eight facets are mapped to the eight vertices of the unit cube,respectively. The (spherical) Gaussian map method is not free of degenerate conditions. Forexample, a facet mapped to a point that lies on the identification arc, or worse yet, coincideswith one of the two poles.The dual representation is extended further, in order to handle all these degeneracies andperform all the necessary operations as efficiently as possible. Each planar map is initializedwith four edges and four vertices that define the unit square — the parallel-axis squarecircumscribing the unit circle. During construction, some of these edges or portions of themalong with some of these vertices may turn into real elements of the Cgm . The introduction Each planar map that corresponds to one of the six unit-cube faces in the Cgm representation alsopartitions the entire plane, but only the [ − , − × [1 , 1] square is relevant. The unbounded face, whichcomprises all the rest, is irrelevant. Other conditions translate to overlapping segments in case of the Cgm or geodesic arcs in case of the(spherical) Gaussian map method as well. .3. The Cubical Gaussian-Map Method XY YX XY XYZ 21 03 32 10 01 32 Y XY XY X X Y Z 13 20 Y X XY XY XY Z 01 32 XY X YY X X YZ 13 20 Figure 3.6: The data structure. Large-font numbersindicate plane ids. Small-font numbers indicate cornerids. X and Y axes in different 2D coordinate systemsare rendered in different colors. The exact mapping from a facet nor-mal in the 3D coordinate-system to a pairthat consists of a planar map and a planarpoint in the 2D coordinate-system is de-fined precisely through the indexing andordering system, illustrated in Figure 3.6.Now before your eyes cross permanently,we advise you to keep reading the nextfew lines, as they reveal the meaning ofsome of the enigmatic numbers that ap-pear in the figure. The six planar mapsare given unique ids from 0 through 5. Ids0, 1, and 2 are associated with the planes x = − y = − 1, and z = − 1, respectively,and ids 3, 4, and 5 are associated with theplanes x = 1, y = 1, and z = 1, respec-tively. The major axes in the 2D Cartesiancoordinate-system of each planar map aredetermined by the 3D coordinate-system.The four corner vertices of each planarmap are also given unique ids from 0 through 3 in lexicographic order in their respective 2Dcoordinate-system, see Table 3.1 columns titled Underlying Plane and 2D Axes .Table 3.1: The coordinate systems and the cyclic chains of corner vertices. PM stands for PlanarMap , and Cr stands for Corner . Underlying 2D Axes CornerPlane 0 (0,0) (0,1) (1,0) (1,1) Id Eq X Y PM Cr PM Cr PM Cr PM Cr x = − Z Y y = − X Z z = − Y X x = 1 Y Z y = 1 Z X z = 1 X Y Dcel used to maintain the incidence relations of the vertices,halfedges, and faces of the Arrangement 2 data structure (see Section 2.3) is extended tohold additional attributes. Some of the attributes are introduced only in order to expeditethe computation of certain operations, but most of them are necessary to handle degeneratecases such as a planar vertex lying on the unit-square boundary. Each planar-map vertex4 Chapter 3. Minkowski Sum Construction v is extended with (i) the coefficients of the plane containing the polygonal facet C − ( v )(see Section 3.1 for the definition of C and C − ), (ii) the location of the vertex — anenumeration indicating whether the vertex coincides with a cube corner, or lies on a cubeedge, or contained in a cube face, (iii) a Boolean flag indicating whether it is non-artificial(there exists a facet that maps to it), and (iv) a pointer to a vertex of a planar map associatedwith an adjacent cube-face that represents the same central projection for vertices thatcoincide with a cube corner or lie on a cube edge. Each planar-map halfedge e is extendedwith a Boolean flag indicating whether it is non-artificial (there exists a polytope edge thatmaps to it). Each planar-map face f is extended with the polytope vertex that maps to it v = C − ( f ). (a) (b) (c)(d) (e) (f)(g) (h) (i)Figure 3.7: (a) An octahedron, (d) a dioctagonal pyramid, (g) an ellipsoid-like polyhedron level 16,(b,e,h) the Cgm of the respective polytope, and (c,f,i) the Cgm unfolded. Each vertex that coincides with a unit-cube corner or lies on a unit-cube edge containsa pointer to a vertex of a planar map associated with an adjacent cube face that representsthe same central projection. Vertices that lie on a unit-cube edge (but do not coincide with .3. The Cubical Gaussian-Map Method Polyhedron 3 and Arrangement 2 data structures for example both use a Dcel datastructure that follows the convention above.We provide a fast clockwise traversal of the faces incident to any givenvertex v . Clockwise traversals around internal vertices are immediatelyavailable by the Dcel . Clockwise traversals around boundary verticesare enabled by the cyclic chains above. This traversal is used to calculatethe normal to the (primary) polytope-facet f = C − ( v ) and to renderthe facet. Fortunately, rendering systems are capable of handling asequence of vertices that define a polygon in clockwise order as well, anorder opposite to the conventional ordering above.The data structure also supports a fast traversal over the planar-maphalfedges that form each one of the four unit-square edges. This traversalis used during construction to quickly locate a vertex that coincides with acube corner or lies on a cube edge. It is also used to update the cyclic chainsof pointers mentioned above; see Section 3.3.2.We maintain a flag that indicates whether a planar vertex coincides with a cube corner,a cube edge, or a cube face. At first glance this looks redundant. After all, this informationcould be derived by comparing the x and y coordinates to − The number of features ofthe six planar maps of the Cgm of thedioctagonal pyramid. Planar map V HE F 0, ( x = − 1) 12 32 61, ( y = − 1) 28 80 142, ( z = − 1) 12 32 63, ( x = 1) 12 32 64, ( y = 1) 21 72 175, ( z = 1) 12 32 6 Total 97 280 55The table to the right shows the number of ver-tices ( V ), halfedges ( HE ), and faces ( F ) of the sixplanar maps that comprise the Cgm of the dioctago-nal pyramid shown in Figure 3.7 (d,e,f). The numberof faces of each planar map include the unboundedface. Table 3.3 shows the number of features in theprimal and dual representations of a small subset ofour polytopes collection, on which we report the re-sults of experiments below. The number of planarfeatures is the total number of features of the six pla-nar maps.6 Chapter 3. Minkowski Sum Construction A similar argument regarding the representation of Minkowski sums using Gaussian mapsmentioned in Section 3.2 holds for the cubical Gaussian maps with the unit cube replacingthe unit sphere. More precisely, a single map that subdivides the unit sphere is replacedby six planar maps, and the computation of a single overlay is replaced by the computationof six overlays of corresponding pairs of planar maps. Recall that each (primal) vertex isassociated with a planar-map face, and is the sum of two vertices associated with the twooverlapping faces of the two Cgm ’s of the two input polytopes, respectively.Each planar map in a Cgm is a convex subdivision. Finke and Hinrichs [FH95] describehow to compute the overlay of such special subdivisions optimally in linear time. However,a preliminary investigation shows that a large constant governs the linear complexity, whichrenders this choice less attractive. Instead, we resort to a sweep-line based algorithm thatexhibits good practical performance, and incurs a mere logarithmic factor over the optimalcomputing time. In particular we use the overlay operation supported by the Arrangementpackage. It requires the provision of a complementary component that is responsible for up-dating the attributes of the Dcel features of the resulting six planar maps; see Section 2.4.2.The overlay operates on two instances of Arrangement 2 . In the description below v , e , and f denote a vertex, a halfedge, and a face of the first operand respectively, and v , e , and f denote the same feature types of the second operand, respectively. When theoverlay operation progresses, new vertices, halfedges, and faces of the resulting planar mapare created based on features of the two operands. Exactly ten cases described below ariseand must be handled. When a new feature is created its attributes are updated. The updatesperformed in all cases except for case (1) are simple and require constant time.1. A new vertex v is induced by coinciding vertices v and v .The location of the vertex v is set to be the same as the location of the vertex v (the locations of v and v must be identical). The induced vertex is not artificial ifand only if (i) at least one of the vertices v or v is not artificial, or (ii) the vertexlies on a cube edge or coincides with a cube corner, and both vertices v and v haveTable 3.3: Complexities of the primal and dual representations. DP — Dioctagonal Pyramid, PH —Pentagonal Hexecontahedron, TI — Truncated Icosidodecahedron, GS4 — Geodesic Sphere level 4, El16— Ellipsoid-like polyhedron made of 16 latitudes and 32 longitudes, t - time consumption in seconds. Object Type Primal SGM CGMV E F V HE F t V HE F t Tetrahedron 4 6 4 4 12 4 0.01 42 102 21 0.01Octahedron 6 12 8 10 28 6 0.01 24 48 12 0.01Icosahedron 12 30 20 21 62 12 0.01 72 192 36 0.01DP 17 32 17 25 80 17 0.01 97 280 55 0.01PH 60 150 92 101 318 60 0.03 200 600 112 0.02TI 120 180 62 77 390 120 0.05 230 840 202 0.03GS4 252 750 500 506 1512 252 0.08 708 2124 366 0.07El16 482 992 512 528 2016 482 0.11 776 2752 612 0.06 .4. Exact Collision Detection v is induced by a vertex v that lies on an edge e .The location of the vertex v is set to be the same as the location of the vertex v . v isnot artificial if and only if v is not artificial or e is not artificial.3. An analogous case of a vertex v that lies on an edge e .4. A new vertex v is induced by a vertex v that is contained in a face f .The attributes of the vertex v are set to be the same as the attributes of the vertex v .5. An analogous case of a vertex v contained in a face f .6. A new vertex v is induced by the intersection of two edges e and e .The vertex v cannot lie on a cube edge and cannot coincide with a cube corner. Thus,it is necessarily not artificial.7. A new edge e is induced by the overlap of two edges e and e .The edge e is not artificial if at least one of e or e is not artificial.8. A new edge e is induced by an edge e that is contained in a face f .The edge e is not artificial if e is not artificial.9. An analogous case of an edge e contained in a face f .10. A new face f is induced by the overlap of two faces f and f .The primal vertex associated with f is set to be the sum of the primal vertices associ-ated with f and f , respectively.After the six map overlays are computed, some maintenance operations must be per-formed to obtain a valid Cgm representation. As mentioned above, the global data consistsof the six planar maps and 24 references to vertices that coincide with the unit-cube cor-ners. For each planar map we traverse its vertices, obtain the four vertices that coincidewith the unit-cube corners, and initialize the global data. We also update the cyclic chainsof pointers to vertices that represent identical central projections. To this end, we exploitthe fast traversal over the halfedges that coincide with the unit-cube edges mentioned inSection 3.3.1.The complexity of a single overlay operation is O ( k log n ), where n is the total numberof vertices in the input planar maps, and k is the number of vertices in the resulting planarmap. The total number of vertices in all the six planar maps in a Cgm that represents apolytope P is of the same order as the number of facets in the primary polytope P . Thus,the complexity of the entire overlay operation is O ( F log( F + F )), where F and F are thenumber of facets in the input polytopes respectively, and F is the number of facets in theMinkowski sum. Computing the separation distance between two polytopes with m and n features respec-tively can be done in O (log m log n ) time, after an investment of at most linear time inpreprocessing [DK90]. Many practical algorithms that exploit spatial and temporal coher-ence between successive queries have been developed, some of which became classic, such asthe GJK algorithm [GJK88] and its improvement [Cam97], and the LC algorithm [LC91] andits optimized variations [EL00, GHZ99, Mir98]. Several general-purpose software libraries8 Chapter 3. Minkowski Sum Construction that offer practical solutions are available today, such as the SOLID library [24] based onthe improved GJK algorithm, the SWIFT library [26] based on an advanced version of theLC algorithm, the QuickCD library [21], and more. For an extensive review of methods andlibraries see the survey [LM04].Given two polytopes P and Q , detecting collision between them and computing their rel-ative placement can be conveniently done in the configuration space, where their Minkowskisum M = P ⊕ ( − Q ) resides. These problems can be solved in many ways, and not all requirethe explicit representation of the Minkowski sum M . However, having it available is attrac-tive, especially when the polytopes are restricted to translations only, as the combinatorialstructure of the Minkowski sum M is invariant to translations of P or Q . The algorithmsdescribed below are based on the following well known observations (see Chapter 1 for defi-nitions): P u ∩ Q w = ∅ ⇔ w − u ∈ M = P ⊕ ( − Q ) ,π ( P u , Q w ) = min {k t k | ( w − u + t ) ∈ M, t ∈ R } ,δ r ( P u , Q w ) = inf { α | ( w − u + α~r ) / ∈ M, α ∈ R } . Given two polytopes P and Q in either (spherical) Gaussian map or Cgm representationrespectively, we reflect Q through the origin to obtain − Q , compute the Minkowski sum M = P ⊕ ( − Q ), and retain it in the respective Gaussian-map representation G ( M ). Then,each time P or Q or both translate by two vectors u and w in R respectively, we apply aprocedure that determines whether the query point s = w − u is inside, on the boundary of, oroutside M . In addition to an enumeration of one of the three conditions above, the procedurereturns a witness of the respective relative placement. Let r be a ray emanating from aninternal point c ∈ M and going through s . If the (spherical) Gaussian map representationis used, the witness data is the vertex v = G ( f ) — a mapping of a facet f of M embeddedon the sphere and stabbed by the ray r . If the Cgm representation is used, the witnessdata is a pair that consists of a vertex v = G ( f ) — a mapping of a facet f of M embeddedin a unit cube face, and the planar map P containing v . This information is used as ahint in consecutive invocations. The internal point c could be the average of all verticesof M computed once and retained along M , or just the midpoint of two vertices that havesupporting planes with opposite normals easily extracted from either map representation.Once f is obtained, determining whether P u and Q w collide is trivial, according to the firstformula (of the three) above. The query point s is contained in the open half-space definedby the supporting plane to f if and only if s is outside of M , this occurs if and only if P u does not collide with Q w . Figure 3.8: Simulation of motion. The collision-detection procedure applies a local walkon the respective Gaussian map faces. It starts with somevertex v , and then performs a loop moving from the cur-rent vertex to a neighboring vertex, until it reaches thefinal vertex. If the Cgm representation is used, the pro-cedure may jump from a planar map associated with onecube-face to a different one associated with an adjacentcube-face. The first time the procedure is invoked, v ischosen to be a vertex that lies on the central projection of .5. Minkowski Sum Complexity r . In consecutive calls, v is chosento be the final vertex of the previous call exploiting spatial and temporal coherence. Thefigure above is a snapshot of a simulation program that detects collision between a staticobstacle and a moving robot, and draws the obstacle and the trail of the robot; instructionsto obtain, install, and execute the program appear in the Appendix. The Minkowski sum isrecomputed only when the robot is rotated, which occurs every other frame. The programhas the distinctive feature of being able to identify the case where the robot grazes the ob-stacle, but does not penetrate it, since it produces exact results. The computation takesjust a fraction of a second on a Pentium PC clocked at 1.7 GHz using either representation.Similar procedures that compute the directional penetration-depth and minimum distanceare available as well. (a) (b) (c)Figure 3.9: (a) The Minkowski sum (of two polytopes) the complexity of which is maximal, (b) the Cgm of the Minkowski sum, and (c) the Cgm unfolded. Red lines are graphs of edges that originatefrom one polytope and blue lines are graphs of edges that originate from the other. In Chapter 4 we show that the exact maximum number of facets of Minkowski sums oftwo polytopes is 4 mn − m − n + 26, where m and n are the number of facets of the twosummands, respectively. This bound is tight [FHW07]. The example depicted in Figure 3.9shows a Minkowski sum that reaches this maximum complexity. It is the sum of two identicalpolytopes, each containing n faces ( n = 11 in Figure 3.9), but one is rotated about thevertical axis approximately ◦ relative to the other. The polytopes are specifically shapedto ensure that the number of intersections between dual edges, which are the mappings ofthe polytope edges, is maximal. A careful counting reveals that the number of vertices in thedual representation excluding the artificial vertices reaches 4 · · − · − · 11 + 26 = 312,which is the number of facets of the Minkowski sum.Not every pair of polytopes yields a Minkowski sum proportional to mn . As a matterof fact, it can be as low as n in the extremely-degenerate case of two identical polytopesvariant under scaling. Even if no degeneracies exist, the complexity can be proportional to The results of all rotations are approximate, as we have not yet dealt with exact rotation. One of ourfuture goals is the handling of exact rotations. Chapter 3. Minkowski Sum Construction only m + n , as in the case of two geodesic spheres level l = 2 slightly rotated with respectto each other, depicted in Figure 3.10. Naturally, an algorithm that accounts for all pairsof vertices, one from each polytope, is rendered inferior compared to an output-sensitivealgorithm like ours in such cases, as we demonstrate in the next section.(a) (b) (c)Figure 3.10: (a) The Minkowski sum of two geodesic spheres level 2 slightly rotated with respect toeach other, (b) the Cgm of the Minkowski sum, and (c) the Cgm unfolded. (a) (b) (c)Figure 3.11: (a) The Minkowski sum of two approximately orthogonal squashed dioctagonal pyramids,(b) the Cgm , and (c) the Cgm unfolded. As mentioned above, the Minkowski sum of two polytopes is the convex hull of the pair-wise sum of the vertices of the two polytopes. We have implemented this straightforwardmethod using the Cgal convex hull 3 function, which uses the Polyhedron 3 data struc-ture to represent the resulting polytope, and used it to verify the correctness of our twomethods. We compared the time it took to compute exact Minkowski sums using our twomethods, a third method implemented by Hachenberger based on Nef polyhedra embed-ded on the sphere [HKM07], a fourth method implemented by Weibel [28], based on anoutput-sensitive algorithm designed by Fukuda [Fuk04], and the naive convex-hull method. An icosahedron, every triangle of which is divided into 4 l triangles using class I aperture 4 partitionmethod, whose vertices are elevated to the circumscribing sphere [SWK03]. .6. Experimental Results O ( δLP (3 , δ ) V ), where δ = δ + δ is the sum of the maximal degrees of vertices, δ and δ , in the two input polytopesrespectively, V is the number of vertices of the resulting Minkowski sum, and LP ( d, m ) isthe time required to solve a linear programming in d variables and m inequalities. Notethat Fukuda’s algorithm is more general, as it can be used to compute the Minkowski sumof polytopes in an arbitrary dimension d , and as far as we know, it has not been optimizedspecifically for d = 3.Table 3.4: Complexities of primal and dual Minkowski-sum representations. DP — Dioctagonal Pyra-mid, ODP — Dioctagonal Pyramid orthogonal to DP, PH — Pentagonal Hexecontahedron, TI —Truncated Icosidodecahedron, GS4 — Geodesic Sphere level 4, RGS4 — Rotated Geodesic Spherelevel 4, El16 — Ellipsoid-like polyhedron made of 16 latitudes and 32 longitudes, OEl16 — Ellipsoid-likepolyhedron made of 16 latitudes and 32 longitudes orthogonal to El16. Summand 1 Summand 2 Minkowski SumPrimal SGM CGMV E F V HE F V HE F Icosahedron Icosahedron 12 30 20 21 62 12 72 192 36DP ODP 131 261 132 141 540 131 242 832 186PH TI 248 586 340 429 1712 429 514 1670 333GS4 RGS4 1048 2582 1536 1564 5220 1048 1906 6288 1250El16 OEl16 2260 4580 2322 2354 9224 2260 2826 10648 2510Table 3.5: Time consumption (in seconds) of the Minkowski-sum computation. CH — the ConvexHull method, SGM — the (spherical) Gaussian map based method, CGM — the cubical Gaussian-mapbased method, NGM — the Nef based method, Fuk — Fukuda’s Linear-Programming based algorithm, F F F — the ratio between the product of the number of input facets and the number of output facets. Summand 1 Summand 2 SGM CGM NGM Fuk CH F F F Icosahedron Icosahedron 0.01 Cgm method in particular is more than fifty times faster than the convex-hull method inone case. The number of models used as summands in the listed experiments is just a smallfraction of the total number of models in our collection, which contains hundreds of modelsof polytopes; see Section A.2 for instructions how to download the corresponding files. The2 Chapter 3. Minkowski Sum Construction listed experiments are just a small sample of all the experiments we have conducted. Thelast column of the table indicates the ratio F F F , where F and F are the number of facetsof the input polytopes respectively, and F is the number of facets of the Minkowski sum. Asthis ratio increases, the relative performance of the output-sensitive algorithms compared tothe convex-hull method improves as expected. Figure 3.12 illustrates some of the resultingMinkowski sums listed in Table 3.5.The SGM , CGM , NGM , and CH methods are all based on Cgal . As the implementa-tion of these methods and of Cgal is generic, it is possible to instantiate specific components,such as the number type, the geometric kernel, and the segment-handling traits class withthe model of the respective concept that achieves the best result. The code of the programsused to obtain the results listed in Table 3.5 are instantiated with the CGAL::Gmpq exactrational number-type, the Simple cartesian geometric kernel, and the Lazy kernel kerneladaptor, which performs lazy exact computations [PF06]. The computation is delayed untila point where the approximation with interval arithmetic is not precise enough to performsafe comparisons. In other words, these programs need only to compute to sufficient preci-sion to evaluate predicates correctly, exploiting a significant relaxation of the naive conceptof numerical exactness.We experimented with two different exact number types: One provided by Leda leda::rational , and another based on GMP CGAL::Gmpq . Theformer does not normalize the rational numbers automatically. Therefore, we had to initiatenormalization operations to contain their bit-length growth. In case of the Cgm method,for example, we chose to do it right after the central projections of the facet-normals arecalculated, and before the chains of segments, which are the mapping of facet-edges, areinserted into the planar maps. Our experience shows that indiscriminate normalizationconsiderably slows down the planar-map construction, and the choice of number type mayhave a drastic impact on the performance of the code overall. .6. Experimental Results (a) The Minkowski sum of two approximately orthogonal dioctagonal pyramids, (d)the Minkowski sum of a Pentagonal Hexecontahedron and a Truncated Icosidodecahedron, (g) theMinkowski sum of two approximately orthogonal ellipsoid-like polyhedra, (b,e,h) the Cgm of the re-spective polytope, and (c,f,i) the Cgm unfolded. Chapter 3. Minkowski Sum Construction was working on the proofof one of my poems all themorning, and took out acomma. In the afternoon Iput it back again. Oscar Wilde Exact Complexity of Minkowski Sums We present a tight bound on the exact maximum complexity of Minkowski sums of polytopesin R . In particular, we prove that the maximum number of facets of the Minkowski sum of k polytopes with m , m , . . . , m k facets respectively is bounded from above by P ≤ i 5) + P ≤ i ≤ k m i + (cid:0) k (cid:1) . Given k positive integers m , m , . . . , m k , we describe how toconstruct two polytopes with corresponding number of facets, such that the number of facetsof their Minkowski sum is exactly P ≤ i 5) + P ≤ i ≤ k m i + (cid:0) k (cid:1) . When k = 2, for example, the expression above reduces to 4 m m − m − m + 26. Figure 4.1illustrates some polytopes that when rotated properly and used as summand, the number offacets of the resulting Minkowski sums is maximal.Various methods to compute the Minkowski sum of two polyhedra in R have beenproposed; see Section 1.2 for details about these methods and about the combinatorialcomplexity of the sum. Recall, that (i) a common approach for computing Minkowski sums(a) (b) (c) (d)Figure 4.1: Gaussian maps of summands of Minkowski sums with maximal number of facets. Thepolytopes represented by the Gaussian maps (a), (b), (c), and (d) consist of 4, 5, 11, and 101 facets,respectively. Chapter 4. Exact Complexity of Minkowski Sums of non-convex polyhedra decomposes each polyhedron into convex pieces, and computespairwise Minkowski sums of the convex pieces, and (ii) all the efficient methods are outputsensitive. Thus, the exact maximum complexity of the Minkowski sum structure has practicalimplications.One method to compute the Minkowski sum of two polytopes is to compute the convexhull of the pairwise sum of the vertices of the two polytopes. While being simple and easyto implement, the time complexity of this method is Ω( mn ) regardless of the size of theresulting sum, which can be as low as ( m + n ) (counting facets) for degenerate cases. InChapter 3 we describe several complete implementations of output-sensitive methods forcomputing exact Minkowski sums (beyond the naive method mentioned above), includingtwo methods that we introduce. These methods exploit efficient innovative techniques inthe area of exact geometric computing to minimize the time it takes to ensure exact results.However, even with the use of these techniques, the amortized time of a single arithmeticoperation is larger than the time it takes to carry out a single arithmetic operation onnative number types, such as floating point. Thus, the constant that scales the dominantelement in the expression of the time complexity of these algorithms increases, which makesthe question this chapter attempts to answer, “What is the exact maximum complexity ofMinkowski sums of polytopes in R ?”, even more relevant.Gritzmann and Sturmfels [GS93] formulated an upper bound on the number of features f di of any given dimension i of the Minkowski sum of many polytopes in d dimensions: f di ( P ⊕ P ⊕ . . . ⊕ P k ) ≤ (cid:0) ji (cid:1) d − i − P h =0 (cid:0) j − i − h (cid:1) for i = 0 , , . . . d − 1, where j denotes the numberof non-parallel edges of P , P , . . . , P k . According to this expression, the number of facets f of the Minkowski sum of two polytopes in R is bounded from above by j ( j − R in terms of the number of vertices of the summands: f ( P ⊕ P ) ≤ f ( P ) f ( P ) + f ( P ) + f ( P ) − 6. They also studied the properties of Minkowski sums ofperfectly centered polytopes and their polars, and provided a tight bound on the number ofvertices of the sum of polytopes in any given dimension.The main result presented in this chapter concerning two polytopes follows. Theorem 4.1. Let P , P , . . . , P k be a set of k polytopes in R , such that the number of facetsof P i is m i for i = 1 , , . . . , k . The number of facets of the Minkowski sum P ⊕ P ⊕ . . . ⊕ P k cannot exceed P ≤ i The maximum exact number of edges in a Gaussian map G ( P ) of apolytope P with m facets is m − . The maximum exact number of faces is m − . Bothmaxima occur at the same Gaussian maps. The above can be obtained by a simple application of Euler’s formula for planar graphsto the Gaussian map G ( P ). It can be used to trivially bound the exact maximum number offacets of the Minkowski sum of two polytopes defined as f ( m, n ) = max { f ( P ⊕ Q ) | f ( P ) = m, f ( Q ) = n } , where f ( P ) is the number of facets of a polytope P . First, we can usethe bound on the number of edges to obtain: f ( m, n ) ≤ m + n + (3 m − · (3 n − 6) =9 mn − m − n + 36. Better yet, we can plug the bound on the number of dual faces,which is the number of primal vertices, in the expression introduced by Fukuda and Weibel,see above, to obtain: f ( m, n ) ≤ (2 m − · (2 n − m − n − − mn − m − n +2.Still, we can improve the bound even further, but first we need to bound the number of facesin G ( M ). Lemma 4.3. Let G and G be two Gaussian maps of convex polytopes, and let G be theiroverlay. Let f , f , and f denote the number of faces of G , G , and G , respectively. Then, f ≤ f · f . Each face in the overlay is an intersection of a face of each map. Since these faces arespherically convex (and smaller than hemispheres), the intersection is also spherically convex(and thus connected). This lemma is similar to the one where convex planar maps replacethe Gaussian maps. Nevertheless, we provide a formal proof directly applied to the sphericalcase. Proof. We label each face in G by a pair of indices of the originating overlaid faces in G and G respectively, and argue that no two faces in G can have the same label. Assume to The number of facets is strictly equal to the given expression, only when no degeneracies occur. Chapter 4. Exact Complexity of Minkowski Sums the contrary that there exist two faces h a and h b in G that have the same label, say h i, j i .That is, the faces h i in G and h j in G induce the two distinct faces h a and h b . Pick twopoints a ∈ h a and b ∈ h b . There must be a geodesic segment between a and b that is entirelycontained in h i and also in h j , as both maps are spherically convex. This implies that noneof the edges in G and G split this geodesic segment, contradicting the fact that they residein two different faces of G .We are ready to tackle the upper bound of Theorem 4.1 for the special case k = 2, thatis, prove that the number of facets of the Minkowski sum P ⊕ Q of two polytopes P and Q with m and n facets respectively cannot exceed 4 mn − m − n + 26; see Page 66. Proof. Let v , e , f and v , e , f denote the number of vertices, edges, and faces of G ( P ) and G ( Q ), respectively. The number of vertices, edges, and faces of G ( M ) is denoted as v , e , and f , respectively. Assume that P and Q are two polytopes, such that the number of facets oftheir Minkowski sum is maximal. Recall that the number of facets of a polytope is equal tothe number of vertices of its Gaussian map. Thus, we have v = m , v = n , and v = f ( m, n ).First, we need to show that vertices of G ( P ), vertices of G ( Q ), and intersections betweenedges of G ( P ) and edges of G ( Q ) do not coincide. Assume to the contrary that some do.Then, one of the polytopes P or Q or both can be slightly rotated to escape this degeneracy,but this would increase the number of vertices v = f ( m, n ), contradicting the fact that f ( m, n ) is maximal. Therefore, the number of vertices v is exactly equal to v + v + v x ,where v x denotes the number of intersections of edges of G ( P ) and edges of G ( Q ) in G ( M ).Counting the degrees of all vertices in G ( M ) implies that 2 e + 2 e + 4 v x = 2 e . Using Euler’sformula, we get e + e + 2 v x = f + v + v + v x − 2. Applying Lemma 4.3, we can bound v x from above v x ≤ f f + v + v − − e − e .Observation 4.2 sets an upper bound on the number of edges e . Thus, e can be expressedin terms of ℓ , a non-negative integer, as follows: e = 3 v − − ℓ . Applying Euler’s formula,the number of facets can be expressed in terms of ℓ as well: f = e + 2 − v = 2 v − − ℓ .Similarly, we have e = 3 v − − ℓ and f = 2 v − − ℓ for some non-negative integer ℓ . v x ≤ (2 v − − ℓ )(2 v − − ℓ ) + v + v − − (3 v − − ℓ ) − (3 v − − ℓ ) ≤ v v − v − v + 26 + h ( ℓ , ℓ ) , (4.1)where h ( ℓ , ℓ ) = ℓ ℓ + 5 ℓ + 5 ℓ − v ℓ − v ℓ . G ( P ) consists of a single connected component. Therefore, the number of edges e mustbe at least v − 1. This is used to obtain an upper bound on ℓ as follows: v − ≤ e =3 v − − ℓ , which implies ℓ ≤ v − 5, and similarly ℓ ≤ v − 5. Thus, we have: h ( ℓ , ℓ ) = ℓ ℓ + 5 ℓ + 5 ℓ − v ℓ − v ℓ = ℓ ( ℓ − (2 v − ℓ ( ℓ − (2 v − ≤ . From Equation (4.1) we get that v x ≤ v v − v − v + 26, and since f ( m, n ) = v + v + v x , we conclude that f ( m, n ) ≤ v v − v − v + 26. The maximum number offacets can be reached when h ( ℓ , ℓ ) vanishes. This occurs when ℓ = ℓ = 0. That is, when .2. The Lower Bound for k = 2 69the number of edges of G ( P ) and G ( Q ) is maximal. This concludes the proof of the upperbound of Theorem 4.1 for the special case k = 2. Corollary 4.4. The maximum number of facets can be attained only when the number ofedges of each of P and Q is maximal for the given number of facets. k = 2 Given two integers m ≥ n ≥ 4, we describe how to construct two polytopes in R with m and n facets respectively, such that the number of facets of their Minkowski sum isexactly 4 mn − m − n + 26, establishing the lower bound of Theorem 4.1 for the specialcase k = 2. More precisely, given i , we describe how to construct a skeleton of a polytope P i with i facets, 3 i − i − P m and P n , properly adjusted and oriented, is exactly 4 mn − m − n + 26.As in the previous sections we mainly operate in the dual space of Gaussian maps. However,the construction of the desired Gaussian maps described below is an involved task, since notevery arrangement of arcs of great circles embedded on the unit sphere, the faces of whichare convex and the edges of which are strictly less than π long constitutes a valid Gaussianmap. uv w Y We defer the treatment of the special case i = 4 to the sequel, andstart with the general case i ≥ 5. The figure to the right depicts theGaussian map of P . We use the subscript letter i in all notations X i to identify some object X with the polytope P i . For example,we give the Gaussian map G ( P i ) of P i a shorter notation G i , but inthis paragraph we omit the subscript letter in all notations for clarity.First, we examine the structure of the Gaussian map G of P to betterunderstand the structure of P . Let V and E denote the set of verticesand edges of G , respectively. Recall that the number of vertices, edges, and faces of G is i ,3 i − 6, and 2 i − 4, respectively. The unit sphere, where G is embedded on, is divided bythe plane y = 0 into two hemispheres H − ⊂ { ( x, y, z ) | y ≤ } and H + ⊂ { ( x, y, z ) | y > } .Three vertices, namely u , v , and w , lie in the plane x = 0. u is located very close to thepole (0 , , − i vertices) that lies in H + . v is locatedexactly at the opposite pole (0 , , w lies in H − very close to v . None of the remaining i − V \ { u, v, w } lie in the plane x = 0; they are all concentrated near the pole(0 , , − 1) and lie in H − . The edge uv , which is contained in the plane x = 0, is the onlyedge whose interior is entirely contained in H + . Every vertex in V \ { u, v, w } is connectedby two edges to v and w , respectively. These edges together with the edge uw , contained inthe plane x = 0, form a set of 2 i − E ′ = E \ { uv } . The length of each ofthe edges in E ′ is almost π , due to the near proximity of u , v , and w to the respective poles.It is easy to verify that if the polytope P is not degenerate, namely, its affine hull is3-space, then any edge of G ( P ) is strictly less than π long. Bearing this in mind, the maindifficulty in arriving at a tight-bound construction is forcing sufficient edges of the set E ′ ofthe Gaussian map of one polytope to intersect sufficient edges of the set E ′ of the Gaussianmap of the other polytope. The remaining pair of edges, one from each Gaussian map,0 Chapter 4. Exact Complexity of Minkowski Sums contributes a single intersection to the total count. As shown below, this is the best one cando in terms of intersections.Figure 4.2: The over-lay of G and G ′ ,where G ′ is G rotated ◦ about the Y axis. The number of facets of the Minkowski sum of P m and P n is max-imal, when the number of vertices in the overlay of G m and G n ismaximal. This occurs, for example, when one of them is rotated 90 ◦ about the Y axis, as depicted on the left for the case of m = n = 5.In this configuration, all the 2 m − E ′ m intersect all the2 n − E ′ n . All intersections occur in H − . In addition, theedge uv m intersects the edge uv n . The intersection point lies in H + exactly at the pole (0 , , m − n − 5) + 1 = 4 mn − m − n + 26. Adding theoriginal vertices of G ( P m ) and G ( P n ), yields the desired result.Next, we explain how P i , i ≥ G i above. The construction of P i is guided by a cylinder.All the vertices of P i lie on the boundary of a cylinder the axis ofwhich coincides with the Z axis. We start with the case i = 5, and show how to generalizethe construction for i > 5. The special case i = 4 is explained last. P Figure 4.3 shows various views of P . Recall that P has 6 vertices, denoted as v , v , . . . , v ,and 9 edges. We omit the subscript digit 5 in all the notations through the rest of thissubsection for clarity. Let v v . . . v n denote the face defined by the sequence of vertices v , v , . . . , v n on the face boundary. The projection of all vertices onto the plane z = 0 lieon the unit circle. As a matter of fact, the entire face f v = v v v v lies in the plane z = 0.It is mapped under G to the vertex v = G ( f v ). Similarly, the faces f u = v v v v and f w = v v v v are mapped under G to the vertices u = G ( f u ) and w = G ( f w ), respectively.Consider the projection of the vertices onto the plane z = 0 best seen in Figure 4.3(b).Once the projection v ′ of v is determined as explained below, v is placed exactly on thebisector of ∠ v ′ ov . The vertices v , v , and v are the reflection of the vertices v , v , and v respectively through the plane x = 0.Two parameters govern the exact placement of v (and v ). One is the size of the exterior-dihedral angle at the edge v v , denoted as α , that is, the length of the geodesic-segmentthat is the mapping of the edge vw of G . This angle is best seen in Figure 4.3(c). Notice,that the Z axis is scaled up for clarity, and the angle in practice is much smaller. The otherparameter is the size of the angle β = ∠ v ′ ov ′ , where v ′ and v ′ are the projections of v and v respectively onto the plane z = 0. This is best seen in Figures 4.3(b) and (e). Given m and n , these angles for each of P m and P n depend on both m and n . For large values of m and n the values of α and β should be small. For example, setting α = β = 10 ◦ is sufficientfor the case m = n = 5 depicted in Figure 4.2. The actual setting is discussed below afterthe description of the general case i > .2. The Lower Bound for k = 2 71 v v v v v v XYZ o oα X YZ (a) (b) (c) o v v v v v v XYZ v v v v v v o X Y Z (d) (e) (f) Figure 4.3: Different views of P . (a) and (d) are perspective views, while (b), (c), (e), and (f) areorthogonal views. Notice that the Z axis is scaled up for clarity. P i , i ≥ o v j v j − v j v j v j v j v j XYZ We construct a polytope, such that two facets are visiblewhen looked at from z = ∞ , and i − z = −∞ . First, we place the projection ofall vertices onto the plane z = 0 along the unit circle, anddenote the projection of a vertex v as v ′ . The projection ofthe vertices v j , v j , v j , v j , v j , and v j , where j = 0, j = ⌊ ( i − / ⌋ , j = ⌊ ( i − / ⌋ + 1, j = i − j = ⌊ (3 i − / ⌋ ,and j = ⌊ (3 i − / ⌋ + 1, are placed at the same locationsas those of the corresponding vertices of P , as depicted onthe right. The projection of the remaining vertices are placed on the arcs \ v ′ j , v j , \ v j , v j , \ v j , v j , and \ v j , v ′ j in cyclic order.The angle γ = ∠ v j ov j − is another parameter that governs the final configuration of P i . Once the placement of the projection of v j − is determined, the projections of thevertices v j +1 , v j +2 , . . . , v j − are arbitrarily spread along the open arc \ v j , v j − . The vertexplacement along the arc \ v ′ j , v j must be a symmetric reflection of the vertex placement alongthe arc \ v j , v j . This guarantees that all the quadrilateral facets are planar. Similarly, the2 Chapter 4. Exact Complexity of Minkowski Sums o v v v v v v v v v v v v v v v v XYZ o v v v v v v v v v v v v v v v v v v XYZ (a) (b)Figure 4.4: (a) An orthogonal view of P . (b) An orthogonal view of P . vertex placement along the arc \ v j , v j is a symmetric reflection of the vertex placement alongthe arc \ v j , v ′ j . For large values of m and n the angle γ should be small as explained below,implying that the projection of the vertices are concentrated near v j and v j , (which lie onthe bisectors of ∠ v j ov ′ j and ∠ v j ov ′ j , respectively). Figure 4.4 depicts the cases i = 10, and i = 11. In these examples we force a regular placement, which is sufficient in many cases.As in the case of i = 5, the face f vi = v , v , . . . , v i − , represented by the vertex v i of G i ,lies in the plane z = 0. The exterior-dihedral angle α at the edge v j v j is made small, sothat the vertex w i of G i representing the adjacent face f wi = v j , v j +1 , . . . , v i − , v , is keptin close proximity to v i .Given m and n , three parameters per polytope listed below govern the final configurationsof P m and P n .1. the exterior-dihedral angle α at the edge v j v j ,2. the angle β = ∠ v ′ j ov ′ j , and3. the angle γ = ∠ v j ov j − .The settings of these angles must satisfy certain conditions, which in turn enable all thenecessary intersections of edges in the Gaussian map of the Minkowski sum. We denote theface v j +1 v j v j v j − adjacent to f u by f x . The vertex x = G ( f x ) is the nearest vertex to u .The y -coordinate of the vertex w n must be greater than the y -coordinate of the edge xv m at z = 0 in P m ’s coordinate system. Similarly, the y -coordinate of the vertex w m must begreater than the y -coordinate of the edge xv n at z = 0 in P n ’s coordinate system. This isbest seen in Figure 4.5(c). The values of the y -coordinates of w n and w m are simply sin( α n )and sin( α m ), respectively. The value of the y -coordinate of the edge xv m at z = 0 howeverdepends on all the three parameters α m , β m , and γ m . Similarly, the y -coordinate of the edge xv n at z = 0 in the respective coordinate system depends on α n , β n , and γ n . Instead ofderiving an expression that directly evaluates these y -coordinates, we suggest an iterativeprocedure that decreases the angles at every iteration until the conditions above are met,and argue that this procedure eventually terminates, because at the limit, we are back atthe case where m = n = 5, for which valid settings exist. The rotation of, say P n , is performed about the Y axis. Thus, it has no bearing on y -coordinates. .2. The Lower Bound for k = 2 73(a) (b) (c) (d)Figure 4.5: (a) The Minkowski sum M , = P ⊕ P ′ , where P ′ is P rotated ◦ about the Y axis. (b) The Gaussian map of M , looked at from z = ∞ . (c) A scaled up view of the Gaussianmap of M , looked at from z = ∞ . (d) The Gaussian map of M , looked at from y = −∞ . P Recall that P has 4 facets, 6 edges, and 4 vertices. Therefore, it cannotbe constructed according to the prescription provided in the previoussection. Applying the same principles though, we place two verticesof G near the pole (0 , , − , , H + . The other three edges that connect vertices nearopposite poles mostly lie in H − . They form a set of 2 i − E ′ . Thelength of every edge in E ′ is almost π . In contrast to the case i ≥ 5, two out of the threeedges in E ′ cross the plane y = 0. Namely, small sections of them lie in H + . As in the case i > 4, one edge, the lightly shaded one, is entirely contained in H + . o v v v v XYZ We construct P , such that the two facets f = v v v and f = v v v are visible when looked at from z = ∞ , and whenlooked at from z = −∞ , the remaining two facets f = v v v and f = v v v are visible. As depicted on the right, theprojection of all four vertices onto the plane z = 0 lie on theunit circle. The vertices v and v lie on the plane z = 0, andthe vertices v and v lie in a parallel plane. The distancebetween the planes is small to form small exterior-dihedralangles at the edges v v and v v .As in the general case, two parameters govern the exact placement of v , v , and v . Oneis the size of the exterior-dihedral angle at the edge v v . The other parameter is the size ofthe angle ∠ v ov ′ , where v ′ is the projection of v onto the plane z = 0. The sizes of theseangles are determined by the same rationale as in the general case.This concludes the proof of the lower bound of Theorem 4.1 for the special case k = 2.4 Chapter 4. Exact Complexity of Minkowski Sums Let P , P , . . . , P k be a set of k polytopes in R , such that the number of facets of P i is m i for i = 1 , , . . . , k . In this section we present a tight bound on the number of facets of theMinkowski sum M = P ⊕ P ⊕ . . . ⊕ P k generalizing the arguments presented above for k = 2. Figure 4.6: Theoverlay of the Gaus-sian maps of threetetrahedra rotatedabout the Y axis ◦ , ◦ , and ◦ ,respectively. Given k positive integers m , m , . . . , m k , such that m i ≥ 4, we de-scribe how to construct k polytopes in R with corresponding numberof facets, such that the number of facets of their Minkowski sum isexactly P ≤ i 5) + (cid:0) k (cid:1) + P ki =1 m i . More precisely,given i , we describe how to construct a skeleton of a polytope P i with i facets, 3 i − i − M = P ⊕ P ⊕ . . . ⊕ P k of the k polytopesproperly adjusted and oriented is exactly the expression above. We usethe same construction described in Section 4.2.The number of facets in the Minkowski sum of P i , i = 1 , , . . . , k is maximal, when the number of vertices in the overlay of G i , i =1 , , . . . , k is maximal. This occurs, for example, when G i is rotated180 ◦ ( i − /k about the Y axis for i = 1 , , . . . , k , as depicted on theleft for the case of m = m = m = 4. (Recall, that E ′ i refers tothe set of edges that span the lowest hemispheres, and its cardinalityis smaller than cardinality of E by one.) In this configuration, all the2 m i − E ′ i intersect all the 2 m j − E ′ j , for 1 ≤ i < j ≤ k . All intersectionsoccur in H − . In addition, the edge uv m i intersects the edge uv m j for 1 ≤ i < j ≤ k . Theseintersection points lie in H + near the pole (0 , , P ≤ i 5) + (cid:0) k (cid:1) . Adding the original vertices of G ( P i ) , i = 1 , , . . . , k ,yields the bound asserted in Theorem 4.1. We can apply the special case k = 2 of Theorem 4.1 to obtain f ( m , m , . . . , m k ) ≤ f ( m , f ( m , m , . . . , m k )) ≤ m f ( m , m , . . . , m k ) − m − f ( m , m , . . . , m k ) + 26 ≤ k k Y i =1 m i + . . . However, we can apply a technique similar to the one used in Section 4.1 and improve thisupper bound, but first we must extend Lemma 4.3. .3. Maximum Complexity of Minkowski Sums of Many Polytopes Lemma 4.5. Let G , G , . . . , G k be a set of k Gaussian maps of convex polytopes, and let G be their overlay. Let f i denote the number of faces of G i , and let f denote the number offaces of G . Then, f ≤ P ≤ i 2, and the maximal number of western-most corners in theoverlay of some G i and G j is f i f j − f ≤ X ≤ i 2) + 2 . This corresponds to summing the western-most corners appearing in the overlay of all pairsof Gaussian maps, and subtracting ( k − 2) times the western-most corners appearing in alloriginal Gaussian maps, since each of them appeared ( k − 1) times in the first sum. Finally,we have: X ≤ i 2) + 2 = X ≤ i 2. Applying Lemma 4.5 and respecting v = P ≤ i ≤ k v i + v x , we can bound v x from above: v x ≤ X ≤ i 2) + k X i =1 ( v i − e i ) − . (4.2)6 Chapter 4. Exact Complexity of Minkowski Sums According to Corollary 4.4 the maximum number of facets of the Minkowski sum of twopolytopes is attained when the number of edges of each summand is maximal. We need toestablish a similar property for the general case. Generalizing the derivation procedure inSection 4.1, we introduce k non-negative integers ℓ i , i = 1 , , . . . , k , such that e i = 3 v i − − ℓ i and f i = 2 v i − − ℓ i . Substituting e i in (4.2) we get: v x ≤ X ≤ i 2) + k X i =1 ( v i − v i + 6 + ℓ i ) − ≤ X ≤ i 5) + k X i =1 ℓ i + k − k . (4.3)Substituting f i in (4.3) we get: v x ≤ X ≤ i 5) + (cid:18) k (cid:19) + h ( ℓ , ℓ , . . . , ℓ k ) , where h ( ℓ , ℓ , . . . , ℓ k ) = X ≤ i 5, which in turn implies that h ( ℓ , ℓ , . . . , ℓ k ) ≤ 0. Thus, we have: v x ≤ X ≤ i 5) + (cid:18) k (cid:19) . (4.4)We conclude that the exact maximum number of facets of the Minkowski sum of k polytopes cannot exceed P ≤ i 5) + P ≤ i ≤ k m i + (cid:0) k (cid:1) , which completes theproof of Theorem 4.1. For example, the exact maximum number of facets of the Minkowskisum of k tetrahedra is 5 k − k . For k = 2 the expression evaluates to 18. ivide each difficulty intoas many parts as is feasibleand necessary to resolve it. Rene Descartes Assembly Planning Assembly partitioning with an infinite translation is the application of an infinite translationto partition an assembled product into two complementing subsets of parts, referred to assubassemblies, each treated as a rigid body. We present an exact implementation of anefficient algorithm based on the general framework described in [HLW00, WL94] to obtainsuch a motion and subassemblies given an assembly of polyhedra in R . As usual, we do notassume general position. Namely, we handle degenerate input, and produce exact results. Asoften occurs, motions that partition a given assembly or subassembly might be isolated in theinfinite space of motions. Any perturbation of the input or of intermediate results, causedby, for example, imprecision, might result with dismissal of valid partitioning-motions. Inthe extreme case, there is only a finite number of valid partitioning motions, as occurs in theassembly shown in Figure 5.1, no motion may be found, even though such exists. Properhandling requires significant enhancements applied to the original algorithmic framework.The implementation is based on software components introduced in Chapters 2 and 3. They(a) (b) (c) (d)Figure 5.1: (a) The Split Star assembly, and (b),(c), and (d) the Split Star six parts divided into threepairs of symmetric parts. The six parts are named according to their color R (ed), G (reen), B (lue), P (urple), Y (ellow), and T (urquoise). Chapter 5. Assembly Planning paved the way to a complete, efficient, and concise implementation. Assembly planning is the problem of finding a sequence of motions that transform the ini-tially separated parts of an assembly to form the assembled product. The reverse order ofsequenced motions separates an assembled product to its parts. Thus, for rigid parts, assem-bly planning and disassembly planning refer to the same problem, and used interchangeably.In this chapter we concentrate on the case where the assembly consists of polyhedra in R and the motions are infinite translations.The motion space is the space of possible motions that assembly parts may undergo. Foreach motion in a motion space, a subset of parts of a given assembly may collide with adifferent subassembly, when transformed as a rigid body according to the motion. Pairs ofsubassemblies that collide constitute constraints. The motion space approach dictates theprecomputation of a decomposition of a motion space into regions, such that the constraintsamong the parts in the assembly are the same for all the motions in the same region. Allconstraints over a region are represented by a graph, called the directional blocking graph (DBG) [WL94]. The collection of all regions in a motion space together with their associatedDBGs, collectively referred to as the nondirectional blocking graph (NDBG), can be used toobtain assembly (or disassembly) sequences.The general framework and some of the techniques presented here have already beendescribed in a series of papers and reports published in the past mainly during the late 90’s.Halperin, Latombe, and Wilson made the connection between previously presented tech-niques that had used the motion space approach, and introduced a unified general frame-work [HLW00] at the end of the previous millennium. Only few publications related to thistopic appeared ever since, to the best of our knowledge, which creates a long gap in the timeline of the respective research. We certainly hope that the tools exposed in this chapter willhelp revive the research on algorithmic assembly planing, a research subject of considerableimportance. Moreover, we believe that the machinery presented here, together with otherrecent advances in the practice of computational-geometry algorithms, can more generallysupport the development of new and improved techniques in algorithmic automation [2].Solution to the assembly planning problem enables better feedback to designers. It pro-vides a design team with additional tools to asses a design, prior to the construction ofphysical mock-ups, and helps creating products that are more cost-effective to manufactureand maintain. This is emphasized in light of the strategy to “plan anywhere, build any-where” many Computer Aided Design (CAD) and Computer Aided Manufacturing (CAM)companies are trying to adopt. Assembly sequences are also useful for selecting assemblytools and equipment, and for laying out manufacturing facilities.We restrict ourselves to two-handed partitioning steps, meaning that we partition thegiven assembly into two complementing subsets each treated as a rigid body. Even for two-handed partitioning, if we allow arbitrary translational motions (and not restrict ourselvesto infinite translations) the problem is NP-hard [KK95]. The more general problem ofassembly sequencing, namely planning a total ordering of assembly operations that merge .1. Introduction infinitetranslations , is technically considerably more complex than partitioning with infinitesimalmotions . Although the latter may sound more general, as it handles infinitesimal transla-tions and rotations , it is far simpler to implement, since it deals only with constraints thatcan be described linearly. Thus, the problem can be reduced to solving linear programs.Indeed, there are several implementations for partitioning with infinitesimal motions (see,e.g., [GHH + 98, SS94]), but none that we are aware of dealing robustly with infinite trans-lations. The shortcoming of using infinitesimal motions only is that suggested disassemblymoves may not be extendible to practical finite-length separation motions.Infinite-translation partitioning was not fully robustly implemented until recently, in spiteof being more useful than infinitesimal partitioning, most probably due to the hardship ofaccurately constructing the underlying geometric primitives. What enables the solution thatwe present here, is the significant headway in the development of computational-geometrysoftware over the past decade, the availability of stable code in the form of the Computa-tional Geometry Algorithms Library ( Cgal ) in general and code for Minkowski sums andarrangements in particular [WFZH07b].The implementation presented in this chapter is based on the Arrangement on surface 2 package of Cgal ; see Chapter 2 for more details. The implementation uses in particulararrangements of geodesic arcs embedded on the sphere; see Section 2.6. The ability to ro-bustly construct such arrangements, and carry out exact operations on them using only(exact) rational arithmetic is a key property that enables an efficient implementation. Theimplementation exploits supported operations, and requires additional operations, e.g., cen-tral projection of polyhedra, which we implemented. We plan to make these new componentsavailable as part of a future public release of Cgal as well. We use the assembly depicted in Figure 5.1 as a running example throughout the chapter.The name “Split Star” was given to this shape by Stewart Coffin in one of his Puzzle Craftbooklet editions back in 1985. He uses the term puzzle to refer to any sort of geometricrecreation having pieces that come apart and fit back together. We use it as an assembly.He describes how to produce the actual pieces out of wood [Cof06], and suggests that they aremade very accurately. He observes that finding the solution requires dexterity and patience,when the pieces are accurately made with a tight fit. Even though the assembly seemsrelatively simple, this should come as little surprise, since the first partitioning motion isone out of only eight possible translations of four symmetric pairs of motions in oppositedirections associated with two complementing subassemblies of three parts each. Evidently,any automatic process that introduces even the slightest error along the way, will most likelyfail to compute the correct first motion in the sequence, and falsely claim that the assemblyis interlocked.0 Chapter 5. Assembly Planning The Split Star assembly has the assembled shape of thefirst stellation of the rhombic dodecahedron [Luk57], illus-trated atop the right pedestal in M. C. Escher’s Waterfallwoodcut [Boo82]. Its orthogonal projection along one ofits fourfold axes of symmetry is a square, while the Starof David is obtained when it is projected along one of itsthreefold axes of symmetry, as seen on the left; for more de-tails see [Cof06]. The assembly is a space-filling solid whenassembled. It consists of six identical concave parts. Eachpart can be decomposed into eight identical tetrahedra yielding 48 tetrahedra in total. Asmanufacturing the pieces requires extreme precision, it is suggested to produce the 48 iden-tical pieces and glue them as necessary. Each part can also be decomposed into three convexpolytopes — two square pyramids and one octahedron, yielding 18 polytopes in total. Thepartitioning described in this chapter requires the decomposition of the parts into convexpieces. The choice of decomposition may have a drastic impact on the time consumption ofthe entire process, as observed in a different study in R [AFH02], and shown by experimentsin Section 5.5. The rest of this chapter is organized as follows. The partitioning algorithm is describedin Section 5.2 along with the necessary terms and definitions. In Section 5.3 we provideimplementation details. Section 5.4 presents optimizations that are not discussed in thepreceding sections, some of which we have already implemented and proved to be useful. Wereport on experimental results in Section 5.5. The main problem we address in this chapter, namely, polyhedral assembly partitioning withinfinite translations , is formally defined as follows: Let A = { P , P , . . . , P n } be a set of n pairwise interior disjoint polyhedra in R . A denotes the assembly that we aim to partition.We look for a proper subset S ⊂ A and a direction ~d in R , such that S can be moved asa rigid body to infinity along ~d without colliding with the rest of the parts of the assembly A \ S . (We allow sliding motion of one part over the other. We disallow the intersection ofthe interior of two polyhedra.)We follow the NDBG approach [WL94], and describe it here using the general formulationand notation of [HLW00]. The motion space in our case, namely the space of all possiblepartitioning directions, is represented by the unit sphere S . Every point p on S defines thedirection from the center of S towards p . Every direction ~d is associated with the directedgraph DBG ( ~d ) = ( V, E ) that encodes the blocking relations between the parts in A whenmoved along ~d as follows: The nodes in V correspond to polyhedra in A ; we denote a nodecorresponding to the polyhedron P i by v ( P i ) ∈ V . There is an edge directed from v ( P i )to v ( P j ), denoted e ( P i , P j ) ∈ E , if and only if the interior of the polyhedron P i intersects .2. The Partitioning Algorithm P j when P i is moved to infinity along the direction ~d , and P j remains stationary.The key idea behind the NDBG approach is that in problems such as ours, where thenumber of parts is finite, and any allowable partitioning motion can be described by a smallnumber of parameters, there is only a relatively small (polynomial) number of distinct DBGsthat need to be constructed in order to detect a possible partitioning direction. Stateddifferently, the motion space can be represented by an arrangement comprising a finitenumber of cells each assigned with a fixed DBG. Once this arrangement is constructed, weconstruct the DBG over each cell of the arrangement, and check it for strong connectivity.A DBG that is not strongly connected is associated with a direction, or a set of directions incase the cell is not a singular point, that partition the given assembly. The desired movablesubset S ⊂ A is a byproduct of the algorithm that checks for strong connectivity. If all theDBGs over all the cells of the arrangement are strongly connected, we conclude that theassembly is interlocked , as a subset of the parts in A that can be separated from the rest ofthe assembly by an infinite translation does not exist.Next we show how to construct the motion-space arrangement and compute the DBGover each one of the arrangement cells. Each ordered pair of distinct polyhedra < P i , P j > defines a region Q ij on S , which is the union of all the directions ~d , such that when P i ismoved along ~d its interior will intersect the interior of P j . How can we effectively computethis region? Let M ij denote the Minkowski sum P j ⊕ ( − P i ) = { b − a | a ∈ P i , b ∈ P j } . Weclaim that the central projection of M ij onto S is exactly Q ij . Lemma 5.1. A direction ~d is in the interior of the central projection of M ij onto S if andonly if when P i is moved along ~d its interior will intersect the interior of P j .Proof. Let ~d be some direction in the central projection of M ij onto S . In other words, thereexists a point m ∈ M ij , such that m = s · ~d , for some positive scalar s . As m is in M ij , thereexist two points p i ∈ P i and p j ∈ P j , such that m = p j − p i . Thus, p j = p i + s · ~d , meaningthat the point p i intersects p j when moved along ~d . A similar argument can be used to showthe converse.Next, we describe how, given two polyhedra P i and P j , we compute the region Q ij , usingrobust and efficient hierarchy of building blocks, which we have developed in recent years.The existing tools that we use are (i) computing the arrangement of spherical polygons[BFH + 07, FSH08b, FSH08a, WFZH07b], and (ii) construction of Minkowski sums of convexpolytopes [BFH + 07, FH07, FSH08b, FSH08a]. We also need some extra machinery, asexplained below.Assume P i is given as the union of a collection of (not necessarily disjoint) convexpolytopes P i , P i , . . . , P im i , and similarly P j is given as the collection of convex polytopes P j , P j , . . . , P jm j . It is easily verified that M ij = S k =1 ,...,m i ,ℓ =1 ,...,m j P jℓ ⊕ ( − P ik ). So we com-pute the Minkowski sum of each pair P jℓ ⊕ ( − P ik ), and centrally project it onto S . Finally,we take the union of all these projections to yield Q ij .There are several ways to effectively compute the central projection of a convex polyhe-dron C (one of the polytopes M ijkℓ = P jℓ ⊕ ( − P ik )) from the origin onto S . We opted for the2 Chapter 5. Assembly Planning following. An edge e of C is called a silhouette edge , if the plane π through the origin and e is tangent to C at e ; namely, π intersects C in e only. We assume for now that no tangentplane contains a facet of C ; we relax this assumption in Section 5.3.5, where we provide adetailed description of the procedure. We traverse the edges of C till we find a silhouetteedge e . One can verify that the silhouette edges form a cycle on C . We start with e , andsearch for a silhouette edge adjacent to e . We proceed in the same manner, till we end updiscovering e again. Projecting this cycle of edges onto S is straightforward.All the boundaries of the regions Q ij form an arrangement of geodesic arcs on the sphere.We traverse the motion-space arrangement in say a breadth-first fashion. For the first facewe check which ones of the regions Q ij contain it. We construct the corresponding DBGand check it for strong connectivity. If it is not strongly connected, we stop and announcea solution as described above. Otherwise we move to an adjacent feature of the currentface. During this move we may have stepped out from a set of regions Q ij , and may havestepped into a new set of regions Q ij . We update the current DBG according to the regionswe left or entered, test the new DBG for strong connectivity, and so on till the traversal ofall the arrangement cells is completed. Notice that it is important to visit also vertices andedges of the arrangement, since the solution may not lie in the interior of a face. Indeed, inour Split Star example, solutions are on vertices of the arrangement. Without careful exactconstructions, such solutions could easily be missed. The implementation of the assembly-partitioning operation consists of eight phases listedbelow. They all exploit arrangements of geodesic arcs embedded on the sphere [BFH + Arrangement on surface 2 package of Cgal already sup-ports the construction and maintenance of such arrangements, the computation of union offaces of such arrangements, the construction of Gaussian maps of polyhedra represented bysuch arrangements, and the computation of their Minkowski sums. It provides the groundfor efficient and generic implementation of the remaining required operations, such as centralprojection.1. Convex Decomposition Sub-part Gaussian map construction Sub-part Gaussian map reflection Pairwise sub-part Minkowski sum construction Pairwise sub-part Minkowski sum projection Pairwise Minkowski sum projection Motion-space construction Motion-space processing .3. Implementation Details The Split Star six parts decomposed into three convex sub-parts each. The partitioning process is implemented as a free function that accepts as input anordered list of polyhedra in R , which are the parts of the assembly. Each part is representedas a polyhedral mesh in R ; see Section 3.2.1 for a definition. We proceed with a detaileddiscussion of the implementation of each phase.We deal below with various details that are typically ignored in reports on geometricalgorithms (for example, under the general position assumption). However, in assemblyplanning, or more generally in movable-separability problems [Tou85] in tight scenarios,much of the difficulty shifts exactly to these technical details and in particular to handlingdegeneracies. This is especially emphasized in Phases 5 and 6 (Subsections 5.3.5 and 5.3.6respectively), but prevails throughout the entire section. We decompose each concave part into convex polyhedra referred to as sub-parts. The outputof this phase is an ordered list of parts, where each part is an ordered list of convex sub-parts represented as polyhedral surfaces. Each polyhedral surface is maintained as a Cgal Polyhedron 3 [Ket07a] data structure, which consists of vertices, edges, and facets andincidence relations on them [Ket99]. A part that is convex to start with is simply convertedinto an object of type Polyhedron 3 .A new package of Cgal that supports convex decomposition of polyhedra has beenrecently introduced [Hac07], but has not become publicly available yet. As we aim for afully automatic process, we intend to exploit such components, once they become available,and study their impact. For the time being we resorted to a manual procedure. A simpledecomposition of the Split Star parts used in the running example is illustrated in Figure 5.2. We convert each sub-part represented as a polyhedral surface into a Gaussian map repre-sented as an arrangement of geodesic arcs embedded on the sphere, where each face f of thearrangement is extended with the coordinates of its associated primal vertex v = G − ( f ),resulting with a unique representation; see Section 3.1 for the exact definition of Gaussian4 Chapter 5. Assembly Planning R R R B B B − R − R − R − B − B − B Figure 5.3: Samples of the Gaussian maps of sub-parts of the Split Star assembly. The bottom rowcontains the reflections of the Gaussian maps at the top row. maps and see Section 3.2 the exact procedure to construct one from a polytope.The output of this phase is an ordered list of parts, where each part is an ordered list ofthe Gaussian maps of the convex sub-parts. Figure 5.3 depicts the Gaussian maps of six ofthe 18 polytopes that comprise the set of sub-parts of our Split Star assembly. We reflect each sub-part P ik through the origin to obtain − P ik . This operation can be per-formed directly on the output of the previous phase, reflecting the underlying arrangementsof geodesic arcs embedded on the sphere, which represent the Gaussian maps, while handlingthe additional data attached to the arrangement faces. As a matter of fact, this phase canbe reduced as part of an optimization discussed in Section 5.4.The output of this phase is an ordered list of parts, where each part is an ordered list ofGaussian maps of the reflected convex sub-parts. Figure 5.3 depicts the Gaussian maps ofsix of the 18 polytopes that comprise the set of reflected sub-parts of the Split Star example. Construct Pairwise Sub-partMinkowski Sums for i in { , , . . . , n } for j in { , , . . . , n } if i == j continuefor k in { , , . . . , m i } for ℓ in { , , . . . , m j } M ijkℓ = P jℓ ⊕ ( − P ik ) We compute the Minkowski sums of the pairwisesub-parts and reflected sub-parts. Aiming for an ef-ficient output sensitive algorithm, the constructionof an individual Minkowski sum from two Gaussianmaps represented as two arrangements respectivelyis performed by overlaying the two arrangements.When the overlay operation progresses, new vertices,edges, and faces of the resulting arrangement are cre-ated based on features of the two operands. When anew feature is created its attributes are updated. There are ten cases that must be handled;see Sections 3.2.2 for details. The Arrangement on surface 2 package conveniently supportsthe overlay operation allowing users to provide their own version of these ten operations. Theoverlay operation is exploited below in several different variants of arrangements of geodesic .3. Implementation Details R ⊕ ( − G ) R ⊕ ( − B ) G ⊕ ( − R ) G ⊕ ( − B ) B ⊕ ( − R ) B ⊕ ( − G )Figure 5.4: Samples of the pairwise Minkowski sums of sub-parts of the Split Star assembly. Themiddle row contains six Minkowski sums. The top row contains the corresponding Gaussian maps. Thebottom row contains the corresponding central projection of the Minkowski sums on S . arcs embedded on the sphere. Each application requires the provision of a different set ofthose ten operations.The output of this phase is a map from ordered pairs of distinct indices into lists ofMinkowski sums represented as Gaussian maps. Each ordered pair < i, j >, i = j is associatedwith the list of Minkowski sums { M ijkℓ | k = 1 , , . . . , m i , ℓ = 1 , , . . . , m j } . In case of our SplitStar the map consists of 30 entries that correspond to all configurations of ordered distinctpairs of parts. Each entry consists of a list of nine Minkowski sums, that is, 270 Minkowskisums in total. Project Pairwise Sub-partMinkowski sums for i in { , , . . . , n } for j in { , , . . . , n } if i == j continuefor k in { , , . . . , m i } for ℓ in { , , . . . , m j } Q ijkℓ = project( M ijkℓ ) We centrally project all pairwise sub-part Minkowskisums computed in the previous phase onto thesphere. Each projection is represented as an arrange-ment of geodesic arcs on the sphere, where each cell c of the arrangement is extended with a Boolean flagthat indicates whether all infinite rays emanatingfrom the origin in all directions ~d ∈ c pierce the cor-responding Minkowski sum. As the Minkowski sumsare convex, their spherical projections are sphericallyconvex.Given a convex Minkowski sum C , we distinguish between four different cases as follows:1. The origin is contained in the interior of a facet of C .2. The origin lies in the interior of an edge of C .3. The origin coincides with a vertex of C .6 Chapter 5. Assembly Planning 4. The origin is separated from C .Computing the projection of a convex polytope C can be done efficiently using dedicatedprocedures that handle the four cases, respectively. Recall that C is represented as a Gaussianmap, which is internally represented as an arrangement of geodesic arcs embedded on thesphere. We traverse the vertices of the arrangement. For each vertex v we extract itsassociated primal facet f = G − ( v ). We dispatch the appropriate computation based on therelative position of the origin with respect to the supporting plane to f , and the supportingplane to adjacent facets of f . If the origin is contained in the interior of a facet f of C , the projection ofthe silhouette of C is a great (full) circle that divides the sphere into two hemispheres. Thenormal to the plane that contains the great circle is identical to the normal to the supportingplane to f , easily extracted from the arrangement representing the Gaussian map of C . The Arrangement on surface 2 package conveniently supports the insertion of a great circle,provided by the normal to the plane that contains it, into an arrangement of geodesic arcsembedded on the sphere.We omit the implementation details of the following two cases, and proceed to the generalcase. If the origin is separated from C , we traverse all edges of C until we find a silhouetteedge characterized as follows: Let v s and v t be the source and target vertices of some edge e inthe arrangement representing the Gaussian map of C , and let f s = G − ( v s ) and f t = G − ( v t )be their associated primal facets, respectively. e is a silhouette edge, if and only if, the originis not in the negative side of the supporting plane to f s and not in the positive side of thesupporting plane to f t . We start with the first silhouette edge we find, and search for anadjacent silhouette edge in a loop, until we rediscover the first one. Weproject only the target vertices of significant silhouette edges, and con-nect consecutive projections using arcs of great circle. Let e and e ′ beadjacent silhouette edges. We skip e , if the projections of e and e ′ lie onthe same great circle. For example, all but the last adjacent silhouetteedges incident to a facet supported by a plane that contains the origin areredundant, as illustrated in the figure above. Here we skip e , e , and e ,and project the target vertex of e .The output of this phase is a map from ordered pairs of distinct indices into lists ofarrangements as described above. Each ordered pair < i, j >, i = j is associated with the listof central projections of the pairwise Minkowski sums of P j ’s sub-parts and the reflectionthrough the origin of P i ’s sub-parts. For each pair of distinct parts P i and P j we compute the union of projections of the pairwiseMinkowski sums of all sub-parts of part P j and reflections of all sub-parts of part P i .The output of this phase is a map from ordered pairs of distinct indices into arrangements.Each ordered pair < i, j >, i = j is associated with a single arrangement extended as describedabove, that represents the central projection Q ij of M ij = P j ⊕ ( − P i ). .3. Implementation Details for i in { , , . . . , n } for j in { , , . . . , n } if i == j continue Q ij = ∅ for k in { , , . . . , m i } for ℓ in { , , . . . , m j } Q ij = Q ij ∪ Q ijkℓ We exploit the overlay operation in this phase thesecond time throughout this process, this time in aloop. Given two distinct parts P i and P j we traverseall projections in the set { Q ijkℓ | k = 1 , , . . . , m i , ℓ =1 , , . . . , m j } , and accumulate the result in the ar-rangement Q ij . As mentioned in Section 5.3.4, whenthe overlay operation progresses, new vertices, edges,and faces of the resulting arrangement are created.When a new face f is created as a result of the over-lay of a face g in some projection Q ijkℓ , and a face in the accumulating arrangement, theBoolean flag associated with f , which indicates whether all directions ~d ∈ f pierce M ij , isturned on, if ~d pierces M ijkℓ , that is, if the flag associated with the face of g is on.The intermediate result of this step are arrangements with potentially redundant edgesand vertices. It is desired (but not necessary) to remove these cells, as it reduces the timeconsumption of the succeeding operations, which is directly related to the complexity of thearrangements. It has even a larger impact when the optimization described in Section 5.4 isapplied, as the optimization decreases the number of preceding operation at the account ofslightly increasing the number of succeeding operations. We remove all edges and verticesthat are in the interior of the projection, that is all marked edges and vertices. We alsoremove spherically collinear vertices on the boundary of the projection, the degree of whichdecreased below three, as a result of the redundant-edge removal.(a) (b) (c) (d) (e) (f)Figure 5.5: Peg-in-the-hole Minkowski sum projections. (a), (b), (c), (d), and (e) are the sub-partprojection. (f) is the union of the former. Cgal also supports Boolean operations applied to general polygons and in particularthe union operation. However, it consumes and produces regularized general polygons; seeSection 2.7.1. This regularization operation is harmful in the realm of assembly planning.Therefore, we work directly on the cells of the arrangements Q ij .Quite often the projection contains isolated vertices and edges,as occurs in the example depicted on the right, referred to as“peg-in-the-hole”. Here the assembled product is translucentlyviewed from two opposite directions. The blue part is stationaryand is decomposed into five sub-parts. Figure 5.5 illustrates the corresponding five pairwiseMinkowski sum projections, and their union. The complement of the union consists of asingle isolated vertex. The generic code supports point sets bounded by algebraic curves embedded on parametric surfacesreferred to as general polygons. Chapter 5. Assembly Planning B ⊕ ( − R ) R ⊕ ( − B ) Recalling our Split Star assembly, the projection of the Minkowskisum of the red part and the reflection of the blue part, and itsreflection, that is, the projection of the Minkowski sum of the bluepart and the reflection of the red part are depicted on the left. We compute a single arrangement that represents the motion space,where each cell c of the arrangement is extended with a DBG. Weuse the adjacency-matrix storage format provided by Boost [3] torepresent each DBG. Recall that for a graph with n vertices suchas ours, an n × n matrix is used, where each element a cij of a DBG associated with cell c isa Boolean flag that indicates whether part P i collides with part P j when moved along anydirection ~d ∈ c . In particular we use the adjacency matrix class. It implements the BoostGraph Library ( BGL ) [SLL02] interface, which supports, among the other, easy insertionsof new edges into existing graphs. Handling large assemblies with sparse blocking relationsmay require different representations of DBGs to reduce memory consumption.We exploit the overlay operation in this phase the third time similar to its applicationin Section 5.3.6. We traverse all central projections in the set { Q ij | ≤ i < j ≤ n } ,and accumulate the result in the final motion-space arrangement. As mentioned above inSection 5.3.4, when the overlay operation progresses, new vertices, edges, and faces of theresulting arrangement are created. When a new cell c is created as aresult of the overlay of a face g in some projection Q ij , and a cell inthe accumulating arrangement, the DBG associated with c is updated.That is, if the flag associated with g is turned on, we insert an edgebetween vertex i and vertex j into the DBG associated with c .Depicted on the left is the motion-space arrangement computed byour program for the Split Star assembly. Table 5.1: The Split Star validpartitioning directions and cor-responding subassemblies. Direction Subset − − − GBT − − 1, 1 RBT − 1, 1, − GP T − 1, 1, 1 RP T 5. 1, − − GBY 6. 1, − 1, 1 RBY 7. 1, 1, − GP Y 8. 1, 1, 1 RP Y We traverse all vertices, edges, and faces of the motion-space arrangement in this order, and test the DBG associatedwith each cell for strong connectivity using the Boost globalfunction strong components() . This function computes thestrongly connected components of a directed graph using Tar-jan’s algorithm based on depth-first search (DFS) [Tar72]. Theset of constraints associated with a vertex v is a proper sub-set of the constraints associated with the edges incident to v .Similarly, the set of constraints associated with an edge e is aproper subset of the constraints associated with the two facesincident to e . Therefore, if the DBGs of all vertices are stronglyconnected, we terminate with the conclusion that the .4. Additional Optimization The reflection of the sub-parts through the origin as described in Section 5.3.3 has beennaively implemented. The computation is applied to the polyhedral-mesh representation ofeach sub-part. An immediate optimization calls for an application of the reflection operationdirectly on the arrangements that represent the Gaussian map. We are planning to intro-duce a generic implementation of the reflection operation that operates on any applicablearrangement. This operation alters the incidence relations between the arrangement featuresand their geometric embeddings. For each vertex, it negates its associated point, and invertsthe order of the halfedges incident to it. For each edge, it negates the associated curve.For each face, it inverts the order of the halfedges along its outer boundary. Similar to theoverlay operation (see Section 2.4.2), where the user can provide a set of ten functions, whichare invoked when new vertices, edges, and faces of the resulting arrangement are created,while the overlay operation progresses, the user can provide a set of three functions that areinvoked when a new vertex, halfedge, and face are created, while the reflection operationprogresses. Extended data associated with these types, such as a primal vertex associatedwith an arrangement face as in the case of an arrangement representing a Gaussian map,can easily be updated with the provision of an appropriate function.The trivial observation that P ⊕ ( − Q ) = − (( − P ) ⊕ Q ) leads to another optimization.Instead of reflecting all sub-parts in the set { P ik | i = 1 , , . . . , n, k = 1 , , . . . , m i } , we reflectonly the sub-parts in the set { P ik | i = 2 , . . . , n, k = 1 , , . . . , m i } , and compute only thepairwise sub-parts Minkowski sums in the set { M ijkℓ | ≤ i < j ≤ n, k = 1 , , . . . , m j , ℓ =1 , , . . . , m j } , their central projection, and the union of the appropriate projections to yieldthe set { Q ij | ≤ i < j ≤ n } . Then, we apply the reflection operation described aboveon each member of this set, and obtain the full set of projections { Q ij | i = 1 , , . . . , n, j =1 , , . . . , n } . The Boolean flag associated with a face of an arrangement that represents acentral projection is equal to the flag associated with its reflection. In other words, a faceof Q ij consists of directions that pierce M ij , if and only if, its reflection in Q ji consists ofdirections that pierce M ji .Phase 8 is purely topological. Thus, we do not expect the time consumption of this phaseto dominate the time consumption of the entire process for any input. Nevertheless, it mightbe possible to reduce its contribution to the total time consumption through efficient testingfor strong connectivity applied to all the DBGs [KMW98], exploiting the similarity betweenDBGs associated with incident cells. Recall, that the set of arcs in a DBG associated witha vertex v is a subset of the set of arcs associated with an edge incident to v . Similarly, theset of arcs in a DBG associated with an edge e is a subset of the set of arcs associated with0 Chapter 5. Assembly Planning Table 5.2: Time consumption (in seconds) of the execution of the eight phases applied to the Split Starassembly as input. Each one of the three rows refers to a different decomposition of the assembly. A — number of convex sub-parts per part. B — number of sub-part vertices per part. C — total numberof convex sub-parts. D — total number of Minkowski sums. E — total number of arrangements ofgeodesic arcs embedded on the sphere constructed throughout the process. A B C D E 1 2 3 4 5 6 7 8 e . The proposed technique reduces the cost from O ( n ) per DBG to anamortized cost of O ( n . ), where n is the maximum number of arcs in any blocking graph. Our program can handle all inputs. However, we limit ourselves to a representative set oftest cases, where we compare the impact of different decompositions on the process timeconsumption. The results listed in Table 5.2 were produced by experiments conducted ona Pentium PC clocked at 1.7 GHz. In all three test cases we use the Split Star assemblyas input. Naturally, in all three cases identical projections are obtained as the intermediateresults of Phase 6, hence the identical time consumption of the succeeding last two phases.Evidently, it is desired to decompose each part into as few as possible sub-parts with as smallas possible number of features. However, an automatic decomposition operation may requirelarge amount of resources to arrive at optimal or near optimal decompositions. Notice thatPhases 4 and 6 dominate the time complexity. This is due to the large number of geometricpredicates that must be evaluated during the execution of the overlay operation. eriousness is the onlyrefuge of the shallow. Oscar Wilde Conclusion and Future Work In this thesis we show how a complete implementation of extendible arrangements with a richset of operations enables a broad spectrum of robust applications that solve problems arisingin domains such as motion planning, assembly planning, and solid modeling. For example,we describe how arrangements embedded on two-dimensional surfaces can be efficiently usedto compute Minkowski sums of two polytopes in R , which in turn, and in conjunction withseveral other operations based on such arrangements, can be used to partition an assemblywith an infinite translation motion. The rest of this chapter is devoted to future prospectsrelated to our research topics. Constructing Minkowski Sums of polytopes in R has been successfully attempted in the past.We introduce a robust, yet efficient method. Table 3.5 shows that both our exact methodsoutperform the other exact methods. However, we believe that both of our methods, and infact all Cgal based methods have great potential for further improvements through futureoptimizations applied to the infrastructure of Cgal , as Cgal is an evolving project. Whilethe space consumption of the Cgm method is greater than the space consumption of thespherical Gaussian-map method, the table also reveals that the Cgm method is currentlymore efficient than its rival. We estimate that the gap will decrease, if not vanish, once alloptimizations for the Arrangement on surface 2 data-structure that are still pending areimplemented and enabled.We are constantly striving to improve the quality of our infrastructure, that is the Arrangement on surface 2 package. We have already identified few weak spots. Elimi-nating them will increase the genericity, extendibility, efficiency, and functionality of thepackage. We provide one example in each category.912 Chapter 6. Conclusion and Future Work There is a certain similarity between observers (see Section 2.4.4) and visitors (see Sec-tion 2.4.1), as typically each of their methods is triggered as a response to a certain event— a member of a pre-determined list of events. Technically, the main difference betweenthem is that observers define a one-to-many mapping between objects, while visitors definea one-to-one mapping. Recall, for example, that a single arrangement may register manyobservers, but it is only natural to relate a single visitor to a specific algorithmic frameworkin order to realize a certain concrete algorithm. Consequently, arrangement observers arederived from a common base class, and their methods must be virtual. This is how mod-ules, which are closed for modification, are extended using object-oriented programming.However, composability of such modules is limited, since independently produced modulesgenerally do not agree on common abstract interfaces from which supplied types must in-herit. In addition, when techniques from the object-oriented programming and the genericprogramming paradigms are mixed, they often clash. There are known methods to replacelists of objects, derived from a common base class, and linked during run time, with a listof syntactically unrelated objects concatenated during compile time (coded using a genericprogramming technique) [Ale01, Chapter 3]. Nevertheless, we would like to simultaneouslyenjoy the benefits of both the object-oriented and the generic programming paradigms, thatis, to enable the immediate production of composable modules that support dynamic poly-morphism. A very important research direction in our context is to explore these possibilities,perhaps pushing the limits of the C ++ programming language along the way. In many cases we need to associate values (called “properties”) with the vertices, thehalfedges, and the faces of the arrangement. In addition, it is often necessary to associatemultiple properties with each vertex, edge, or face; this is what Boost literature refers toas multi-parameterization. BGL [SLL02] graph classes have template parameters for vertexand edge “properties”. A property specifies the parameterized type of the property and alsoassigns an identifying tag to the property.There are various ways to associate properties with arrangement cells. One option isto extend the geometric types of the kernel, as the kernel is fully adaptable and exten-sible [HHK + Another option,is to extend the vertex, halfedge, or face records of the Dcel ; see Section 2.3.3. This may alsolead to excessive space-consumption, for example, when the data associated with a halfedgeis in fact tied to the embedded geometric curve. In this case the data, or at least a handleto the data, must be stored twice in both twin halfedges. A third option is to extend thecurve (or point) types defined by the geometry-traits class; see Section 2.3.1. But even thisoption leads to unjustified space-consumption, when only a limited number of arrangement They also differ in essence. While an observer typically implements a notifier, a visitor is usually acoherent part of an algorithm based on a fundamental and flexible framework. [GHJV95] It also requires nontrivial knowledge about the kernel structure and the techniques to extend it. .1. Arrangements on Two-Dimensional Surfaces ++ programming language. Point location is one of the most fundamental operations applied to arrangements; see Sec-tion 2.4.5. The contest between the different point-location strategies for arrangement em-bedded in the plane was settled in favor of the “landmark” variants for many types ofarrangements [HH08]. Unfortunately, at this point, these strategies cannot be applied on ar-rangements embedded on surfaces other than the plane. This is due to limitations that arisefrom the specific implementations of geometry traits that support these types of embeddedsurfaces, e.g., the geometry traits that handles geodesic arcs embedded on the sphere; seeSection 2.6. The problem should be attacked from both its ends. That is, we can try toenhance the geometry traits implementations, and add all ingredients required by the con-cept ArrangementLandmarksTraits 2 as defined today, and at the same time, try to comeup with alternative, perhaps similar, strategies that induce different requirements that areeasier to satisfy by the geometry-traits classes. Continuous and steady effort is made to further extend the arsenal of geometry-traits models,or simply to improve the existing ones; see Section 2.5.2. Supporting arrangements inducedby rich families of curves opens the door for numerous applications.The dominant bottleneck of all applications mentioned in this thesis is the applicationof the geometry operations implemented in the geometry-traits classes. Expediting theirperformance, while containing the growth of their memory footprint is always desired. Forexample, the arrangement package provides two traits classes that handle line segments.Both are parameterized by a geometric kernel; see Section 1.3.3. Segments defined by most Cgal kernels are represented only by their two endpoints. When a segment is split severaltimes, the bit-length needed to represent the coordinates of its endpoints may grow exponen-tially (see [FWH04] for a discussion), which may significantly slow down the computation.Therefore, one of the two traits classes represents a segment by its supporting line in addi-tion to its two endpoints. When the traits class computes an intersection point of two linesegments, it uses the coefficients of their supporting lines. When a segment is split at anintersection point, the underlying line of the two resulting sub-segments remains the same,and only its endpoints are updated. This traits class thus overcomes the undesired effect ofcascading intersection-point representation, as described above, at the account of a largermemory footprint. A similar idea can be applied to the traits class that handles geodesic arcsembedded on the sphere. An implementation of a geometry-traits that handles such arcsthat stores the projections of all the arrangement geometric features once calculated, and4 Chapter 6. Conclusion and Future Work retrieves them when subsequently needed, has great potential to reduce time consumptionat the price of growth in space consumption.Another direction is to expand existing implementations to meet the requirements ofthe various concepts in the hierarchy described in Chapter 2. Consider, for example, thegeometry traits class that handles geodesic arcs embedded on the sphere. Currently, itsupports the basic operations required to construct and maintain arrangements induced bysuch arcs. As mentioned in the previous section, it might be possible to enhance it to enablethe use of one or more of the variants of the “landmark” point-location strategies. Similarly,enabling Boolean set operations of spherical patches bounded by geodesic arcs embedded onthe sphere requires the provision of few additional operations by the traits class. Consider the following task: Given a set S = { S , S , . . . , S n } of two-dimensional surfacesin R , construct the three-dimensional arrangement A ( S ) induced by S . Fulfilling this taskin an efficient, complete, and robust manner has not been attempted yet, and is consideredchallenging. Implementing various strategies of point-location that operate on arrangementsin R and a plane-sweep and zone-construction frameworks for such a data structure isgreatly desired, but extremely ambitious. In analogy to two-dimensional arrangements, ageneric implementation of a plane-sweep framework can enable various operations, such asthe overlay of spatial subdivisions and ordinary and regularized Boolean set operations ofpoint sets bounded by general algebraic surfaces. These operations, in turn, can enable theimplementation of a multitude of applications.Arrangements embedded on two-dimensional surfaces can be used as building blocksin the implementation of a data structure that represents a three-dimensional arrange-ment [Wei07]. We can consider each surface S k separately, and construct the arrangement A k = A k ( S ) induced by intersection curves between S k and S \ S k embedded on S . Thearrangements A , A ,. . . , A n can subsequently be connected together to properly representthe spatial subdivision A ( S ). In some sense and to some extent this thesis attempts to close gaps between theoreticalresults and practical needs. It is not accidental that great parts of the thesis are closelyrelated to Cgal , as one of the goals of Cgal stated in the Introduction chapter of the thesisis to translate (theoretical) results into useful, reliable, and efficient programs for industrialand academic applications. Evidently, the Boolean set operations package of Cgal , whichis based on the Arrangement on surface 2 package (see Section 2.7.1) is one of the mostpopular packages among Cgal packages in the commercial market. Naturally, we would liketo continue improving this package. The problems addressed in the next two subsectionswere raised during the 3 rd Cgal User Workshop [dW08] [1]. .3. Boolean Set-Operations Input data of Boolean set operations, namely, a set of one or more polygons, used in real-world applications is occasionally corrupted, as it originates from measuring devices thatare susceptive to noise and physical disturbances. In some other cases, it contains manydegeneracies, which either disable computations based on fixed-precision arithmetic, or slowdown further computation using exact geometric computation. Invalid Data a b cde Figure 6.1: A relativelysimple polygon that isnot simple, given by itsboundary { a, b, c, d, b, e } . A polygon P is said to be simple (or Jordan), if the only points ofthe plane belonging to two edges of P are vertices of consecutiveedges P [23]. Namely, no two edges intersect, except for every twoconsecutive edges, which share one endpoint. A simple polygon istopologically equivalent to a disk. A polygon P is said to be weaklysimple , if the chain of the edges of P does not cross itself. A polygon P is said to be relatively simple , if it is weakly simple and the edgesof P do not intersect in their relative interior. Observe, that arelatively simple polygon, the vertices of which appear only once inthe boundary (the degree of each vertex is two), is simple. a b cdefgh Figure 6.2: A poly-gon with a hole givenby its outer boundary { a, b, c, d, e, f, g, h } andits hole { h, f, d, b } . Input data for any Boolean set operation represents points setthat may be bounded or unbounded, and may have holes. Items ofsuch input take the form of general polygons or general polygonswith holes with well-defined interiors and exteriors. A valid polygonmust be weakly simple (but not necessarily simple) and its verticesmust be ordered in counterclockwise direction around its interior. Avalid polygon with holes that represents a bounded point set, has anouter boundary represented as a weakly simple (but not necessarilysimple) polygon, the vertices of which are oriented in counterclock-wise order around its interior. In addition, the set may containholes, where each hole is represented as a simple polygon, the ver-tices of which are oriented in clockwise order around the interiorof the hole. Note that an unbounded polygon without holes spansthe entire plane. Vertices of holes may coincide with vertices of theboundary.As mentioned above, real-world data is often corrupted. Naturally, passing invalid poly-gons (polygons with holes, respectively) as input to a Boolean set operation must be avoided.Apparently, automatically “fixing” corrupted data, that is, converting invalid input polygonsor polygons with holes to valid ones, is not a simple task. Consider, for example, the selfintersecting star depicted in Figure 6.3(a). A point is considered inside the point set, if andonly if the number of counterclockwise turns the oriented boundary makes around the point,also called the winding number, is greater than zero. It can be efficiently calculated usingan Arrangement 2 data structure as follows. We extend each halfedge h with a Boolean flagthat indicates whether the winding number increases or decreases when we cross h , that is,6 Chapter 6. Conclusion and Future Work adb e c ace g ibdf h j ace g ibdf h j ace g ibdf h j (a) (b) (c) (d)Figure 6.3: (a) A self crossing polygon given by { a, b, c, d, e } . (b) The Arrangement 2 data structureconstructed from the polygon edges. (c) The Arrangement 2 data structure with updated face winding-numbers. (d) The Arrangement 2 data structure with internal edges removed. when we move from the face incident to h to the face incident to the twin of h . Observethat twin halfedges always have opposite flag values. We extend each face with an integerthat counts the winding number of every point in the face. Finally, we apply a Breath-FirstSearch (BFS) on all the arrangement faces starting from the unbounded face and updatingthe face counters as we cross halfedges. The BGL , for example, can be employed for thistask. Figure 6.3 illustrates the process. Once the winding number of every face is updated,we remove internal edges, as they are redundant, and convert the arrangement back into avalid polygon or polygon with holes.The problem becomes more complicated when the input polygons or polygons with holesviolate several validity properties at the same time. Naturally, converting corrupted datainto valid data consumes time. The challenge is to perform this task flawlessly and efficiently,while presenting a convenient interface to the user. Degenerate Data In computational geometry there are two main techniques to eliminate degeneracies and near-degeneracies. One is snap rounding [GM95, GGHT97, Hob99] and the other is controlledperturbation [HL04, MO06]. Both techniques aim at processing geometric data, e.g., curvesof the boundaries of general polygons, to yield new data that can be further robustly andmore efficiently processed, perhaps using only limited precision. Traditionally, snap roundinghas been applied to linear objects embedded in the plane [HP01, HP02], where it replacessets of linear segments with sets of polylines. It can be extended though to other types ofcurves in the plane, such as B´ezier curves [EKW07], or even other curve types embeddedon other surfaces that have a well defined grid (and perhaps other properties), such asgeodesic arcs embedded on the sphere. Controlled perturbation has even a larger spectrumof applicable platforms. With respect to our polygons, applying any one of the two techniquesabove may result with (partially) overlapping curves, which belong to two different polygons,respectively. In this case, we need to merge the incremental winding contributions of theoriginal curves mentioned above. .3. Boolean Set-Operations The Boolean set operations 2 package provides efficient operations that compute the reg-ularized union or the regularized intersection of a set of input polygons. There is no restric-tion on the polygons in the set; naturally, they may intersect each other. The package alsoprovides an efficient predicate that determines whether all polygons in a given set intersect.Figure 6.4: The unionof eight discs. There are at least three different methods to compute the unionof a set of polygons P , . . . P m . We can do it incrementally as fol-lows. At each step we compute the union of S k − = S k − i =1 P i with P k and obtain S k . A second option is to use a divide-and-conquer ap-proach. First, we divide the set of polygons into two subsets. Then,we compute the union of each subset recursively, and obtain thepartial results in S and S , respectively. Finally, we compute theunion S ∪ S . A third option aggregately computes the union of allpolygons. We construct an arrangement inserting the polygon edgesat once, utilizing the sweep-line framework, (see Section 2.4.1) andextract the result from the arrangement. Similarly, it is also possi-ble to aggregately compute the intersection T mi =1 P i of a set of inputpolygons.The incremental method is more efficient for a small (constant) size of input polygons,and the aggregate method is more efficient for sparse polygons with a relatively small numberof intersections. It is also possible to mix between the three methods, reaping the benefitsof them all. We would like to figure out what are the exact conditions that should be usedto determine when to use each method, or when to switch from one to another. The Cgal package Nef 2 supports ordinary set-operations on point sets in R [See07]. Thepoint-set operands and results are rectilinear polygonal model. Such a point set can bedefined by a finite set of open halfspaces, or obtained by set complement and set intersectionoperations. The package supports operations that consume and produce linear polygonsdefined by linear edges. The Boolean set operations package, on the other hand, supportsonly regularized set-operations, but the operations consume and produce general polygons.Recall, that a general polygon is a point set in R that has a topology of a polygon, butits boundary edges map to arcs of curves, which are not necessarily linear. Extending thepackage to support not only regularized operations, but also ordinary ones, will make ituseful for more applications; see, for example, Section 5.3.6. Boolean set operations are intuitive and therefore popular in many fields. For example, CSGis a representation model for solids based on Boolean set operations. Solids represented usingCSG result from Boolean set operations applied to elementary solids called primitives, e.g.,cubes, spheres, cones, and cylinders. A CSG solid is represented in a tree structure, where8 Chapter 6. Conclusion and Future Work the leaves represent primitives, and internal nodes represent Boolean operations.The Cgal package Nef 3 supports ordinary set-operations on point sets in R [HK07].Similar to the planar case, the package supports operations that consume and produce (lin-ear) polyhedra defined by (linear) halfspaces. Having the ability to construct and maintainarrangements in R , will enable the development of a new package, that will support eitherregularized, or even non-regularized, robust Boolean set-operations that consume and pro-duce general (curved) polyhedra in R , the boundaries of which are general surfaces. Manyfundamental problems in solid modeling, motion planning, and other domains, can benefitfrom such a package. One possible progression of the collision detection algorithm and its implementation de-scribed in Section 3.4 is a complete integrated framework that answers proximity queriesabout the relative placement of polytopes that undergo rigid motions including rotation .The framework may use either the spherical or the cubical Gaussian-map to represent poly-topes. The interface of these two data structures should be consolidated to allow rapidinterchanging.Some of the methods we foresee compute only those portions of the Minkowski sum thatare absolutely necessary, making our approach even more competitive. Briefly, instead ofcomputing the Minkowski sum of P and − Q , we walk simultaneously on the two respective Cgm ’s, producing one feature of the Minkowski sum at each step of the walk. Such a strategycould be adapted to the case of rotation by rotating the trajectory of the walk, keeping the Cgm of − Q intact, instead of rotating the Cgm itself. Figure 6.5: Environmentof the St. Peters Cathedralmapped on a teapot with a sil-ver material applied. Taken inRTHDRIBLE [22]. We have developed a new data structure that can be usedto construct and maintain cubical Gaussian-maps and com-pute Minkowski sums of pairs of polytopes represented by thenew data structure; see Chapter 3. The name of the datastructure is Cubical gaussian map 3 , and we are consideringintroducing a package by the same name to a prospective re-lease of Cgal . The implementation is generic and can beused for other purposes, where six planar subdivisions embed-ded on a unit cube and stitched properly at the edges of thecube is useful, for example Cubical Environment Mapping; see,e.g., [FvDFH95, Section 16.6].In computer graphics, reflection mapping is an efficientmethod of simulating a complex mirroring surface by meansof a precomputed texture image. The texture is used to storethe image of the environment surrounding the rendered ob- .6. Exact Complexity of Minkowski Sums Dimension Exact Maximum Complexity d = 2 m + m + . . . + m k d = 3 P ≤ i 5) + P ≤ i ≤ k m i + (cid:0) k (cid:1) The table to the rightsummarizes the knownexact bounds on themaximum complexity ofMinkowski sums of polytopes in terms of number of facets (( d − m and n facets can be as low as m , when the two polytopes have the samenumber of facets m and parallel features, but it is unknown what is the minimum exactcomplexity of Minkowski sums of polytopes that have only a limited number of parallelfeatures, or none at all.00 Chapter 6. Conclusion and Future Work Software Components, Libraries, and Packages A.1 Visual Simulation We have developed a toolkit, called Sgal (Scene Graph Algorithm Library), that supportsthe construction and maintenance of directed acyclic graphs that represent scenes and modelsin R . The toolkit includes, among the other, two interactive 3D applications. The firstdetects collisions and answers proximity queries for polytopes that undergo translation androtation. The second enables users to visualize a scene in an interactive manner. It parsesinput files that describe the scene in a degenerate yet extended VRML format [27]. Theformat is degenerate, as not all VRML features are supported (yet). However, it has beensignificantly extended as described below.Both applications are linked with (i) Cgal , (ii) a library that provides the exact rationalnumber-type, (iii) several Boost libraries, (iv) imagemagick [15] — a library that creates,edits, and composes bitmap images, and (v) internal libraries that construct and maintain3D scene-graphs, written in C++, and built on top of OpenGL . The internal code is dividedinto two libraries; SGAL — The main 3D scene-graph library and SCGAL — Extensions thatdepend on Cgal .We added several geometry nodes that represent polytopes using various sub-representations,such as Gaussian maps, a few other related nodes that handle coordinates, and many othermiscellaneous nodes that provide services, such as antialiasing and snapshoting. The de-scriptions of some of the geometry nodes follows. ArrangementOnSphere This node represents models as arrangements of geodesic arcs em-bedded on the sphere. ExactPolyhedron This node represents polyhedra as meshes using the Cgal Polyhedron 3 We plan to offer Sgal with an open-source license in the future, making it available to the public. Appendix A. Software Components, Libraries, and Packages data structure. SphericalGaussianMap This node represents polytopes as spherical Gaussian maps usingthe Arrangement on surface 2 data structure instantiated as an arrangement embed-ded on the sphere. CubicalGaussianMap This node represents polytopes as cubical Gaussian maps using the Cubical gaussian map 3 data structure.The implementation relies on inputing exact coordinates. To this end, the format wasfurther extended with a node called ExactCoordinate that represents exact coordinates,and enables the provision of exact rational coordinates as input.The entire partitioning process described in Chapter 5 is realized within Sgal . Thelibrary contains all the necessary ingredients required to represent and visualize the inputand the output, and to simulate the process. In particular, it has been extended with twogeometry node types: the Assembly node type represents assemblies or subassemblies, andthe AssemblyPart node type represents parts of assemblies. Notice, that each node ob-ject of the three types AssemblyPart , SphericalGaussianMap , and ArrangementOnSphere internally maintains the Cgal data structure that represents an arrangement of geodesicarcs embedded on the sphere [BFH + Arrangement on surface 2 classtemplate. A.2 Software Availability Compiling and executing the programs require the latest internal release of Cgal (postversion 3.31) and the components listed above including an internal package of Cgal thatsupports the Cubical gaussian map 3 data structure. The source code is available uponrequest. Precompiled executables compiled either with g++ VC 8 onWindows, data sets, and documentation can be downloaded from . Send email to [email protected] . ibliography [AFH02] Pankaj Kumar Agarwal, Eyal Flato, and Dan Halperin. Polygon decomposi-tion for efficient construction of Minkowski sums. Computational Geometry:Theory and Applications , 21:39–61, 2002.[AK00] Franz Aurenhammer and Rolf Klein. Voronoi diagrams. In J¨org-R¨udiger Sackand Jorge Urrutia, editors, Handbook of Computational Geometry , chapter 5,pages 201–290. Elsevier Science Publishers, B.V. North-Holland, Amsterdam,North-Holland, 2000.[Ale01] Andrei Alexandrescu. Modern C ++ Design . Addison-Wesley, 2001.[AS97] Boris Aronov and Micha Sharir. On translational motion planning of a convexpolyhedron in 3-space. SIAM Journal on Computing , 26(6):1785–1803, 1997.[AS00] Pankaj Kumar Agarwal and Micha Sharir. Arrangements and their applica-tions. In J¨org-R¨udiger Sack and Jorge Urrutia, editors, Handbook of Compu-tational Geometry , chapter 2, pages 49–119. Elsevier Science Publishers, B.V.North-Holland, Amsterdam, North-Holland, 2000.[AS01] Marcus Vin´ıcius Alvim Andrade and Jorge Stolfi. Exact algorithms for cir-cles on the sphere. International Journal of Computational Geometry andApplications , 11(3):267–290, jun 2001.[Aus99] Matthew H. Austern. Generic Programming and the STL . Addison-Wesley,1999.[BdLT97] Jean-Daniel Boissonnat, Eduard E. de Lange, and Monique Teillaud. Min-kowski operations for satellite antenna layout. In Proceedings of 13th AnnualACM Symposium on Computational Geometry (SoCG) , pages 67–76. Associ-ation for Computing Machinery (ACM) Press, 1997.[BDTY00] Jean-Daniel Boissonnat, Olivier Devillers, Monique Teillaud, and MarietteYvinec. Triangulations in C gal . In Proceedings of 16th Annual ACM Sym-posium on Computational Geometry (SoCG) , pages 11–18. Association forComputing Machinery (ACM) Press, 2000.[BEH + 02] Eric Berberich, Arno Eigenwillig, Michael Hemmer, Susan Hert, KurtMehlhorn, and Elmar Sch¨omer. A computational basis for conic arcs and10304 BIBLIOGRAPHY Boolean operations on conic polygons. In Proceedings of 10th Annual Euro-pean Symposium on Algorithms (ESA) , volume 2461 of LNCS , pages 174–186.Springer-Verlag, 2002.[BEH + 05] Eric Berberich, Arno Eigenwillig, Michael Hemmer, Susan Hert, Lutz Ket-tner, Kurt Mehlhorn, Joachim Reichel, Susanne Schmitt, Elmar Sch¨omer,and Nicola Wolpert. E xacus : Efficient and exact algorithms for curves andsurfaces. In Proceedings of 13th Annual European Symposium on Algorithms(ESA) , volume 3669 of LNCS , pages 155–166. Springer-Verlag, 2005.[BFH + 07] Eric Berberich, Efi Fogel, Dan Halperin, Kurt Melhorn, and Ron Wein. Sweep-ing and maintaining two-dimensional arrangements on surfaces: A first step.In Proceedings of 15th Annual European Symposium on Algorithms (ESA) ,volume 4698 of LNCS , pages 645–656. Springer-Verlag, 2007.[BFH + + Abstracts of23rd European Workshop on Computational Geometry , pages 223–226, 2007.[BGRR96] Julien Basch, Leonidas J. Guibas, G. D. Ramkumar, and L. Ramshaw. Poly-hedral tracings and their convolution. In Proceedings of 2nd Workshop onAlgorithmic Foundations of Robotics , 1996.[BHK + 05] Eric Berberich, Michael Hemmer, Lutz Kettner, Elmar Sch¨omer, and NicolaWolpert. An exact, complete and efficient implementation for computing pla-nar maps of quadric intersection curves. In Proceedings of 21st Annual ACMSymposium on Computational Geometry (SoCG) , pages 99–106. Associationfor Computing Machinery (ACM) Press, 2005.[BK08] Eric Berberich and Michael Kerber. Exact arrangements on tori and Dupincyclides. In Proceedings of ACM Symposium on Solid and Physical Modeling(SPM) , pages 59–66. Association for Computing Machinery (ACM) Press,2008.[BM07] Eric Berberich and Michal Meyerovitch. Computing envelopes of quadrics.In Abstracts of 23rd European Workshop on Computational Geometry , pages235–238, 2007.[BO79] Jon Louis Bentley and Thomas Ottmann. Algorithms for reporting and count-ing geometric intersections. IEEE Transactions on Computers , 28(9):643–647,1979. IBLIOGRAPHY M. C. Escher: His Life and Complete Graphic Work . HarryN. Abrams, Inc., 1982.[BR01] Henk Bekker and Jos B. T. M. Roerdink. An efficient algorithm to calculatethe Minkowski sum of convex 3D polyhedra. In Proceedings of InternationalConference on Computational Science Part I , volume 2073 of LNCS , pages619–628. Springer-Verlag, 2001.[Cam97] Stephen A. Cameron. Enhancing GJK: Computing minimum and penetrationdistances between convex polyhedra. In Proceedings of IEEE InternationalConference on Robotics and Automation , pages 3112–3117, 1997.[CC86] Stephen A. Cameron and R. K. Culley. Determining the minimum transla-tional distance between two convex polyhedra. In Proceedings of IEEE Inter-national Conference on Robotics and Automation , pages 591–596, 1986.[CDR92] John Canny, Bruce Donald, and Eugene K. Ressler. A rational rotationmethod for robust geometric algorithms. In Proceedings of 8th Annual ACMSymposium on Computational Geometry (SoCG) , pages 251–260. Associationfor Computing Machinery (ACM) Press, 1992.[cga07] Cgal User and Reference Manual , 3.3 edition, 2007. .[CKF + 04] Tim Culver, John Keyser, Mark Foskey, Shankar Krishnan, and DineshManocha. E solid — a system for exact boundary evaluation. Computer-Aided Design , 36(2):175–193, 2004.[CL06] Frederic Cazals and Sebastien Loriot. Computing the exact arrangement ofcircles on a sphere, with applications in structural biology. Technical Report6049, Inria Sophia-Antipolis, 2006.[Cof06] Stewart T. Coffin. Geometric Puzzle Design . A.K. Peters, Ltd., 2nd edition,2006.[COLYKT03] Daniel Cohen-Or, Shuly Lev-Yehudi, Adi Karol, and Ayellet Tal. Inner-coverof non-convex shapes. International Journal on Shape Modeling , 9(2):223–238,Dec 2003.[dBvKOS00] Mark de Berg, Mark van Kreveld, Mark Overmars, and Otfried Schwarzkopf. Computational Geometry: Algorithms and Applications . Springer-Verlag,Berlin, Germany, 2nd edition, 2000.[dCPT07] Pedro M. M. de Castro, Sylvain Pion, and Monique Teillaud. Exact andefficient computations on circles in C gal . In Abstracts of 23rd EuropeanWorkshop on Computational Geometry , pages 219–222, 2007.06 BIBLIOGRAPHY [DHH01] Duong Anh Duc, Nguyen Dong Ha, and Lethi Thuy Hang. Proposing amodel to store and a method to edit spatial data in topological maps. Tech-nical report, Ho Chi Minh University of Natural Sciences, Ho Chi Minh City,Vietnam, 2001.[DK90] David P. Dobkin and David G. Kirkpatrick. Determining the separation ofpreprocessed polyhedra — a unified approach. In Proceedings of 17th Inter-national Colloquium on Automata, Languages and Programming , volume 443of LNCS , pages 400–413. Springer-Verlag, 1990.[dW08] Michiel de Wilde. Using C gal for robust planar geometry processing inagilent ADS. 3rd Cgal User Workshop, 2008. .[EK08] Arno Eigenwillig and Michael Kerber. Exact and efficient 2D-arrangements ofarbitrary algebraic curves. In Proceedings of 19th Annual ACM-SIAM Sympo-sium on Discrete Algorithms (SODA) , pages 122–131, Philadelphia, PA, USA,2008. Society for Industrial and Applied Mathematics (SIAM).[EKP + 04] Ioannis Z. Emiris, Athanasios Kakargias, Sylvain Pion, Monique Teillaud, andElias P. Tsigaridas. Towards and open curved kernel. In Proceedings of 20thAnnual ACM Symposium on Computational Geometry (SoCG) , pages 438–446. Association for Computing Machinery (ACM) Press, 2004.[EKSW04] Arno Eigenwillig, Lutz Kettner, Elmar Sch¨omer, and Nicola Wolpert. Com-plete, exact and efficient computations with cubic curves. In Proceedings of20th Annual ACM Symposium on Computational Geometry (SoCG) , pages409–418. Association for Computing Machinery (ACM) Press, 2004.[EKW07] Arno Eigenwillig, Lutz Kettner, and Nicola Wolpert. Snap rounding of B´eziercurves. In Proceedings of 23rd Annual ACM Symposium on Computational Ge-ometry (SoCG) , pages 158–167. Association for Computing Machinery (ACM)Press, 2007.[EL00] Stephen A. Ehmann and Ming C. Lin. Accelerated proximity queries betweenconvex polyhedra by multi-level Voronoi marching. In Proceedings of IEEEConference on Intelligent Robots and Systems , pages 2101–2106, 2000.[EOR92] Roger C. Evans, Michael A. O’Connor, and Jarek R. Rossignac. Constructionof Minkowski sums and derivatives morphological combinations of arbitrarypolyhedra in CAD/CAM systems, 1992. US Patent 5159512.[ES86] Herbert Edelsbrunner and Raimund Seidel. Voronoi diagrams and arrange-ments. Discrete & Computational Geometry , 1:25–44, 1986.[FFHL02] Eyal Flato, Efi Fogel, Dan Halperin, and Eyal Leiserowitz. Movie: ExactMinkowski sums and applications. In Proceedings of 18th Annual ACM Sym-posium on Computational Geometry (SoCG) , pages 273–274. Association forComputing Machinery (ACM) Press, 2002. IBLIOGRAPHY + 00] Andreas Fabri, Geert-Jan Giezeman, Lutz Kettner, Stefan Schirra, and SvenSch¨onherr. On the design of C gal a computational geometry algorithmslibrary. Software — Practice and Experience , 30(11):1167–1202, 2000.[FH95] Ulrich Finke and Klaus H. Hinrichs. Overlaying simply connected planarsubdivisions in linear time. In Proceedings of 11th Annual ACM Symposium onComputational Geometry (SoCG) , pages 119–126. Association for ComputingMachinery (ACM) Press, 1995.[FH05] Efi Fogel and Dan Halperin. Movie: Exact Minkowski sums of convex polyhe-dra. In Proceedings of 21st Annual ACM Symposium on Computational Ge-ometry (SoCG) , pages 382–383. Association for Computing Machinery (ACM)Press, 2005.[FH06] Efi Fogel and Dan Halperin. Exact and efficient construction of Minkowskisums of convex polyhedra with applications. In Proceedings of 8th Workshopon Algorithm Engineering and Experiments , 2006.[FH07] Efi Fogel and Dan Halperin. Exact and efficient construction of Min-kowski sums of convex polyhedra with applications. Computer-Aided Design ,39(11):929–940, 2007.[FH08] Efi Fogel and Dan Halperin. Polyhedral assembly partitioning with infinitetranslations or the importance of being exact. 2008. Proceedings of 8thWorkshop on Algorithmic Foundations of Robotics.[FHH + 00] Eyal Flato, Dan Halperin, Iddo Hanniel, Oren Nechushtan, and Ester Ezra.The design and implementation of planar maps in C gal . The ACM Journalof Experimental Algorithmics , 5:1–23, 2000.[FHK + 07] Efi Fogel, Dan Halperin, Lutz Kettner, Monique Teillaud, Ron Wein, andNicola Wolpert. Arrangements. In J.-D. Boissonnat and M. Teillaud, editors, Effective Computational Geometry for Curves and Surfaces , chapter 1, pages1–66. Springer-Verlag, 2007.[FHW] Efi Fogel, Dan Halperin, and Christophe Weibel. On the exact maximumcomplexity of Minkowski sums of convex polyhedra. Discrete & ComputationalGeometry . Accepted for publication.[FHW + 04] Efi Fogel, Dan Halperin, Ron Wein, Sylvain Pion, Monique Teillaud, IoannisEmiris, Athanasios Kakargias, Elias Tsigaridas, Eric Berberich, Arno Eigen-willig, Michael Hemmer, Lutz Kettner, Kurt Mehlhorn, Elmar Schomer, andNicola Wolpert. An empirical comparison of software for constructing ar-rangements of curved arcs (preliminary version). Technical Report ECG-TR-361200-01, Tel-Aviv University, INRIA Sophia-Antipolis, MPI Saarbr¨ucken,2004.08 BIBLIOGRAPHY [FHW07] Efi Fogel, Dan Halperin, and Christophe Weibel. On the exact maximumcomplexity of Minkowski sums of convex polyhedra. In Proceedings of 23rdAnnual ACM Symposium on Computational Geometry (SoCG) , pages 319–326. Association for Computing Machinery (ACM) Press, 2007.[FSH08a] Efi Fogel, Ophir Setter, and Dan Halperin. Exact implementation of arrange-ments of geodesic arcs on the sphere with applications. In Abstracts of 24thEuropean Workshop on Computational Geometry , pages 83–86, 2008.[FSH08b] Efi Fogel, Ophir Setter, and Dan Halperin. Movie: Arrangements of geodesicarcs on the sphere. In Proceedings of 24th Annual ACM Symposium on Com-putational Geometry (SoCG) , pages 218–219. Association for Computing Ma-chinery (ACM) Press, 2008.[FT07] Efi Fogel and Monique Teillaud. Generic programming and the C gal library.In J.-D. Boissonnat and M. Teillaud, editors, Effective Computational Ge-ometry for Curves and Surfaces , chapter 8, pages 313–320. Springer-Verlag,2007.[Fuk04] Komei Fukuda. From the zonotope construction to the Minkowski addition ofconvex polytopes. Journal of Symbolic Computation , 38(4):1261–1272, 2004.[FvDFH95] James D. Foley, Andries van Dam, Steven K. Feiner, and John F. Hughes. Computer Graphics: Principles and Practice in C . Addison-Wesley, 2nd edi-tion, 1995.[FW07] Komei Fukuda and Christophe Weibel. f-vectors of Minkowski additions ofconvex polytopes. Discrete & Computational Geometry , 37:503–516, 2007.[FWH04] Efi Fogel, Ron Wein, and Dan Halperin. Code flexibility and program effi-ciency by genericity: Improving C gal ’s arrangements. In Proceedings of 12thAnnual European Symposium on Algorithms (ESA) , volume 3221 of LNCS ,pages 664–676. Springer-Verlag, 2004.[FWZH07] Efi Fogel, Ron Wein, Baruch Zukerman, and Dan Halperin. 2d regularizedboolean set-operations. In Cgal Editorial Board, editor, Cgal User andReference Manual . 3.3 edition, 2007.[GGHT97] Michael T. Goodrich, Leonidas J. Guibas, John Hershberger, and Paul J.Tanenbaum. Snap rounding line segments efficiently in two and three dimen-sions. In Proceedings of 13th Annual ACM Symposium on Computational Ge-ometry (SoCG) , pages 284–293. Association for Computing Machinery (ACM)Press, 1997.[GHH + 98] Leonidas J. Guibas, Dan Halperin, Hirohisa Hirukawa, Jean-Claude Latombe,and Randall H. Wilson. Polyhedral assembly partitioning using maximallycovered cells in arrangements of convex polytopes. International Journal ofComputational Geometry and Applications , 8(2):179–200, 1998. IBLIOGRAPHY DesignPatterns — Elements of Reusable Object-Oriented Software . Addison-Wesley,1995.[Gho93] Pijush K. Ghosh. A unified computational framework for Minkowski opera-tions. Computers & Graphics , 17(4):357–378, 1993.[GHZ99] Leonidas J. Guibas, David Hsu, and Li Zhang. H-walk: Hierarchical distancecomputation for moving convex bodies. In Proceedings of 15th Annual ACMSymposium on Computational Geometry (SoCG) , pages 265–273. Associationfor Computing Machinery (ACM) Press, 1999.[GJK88] Elmer G. Gilbert, Daniel W. Johnson, and Sathiya S. Keerthi. A fast pro-cedure for computing the distance between complex objects. 4(2):193–203,1988.[GM95] Leonidas J. Guibas and David H. Marimont. Rounding arrangements dynam-ically. In Proceedings of 11th Annual ACM Symposium on Computational Ge-ometry (SoCG) , pages 190–199. Association for Computing Machinery (ACM)Press, 1995.[GRS83] Leonidas J. Guibas, Leo Ramshaw, and Jorge Stolfi. A kinetic framework forcomputational geometry. In Proceedings of 24th Annual IEEE Symposium onthe Foundations of Computer Science , pages 100–111, 1983.[GS87] Leonidas J. Guibas and Raimund Seidel. Computing convolutions by recipro-cal search. Discrete & Computational Geometry , 2:175–193, 1987.[GS93] Peter Gritzmann and Bernd Sturmfels. Minkowski addition of polytopes:Computational complexity and applications to Gr¨obner bases. SIAM Journalon Discrete Math , 6(2):246–269, 1993.[Hac07] Peter Hachenberger. Exact Minkowski sums of polyhedra and exact and ef-ficient decomposition of polyhedra into convex pieces. In Proceedings of 15thAnnual European Symposium on Algorithms (ESA) , volume 4698 of LNCS ,pages 669–680. Springer-Verlag, 2007.[Hal04] Dan Halperin. Arrangements. In Jacob E. Goodman and Joseph O’Rourke,editors, Handbook of Discrete and Computational Geometry , chapter 24, pages529–562. Chapman & Hall/CRC, 2nd edition, 2004.[Han00] Iddo Hanniel. The design and implementation of planar arrangements ofcurves in C gal . M.Sc. thesis, School of Computer Science, Tel Aviv Univer-sity, 2000.[HH00] Iddo Hanniel and Dan Halperin. Two-dimensional arrangements in C gal andadaptive point location for parametric curves. In Proceedings of InternationalWorkshop on Algorithm Engineering (WAE) , volume 1982 of LNCS , pages171–182. Springer-Verlag, 2000.10 BIBLIOGRAPHY [HH03] Shai Hirsch and Dan Halperin. Hybrid motion planning: Coordinating twodiscs moving among polygonal obstacles in the plane. In Jean-Daniel Boisson-nat, Joel Burdick, Ken Goldberg, and Seth Hutchinson, editors, AlgorithmicFoundations of Robotics V , pages 239–255. Springer-Verlag, 2003.[HH08] Idit Haran and Dan Halperin. An experimental study of point location in pla-nar arrangements in C gal . The ACM Journal of Experimental Algorithmics ,13, 2008.[HHK + 07] Susan Hert, Michael Hoffmann, Lutz Kettner, Sylvain Pion, and Michael Seel.An adaptable and extensible geometry kernel. Computational Geometry: The-ory and Applications , 38(1-2):16–36, 2007.[HK07] Peter Hachenberger and Lutz Kettner. 3D Boolean operations on nef polyhe-dra. In Cgal Editorial Board, editor, Cgal User and Reference Manual . 3.3edition, 2007.[HKL + 99] Kenneth E. Hoff III, John Keyser, Ming Lin, Dinesh Manocha, and Tim Cul-ver. Fast computation of generalized voronoi diagrams using graphics hard-ware. In Proceedings of 26th Annual International Conference on ComputerGraphics and Interactive Techniques , pages 277–286. Association for Comput-ing Machinery (ACM) Press, 1999.[HKL04] Dan Halperin, Lydia E. Kavraki, and Jean-Claude Latombe. Robotics. InJacob E. Goodman and Joseph O’Rourke, editors, Handbook of Discreteand Computational Geometry , chapter 48, pages 1065–1093. Chapman &Hall/CRC, 2nd edition, 2004.[HKM07] Peter Hachenberger, Lutz Kettner, and Kurt Mehlhorn. Boolean operations on3D selective Nef complexes: Data structure, algorithms, optimized implemen-tation and experiments. Computational Geometry: Theory and Applications ,38(1-2):64–99, 2007. Special issue on C gal .[HL04] Dan Halperin and Eran Leiserowitz. Controlled perturbation for arrangementsof circles. volume 14, pages 277–310, 2004.[HLW00] Dan Halperin, Jean-Claude Latombe, and Randall H. Wilson. A generalframework for assembly planning: The motion space approach. Algorithmica ,26:577–601, 2000.[Hob99] John D. Hobby. Practical segment intersection with finite precision output. Computational Geometry: Theory and Applications , 13(4):199–214, 1999.[Hof04] Christoph M. Hoffmann. Solid modeling. In Jacob E. Goodman and JosephO’Rourke, editors, Handbook of Discrete and Computational Geometry , chap-ter 56, pages 1257–1278. Chapman & Hall/CRC, 2nd edition, 2004. IBLIOGRAPHY Abstracts of 17thEuropean Workshop on Computational Geometry , pages 82–85. Freie Univer-sit¨at Berlin, 2001.[HP02] Dan Halperin and Eli Packer. Iterated snap rounding. Computational Geom-etry: Theory and Applications , 23:209–225, 2002.[HRS92] Craig D. Hodgson, Igor Rivin, and Warren D. Smith. A characterization ofconvex hyperbolic polyhedra and of convex polyhedra inscribed in the sphere. Bull. (New Series) of the AMS , 27:246–251, 1992.[HS98] Dan Halperin and Christian R. Shelton. A perturbation scheme for sphericalarrangements with application to molecular modeling. Computational Geom-etry: Theory and Applications , 10:273–287, 1998.[HW07] Iddo Hanniel and Ron Wein. An exact, complete and efficient computation ofarrangements of B´ezier curves. In Proceedings of ACM Symposium on Solidand Physical Modeling (SPM) , pages 253–263. Association for Computing Ma-chinery (ACM) Press, 2007.[IIM85] Hiroshi Imai, Masao Iri, and Kazuo Murota. Voronoi diagram in the Laguerregeometry and its applications. SIAM Journal on Computing , 14(1):93–105,1985.[KCMK00] John Keyser, Tim Culver, Dinesh Manocha, and Shankar Krishnan. Effi-cient and exact manipulation of algebraic points and curves. Computer-AidedDesign , 32(11):649–662, 2000.[Ket99] Lutz Kettner. Using generic programming for designing a data structurefor polyhedral surfaces. Computational Geometry: Theory and Applications ,13(1):65–90, 1999.[Ket07a] Lutz Kettner. 3D polyhedral surfaces. In Cgal Editorial Board, editor, Cgal User and Reference Manual . 3.3 edition, 2007.[Ket07b] Lutz Kettner. Halfedge data structures. In Cgal Editorial Board, editor, Cgal User and Reference Manual . 3.3 edition, 2007.[KK95] Lydia E. Kavraki and Mihail N. Kolountzakis. Partitioning a planar assemblyinto two connected parts is NP-complete. Information Processing Letters ,55:159–165, 1995.[KLPY99] Vijay Karamcheti, Chen Li, Igor Pechtchanski, and Chee K. Yap. A corelibrary for robust numeric and geometric computation. In Proceedings of 15thAnnual ACM Symposium on Computational Geometry (SoCG) , pages 351–359. Association for Computing Machinery (ACM) Press, 1999.12 BIBLIOGRAPHY [kLsKE98] In kwon Lee, Myung soo Kim, and Gershon Elber. Polynomial/rational ap-proximation of Minkowski sum boundary curves. Graphical Models and ImageProcessing , 60(2):136–165, 1998.[KMP + 08] Lutz Kettner, Kurt Mehlhorn, Sylvain Pion, Stefan Schirra, and Chee Yap.Classroom examples of robustness problems in geometric computations. Com-putational Geometry: Theory and Applications , 40(1):61–78, 2008.[KMW98] Sanjeev Khanna, Rajeev Motwani, and Randall H. Wilson. On certificates andlookahead in dynamic graph problems. Algorithmica , 21(4):377–394, 1998.[KN04] Lutz Kettner and Stefan N¨aher. Two computational geometry libraries: L eda and C gal . In Jacob E. Goodman and Joseph O’Rourke, editors, Handbook ofDiscrete and Computational Geometry , chapter 65, pages 1435–1463. Chap-man & Hall/CRC, Boca Raton, FL, 2nd edition, 2004.[KR91] Anil Kaul and Jarek R. Rossignac. Solid-interpolating deformations: Con-struction and animation of PIPs. In Proceedings of European Computer Graph-ics Conference (Eurographics) , pages 493–505, 1991.[Lat91] Jean-Claude Latombe. Robot Motion Planning . Kluwer Academic Publishers,Boston, 1991.[LC91] Ming C. Lin and John F. Canny. A fast algorithm for incremental distancecalculation. In Proceedings of IEEE International Conference on Robotics andAutomation , pages 1008–1014, 1991.[LM04] Ming C. Lin and Dinesh Manocha. Collision and proximity queries. In Ja-cob E. Goodman and Joseph O’Rourke, editors, Handbook of Discrete andComputational Geometry , chapter 35, pages 787–807. Chapman & Hall/CRC,2nd edition, 2004.[LPT08] Sylvain Lazard, Luis Pe˜naranda, and Elias Tsigaridas. A C gal based alge-braic kernel based on RS and application to arrangements. In Abstracts of24th European Workshop on Computational Geometry , pages 91–94, 2008.[Luk57] Dorman Luke. Stellations of the rhombic dodecahedron. The MathematicalGazette , 41(337):189–194, 1957.[Mey06] Michal Meyerovitch. Robust, generic and efficient construction of envelopesof surfaces in three-dimensional space. In Proceedings of 14th Annual Euro-pean Symposium on Algorithms (ESA) , volume 4168 of LNCS , pages 792–803.Springer-Verlag, 2006.[Mir98] Brian Mirtich. V-clip: Fast and robust polyhedral collision detection. ACMTransactions on Graphics , 17(3):177–208, 1998.[MN00] Kurt Mehlhorn and Stefan N¨aher. Leda : A Platform for Combinatorial andGeometric Computing . Cambridge University Press, Cambridge, UK, 2000. IBLIOGRAPHY Automata, Languages and Programming ,volume 4051 of LNCS , pages 299–310. Springer-Verlag, 2006.[MS88] David A. Musser and Alexander A. Stepanov. Generic programming. In Pro-ceedings of International Conference on Symbolic and Algebraic Computation ,volume 358 of LNCS , pages 13–25. Springer-Verlag, 1988.[MS03] Kurt Mehlhorn and Michael Seel. Infimaximal frames: A technique for makinglines look like segments. International Journal of Computational Geometryand Applications , 13(3):241–255, 2003.[Mul90] Ketan Mulmuley. A fast planar partition algorithm, I. Journal of SymbolicComputation , 10(3-4):253–280, 1990.[MWZ07] Michal Meyerovitch, Ron Wein, and Baruch Zukerman. 3D envelopes. In Cgal Editorial Board, editor, Cgal User and Reference Manual . 3.3 edition,2007.[Mye98] Nathan Myers. A new and useful template technique: “Traits”. In Stanly B.Lippman, editor, C ++ Gems , volume 5 of SIGS Reference Library , pages 451–458. Cambridge University Press, Cambridge, UK, 1998.[Nat88] B. K. Natarajan. On planning assemblies. In Proceedings of 4th Annual ACMSymposium on Computational Geometry (SoCG) , pages 299–308. Associationfor Computing Machinery (ACM) Press, 1988.[NLC02] Hyeon-Suk Na, Chung-Nim Lee, and Otfried Cheong. Voronoi diagrams on thesphere. Computational Geometry: Theory and Applications , 23(2):183–194,2002.[OBSC00] Atsuyuki Okabe, Barry Boots, Kokichi Sugihara, and Sung Nok Chiu. SpatialTessellations: Concepts and Applications of Voronoi Diagrams . John Wiley& Sons, NYC, 2nd edition, 2000.[Ove96] Mark H. Overmars. Designing the computational geometry algorithms libraryC gal . In Proceedings of ACM Workshop on Applied Computational Geome-try, Towards Geometric Engineering , volume 1148, pages 53–58, London, UK,1996. Springer-Verlag.[PF06] Sylvain Pion and Andreas Fabri. A generic lazy evaluation scheme for exactgeometric computations. In ,2006.[PS85] Franco P. Preparata and Michael Ian Shamos. Computational Geometry: AnIntroduction . Springer-Verlag, New York, NY, 1985.[PY07] Sylvain Pion and Mariette Yvinec. 2D triangulation data structure. In Cgal Editorial Board, editor, Cgal User and Reference Manual . 3.3 edition, 2007.14 BIBLIOGRAPHY [Rog03] Vadim Rogol. Maximizing the area of an axially-symmetric polygon inscribedby a simple polygon. Master’s thesis, Technion, Haifa, Israel, 2003.[Sch00] Stefan Schirra. Robustness and precision issues in geometric computation.In J¨org-R¨udiger Sack and Jorge Urrutia, editors, Handbook of ComputationalGeometry , chapter 14, pages 597–632. Elsevier Science Publishers, B.V. North-Holland, Amsterdam, North-Holland, 2000.[See07] Michael Seel. 2D Boolean operations on nef polygons. In Cgal Edito-rial Board, editor, Cgal User and Reference Manual . 2007.[SH75] Michael Ian Shamos and Dan Hoey. Closest-point problems. In Proceedingsof 16th IEEE Symposium on the Foundations of Computer Science , pages151–162, 1975.[SH89] Jack Snoeyink and John Hershberger. Sweeping arrangements of curves.In Proceedings of 5th Annual ACM Symposium on Computational Geometry(SoCG) , pages 354–363. Association for Computing Machinery (ACM) Press,1989.[Sha04] Micha Sharir. Algorithmic motion planning. In Jacob E. Goodman and JosephO’Rourke, editors, Handbook of Discrete and Computational Geometry , chap-ter 47, pages 1037–1064. Chapman & Hall/CRC, 2nd edition, 2004.[SKS02] Joon-Kyung Seong, Myung-Soo Kim, and Kokichi Sugihara. The Minkowskisum of two simple surfaces generated by slope-monotone closed curves. In Geometric Modeling and Processing — Theory and Applications , pages 33–42, 2002.[SLL02] Jeremy G. Siek, Lie-Quan Lee, and Andrew Lumsdaine. The B oost GraphLibrary . Addison-Wesley, 2002.[SS94] Jack Snoeyink and Jorge Stolfi. Objects that cannot be taken apart with twohands. Discrete & Computational Geometry , 12:367–384, 1994.[SSH08] Ophir Setter, Micha Sharir, and Dan Halperin. Construction two-dimensionalVoronoi diagrams via divide and conquer of envelopes in space, 2008.Manuscript.[Sug02] Kokichi Sugihara. Laguerre Voronoi diagram on the sphere. Journal forGeometry and Graphics , 6(1):69–81, 2002.[SWK03] Kevin Sahr, Denis White, and A. Jon Kimerling. Geodesic discrete globalgrid systems. Cartography and Geographic Information Science , 30(2):121–134, 2003.[Tar72] Robert E. Tarjan. Depth first search and linear graph algorithms. SIAMJournal on Computing , 1(2):146–160, 1972. IBLIOGRAPHY Movable separability of sets . North-Holland, 1985.[VKSM05] Gokul Varadhan, Shankar Krishnan, T. V. N. Sriram, and Dinesh Manocha.A simple algorithm for complete motion planning of translating polyhedralrobots. International Journal of Robotics Research , 24(11):983–995, 2005.[VM06] Gokul Varadhan and Dinesh Manocha. Accurate Minkowski sum approxima-tion of polyhedral models. Graphical Models and Image Processing , 68(4):343–355, 2006.[Wei02] Ron Wein. High-level filtering for arrangements of conic arcs. In Proceedingsof 10th Annual European Symposium on Algorithms (ESA) , volume 2461 of LNCS , pages 884–895. Springer-Verlag, 2002.[Wei05] Ron Wein. Efficient implementation of red-black trees with split and catenateoperations. Technical report, Tel-Aviv University, 2005. .[Wei07] Ron Wein. The Integration of Exact Arrangements with Effective MotionPlanning . Ph.D. thesis, The Blavatnik School of Computer Science, Tel AvivUniversity, 2007.[WF05] Ron Wein and Efi Fogel. The new design of C gal ’s arrangement package.Technical report, Tel-Aviv University, 2005. .[WFZH05] Ron Wein, Efi Fogel, Baruch Zukerman, and Dan Halperin. Advanced pro-gramming techniques applied to C gal ’s arrangement package. In , 2005.[WFZH07a] Ron Wein, Efi Fogel, Baruch Zukerman, and Dan Halperin. 2D arrangements.In Cgal Editorial Board, editor, Cgal User and Reference Manual . 3.3 edi-tion, 2007.[WFZH07b] Ron Wein, Efi Fogel, Baruch Zukerman, and Dan Halperin. Advanced pro-gramming techniques applied to C gal ’s arrangement package. ComputationalGeometry: Theory and Applications , 38(1–2):37–63, 2007. Special issue onC gal .[WL94] Randall H. Wilson and Jean-Claude Latombe. Geometric reasoning aboutmechanical assembly. Artificial Intelligence , 71(2):371–396, 1994.[WSD03] Yanyan Wu, Jami J. Shah, and Joseph K. Davidson. Improvements to al-gorithms for computing the Minkowski sum of 3-polytopes. Computer-AidedDesign , 35(13):1181–1192, 2003.[WvH07] Ron Wein, Jur P. van den Berg, and Dan Halperin. The visibility-Voronoicomplex and its applications. Computational Geometry: Theory and Applica-tions , 36(1):66–87, 2007.16 BIBLIOGRAPHY [Yap04] Chee K. Yap. Robust geomtric computation. In Jacob E. Goodman andJoseph O’Rourke, editors, Handbook of Discrete and Computational Geome-try , chapter 41, pages 927–952. Chapman & Hall/CRC, 2nd edition, 2004. inks [1] Using Cgal for robust planar geometry processing in Agilent ADS. http://acg.cs.tau.ac.il/projects/external-projects/agilent-ads/project-page .[2] Algorithmic automation. http://goldberg.berkeley.edu/algorithmic-automation .[3] Boost — portable C ++ libraries. .[4] Boost , Generic Programming Techniques. .[5] Cgal — computational geometry algorithms library. .[6] The Core number library homepage. http://cs.nyu.edu/exact/core_pages .[7] Ecg — effective computational geometry for curves and surfaces. .[8] Esolid — exact boundary evaluation of low-degree curved solids. .[9] Exacus — efficient and exact algorithms for curves and surfaces. .[10] Function object. http://en.wikipedia.org/wiki/Function_object .[11] GeometryFactory . .[12] Gmp — GNU multiple precision arithmetic library. http://gmplib.org .[13] Gnuplot homepage. . 11718 LINKS [14] Google maps. http://maps.google.com .[15] Imagemagick homepage. .[16] Cgal arrangement of Irit free-form curves. .[17] Leda — library for efficient data types and algorithms. .[18] Leda external package: Sphere geometry. .[19] Mapc — efficient and exact manipulation of algebraic points and curves. .[20] Minkowski sum related movies. http://acg.cs.tau.ac.il/movies .[21] The Q uick CD library homepage. .[22] Real-time high dynamic range image-based lighting. .[23] Wolfram Mathworld simple polygon. Simple polygon. http://mathworld.wolfram.com/SimplePolygon.html .[24] The S olid library homepage. .[25] Stl — C ++ standard template library. .[26] S wift++ library homepage. http://gamma.cs.unc.edu/SWIFT++/ .[27] The web3D homepage. .[28] Christophe Weibel. Minkowski sums. http://roso.epfl.ch/cw/poly/public.phphttp://roso.epfl.ch/cw/poly/public.php