Feature/robust predicates

[tinko92]

Introduction

This is the first (and most simple and self-contained) of three branches that I created during the GSoC 2019 that has now concluded. It provides two strategies for 2d cartesian points:

• a side-strategy that is basically a robust, adaptive alternative to side_by_triangle (more details below). Given p1, p2, p3 it calculates the signed area (times two) of the triangle p1, p2, p3 and (based on the sign) answers the question: Does p3 lie on the left or right side or exactly on the segment p1, p2.
• an in-circle-strategy for 4 points p1, p2, p3, p, that answers the question: Does p lie in the circumcircle of the triangle formed by p1, p2, p3?

The strategy classes come with documentation and tests, currently in the extension-directories.

Rationale

The new side-strategy can be used where side_by_triangle is used and may help in situations where computations break down because of incorrect results in side_by_triangle.
The in_circle-strategy will be required by algorithms for the construction of Delaunay triangulations such as the one I will propose in a later pull request. It may also be of use in constructing an adaptive, robust side-strategy for circular segments, which were recently proposed by @meastp in the Gitter room and generally useful on its own for users.

Implementation details

Both strategies are based on evaluating certain determinants (like the existing side_by_triangle), They are robust in the sense that they return correct results even on inputs where direct float arithmetics would fail due to precision issues, such as side of POINT(1, 1) with respect to the segment of POINT( 1e-20, 0) and POINT(1e20, 1e20). They are adaptive in the sense that for inputs where direct float arithmetic can provably produce correct results (like in all non-degenerate inputs), the computation terminates early to save time. This is done using analytically computed error bounds. The algorithm is intended for floating-point parameters but returns correct results for integers that can be translated to the calculation_type without loss of precision. It is not robust with respect to inputs that cause an overflow in the computation of interim results.

All of the math behind this is explained in the papers linked here https://www.cs.cmu.edu/~quake/robust.html. The details are implemented almost directly as explained in the paper. All of the implementation details (the mathematics) are in the file
include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/precise_math.hpp

The robustness-parameter is my own addition and can be used to limit the maximum precision. If 0 is provided, the predicates are comparatively fast.

[vissarion]

Thanks for this PR! I am ok with the math, the computations seem ok too. However, the strategy needs some cleaning/revision.

[vissarion]

You could acknowledge GSoC here if you like. See for example https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/formulas/karney_direct.hpp#L5

[vissarion]

this is not very clear. Do you want to say values from 0 to 3, where 0 is the fastest and less accurate and 3 the slowest and most accurate ?

[tinko92]

I understand, I'll reconsider the wording there. It is a little hard to say it in one sentence.
By default (robustness = 3), the strategy will compute up to four successively more precise approximations of the true side_value (with the fourth being guaranteed to be correct) and check for each approximation whether its sign can already be guaranteed using relative error bounds. If that is the case, the remaining approximations (if any) will not be computed and the current approximation will be returned.

If the robustness parameter is set to 0, it will return the first approximation without any test against the error bounds, if it is set to 1, it will compute the first approximation, test it against the error bound and either return if the precision is sufficient or compute the second approximation. With robustness = 2, up to three approximations will be computed.

[vissarion]

For me, the technical details are not needed here. Hence, saying something like 0 is the fastest and less accurate and 3 the slowest and most accurate, is fine.

[vissarion]

unused template parameter

[vissarion]

maybe this is too generic. What is the case where p1 and p2 are of different types?

[tinko92]

I cannot imagine such a case, however, my choice to be that generic was motivated by side_by_triangle having the same interface and I wanted to be as compatible with it as possible.

[barendgehrels]

Yes, the same interface is OK.

[vissarion]

there is a select_most_precise with 3 arguments

[vissarion]

this is implemented in utilities as select_calculation_type_alt

[vissarion]

coordinate_type here is unused

[barendgehrels]

Tinko, thanks a lot for your contributions!
Can you spell the robustness template parameter as Robustness and maybe make it an std::size_t

[barendgehrels]

Yes, the same interface is OK.

[vissarion]

Tinko thanks for the revision, I am ok with merging now.

[vissarion]

For me, the technical details are not needed here. Hence, saying something like 0 is the fastest and less accurate and 3 the slowest and most accurate, is fine.

[awulkiew]

Thanks for this PR! I'm ok with merging after the issues above are addressed.

I'd like to be explicit regarding one thing though. This PR depends on C++11. After merging it the library will partially depend on newer standard. I'm ok with that since we plan to move to C++11 anyway, but I'd like to get explicit confirmation from you @barendgehrels @vissarion.

[vissarion]

@awulkiew I am OK to merge this since it is currently in extensions. Also, since in the discussion here: #590 we haven't reached to a conclusion on a plan to move to newer standards, adding new features that depend on newer standards in extensions could be a start.

[awulkiew]

Thanks Tinko!