CGAL 6.1.2 - Manual

Loading...
Searching...
No Matches
Hello World

Author
  CGAL Editorial Board

This tutorial is for the CGAL newbie, who knows C++ and has a basic knowledge of geometric algorithms. The first section shows how to define a point and segment class, and how to apply geometric predicates on them. The section further raises the awareness that that there are serious issues when using floating point numbers for coordinates. In the second section, you will meet a typical CGAL function, which computes a 2D convex hull. The third section shows what we mean with a Traits class, and the fourth section explains the notion of concept and model.

# Three Points and One Segment

In this first example, we demonstrate how to construct some points and a segment, and perform some basic operations on them.

All CGAL header files are in the subdirectory `include/CGAL`. All CGAL classes and functions are in the namespace `CGAL`. Classes start with a capital letter, global function with a lowercase letter, and constants are all uppercase. The dimension of an object is expressed with a suffix.

The geometric primitives, like the point type, are defined in a kernel. The kernel we have chosen for this first example uses `double` precision floating point numbers for the Cartesian coordinates of the point.

Besides the types we see predicates like the orientation test for three points, and constructions like the distance and midpoint computation. A predicate has a discrete set of possible results, whereas a construction produces either a number, or another geometric entity.

File Kernel_23/points_and_segment.cpp

#include <iostream>
#include <CGAL/Simple_cartesian.h>
typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Point_2 Point_2;
typedef Kernel::Segment_2 Segment_2;
int main()
{
Point_2 p(1,1), q(10,10);
std::cout << "p = " << p << std::endl;
std::cout << "q = " << q.x() << " " << q.y() << std::endl;
std::cout << "sqdist(p,q) = "
<< CGAL::squared_distance(p,q) << std::endl;
Segment_2 s(p,q);
Point_2 m(5, 9);
std::cout << "m = " << m << std::endl;
std::cout << "sqdist(Segment_2(p,q), m) = "
<< CGAL::squared_distance(s,m) << std::endl;
std::cout << "p, q, and m ";
switch (CGAL::orientation(p,q,m)){
case CGAL::COLLINEAR:
std::cout << "are collinear\n";
break;
case CGAL::LEFT_TURN:
std::cout << "make a left turn\n";
break;
case CGAL::RIGHT_TURN:
std::cout << "make a right turn\n";
break;
}
std::cout << " midpoint(p,q) = " << CGAL::midpoint(p,q) << std::endl;
return 0;
}
Kernel::Point_2
Kernel::Segment_2
CGAL::RIGHT_TURN
const CGAL::Orientation RIGHT_TURN
CGAL::LEFT_TURN
const CGAL::Orientation LEFT_TURN
CGAL::COLLINEAR
const CGAL::Orientation COLLINEAR
CGAL::midpoint
CGAL::Point_2< Kernel > midpoint(const CGAL::Point_2< Kernel > &p, const CGAL::Point_2< Kernel > &q)
CGAL::orientation
Orientation orientation(const CGAL::Point_2< Kernel > &p, const CGAL::Point_2< Kernel > &q, const CGAL::Point_2< Kernel > &r)
CGAL::squared_distance
Kernel::FT squared_distance(Type1< Kernel > obj1, Type2< Kernel > obj2)
Kernel
CGAL::Simple_cartesian

To do geometry with floating point numbers can be surprising as the next example shows.

File Kernel_23/surprising.cpp

#include <iostream>
#include <CGAL/Simple_cartesian.h>
typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Point_2 Point_2;
int main()
{
{
Point_2 p(0, 0.3), q(1, 0.6), r(2, 0.9);
std::cout << (CGAL::collinear(p,q,r) ? "collinear\n" : "not collinear\n");
}
{
Point_2 p(0, 1.0/3.0), q(1, 2.0/3.0), r(2, 1);
std::cout << (CGAL::collinear(p,q,r) ? "collinear\n" : "not collinear\n");
}
{
Point_2 p(0,0), q(1, 1), r(2, 2);
std::cout << (CGAL::collinear(p,q,r) ? "collinear\n" : "not collinear\n");
}
return 0;
}
CGAL::collinear
bool collinear(const CGAL::Point_2< Kernel > &p, const CGAL::Point_2< Kernel > &q, const CGAL::Point_2< Kernel > &r)

Reading the code, we could assume that it would print three times "collinear". However the actual output is the following:

    not collinear
    not collinear
    collinear

This is because these fractions are not representable as double-precision numbers, and the collinearity test will internally compute a determinant of a 3x3 matrix which is close but not equal to zero, and hence the non collinearity for the first two tests.

Something similar can happen with points that perform a left turn, but due to rounding errors during the determinant computation, it seems that the points are collinear, or perform a right turn.

If you must ensure that your numbers get interpreted at their full precision you can use a CGAL kernel that performs exact predicates and exact constructions.

File Kernel_23/exact.cpp

#include <iostream>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <sstream>
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::Point_2 Point_2;
int main()
{
Point_2 p(0, 0.3), q, r(2, 0.9);
{
q = Point_2(1, 0.6);
std::cout << (CGAL::collinear(p,q,r) ? "collinear\n" : "not collinear\n");
}
{
std::istringstream input("0 0.3 1 0.6 2 0.9");
input >> p >> q >> r;
std::cout << (CGAL::collinear(p,q,r) ? "collinear\n" : "not collinear\n");
}
{
q = CGAL::midpoint(p,r);
std::cout << (CGAL::collinear(p,q,r) ? "collinear\n" : "not collinear\n");
}
return 0;
}
CGAL::Exact_predicates_exact_constructions_kernel

Here comes the output and you may still be surprised.

    not collinear
    collinear
    collinear

In the first block the points are still not collinear, for the simple reason that the coordinates you see as text get turned into floating point numbers. When they are then turned into arbitrary precision rationals, they exactly represent the floating point number, but not the text!

This is different in the second block, which corresponds to reading numbers from a file. The arbitrary precision rationals are then directly constructed from a string so that they represent exactly the text.

In the third block you see that constructions as midpoint constructions are exact, just as the name of the kernel type suggests.

In many cases, you will have floating point numbers that are "exact", in the sense that they were computed by some application or obtained from a sensor. They are not the string "0.1" or computed on the fly as "1.0/10.0", but a full precision floating point number. If they are input to an algorithm that makes no constructions, you can use a kernel that provides exact predicates but inexact constructions. One such example is the convex hull algorithm which we will see in the next section. The output is a subset of the input, and the algorithm only compares coordinates and performs orientation tests.

At a first glance the kernel doing exact predicates and constructions seems to be the perfect choice, but performance requirements or limited memory resources make that it is not. Furthermore, for many algorithms it is irrelevant to do exact constructions. For example a surface mesh simplification algorithm that iteratively contracts an edge by collapsing it to the midpoint of the edge.

Most CGAL packages explain which kind of kernel they should use or support.

# The Convex Hull of a Sequence of Points

All examples in this section compute the 2D convex hull of a set of points. We show that algorithms get their input as a begin/end iterator pair denoting a range of points, and that they write the result (in the example the points on the convex hull) into an output iterator.

## The Convex Hull of Points in a Built-in Array

In the first example, we have as input an array of five points. As the convex hull of these points is a subset of the input, it is safe to provide an array for storing the result which has the same size.

File Convex_hull_2/array_convex_hull_2.cpp

#include <iostream>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/convex_hull_2.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
int main()
{
Point_2 points[5] = { Point_2(0,0), Point_2(10,0), Point_2(10,10), Point_2(6,5), Point_2(4,1) };
Point_2 result[5];
Point_2 *ptr = CGAL::convex_hull_2( points, points+5, result );
std::cout << ptr - result << " points on the convex hull:" << std::endl;
for(int i = 0; i < ptr - result; i++){
std::cout << result[i] << std::endl;
}
return 0;
}
CGAL::Exact_predicates_inexact_constructions_kernel
CGAL::convex_hull_2
OutputIterator convex_hull_2(InputIterator first, InputIterator beyond, OutputIterator result, const Traits &ch_traits)

We have seen in the previous section that CGAL comes with several kernels. Since the convex hull algorithm only makes comparisons of coordinates and orientation tests of input points, we can choose a kernel that provides exact predicates, but no exact geometric constructions.

The convex hull function takes three arguments, the start and past-the-end pointer for the input, and the start pointer of the array for the result. The function returns the pointer into the result array just behind the last convex hull point written, so the pointer difference tells us how many points are on the convex hull.

## The Convex Hull of Points in a Vector

In the second example, we replace the built-in array by an `std::vector` of the Standard Template Library.

File Convex_hull_2/vector_convex_hull_2.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/convex_hull_2.h>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
typedef std::vector<Point_2> Points;
int main()
{
Points points, result;
points.push_back(Point_2(0,0));
points.push_back(Point_2(10,0));
points.push_back(Point_2(10,10));
points.push_back(Point_2(6,5));
points.push_back(Point_2(4,1));
CGAL::convex_hull_2( points.begin(), points.end(), std::back_inserter(result) );
std::cout << result.size() << " points on the convex hull" << std::endl;
return 0;
}

We put some points in the vector, calling the `push_back()` method of the `std::vector` class.

We then call the convex hull function. The first two arguments, `points.begin()` and `points.end()` are iterators, which are a generalization of pointers: they can be dereferenced and incremented. The convex hull function is generic in the sense that it takes as input whatever can be dereferenced and incremented.

The third argument is where the result gets written to. In the previous example we provided a pointer to allocated memory. The generalization of such a pointer is the output iterator, which allows to increment and assign a value to the dereferenced iterator. In this example we start with an empty vector which grows as needed. Therefore, we cannot simply pass it `result.begin()`, but an output iterator generated by the helper function `std::back_inserter(result)`. This output iterator does nothing when incremented, and calls `result.push_back(..)` on the assignment.

If you know the STL, the Standard Template Library, the above makes perfect sense, as this is the way the STL decouples algorithms from containers. If you don't know the STL, you maybe better first familiarize yourself with its basic ideas.

# About Kernels and Traits Classes

In this section, we show how we express the requirements that must be fulfilled in order that a function like `convex_hull_2()` can be used with an arbitrary point type.

If you look at the manual page of the function `convex_hull_2()` and the other 2D convex hull algorithms, you see that they come in two versions. In the examples we have seen so far, the function that takes two iterators for the range of input points and an output iterator for writing the result to. The second version has an additional template parameter `Traits`, and an additional parameter of this type.

template<class InputIterator , class OutputIterator , class Traits >
OutputIterator
convex_hull_2(InputIterator first,
InputIterator beyond,
OutputIterator result,
const Traits & ch_traits)
InputIterator
Concept from the C++ standard. See https://en.cppreference.com/w/cpp/named_req/InputIterator.
Definition: General.txt:143
OutputIterator
Concept from the C++ standard. See https://en.cppreference.com/w/cpp/named_req/OutputIterator.
Definition: General.txt:138

What are the geometric primitives a typical convex hull algorithm uses? Of course, this depends on the algorithm, so let us consider what is probably the simplest efficient algorithm, the so-called "Graham/Andrew Scan". This algorithm first sorts the points from left to right, and then builds the convex hull incrementally by adding one point after another from the sorted list. To do this, it must at least know about some point type, it should have some idea how to sort those points, and it must be able to evaluate the orientation of a triple of points.

And that is where the template parameter `Traits` comes in. For `ch_graham_andrew()` it must provide the following nested types:

- `Traits::Point_2`
- `Traits::Less_xy_2`
- `Traits::Left_turn_2`
- `Traits::Equal_2`

As you can guess, `Left_turn_2` is responsible for the orientation test, while `Less_xy_2` is used for sorting the points. The requirements these types have to satisfy are documented in full with the concept `ConvexHullTraits_2`.

The types are regrouped for a simple reason. The alternative would have been a rather lengthy function template, and an even longer function call.

template <class InputIterator, class OutputIterator, class Point_2, class Less_xy_2, class Left_turn_2, class Equal_2>
OutputIterator
ch_graham_andrew( InputIterator first,
InputIterator beyond,
OutputIterator result);
CGAL::ch_graham_andrew
OutputIterator ch_graham_andrew(InputIterator first, InputIterator beyond, OutputIterator result, const Traits &ch_traits=Default_traits)

There are two obvious questions: What can be used as argument for this template parameter? And why do we have template parameters at all?

To answer the first question, any model of the CGAL concept `Kernel` provides what is required by the concept `ConvexHullTraits_2`.

As for the second question, think about an application where we want to compute the convex hull of 3D points projected into the `yz` plane. Using the class `Projection_traits_yz_3` this is a small modification of the previous example.

File Convex_hull_2/convex_hull_yz.cpp

#include <iostream>
#include <iterator>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Projection_traits_yz_3.h>
#include <CGAL/convex_hull_2.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K3;
typedef CGAL::Projection_traits_yz_3<K3> K;
typedef K::Point_2 Point_2;
int main()
{
std::istream_iterator< Point_2 > input_begin( std::cin );
std::istream_iterator< Point_2 > input_end;
std::ostream_iterator< Point_2 > output( std::cout, "\n" );
CGAL::convex_hull_2( input_begin, input_end, output, K() );
return 0;
}
CGAL::Projection_traits_yz_3

Another example would be about a user defined point type, or a point type coming from a third party library other than CGAL. Put the point type together with the required predicates for this point type in the scope of a class, and you can run `convex_hull_2()` with these points.

Finally, let us explain why a traits object that is passed to the convex hull function? It would allow to use a more general projection traits object to store state, for example if the projection plane was given by a direction, which is hardwired in the class `Projection_traits_yz_3`.

# Concepts and Models

In the previous section, we wrote that: any model of the CGAL concept `Kernel` provides what is required by the concept `ConvexHullTraits_2`.

A concept is a set of requirements on a type, namely that it has certain nested types, certain member functions, or comes with certain free functions that take the type as it. A model of a concept is a class that fulfills the requirements of the concept.

Let's have a look at the following function.

template <typename T>
T
duplicate(T t)
{
return t;
}

If you want to instantiate this function with a class `C`, this class must at least provide a copy constructor, and we say that class `C` must be a model of `CopyConstructible`. A singleton class does not fulfill this requirement.

Another example is the function:

template <typename T>
T& std::min(const T& a, const T& b)
{
return (a<b)?a:b;
}

This function only compiles if the `operator<(..)` is defined for the type used as `T`, and we say that the type must be a model of `LessThanComparable`.

An example for a concept with required free functions is the `HalfedgeListGraph` in the CGAL package CGAL and the Boost Graph Library. In order to be a model of `HalfedgeListGraph` a class `G` there must be a global function `halfedges(const G&)`, etc.

An example for a concept with a required traits class is `InputIterator`. For a model of an `InputIterator` a specialization of the class `std::iterator_traits` must exist (or the generic template must be applicable).

# Further Reading

We also recommend the standard text books "The C++ Standard Library, A Tutorial and Reference" by Nicolai M. Josuttis from Addison-Wesley, or "Generic Programming and the STL" by Matthew H. Austern for the STL and its notion of concepts and models.

Other resources for CGAL are the rest of the tutorials and the user support page at https://www.cgal.org/.