GGL and Robustness

Introduction

Floating point coordinates have limited precision. Geometry algorithms have to take this into account.

If differences between points are very small, it may lead to false result of a mathematical calculation performed on such points, what in turn, may cause algorithm result as inadequate to actual geometric situation. For example, a point which is located left from a segment, but very close to it, can be reported on that segment or right from it. Also if differences are a little larger, artifacts can be shown in geometric algorithms.

See for more backgrounds e.g.:

GGL is aware of these issues and provides several approaches to minimize the problems, or avoid them completely using

Example

An example. Consider the elongated triangle and box below.

The left edge of the triangle has a length of about the precision of the floating point grid. It is not possible to do this intersection correctly, using floating point. No library (using floating point) can do that, by nature of float point numerical representation. It is not possible to express the four different coordinates in the resulting intersected polygon. Theoretically distinct points will be assigned to the very same location.

Also if the left edge is longer than that, or using double coordinates, those effects will be there. And of course not only in triangles, but any spiky feature in a polygon can result in non representable coordinates or zero-length edges.

Coordinate types

All geometry types provided by GGL, and types by the user, do have a coordinate type. For example

ggl::point_xy<float> p1;
ggl::point_xy<double> p2;
ggl::point_xy<long double> p3;

describes three points with different coordinate types, a 32 bits float a 64 bits double a long double, not standardized by IEEE and is on some machines 96 bits (using a MSVC compiler it is a synonym for a double).

By default, algorithms select the coordinate type of the input geometries. If there are two input geometries, and they have different coordinate types, the coordinate type with the most precision is selected. This is done by the meta-function select_most_precise.

GGL supports also high precision arithmetic types, by adaption. The numeric_adaptor, used for adaption, is not part of GGL itself but developed by us and sent (as preview) to the Boost List (as it turned out, that functionality might also be provided by Boost.Math bindings, but the mechanism is the same). Types from the following libraries are supported:

GMP (http://gmplib.org)
CLN (http://www.ginac.de/CLN)

Note that the libraries themselves are not included in GGL, they are completely independant of each other.

These numeric types can be used as following:

ggl::point_xy<boost::numeric_adaptor::gmp_value_type> p4;
ggl::point_xy<boost::numeric_adaptor::cln_value_type> p5;

All algorithms using these points will use the GMP resp. CLN library for calculations.

Calculation types

If high precision arithmetic types are used as shown above, coordinates are stored in these points. That is not always necessary. Therefore, GGL provides a second approach. It is possible to specify that only the calculation is done in high precision. This is done by specifying a strategy. For example, in area:

Example:

The code below shows the calculation of the area. Points are stored in double; calculation is done using GMP

{
    typedef ggl::point_xy<double> point_type;
    ggl::linear_ring<point_type> ring;
    ring.push_back(ggl::make<point_type>(0.0, 0.0));
    ring.push_back(ggl::make<point_type>(0.0, 0.0012));
    ring.push_back(ggl::make<point_type>(1234567.89012345, 0.0));
    ring.push_back(ring.front());

    typedef boost::numeric_adaptor::gmp_value_type gmp;

    gmp area = ggl::area(ring, ggl::strategy::area::by_triangles<point_type, gmp>());
    std::cout << area << std::endl;
}

Above shows how this is used to use GMP or CLN for double coordinates. Exactly the same mechanism works (of course) also to do calculation in double, where coordinates are stored in float.

Strategies

In the previous section was shown that strategies have an optional template parameter CalculationType so enhance precision. However, the design of GGL also allows that the user can implement a strategy himself. In that case he can implement the necessary predicates, or use the necessary floating point types, or even go to integer and back. Whatever he prefers.

Examples

We show here some things which can occur in challenging domains.

The image below is drawn in PowerPoint to predict what would happen at an intersection of two triangles using float coordinates in the range 1e-45.

If we perform the intersection using GGL, we get the effect that is presented on the pictures below, using float (left) and using double (right).

This shows how double can solve issues which might be present in float. However, there are also cases which cannot be solved by double or long double. And there are also cases which give more deviations than just a move of the intersection points.

We investigated this and created an artificial case. In this case, there are two stars, they are the same but one of them is rotated over an interval of about 1e-44. When those stars are intersected, the current GGL implementation using float, double or long double will give some artifacts.

Those artifacts are caused by taking the wrong paths on points where distances are zero (because they cannot be expressed in the used coordinate systems).

If using GMP or CLN, the intersection is correct.

Note again, these things happen in differences of about 1e-45. We can investigate if they can be reduced or sometimes even solved. However, they belong to the floating point approach, which is not exact.

For many users, this all is not relevant. Using double they will be able to do all operations without any problems or artefacts. They can occur in real life, where differences are very small, or very large. Those users can use GMP to use GGL without any problem.

Future work

There are several methods to avoid instability and we don't know them all, some of them might be applicable to our algorithms. Therefore is stated that GGL is "not checked on 100% robustness". As pointed out in the discussions on the Boost mailing list in spring '09 (a.o. "The dependent concept should explicitely require unlimited precision since the library specifically uses it to *guarantee* robustness. [...] Users should be left with no choice but to pick some external component fulfilling the unlimited precision requirement"), it seems that it is not possible to solve all these issues using any FP number, that it is necessary to use a library as GMP or CLN for this guarantee.

Therefore we decided to go for supporting high precision numeric types first, and they are currently supported in most algorithms (a.o. area, length, perimeter, distance, centroid, intersection, union). However, we certainly are willing to take other measures as well.

Summary

The GGL approach to support high precision:

library users can use points with float, double, or long double coordinates (the default)
for higher numerical robustness: users can call algorithms using another calculation type (e.g. GMP or CLN, ...)
for higher numerical robustness: users can use points with GMP or CLN coordinates
users can implement their own implementation and provide this as a strategy, the GGL design allows this
other measures can be implemented, described as future work.