Chapter 11
2D and 3D Geometry Kernel

Table of Contents

11.1   Introduction

This part of the reference manual covers the kernel. The kernel contains objects of constant size, such as point, vector, direction, line, ray, segment, triangle, iso-oriented rectangle and tetrahedron. With each type comes a set of functions which can be applied to an object of this type. You will typically find access functions (e.g. to the coordinates of a point), tests of the position of a point relative to the object, a function returning the bounding box, the length, or the area of an object, and so on. The Cgal kernel further contains basic operations such as affine transformations, detection and computation of intersections, and distance computations.

11.1.1   Robustness

In Cgal you can choose the underlying number types and arithmetic. You can use different types of arithmetic simultaneously and the choice can be easily changed, e.g. for testing. So you can choose between implementations with fast but occasionally inexact arithmetic and implementations guaranteeing exact computation and exact results. Of course you have to pay for the exactness in terms of execution time and storage space. See the dedicated chapter for more details on number types and their capabilities and performance.

11.2   Kernel Representations

11.2.1   Genericity through Parameterization

Almost all the kernel objects (and the corresponding functions) are templates with a parameter that allows the user to choose the representation of the kernel objects. A type that is used as an argument for this parameter must fulfill certain requirements on syntax and semantics. The list of requirements defines an abstract kernel concept. For all kernel objects types, the types CGAL::Type<Kernel> and Kernel::Type are identical.

CGAL offers four families of concrete models for the concept Kernel, two based on the Cartesian representation of points and two based on the homogeneous representation of points. The interface of the kernel objects is designed such that it works well with both Cartesian and homogeneous representation. For example, points in 2D have a constructor with three arguments as well (the three homogeneous coordinates of the point). The common interfaces parameterized with a kernel class allow one to develop code independent of the chosen representation. We said ``families'' of models, because both families are parameterized too. A user can choose the number type used to represent the coordinates.

For reasons that will become evident later, a kernel class provides two typenames for number types, namely Kernel::FT and Kernel::RT. The type Kernel::FT must fulfill the requirements on what is called a FieldNumberType in Cgal. This roughly means that Kernel::FT is a type for which operations +, -, * and / are defined with semantics (approximately) corresponding to those of a field in a mathematical sense. Note that, strictly speaking, the built-in type int does not fulfill the requirements on a field type, since ints correspond to elements of a ring rather than a field, especially operation / is not the inverse of *. The requirements on the type Kernel::RT are weaker. This type must fulfill the requirements on what is called a RingNumberType in Cgal. This roughly means that Kernel::RT is a type for which operations +, -, * are defined with semantics (approximately) corresponding to those of a ring in a mathematical sense.

11.2.2   Cartesian Kernels

Cartesian<FieldNumberType> uses reference counting internally to save copying costs. CGAL also provides Simple_cartesian<FieldNumberType>, a kernel that uses Cartesian representation but no reference counting. Debugging is easier with Simple_cartesian<FieldNumberType>, since the coordinates are stored within the class and hence direct access to the coordinates is possible. Depending on the algorithm, it can also be slightly more or less efficient than Cartesian<FieldNumberType>. Again, in Simple_cartesian<FieldNumberType> both Simple_cartesian<FieldNumberType>::FT and Simple_cartesian<FieldNumberType>::RT are mapped to FieldNumberType.

11.2.3   Homogeneous Kernels

Homogeneous<RingNumberType> uses reference counting internally to save copying costs. CGAL also provides Simple_homogeneous<RingNumberType>, a kernel that uses homogeneous representation but no reference counting. Debugging is easier with Simple_homogeneous<RingNumberType>, since the coordinates are stored within the class and hence direct access to the coordinates is possible. Depending on the algorithm, it can also be slightly more or less efficient than Homogeneous<RingNumberType>. Again, in Simple_homogeneous<RingNumberType> the type Simple_homogeneous<RingNumberType>::FT is equal to Quotient<RingNumberType> while Simple_homogeneous<RingNumberType>::RT is equal to RingNumberType.

11.2.4   Naming conventions

The use of kernel classes not only avoids problems, it also makes all Cgal classes very uniform. They always consist of:

1. The capitalized base name of the geometric object, such as Point, Segment, or Triangle.

2. An underscore followed by the dimension of the object, for example _2, _3, or _d.

3. A kernel class as parameter, which itself is parameterized with a number type, such as Cartesian<double> or Homogeneous<leda_integer>.

11.2.5   Kernel as a Traits Class

11.2.6   Choosing a Kernel and Predefined Kernels

If new points are to be constructed, for example the intersection point of two lines, computation of Cartesian coordinates usually involves divisions. Hence, one needs to use a FieldNumberType with Cartesian representation, or alternatively, switch to homogeneous representation. The type double is a - though imprecise - model for FieldNumberType. You can also put any RingNumberType into the Quotient adaptor to get a field type which then can be put into Cartesian. But using homogeneous representation on the RingNumberType is usually the better option. Other valid FieldNumberTypes are leda_rational and leda_real.

If it is crucial for you that the computation is reliable, the right choice is probably a number type that guarantees exact computation. The Filtered_kernel provides a way to apply filtering techniques [BBP01] to achieve a kernel with exact and efficient predicates. Still other people will prefer the built-in type double, because they need speed and can live with approximate results, or even algorithms that, from time to time, crash or compute incorrect results due to accumulated rounding errors.

Predefined kernels. For the user's convenience, Cgal provides 3 typedefs to generally useful kernels.

• They are all Cartesian kernels.
• They all support constructions of points from double Cartesian coordinates.
• All these 3 kernels provide exact geometric predicates.
• They handle geometric constructions differently:
  • Exact_predicates_exact_constructions_kernel: provides exact geometric constructions, in addition to exact geometric predicates.
  • Exact_predicates_exact_constructions_kernel_with_sqrt: same as Exact_predicates_exact_constructions_kernel, but the number type it provides (Exact_predicates_exact_constructions_kernel_with_sqrt::FT) supports the square root operation exactly ¹.
  • Exact_predicates_inexact_constructions_kernel: provides exact geometric predicates, but geometric constructions may be inexact due to round-off errors. It is however enough for most Cgal algorithms, and faster than both Exact_predicates_exact_constructions_kernel and Exact_predicates_exact_constructions_kernel_with_sqrt.

11.3   Kernel Geometry

11.3.1   Points and Vectors

^d

Cgal defines a symbolic constant ORIGIN of type Origin which denotes the point at the origin. This constant is used in the conversion between points and vectors. Subtracting it from a point p results in the locus vector of p.

  Point_2< Cartesian<double> >  p(1.0, 1.0), q;
  Vector_2< Cartesian<double> >  v;
  v = p - ORIGIN;
  q = ORIGIN + v;  
  assert( p == q );

In order to obtain the point corresponding to a vector v you simply have to add v to ORIGIN. If you want to determine the point q in the middle between two points p₁ and p₂, you can write²

  q = p_1 + (p_2 - p_1) / 2.0;

Note that these constructions do not involve any performance overhead for the conversion with the currently available representation classes.

11.3.2   Kernel Objects

Lines (Line_2<Kernel>, Line_3<Kernel>) in Cgal are oriented. In two-dimensional space, they induce a partition of the plane into a positive side and a negative side. Any two points on a line induce an orientation of this line. A ray (Ray_2<Kernel>, Ray_3<Kernel>) is semi-infinite interval on a line, and this line is oriented from the finite endpoint of this interval towards any other point in this interval. A segment (Segment_2<Kernel>, Segment_3<Kernel>) is a bounded interval on a directed line, and the endpoints are ordered so that they induce the same direction as that of the line.

Planes are affine subspaces of dimension two in ³, passing through three points, or a point and a line, ray, or segment. Cgal provides a correspondence between any plane in the ambient space ³ and the embedding of ² in that space. Just like lines, planes are oriented and partition space into a positive side and a negative side. In Cgal, there are no special classes for half-spaces. Half-spaces in 2D and 3D are supposed to be represented by oriented lines and planes, respectively.

Concerning polygons and polyhedra, the kernel provides triangles, iso-oriented rectangles, iso-oriented cuboids and tetrahedra. More complex polygons³ and polyhedra or polyhedral surfaces can be obtained from the basic library (Polygon_2, Polyhedron_3), so they are not part of the kernel. As with any Jordan curves, triangles, iso-oriented rectangles and circles separate the plane into two regions, one bounded and one unbounded.

11.3.3   Orientation and Relative Position

For these objects there is a function oriented_side() that determines whether a test point is on the positive side, the negative side, or on the oriented boundary. These function returns a value of type Oriented_side.

Those objects that split the space in a bounded and an unbounded part, have a member function bounded_side() with return type Bounded_side.

If an object is lower dimensional, e.g. a triangle in three-dimensional space or a segment in two-dimensional space, there is only a test whether a point belongs to the object or not. This member function, which takes a point as an argument and returns a Boolean value, is called has_on()

11.4   Predicates and Constructions

11.4.1   Predicates

Cgal provides predicates for the orientation of point sets (orientation, leftturn, rightturn, collinear, coplanar), for comparing points according to some given order, especially for comparing Cartesian coordinates (e.g. lexicographically_xy_smaller), in-circle and in-sphere tests, and predicates to compare distances.

11.4.2   Constructions

Affine transformations (Aff_transformation_2<Kernel>, Aff_transformation_3<Kernel>) allow to generate new object instances under arbitrary affine transformations. These transformations include translations, rotations (in 2D only) and scaling. Most of the geometric objects in a kernel have a member function transform(Aff_transformation t) which applies the transformation to the object instance.

Cgal also provides a set of functions that detect or compute the intersection between objects of the 2D kernel, and many objects in the 3D kernel, and functions to calculate their squared distance . Moreover, some member functions of kernel objects are constructions.

So there are routines that compute the square of the Euclidean distance, but no routines that compute the distance itself. Why? First of all, the two values can be derived from each other quite easily (by taking the square root or taking the square). So, supplying only the one and not the other is only a minor inconvenience for the user. Second, often either value can be used. This is for example the case when (squared) distances are compared. Third, the library wants to stimulate the use of the squared distance instead of the distance. The squared distance can be computed in more cases and the computation is cheaper. We do this by not providing the perhaps more natural routine, The problem of a distance routine is that it needs the sqrt operation. This has two drawbacks:

• The sqrt operation can be costly. Even if it is not very costly for a specific number type and platform, avoiding it is always cheaper.
• There are number types on which no sqrt operation is defined, especially integer types and rationals.

11.4.3   Polymorphic Return Values

Example

{
    typedef Cartesian<double>  K;
    typedef K::Point_2         Point_2;
    typedef K::Segment_2       Segment_2;

    Segment_2 segment_1, segment_2;

    std::cin >> segment_1 >> segment_2;

    Object obj = intersection(segment_1, segment_2);

    if (const Point_2 *point = object_cast<Point_2>(&obj)) {
        /* do something with *point */
    } else if (const Segment_2 *segment = object_cast<Segment_2>(&obj)) {
        /* do something with *segment*/
    }

    /*  there was no intersection */
}

The intersection routine itself looks roughly as follows:

template < class Kernel >
Object  intersection(Segment_2<Kernel> s1, Segment_2<Kernel> s2)
{

    if (/* intersection in a point */ ) {

       Point_2<Kernel> p = ... ;
       return make_object(p);

    } else if (/* intersection in a segment */ ) {

       Segment_2<Kernel> s = ... ;
       return make_object(s);
    }
    return Object();
}

11.4.4   Constructive Predicates

In Cgal, we provide predicates in which such geometric decisions are made directly with a reference to the input points p, q, r, s, without an intermediary object like a plane. For the above test, the recommended way to get the result is to use orientation(p,q,r,s). For exact number types, the situation is different. If several tests are to be made with the same plane, it pays off to construct the plane and to use oriented_side(p).

11.5   Extensible Kernel

This manual section describe how users can plug user defined geometric classes in existing Cgal kernels. This is best illustrated by an example.

11.5.1   Introduction

Cgal defines the concept of a geometry kernel. Such a kernel provides types, construction objects and generalized predicates. Most implementations of Computational Geometry algorithms and data structures in the basic library of Cgal were done in a way that classes or functions can be parametrized with a geometric traits class.

In most cases this geometric traits class must be a model of the Cgal geometry kernel concept (but there are some exceptions).

11.5.2   An Extensive Example

Assume you have the following point class, where the coordinates are stored in an array of doubles, where we have another data member color, which shows up in the constructor.

class MyPointC2 {

private:
  double vec[2];
  int col;

public:

  MyPointC2()
    : col(0)
  {
    *vec = 0;
    *(vec+1) = 0;
  }

  
  MyPointC2(const double x, const double y, int c)
    : col(c)
  {
    *vec = x;
    *(vec+1) = y;
  }

  const double& x() const  { return *vec; }

  const double& y() const { return *(vec+1); }

  double & x() { return *vec; }

  double& y() { return *(vec+1); }

  int color() const { return col; }

  int& color() { return col; }
  
  
  bool operator==(const MyPointC2 &p) const
  {
    return ( *vec == *(p.vec) )  && ( *(vec+1) == *(p.vec + 1) && ( col == p.col) );
  }

  bool operator!=(const MyPointC2 &p) const
  {
      return !(*this == p);
  }

};

As said earlier the class is pretty minimalistic, for example it has no bbox() method. One might assume that a basic library algorithm which computes a bounding box (e.g, to compute the bounding box of a polygon), will not compile. Luckily it will, because it does not use of member functions of geometric objects, but it makes use of the functor Kernel::Construct_bbox_2.

To make the right thing happen with MyPointC2 we have to provide the following functor.

template <class ConstructBbox_2>
class MyConstruct_bbox_2 : public ConstructBbox_2 {
public:
  using ConstructBbox_2::operator();

  CGAL::Bbox_2 operator()(const typename MyPointC2& p) const {
    return CGAL::Bbox_2(p.x(), p.y(), p.x(), p.y());
  }
};

Things are similar for random access to the Cartesian coordinates of a point. As the coordinates are stored in an array of doubles we can use double* as random access iterator.

class MyConstruct_coord_iterator {
public:
  const double* operator()(const MyPointC2& p)
  {
    return &p.x();
  }

  const double* operator()(const MyPointC2& p, int)
  {
    const double* pyptr = &p.y();
    pyptr++;
    return pyptr;
  }
};

The last functor we have to provide is the one which constructs points. That is you are not forced to add the constructor with the Origin as parameter to your class, nor the constructor with homogeneous coordinates. The functor is a kind of glue layer between the Cgal algorithms and your class.

 template <typename K>
  class MyConstruct_point_2
  {
    typedef typename K::RT         RT;
    typedef typename K::Point_2    Point_2;
  public:
    typedef Point_2          result_type;

    Point_2
    operator()(CGAL::Origin o) const
    { return Point_2(0,0, 0); }

    Point_2
    operator()(const RT& x, const RT& y) const
    { return Point_2(x, y, 0); }

    
    // We need this one, as such a functor is in the Filtered_kernel
    Point_2
    operator()(const RT& x, const RT& y, const RT& w) const
    { 
      if(w != 1){
	return Point_2(x/w, y/w, 0); 
      } else {
	return Point_2(x,y, 0);
      }
    }
  };

Now we are ready to put the puzzle together. We won't explain it in detail, but you see that there are typedefs to the new point class and the functors. All the other types are inherited.

File: examples/Kernel_23/MyKernel.h

#ifndef MYKERNEL_H
#define MYKERNEL_H

#include <CGAL/Cartesian.h>
#include "MyPointC2.h"
#include "MySegmentC2.h"

// K_ is the new kernel, and K_Base is the old kernel
template < typename K_, typename K_Base >
class MyCartesian_base
  : public K_Base::template Base<K_>::Type
{
  typedef typename K_Base::template Base<K_>::Type   OldK;
public:
  typedef K_                                Kernel;
  typedef MyPointC2                         Point_2;
  typedef MySegmentC2<Kernel>               Segment_2;
  typedef MyConstruct_point_2<Kernel, OldK>       Construct_point_2;
  typedef const double*                     Cartesian_const_iterator_2;
  typedef MyConstruct_coord_iterator        Construct_cartesian_const_iterator_2;
  typedef MyConstruct_bbox_2<typename OldK::Construct_bbox_2>
                                            Construct_bbox_2;

  Construct_point_2
  construct_point_2_object() const
  { return Construct_point_2(); }

  Construct_bbox_2
  construct_bbox_2_object() const
  { return Construct_bbox_2(); }

  Construct_cartesian_const_iterator_2
  construct_cartesian_const_iterator_2_object() const
  { return Construct_cartesian_const_iterator_2(); }

  template < typename Kernel2 >
  struct Base { typedef MyCartesian_base<Kernel2, K_Base>  Type; };
};


template < typename FT_ >
struct MyKernel
  : public CGAL::Type_equality_wrapper<
                MyCartesian_base<MyKernel<FT_>, CGAL::Cartesian<FT_> >,
                MyKernel<FT_> >
{};

#endif // MYKERNEL_H

The kernel itself can be made robust by plugging it in the Filtered_kernel.

File: examples/Kernel_23/MyKernel.cpp

#include <CGAL/basic.h>
#include <CGAL/Filtered_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/squared_distance_2.h>
#include <cassert>
#include "MyKernel.h"


typedef MyKernel<double>                   MK;
typedef CGAL::Filtered_kernel_adaptor<MK>  K;
typedef CGAL::Delaunay_triangulation_2<K>  Delaunay_triangulation_2;


typedef K::Point_2         Point;
typedef K::Segment_2       Segment;
typedef K::Ray_2           Ray;
typedef K::Line_2          Line;
typedef K::Triangle_2      Triangle;
typedef K::Iso_rectangle_2 Iso_rectangle;

const int RED= 1;
const int BLACK=2;

int main()
{
  Point a(0,0), b(1,0), c(1,1), d(0,1);
  a.color()=RED;
  b.color()=BLACK;
  d.color()=RED;

  Delaunay_triangulation_2 dt;
  dt.insert(a);

  K::Orientation_2 orientation;
  orientation(a,b,c);

  Point p(1,2), q;
  p.color() = RED;
  q.color() = BLACK;
  std::cout << p << std::endl;

  K::Compute_squared_distance_2 squared_distance;

  std::cout << "squared_distance(a, b) == "
            << squared_distance(a, b) << std::endl;

  Segment s1(p,q), s2(a, c);

  K::Construct_midpoint_2 construct_midpoint_2;

  Point mp = construct_midpoint_2(p,q);

  assert(s1.source().color() == RED);


  K::Intersect_2 intersection;

  CGAL::Object o = intersection(s1, s2);

  K::Construct_cartesian_const_iterator_2 construct_it;
  K::Cartesian_const_iterator_2  cit = construct_it(a);
  assert(*cit == a.x());

  cit = construct_it(a,0);

  cit--;
  assert(*cit == a.y());

  Line l1(a,b), l2(p, q);

  intersection(l1, l2);


  intersection(s1, l1);

  Ray r1(d,b), r2(d,c);
  intersection(r1, r2);

  intersection(r1, l1);

  squared_distance(r1, r2);
  squared_distance(r1, l2);
  squared_distance(r1, s2);

  Triangle t1(a,b,c), t2(a,c,d);
  intersection(t1, t2);
  intersection(t1, l1);

  intersection(t1, s1);

  intersection(t1, r1);

  Iso_rectangle i1(a,c), i2(d,p);
  intersection(i1, i2);
  intersection(i1, s1);

  intersection(i1, r1);
  intersection(i1, l1);

  t1.orientation();

  std::cout << s1.source() << std::endl;

  std::cout << t1.bbox() << std::endl;

  std::cout << "done" << std::endl;
  return 0;
}

11.5.3   Limitations

The point class must have member functions x() and y() (and z() for the 3d point). We will probably introduce function objects that take care of coordinate access.

As we enforce type equality between MyKernel::Point_2 and Point_2<MyKernel>, the constructor with the color as third argument is not available.

11.6   Design and Implementation History

At a meeting at Utrecht University in January 1995, Olivier Devillers, Andreas Fabri, Wolfgang Freiseisen, Geert-Jan Giezeman, Mark Overmars, Stefan Schirra, Otfried Schwarzkopf (now Otfried Cheong), and Sven Schönherr discussed the foundations of the Cgal kernel. Many design and software engineering issues were addressed, e.g. naming conventions, coupling of classes (flat versus deep class hierarchy), memory allocation, programming conventions, mutability of atomic objects, points and vectors, storing additional information, orthogonality of operations on the kernel objects, viewing non-constant-size objects like polygons as dynamic data structures (and hence not as part of the (innermost) kernel).

The people attending the meeting delegated the compilation of a draft specification to Stefan Schirra. The resulting draft specification was intentionally modeled on Cgal's precursors C++gal and Plageo as well as on the geometric part of Leda. The specification already featured coexistence of Cartesian and homogeneous representation of point/vector data and parameterization by number type(s). During the discussion of the draft a kernel design group was formed. The members of this group were Andreas Fabri, Geert-Jan Giezeman, Lutz Kettner, Stefan Schirra, and Sven Schönherr. The work of the kernel design group led to significant changes and improvements of the original design, e.g. the strong separation between points and vectors. Probably the most important enhancement was the design of a common superstructure for the previously uncoupled Cartesian and homogeneous representations. One can say, that the kernel was designed by this group. The kernel was later revised based on suggestions by Hervé Brönnimann, Bernd Gärtner, Michael Hoffmann, and Lutz Kettner.

A first version of the kernel was internally made available at the beginning of the Cgal-project (esprit ltr iv project number 21957). Since then many more people contributed to the evolution of the kernel through discussions on the Cgal mailing lists. The implementation based on Cartesian representation was (initially) provided by Andreas Fabri, the homogeneous representation (initially) by Stefan Schirra. Intersection and distance computations were implemented by Geert-Jan Giezeman. Further work has been done by Susan Hert on the overall maintenance of the kernel. Philippe Guigue has provided efficient intersection tests for 3D triangles. Andreas Fabri, Michael Hoffmann and Sylvain Pion have improved the support for the extensibility and adaptability of the kernel. Pedro Machado Manhães de Castro and Monique Teillaud introduced 3D circles. In 2010, Pierre Alliez, Stéphane Tayeb and Camille Wormser added intersection constructions for 3D triangles and efficient intersection tests for bounding boxes.

11.6.1   Acknowledgment

This work was supported by the Graduiertenkolleg `Algorithmische Diskrete Mathematik', under grant DFG We 1265/2-1, and by ESPRIT IV Long Term Research Projects No. 21957 (CGAL) and No. 28155 (GALIA).

Footnotes