6. CGAL, the Computational Geometry Algorithms Library

• Introduction
• Related Work
• Overview of the Library Structure
• Geometric Kernel
• Basic Library
• Geometric Traits Classes

Introduction

Geometric algorithms arise in various areas of computer science. Computer graphics and virtual reality, computer aided design and manufacturing, solid modeling, robotics, geographical information systems, computer vision, shape reconstruction, molecular modeling, and circuit design are well-known examples.

To a large extent the theory has been developed with asymptotic worst-case complexity analysis and under the assumption of the real RAM model. Computations with real numbers are assumed to be in constant time. For many (perhaps most) algorithms and problems this is a justified assumption that can be perfectly simulated with finite precision numbers if the input has limited precision. In 1996 the Computational Geometry Impact Task Force published a task force report, Application Challenges to Computational Geometry. Although crediting the remarkable success of the field in theory, the report now demanded to address also the applicability in practice. Recommendation number one on the list of four was ``production and distribution of usable (and useful) geometric codes''. Where are the difficulties in doing so? There are four major reasons why implementing geometric algorithms in particular is seen to be more difficult than in other fields:

• Algorithms in geometry are among the most advanced algorithms in algorithm design and they make frequent use of complicated data structures.

  Furthermore, software engineering has ignored the difficulties in algorithm engineering for quite a while. For example, the main focus in object-oriented design is on data abstraction, data encapsulation, relationships among data, its reuse and the design of large-scale systems. Only recently have implementations of algorithms been rediscovered as an active topic in software engineering, for example with the generic programming paradigm [Musser89].

• The asymptotic worst-case complexity analysis does not match the practical need for two reasons. In practice, the constant hidden in the asymptotic analysis can easily outweigh asymptotic factors, such as log-factors, and the worst-case usually depends on a general input model that may be unrealistic. More realistic input models, such as fatness, can help understanding the practicability of algorithms.

• Theoretical papers assume exact arithmetic with real numbers. The correctness proofs of the algorithms rely on exact computation, and replacing exact arithmetic by imprecise built-in floating-point arithmetic does not work in general. Geometric algorithms in particular are sensitive to rounding errors since numerical data and control flow decisions have usually a strong interrelation. The numerical problems may destroy the geometric consistency that an algorithm may rely on. As a result the program may crash, may run into an infinite loop, or -- perhaps worst of all -- may produce unpredictable erroneous output.

  The requirements on the arithmetic vary with the algorithm. Some algorithms require only sign computations of polynomial expressions of bounded degree in the input variables. Others require unbounded degree or algebraic roots. Various packages for exact arithmetic are available for different needs. Another approach is to redesign the algorithm to cope with inexact arithmetic. Usually the output is only an approximation of the exact solution. As a common prerequisite for exact arithmetic the input is rounded, either to convert the input into the format required for the arithmetic (rounding floating point to integer), or to lessen the precision requirements on the arithmetic.

• Often, theoretical papers exclude degenerate configurations in the input. Typically, these degeneracies are specific to the algorithm and the problem, and would involve the treatment of special cases in the algorithm. Simple examples of configurations considered as degenerate are duplicate points in a point set or three lines intersecting in one point. For some problems, it is not difficult to handle the degeneracies, but for other problems the special case treatment distracts from the solution of the general problem and it can amount to a considerable fraction of the coding effort.

  In theory, this approach of excluding degeneracies from consideration is justified by the argument that degenerate cases are very rare in the set of all possible input over the real numbers, i.e., they have zero probability if the input set is randomly chosen over the real numbers. Another argument is that it is first of all important to understand the general case before treating special cases.

  In practice, however, degenerate input occurs frequently. For instance, the coordinates of the geometric objects may not be randomly chosen over the real numbers, but lie on a grid. They may be created by clicking in a window in a graphical user interface. In some applications, what are called degeneracies are even high-valued design criteria. In architecture features of buildings do align on purpose. As a consequence, practical implementations usually must address the handling of degeneracies.

  General approaches in handling degeneracies are symbolic perturbation or randomized perturbation with performance and correctness guarantee.

To remedy this situation the Computational Geometry Algorithms Library, CGAL (http://www.cgal.org/), has been started five years ago in Europe in order to provide correct, efficient, and reusable implementations [Fabri99]. The library is being developed by several universities and research institutes in Europe and Israel.

The major design goals for CGAL include correctness, robustness, flexibility, efficiency, and ease-of-use. One aspect of flexibility is that CGAL algorithms can be easily adapted to work on data types in applications that already exist. The design goals, especially flexibility and efficient robust computation, led us to opt for the generic programming paradigm using templates in C++, and to reject the object-oriented paradigm in C++ (as well as in Java). In several appropriate places, however, we make use of object-oriented solutions and design patterns. Generic programming with templates in C++ also provides us with the help of strong type checking at compile time. Moreover, the C++ abstractions used by us do not cause any runtime overhead.

The birth of the CGAL-library dates back to a meeting in Utrecht in January 1995. Shortly afterwards, the five authors of [Fabri99] started developing the kernel. The CGAL-project has been funded officially since October 1996 and the team of developers has grown considerably among professionals in the field of computational geometry and related areas, in academia especially research assistants, PhD students and postdocs. The CGAL release 1.2 of January 1999 consists of approximately 110,000 lines of C++ source code for the library, plus 50,000 lines for accompanying sources, such as the test suite and example programs, not counting C++ comments or empty lines (the release 2.4 of May 2002 consists of approximately 290,000 lines of code, comments and empty lines included this time, plus test suite and example programs). In terms of the elder Constructive Cost Model (COCOMO) the line counts, the people involved, and the time schedule indicate a large project comparable to operating systems or database management systems. The WWW home-page of CGAL (http://www.cgal.org/) provides a list of publications about CGAL and related research.

Related Work

Collecting isolated implementations, also called the Gems approach according to the successful Graphics Gems series, usually requires some adaption effort to make things work together. Compared to the graphics gems, computational geometry implementations usually use more involved data structures and more advanced algorithms. This makes adaption harder. A good collection provides the Directory of Computational Geometry Software (http://www.geom.umn.edu/software/cglist/).

Integrated applications and workbenches provide a homogeneous environment, for example with animation and interaction capabilities, and all parts work smoothly together. However, they tend to be monolithic, hard to extend, and hard to reuse in other projects. Examples date back to the end of the Eighties, specifically XYZ GeoBench (http://wwwjn.inf.ethz.ch/geobench/XYZGeoBench.html) developed at ETH Zurich, Switzerland, is one of the precursors of CGAL.

Software libraries promise that the components work seamlessly together, that the library is extensible and that the components can be reused in other projects. Examples are the precursors of CGAL developed by members of the CGAL consortium. These precursors are PlaGeo developed at Utrecht University, C++Gal developed at Inria Sophia-Antipolis, and the geometric part of LEDA (http://www.mpi-sb.mpg.de/LEDA/), a library for combinatorial and geometric computing, which has been developed at Max-Planck-Institut für Informatik, Saarbrücken. Another example is GeomLib, a computational geometry library implemented in Java at the Center for Geometric Computing, located at Brown University, Duke University, and Johns Hopkins University in the United States. They state their goal as an effective technology transfer from Computational Geometry to relevant applied fields.

Overview of the Library Structure

The library layers and the support library are further subdivided into smaller modular units. The modular approach has several benefits: The library is easier to learn, the implementation work is more easily spread among the project partners, and the reduction of dependencies facilitates testing and maintenance.

The geometric kernel contains simple geometric objects of constant size such as points, lines, segments, triangles, tetrahedra, circle and more. It provides geometric predicates on those objects, operations to compute intersections of and distances between objects, and affine transformations. The kernel objects are closed under affine transformations, e.g., the existence of circles implies that there are also ellipses in the kernel.

The geometric kernel is split in three parts, one for two-dimensional objects, one for three-dimensional objects, and one for general-dimensional objects. Geometry in two and three dimensions is well studied and has lots of applications which explains their special status. For all dimensions there are Cartesian and homogeneous representations available for the coordinates.

To solve robustness problems, CGAL advocates the use of exact arithmetic instead of floating point arithmetic. An arithmetic is associated with a number type in CGAL and the classes in the geometric kernel are parameterized by number types. CGAL provides own number types and supports number types from other sources, e.g., from LEDA or the Gnu Multiple Precision library. Since the arithmetic operations needed in CGAL are quite basic, every library supplying number types can be adapted easily to work with CGAL.

The basic library contains more complex geometric objects and data structures: polygons, triangulations, planar maps, polyhedra and so on. It also contains algorithms, such as for computing the convex hull of a set of points, the union of two polygons, smallest enclosing ellipse and so on. The figure above indicates the major parts in the basic library. These parts are mostly independent from each other and even independent from the kernel. This independence has been achieved with geometric traits classes, to be discussed later. Default implementations of the traits classes use the CGAL kernel for the types and primitive operations. Other implementations of the traits classes provided in CGAL use the LEDA geometric part. The traits class requirements are simple enough for a user to be able to write a traits class for own geometric data types and operations.

The core library offers basic non-geometric functionality that is needed in the geometric kernel or the basic library, for example support for coping with different C++ compilers which all have their own limitations. The core library contains the support for assertions, preconditions and postconditions. Circulators and random number generators belong here as well.

The support library also contains functionality with non-geometric aspects. In contrast to the core library, this functionality is not needed by the geometric kernel nor the basic library. The support library interfaces the geometric objects with external representations, like visualizations or external file formats. Among the list of supported formats are VRML and PostScript as well as the GeomView program and LEDA windows for 2D and 3D visualization. The support library also contains generators for synthetic test data sets, for example random points uniformly distributed in a certain domain. The adaptation of number types from other libraries is contained in the support library as well. The separation from the kernel and the basic library makes the functionality in the support library orthogonal and therefore open for future extensions.

Geometric Kernel

The current geometric kernel provides two families of geometric objects: One family is based on the representation of points using Cartesian coordinates. The other family is based on the representation of points using homogeneous coordinates. The homogeneous representation extends the representation with Cartesian coordinates by an additional coordinate which is used as a common denominator. More formally, in d-dimensional space, a point with homogeneous coordinates $(x$₀, x₁, ..., x_d-1, x_d), where $x$_d != 0, has Cartesian coordinates $(x$₀/x_d, x₁/x_d, ..., x_d-1/x_d). This avoids divisions and reduces many computations in geometric algorithms to calculations over the integers. The homogeneous representation is used for affine geometry in CGAL, and not projective geometry, where the homogeneous representation is usually known from. Both families are parameterized by the number type used to represent the Cartesian or homogeneous coordinates. The type CGAL::Cartesian<double> specifies the Cartesian representation with coordinates of type double, and the type CGAL::Homogeneous<int> specifies the homogeneous representation with coordinates of type int. These representation types are used as template argument in all geometric kernel types, like a two-dimensional point declared as

template <class R>
CGAL::Point_2;

typedef CGAL::Point_2< CGAL::Homogeneous<int> > Point_2;

CGAL provides clean mathematical concepts to the user without sacrificing efficiency. For example, CGAL strictly distinguishes points and (mathematical) vectors, i.e., it distinguishes affine geometry from the underlying linear algebra. Points and vectors are not the same with regard to illicit computations resulting from identification of points and vectors in geometric computations. In particular, points and vectors behave differently under affine transformations. We do not even provide automatic conversion between points and vectors but use the geometric concept of an origin instead. The symbolic constant CGAL::ORIGIN represents a point and can be used to compute the locus vector as the difference between a point and the origin. Function overloading is used to implement this operation internally as a simple conversion without any overhead. Note that we do not provide the geometrically invalid addition of two points, since this might lead to ambiguous expressions: Assuming three points p, q, and r and an affine transformation A, one can write in CGAL the perfectly legal expression $A(\; p\; +\; (\; q\; -\; r))$. The slightly different expression $A((\; p\; +\; q)\; -\; r)$ contains the illegal addition of two points. However, if we allow this addition, we would expect the same result coordinatewise as in the previous, legal expression. But this is not necessarily intended, since the expression within the affine transformation is meant to evaluate to a vector, and not to a point as in the previous expression. Vectors and points behave differently under affine transformations. To avoid these ambiguities, the automatic conversion between points and vectors is not provided as well. (see Origin.C for an example point/vector/origin implementation)

Class hierarchies are used rarely in CGAL. One example are affine transformations, which maintain distinct internal representations specialized on restricted transformations. The internal representations differ considerably in their space requirements and the efficiency of their member functions. For all but the most general representation we gain performance in terms of space and time. And for the most general representation, the performance penalty caused by the virtual functions is negligible, because the member functions are computationally expensive for this general representation. Alternatively we could have used this general representation for affine transformations only. But the use of a hierarchy is justified, since the specialized representations, namely translation, rotation and scaling, arise frequently in geometric computing.

Another design decision was to make the (constant-size) geometric objects in the kernel non-modifiable (value semantics). For example, there are no member functions to set the Cartesian coordinates of a point. Points are viewed as atomic units, and no assumption is made on how these objects are represented. In particular, there is no assumption that points are represented with Cartesian coordinates. They might use polar coordinates or homogeneous coordinates instead. Then, member functions to set the Cartesian coordinates are expensive. Nevertheless, in current CGAL the types based on the Cartesian representation as well as the types based on the homogeneous representation have both member functions returning Cartesian coordinates and member functions returning homogeneous coordinates. These access functions are provided to make implementing own predicates and operations more convenient.

Like other libraries we use reference counting for the kernel objects. Objects point to a shared representation. Each representation counts the number of objects pointing to it. Copying objects increments the counter of the shared representation, deleting an object decrements the counter of its representation. If the counter reaches zero by the decrement, the representation itself is deleted. The implementation of reference counting is simplified by the non-modifiability of the objects. However, the use of reference counting was not the reason for choosing non-modifiability. Using `copy on write' (a new representation is created for an object whenever its value is changed by a modifying operation), reference counting with modifiable objects is possible and only slightly more involved. Reference counting costs about 15% to 30% runtime for the types double and float, but gains 2% to 11% runtime for the type leda_real.

Basic Library

Following the generic programming paradigm as introduced above, CGAL is made to comply with STL. The interfaces of geometric objects and data structures in the basic library make extensive use of iterators, circulators, and handles (trivial iterator), so that algorithms and data structures can be easily combined with each other and with those provided by STL and other libraries.

An example of a geometric algorithmic problem is the computation of the convex hull. The algorithm takes a set of points and outputs the sequence of extreme points on the boundary of the convex hull. The following program computes the convex hull of 100 random points uniformly distributed in the disc of radius one centered at the origin. The point generator gets as parameter a random source. This random source is initialized with the fixed seed 1 in this example. The result is drawn in a LEDA window, first all points in black, then the hull as polygon in green and finally the vertices in red. The header file hides the usual typedefs and declares all types parameterized with the representation class CGAL::Cartesian<double>.

#include "cartesian_double.h"

int main () {
    Random                   rnd(1);
    Random_points_in_disc_2  rnd_pts( 1.0, rnd);
    list<Point_2>            pts;
    copy_n( rnd_pts, 100, back_inserter( pts));

    Polygon_2  ch;
    CGAL::convex_hull_points_2( pts.begin(), pts.end(), back_inserter(ch));

    Window* window = demo_window();
    Window_iterator_point_2  wout( *window);
    copy( pts.begin(), pts.end(), wout);
    *window << CGAL::GREEN << ch << CGAL::RED;
    copy( ch.vertices_begin(), ch.vertices_end(), wout);

    Point_2  p;
    *window >> p;     // wait for mouse click
    delete window;
    return 0;
}

Triangulations are another example of a container-like data structure in the basic library. Triangulations in CGAL support the incremental construction. The following program is therefore even simpler than the previous one; without using any intermediate container to store all input points, 100 random points are copied into the triangulation data structure.

#include "cartesian_double.h"

int main () {
    Random                    rnd(1);
    Random_points_in_disc_2   rnd_pts( 1.0, rnd);

    Delaunay_triangulation_2  dt;
    copy_n( rnd_pts, 100, back_inserter( dt));

    Window* window = demo_window();
    *window << dt;
    Point_2  p;
    *window >> p;     // wait for mouse click
    delete window;
    return 0;
}

#include "cartesian_double.h"
#include <CGAL/Triangulation_euclidean_traits_xy_3.h>

typedef CGAL::Triangulation_euclidean_traits_xy_3<REP>         Traits;
typedef CGAL::Delaunay_triangulation_2<Traits>       Triangulation_xy;

int main () {
    Triangulation_xy  dt;
    copy( istream_iterator<Point_3>(cin), istream_iterator<Point_3>(),
          back_inserter( dt));
    cout << dt;
    return 0;
}

Geometric Traits Classes

The following example of a convex hull algorithm illustrates the use of a geometric traits class. The algorithm used is Andrew's variant of Graham's scan [Andrew79]. The actual implementation presented here stems from our framework for one sided error predicates [Kettner98b]. The implementation has been modified to use the iterator based interface from above and the traits class.

Andrew's variant of Graham's scan needs only a point type and a leftturn predicate from the geometric kernel given that the input points are already sorted lexicographically. Thus, the geometric traits class is quite short in this example. The leftturn predicate for three points p, q, and r in the plane is true if the points in this order perform a left turn. For points represented in Cartesian coordinates, the predicate is equivalent to the sign of the following determinant:

In the following example, the geometric traits class for the convex hull algorithm is itself a class template parameterized with the number type NT for the point coordinates. NT is the same type as used for the Cartesian point type Point_2. The evaluation of the determinant is implemented as a function object assuming exact arithmetic and the leftturn member function gives access to the function object.

template <class NT> struct Point_2 {
    typedef NT Number_type;
    NT x;
    NT y;
};
template <class NT> struct Convex_hull_traits {
    typedef Point_2<NT> Point;
    struct Leftturn {
        bool operator()( const Point& p, const Point& q, const Point& r) {
            return (q.x-p.x) * (r.y-p.y) > (r.x-p.x) * (q.y-p.y);
        }
    };
    Leftturn leftturn() const { return Leftturn(); }
};

template <class BidirectionalIterator, class OutputIterator, class Traits>
OutputIterator 
convex_hull( BidirectionalIterator first, BidirectionalIterator beyond,
             OutputIterator result, const Traits& traits) {
    typedef typename Traits::Point Point;
    vector<Point> hull;
    hull.push_back( *first); // sentinel
    hull.push_back( *first);

    // lower convex hull (left to right)
    BidirectionalIterator i = first;
    for ( ++i; i != beyond; ++i) {
        while ( traits.leftturn()( hull.end()[-2], *i, hull.back()))
            hull.pop_back();
        hull.push_back( *i);
    }
    // upper convex hull (right to left)
    i = beyond;
    for ( --i; i != first; ) {
        --i;
        while ( traits.leftturn()( hull.end()[-2], *i, hull.back()))
            hull.pop_back();
        hull.push_back( *i);
    }
    // clean up and copy hull to output iterator
    hull.pop_back();
    hull.front() = hull.back();
    hull.pop_back();
    return copy( hull.begin(), hull.end(), result);
}

template <class BidirectionalIterator, class OutputIterator>
OutputIterator 
convex_hull( BidirectionalIterator first, BidirectionalIterator beyond,
             OutputIterator result) {
    typedef typename iterator_traits<BidirectionalIterator>::value_type P;
    typedef typename P::Number_type Number_type;
    typedef Convex_hull_traits<Number_type> Traits;
    return convex_hull( first, beyond, result, Traits());
}

In our framework on one sided error predicates [Kettner98b] we introduced the notion of a conservative implementation of a predicate. If a conservative implementation of the leftturn predicate returns true, the three points perform a left turn, but if it returns false, we do not know the orientation of the points. Thus, the decision errors due to rounding errors in inexact arithmetic are limited to one side of the two possible answers. Useful applications for such predicates will assume that these false answers occur rarely, but in principle an implementation saying always false is a legal implementation. Examples for convex hull algorithms and a triangulation of point sets are given in the paper that compute a well defined output if the predicate is not exact but a conservative implementation. Andrew's variant of Graham's scan as presented above computes a sequence of points of which the points on the convex hull are a subsequence if it is used with a conservative implementation of the predicate. The output can be easily postprocessed with the same algorithm but with an exact implementation of the predicate. For more properties of the computed output and other examples see [Kettner98b].

For floating point arithmetic an error bound can be computed such that if the expression computing the determinant is larger than the error bound, the exact value of the determinant is greater than zero. The expression to compute the determinant and its error bound give us the conservative implementation

(q.x-p.x)*(r.y-p.y) - (q.y-p.y)*(r.x-p.x) > $$8(3u + 6u² + 4u³ + u⁴) B² ,

where u is the unit roundoff of the floating-point number system and B is the absolute value of the maximal coordinate value of the points. We have $$u = 2^-53 for IEEE double precision, and $$u=2^-24 for IEEE single precision floating-point numbers. In the following traits class we assume a built-in type double following the IEEE standard and we make the value B a state value of the traits class. In order to make the error bound representable as double, we round it to $$(3 2^-50 + 2^-100) B² and require B to be a power of two.

class Convex_hull_traits_2 {
    double B;
public:
    typedef Point_2<double> Point;
    Convex_hull_traits_2( double b) : B(b) {}
    struct Leftturn {
        double B;
        Leftturn( double b) : B(b) {}
        bool operator()( const Point& p, const Point& q, const Point& r) {
            const double C = 1.0 / 1024.0 / 1048576.0 / 1048576.0; //2^-50
            return   (q.x-p.x) * (r.y-p.y) - (r.x-p.x) * (q.y-p.y)
                   > (3.0 * C + C * C) * B * B;
        }
    };
    Leftturn leftturn() const { return Leftturn(B); }
};

template <>
struct Convex_hull_traits<double> : public Convex_hull_traits_2 {
    Convex_hull_traits() : Convex_hull_traits_2(1.0) {}
};

Another example of the flexibility of the geometric traits classes is the reconstruction of a terrain from a set of three-dimensional sample points. A common approach is to triangulate the sample points using a Delaunay triangulation in the xy-projection, just ignoring the elevation in the z-coordinate. Similar to the convex hull algorithm, a geometric class can be used to parameterize the two-dimensional triangulation algorithm to work on the three-dimensional data set. See the example for CGAL in the previous Section.