|
CGAL 5.3 - dD Geometry Kernel
|
|
CGAL 5.3 - dD Geometry Kernel
|
This part of the reference manual covers the higher-dimensional kernel. The kernel contains objects of constant size, such as point, vector, direction, line, ray, segment, circle. 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. Note that this section partly recapitulates facts already mentioned for the lower-dimensional kernel.
The correctness proof of nearly all geometric algorithms presented in theory papers assumes exact computation with real numbers. This leads to a fundamental problem with the implementation of geometric algorithms. Naively, often the exact real arithmetic is replaced by inexact floating-point arithmetic in the implementation. This often leads to acceptable results for many input data. However, even for the implementation of the simplest geometric algorithms this simplification occasionally does not work. Rounding errors introduced by inaccurate arithmetic may lead to inconsistent decisions, causing unexpected failures for some correct input data. There are many approaches to this problem, one of them is to compute exactly (compute so accurate that all decisions made by the algorithm are exact) which is possible in many cases but more expensive than standard floating-point arithmetic. C. M. Hoffmann [2], [1] illustrates some of the problems arising in the implementation of geometric algorithms and discusses some approaches to solve them. A more recent overview is given in [3]. The exact computation paradigm is discussed by Yap and Dubé [4] and Yap [5].
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.
To increase generic usage of objects and predicates the higher-dimensional kernel makes heavy use of iterator ranges as defined in the STL for modeling tuples. Iterators conceptualize C++ pointers.
For an iterator range [first,last) we define T = tuple [first,last) as the ordered tuple \( (T[0],T[1], \ldots T[d-1])\) where \( S[i] = *++^{(i)}\mathit{first}\) (the element obtained by \( i\) times forwarding the iterator by operator ++ and then dereferencing it to get the value to which it points). We write d = size [first,last) and S = set [first,last) to denote the unordered set of elements of the corresponding tuple.
This extends the syntax of random access iterators to input iterators. If we index the tuple as above then we require that \( ++^{(d)}\mathit{first} = \mathit{last}\).
Our object of study is the \( d\)-dimensional affine Euclidean space, where \( d\) is a parameter of our geometry. Objects in that space are sets of points. A common way to represent the points is the use of Cartesian coordinates, which assumes a reference frame (an origin and \( d\) orthogonal axes). In that framework, a point is represented by a \( d\)-tuple \( (c_0,c_1,\ldots,c_{d-1})\), and so are vectors in the underlying linear space. Each point is represented uniquely by such Cartesian coordinates.
Another way to represent points is by homogeneous coordinates. In that framework, a point is represented by a \( (d+1)\)-tuple \( (h_0,h_1,\ldots,h_d)\). Via the formulae \( c_i = h_i/h_d\), the corresponding point with Cartesian coordinates \( (c_0,c_1,\ldots,c_{d-1})\) can be computed. Note that homogeneous coordinates are not unique. For \( \lambda\ne 0\), the tuples \( (h_0,h_1,\ldots,h_d)\) and \( (\lambda\cdot h_0,\lambda\cdot h_1,\ldots,\lambda\cdot h_d)\) represent the same point. For a point with Cartesian coordinates \( (c_0,c_1,\ldots,c_{d-1})\) a possible homogeneous representation is \( (c_0,c_1,\ldots,c_{d-1},1)\). Homogeneous coordinates in fact allow to represent objects in a more general space, the projective space \( \mathbb{P}^d\). In CGAL, we do not compute in projective geometry. Rather, we use homogeneous coordinates to avoid division operations, since the additional coordinate can serve as a common denominator.
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. In CGAL such a kernel concept is often also called a representation class and denoted by R. A representation class provides the actual implementations of the kernel objects. For all kernel objects Kernel_object of a representation class R based on Cartesian_d or Homogeneous_d, the types CGAL::Kernel_object<R> and R::Kernel_object are identical.
CGAL offers three families of concrete models for the concept representation class, two based on the Cartesian representation of points and one 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 have a constructor with a range of coordinates plus a common denominator (the \( d+1\) homogeneous coordinates of the point). The common interfaces parameterized with a representation 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 and the linear algebra module used to calculate the result of predicates and constructions.
For reasons that will become evident later, a representation class provides two typenames for number types, namely R::FT and R::RT. The type R::FT must fulfill the requirements on what is called a field type in CGAL. This roughly means that R::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 R::RT are weaker. This type must fulfill the requirements on what is called a Euclidean ring type in CGAL. This roughly means that R::RT is a type for which operations \( +\), \( -\), \( *\) are defined with semantics (approximately) corresponding to those of a ring in a mathematical sense. A very limited division operation \( /\) must be available as well. It must work for exact (i.e., no remainder) integer divisions only. Furthermore, both number types should fulfill CGAL's requirements on a number type.
With Cartesian_d<FieldNumberType,LinearAlgebra> you can choose Cartesian representation of coordinates. The type LinearAlgebra must me a linear algebra module working on numbers of type FieldNumberType. The second parameter defaults to module delivered with the kernel so for short a user can just write Cartesian_d<FieldNumberType> when not providing her own linear algebra.
When you choose Cartesian representation you have to declare at least the type of the coordinates. A number type used with the Cartesian_d representation class should be a field type as described above. Both Cartesian_d<FieldNumberType>::FT and Cartesian_d<FieldNumberType>::RT are mapped to number type FieldNumberType. Cartesian_d<FieldNumberType,LinearAlgebra>::LA is mapped to the type LinearAlgebra. Cartesian_d<FieldNumberType> uses reference counting internally to save copying costs.
As we mentioned before, homogeneous coordinates permit to avoid division operations in numerical computations, since the additional coordinate can serve as a common denominator. Avoiding divisions can be useful for exact geometric computation. With Homogeneous_d<RingNumberType,LinearAlgebra> you can choose homogeneous representation of coordinates with the kernel objects. As for Cartesian representation you have to declare at the same time the type used to store the homogeneous coordinates. Since the homogeneous representation allows one to avoid the divisions, the number type associated with a homogeneous representation class must be a model for the weaker concept Euclidean ring type only.
The type LinearAlgebra must me a linear algebra module working on numbers of type RingNumberType. Again the second parameter defaults to module delivered with the kernel so for short one can just write Homogeneous_d<RingNumberType> when replacing the default is no issue.
However, some operations provided by this kernel involve division operations, for example computing squared distances or returning a Cartesian coordinate. To keep the requirements on the number type parameter of Homogeneous low, the number type Quotient<RingNumberType> is used instead. This number type turns a ring type into a field type. It maintains numbers as quotients, i.e. a numerator and a denominator. Thereby, divisions are circumvented. With Homogeneous_d<RingNumberType>, Homogeneous_d<RingNumberType>::FT is equal to Quotient<RingNumberType> while Homogeneous_d<RingNumberType>::RT is equal to RingNumberType. Homogeneous_d<RingNumberType,LinearAlgebra>::LA is mapped to the type LinearAlgebra.
The kernel Epick_d<DimensionTag>, short for Exact Predicates Inexact Constructions Kernel is a kernel particularly useful when the dimension of the space is known at compile-time; The template parameter DimensionTag is then Dimension_tag<d> where d is an integer representing the dimension. It may also be used with parameter Dynamic_dimension_tag, in which case the dimension does not need to be known at compile-time. It uses a Cartesian representation and supports construction of points from double coordinates. It provides exact geometric predicates, but the geometric constructions are not guaranteed to be exact.
Note that it provides few interfaces in addition to those documented in the Kernel_d concept. In particular, the type of a point is only available as Epick_d<DimensionTag>::Point_d, not Point_d<Epick_d<DimensionTag>>.
The kernel Epeck_d<DimensionTag>, short for Exact Predicates Exact Constructions Kernel is a kernel particularly useful when the dimension of the space is known at compile-time; The template parameter DimensionTag is then Dimension_tag<d> where d is an integer representing the dimension. It may also be used with parameter Dynamic_dimension_tag, in which case the dimension does not need to be known at compile-time. It uses a Cartesian representation and supports construction of points from double coordinates. It provides exact geometric predicates and constructions.
Note that it provides few interfaces in addition to those documented in the Kernel_d concept. In particular, the type of a point is only available as Epeck_d<DimensionTag>::Point_d, not Point_d<Epeck_d<DimensionTag>>.
The use of representation classes does not only avoid problems, it also makes all CGAL classes very uniform. Like Cartesian_d<double>::Point_d, they always consist of:
The capitalized base name of the geometric object, such as Point, Segment, Triangle.
Followed by _d.
Cartesian_d<double> or Homogeneous_d<leda_integer>, where the type can be found, except for Epick_d<DimensionTag> where the number type is implicitly double. Algorithms and data structures in the basic library of CGAL are parameterized by a traits class that subsumes the objects on which the algorithm or data structure operates as well as the operations to do so. For most of the algorithms and data structures in the basic library you can use a kernel as a traits class. For some algorithms you even do not have to specify the kernel; it is detected automatically using the types of the geometric objects passed to the algorithm. In some other cases, the algorithms or data structures need more than is provided by a kernel. In these cases, a kernel can not be used as a traits class.
If you start with integral Cartesian coordinates, many geometric computations will involve integral numerical values only. Especially, this is true for geometric computations that evaluate only predicates, which are tantamount to determinant computations. Examples are triangulation of point sets and convex hull computation.
The dimension \( d\) of our affine space determines the dimension of the matrix computations in the mathematical evaluation of predicates. As rounding errors accumulate fast the homogeneous representation used with multi-precision integers is the kernel of choice for well-behaved algorithms. Note, that unless you use an arbitrary precision integer type, incorrect results might arise due to overflow.
If new points are to be constructed, for example the intersection point of two lines, computation of Cartesian coordinates usually involves divisions, so you need to use a field type with Cartesian representation or have to switch to homogeneous representation. double is a possible, but imprecise field type. You can also put any ring type into Quotient to get a field type and put it into Cartesian_d, but you better put the ring type into Homogeneous. leda_rational and leda_real are valid field types, too.
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.
The Epick_d kernel provides a compromise using double coordinates. It evaluates predicates exactly, which is slower than plain double computations, but still faster than using an exact number type thanks to filtering techniques. Constructions are inexact, computed with double.
You need just to include a representation class to obtain the geometric objects of the kernel that you would like to use with the representation class, i.e., CGAL/Cartesian_d.h or CGAL/Homogeneous_d.h
In CGAL, we strictly distinguish between points, vectors and directions. A point is a point in the Euclidean space \( \E^d\), a vector is the difference of two points \( p_2\), \( p_1\) and denotes the direction and the distance from \( p_1\) to \( p_2\) in the vector space \( \mathbb{R}^d\), and a direction is a vector where we forget about its length. They are different mathematical concepts. For example, they behave different under affine transformations and an addition of two points is meaningless in affine geometry. By putting them in different classes we not only get cleaner code, but also type checking by the compiler which avoids ambiguous expressions. Hence, it pays twice to make this distinction.
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\).
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_1\) and \( p_2\), you can writeyou might call midpoint(p_1,p_2) instead
Note that these constructions do not involve any performance overhead for the conversion with the currently available representation classes.
Besides points (R::Point_d), vectors (R::Vector_d), and directions (R::Direction_d), CGAL provides lines, rays, segments, hyperplanes, and spheres.
Lines (R::Line_d) in CGAL are oriented. A ray (R::Ray_d) is a 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 (R::Segment_d) 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.
Hyperplanes are affine subspaces of dimension \( d-1\) in \( \E^d\), passing through \( d\) points. Hyperplanes are oriented and partition space into a positive side and a negative side. In CGAL, there are no special classes for halfspaces. Halfspaces are supposed to be represented by oriented hyperplanes. All kernel objects are equality comparable via operator== and operator!=. For those oriented objects whose orientation can be reversed (segments, lines, hyperplanes, spheres) there is also a global function weak_equality() that allows to test for point set equality disregarding the orientation.
Geometric objects in CGAL have member functions that test the position of a point relative to the object. Full dimensional objects and their boundaries are represented by the same type, e.g. halfspaces and hyperplanes are not distinguished, neither are balls and spheres. Such objects split the ambient space into two full-dimensional parts, a bounded part and an unbounded part (e.g. spheres), or two unbounded parts (e.g. hyperplanes). By default these objects are oriented, i.e., one of the resulting parts is called the positive side, the other one is called the negative side. Both of these may be unbounded.
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 segment in \( d\)-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()
Predicates are at the heart of a geometry kernel. They are basic units for the composition of geometric algorithms and encapsulate decisions. Hence their correctness is crucial for the control flow and hence for the correctness of an implementation of a geometric algorithm. CGAL uses the term predicate in a generalized sense. Not only components returning a Boolean value are called predicates but also components returning an enumeration type like a Comparison_result or an Orientation. We say components, because predicates are implemented both as functions and function objects (also called functors and provided by a kernel class).
CGAL provides predicates for the orientation of point sets (orientation), for comparing points according to some given order, especially for comparing Cartesian coordinates (e.g. lexicographically_xy_smaller), in-sphere tests, and predicates to compare distances.
Functions and function objects that generate objects that are neither of type bool nor enum types are called constructions. Constructions involve computation of new numerical values and may be imprecise due to rounding errors unless a kernel with an exact number type is used.
Affine transformations (R::Aff_transformation_d) allow to generate new object instances under arbitrary affine transformations. These transformations include translations, rotations (within planes) and scaling. Most of the geometric objects in a kernel have a member function transform(Aff_transformation_d t) which applies the transformation to the object instance.
CGAL also provides a set of functions that detect or compute the intersection between objects 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:
sqrt operation can be costly. Even if it is not very costly for a specific number type and platform, avoiding it is always cheaper. sqrt operation is defined, especially integer types and rationals. Intersections on kernel objects currently cover only those objects that are part of flats (R::Segment_d, R::Ray_d, R::Line_d, and R::Hyperplane_d). For any pair of objects \( o1\), \( o2\) of these types the operation intersection(o1,o2) returns a boost::optional< boost::variant< T... > > where T... is a list of all possible resulting geometric objects.
The return type of intersecting two objects of the types Type1 and Type2 can be specified through the placeholder type specifier auto.
In the following example, the auto is used for the intersection computation, as there are possibly different return values.
For testing where a point \( p\) lies with respect to a hyperplane defined by an array \( P\) of points \( p_1\), ... , \( p_d\), one may be tempted to construct the hyperplane R::Hyperplane_d(d,P,P+d) and use the method oriented_side(p). This may pay off if many tests with respect to the plane are made. Nevertheless, unless the number type is exact, the constructed plane is only approximated, and round-off errors may lead oriented_side(p) to return an orientation which is different from the orientation of \( p_1\), ... , \( p_d\), \( p\).
In CGAL, we provide predicates in which such geometric decisions are made directly with a reference to the input points in \( P\) without an intermediary object like a plane. For the above test, the recommended way to get the result is to use \(\mathrm{orientation}(P',P'+d)\), where \( P'\) is an array containing the points \( p_1\), ... , \( p_d\), \( p\).
For exact number types like leda_real, 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).
This higher-dimensional kernel is the result of a long evolving development. A first version of the kernel was offered as a LEDA extension package ddgeo by Kurt Mehlhorn and Michael Seel. The original design was driven by the realization of a d-dimensional convex hull data type developed at the Max-Planck Institut für Informatik.
The code base was discussed and reviewed within the CGAL kernel group (of the low-dimensional kernel). This led to the identification of the concept interfaces and in parallel to adaptations according to the evolution of the low-dimensional kernel. The kernel was revised based on suggestions by Hervé Brönnimann, Michael Hoffmann, and Stefan Schirra.
Marc Glisse added Epick_d in 2014 and Epeck_d in 2019.
This work was supported by ESPRIT IV Long Term Research Projects No. 21957 (CGAL) and No. 28155 (GALIA).
The Epick_d kernel was partially supported by the IST Programme of the EU (FET Open) Project under Contract No IST-25582 – (CGL - Computational Geometric Learning).
1.8.13