1 // Boost.Geometry (aka Boost.Geometry, Generic Geometry Library) 2 // 3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. 4 // Use, modification and distribution is subject to the Boost Software License, 5 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at 6 // http://www.boost.org/LICENSE_1_0.txt) 7 8 #ifndef _DOXYGEN_ROBUSTNESS_HPP 9 #define _DOXYGEN_ROBUSTNESS_HPP 10 11 /--------------------------------------------------------------------------------------------------- 12 13 /*! 14 \page robustness Boost.Geometry and Robustness 15 16 \section robustness_par1 Introduction 17 18 Floating point coordinates have limited precision. Geometry algorithms have to take this into account. 19 20 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 \b 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. 21 22 See for more backgrounds e.g.: 23 24 - <a href="http://www.mpi-inf.mpg.de/~kettner/pub/nonrobust_cgta_06.pdf">Classroom Examples of Robustness Problems in Geometric Computations</a> 25 - <a href="http://groups.csail.mit.edu/graphics/classes/6.838/S98/meetings/m12/pred/m12.html">Robust Predicates and Degeneracy</a> 26 27 Boost.Geometry is aware of these issues and provides several approaches to minimize the problems, or avoid them completely using 28 29 - <a href="http://en.wikipedia.org/wiki/GNU_Multi-Precision_Library">GMP</a> 30 - <a href="http://en.wikipedia.org/wiki/Class_Library_for_Numbers">CLN</a> 31 32 33 \section robustness_par2 Example 34 35 An example. Consider the elongated triangle and box below. 36 37 \image html robust_triangle_box.png 38 39 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. 40 41 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. 42 43 \section robustness_par3 Coordinate types 44 45 All geometry types provided by Boost.Geometry, and types by the user, do have a coordinate type. For example 46 47 \code 48 boost::geometry::point_xy<float> p1; 49 boost::geometry::point_xy<double> p2; 50 boost::geometry::point_xy<long double> p3; 51 \endcode 52 53 describes three points with different coordinate types, 54 a 32 bits <a href="http://en.wikipedia.org/wiki/Single_precision">float</a>, 55 a 64 bits <a href="http://en.wikipedia.org/wiki/Double_precision_floating-point_format">double</a>, 56 a <a href="http://en.wikipedia.org/wiki/Long_double">long double</a>, not standardized by IEEE and is on some machines 96 bits (using a MSVC compiler it is a synonym for a double). 57 58 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 \b select_most_precise. 59 60 Boost.Geometry supports also high precision arithmetic types, by adaption. The numeric_adaptor, used for adaption, is not part of Boost.Geometry 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: 61 62 - GMP (http://gmplib.org) 63 - CLN (http://www.ginac.de/CLN) 64 65 Note that the libraries themselves are not included in Boost.Geometry, they are completely independant of each other. 66 67 These numeric types can be used as following: 68 \code 69 boost::geometry::point_xy<boost::numeric_adaptor::gmp_value_type> p4; 70 boost::geometry::point_xy<boost::numeric_adaptor::cln_value_type> p5; 71 \endcode 72 73 All algorithms using these points will use the \b GMP resp. \b CLN library for calculations. 74 75 \section robustness_par4 Calculation types 76 77 If high precision arithmetic types are used as shown above, coordinates are stored in these points. That is not always necessary. Therefore, Boost.Geometry 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: 78 79 Example: 80 The code below shows the calculation of the area. Points are stored in double; calculation is done using \b GMP 81 82 \code 83 { 84 typedef boost::geometry::point_xy<double> point_type; 85 boost::geometry::linear_ring<point_type> ring; 86 ring.push_back(boost::geometry::make<point_type>(0.0, 0.0)); 87 ring.push_back(boost::geometry::make<point_type>(0.0, 0.0012)); 88 ring.push_back(boost::geometry::make<point_type>(1234567.89012345, 0.0)); 89 ring.push_back(ring.front()); 90 91 typedef boost::numeric_adaptor::gmp_value_type gmp; 92 93 gmp area = boost::geometry::area(ring, boost::geometry::strategy::area::by_triangles<point_type, gmp>()); 94 std::cout << area << std::endl; 95 } 96 \endcode 97 98 Above shows how this is used to use \b GMP or \b CLN for double coordinates. Exactly the same mechanism works (of course) also to do calculation in double, where coordinates are stored in float. 99 100 \section robustness_par5 Strategies 101 102 In the previous section was shown that strategies have an optional template parameter \b CalculationType so enhance precision. However, the design of Boost.Geometry 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. 103 104 \section robustness_par6 Examples 105 106 We show here some things which can occur in challenging domains. 107 108 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. 109 110 \image html robust_float.png 111 112 If we perform the intersection using Boost.Geometry, we get the effect that is presented on the pictures below, using float (left) and using double (right). 113 114 \image html robust_triangles.png 115 116 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. 117 118 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 Boost.Geometry implementation using float, double or long double will give some artifacts. 119 120 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). 121 122 If using \b GMP or \b CLN, the intersection is correct. 123 124 \image html robust_stars.png 125 126 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. 127 128 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 \b GMP to use Boost.Geometry without any problem. 129 130 \section robustness_par7 Future work 131 132 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 Boost.Geometry is "not checked on 100% robustness". As pointed out in the discussions on the Boost mailing list in spring '09 (a.o. <i>"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"</i>), 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. 133 134 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. 135 136 \section robustness_par8 Summary 137 138 Boost.Geometry approach to support high precision: 139 140 - library users can use points with float, double, or long double coordinates (the default) 141 - for higher numerical robustness: users can call algorithms using another calculation type (e.g. \b GMP or \b CLN, ...) 142 - for higher numerical robustness: users can use points with \b GMP or \b CLN coordinates 143 - users can implement their own implementation and provide this as a strategy, the Boost.Geometry design allows this 144 - other measures can be implemented, described as future work. 145 146 147 148 */ 149 150 #endif // _DOXYGEN_ROBUSTNESS_HPP 151