Difference or intersection fail depending on BOOST_GEOMETRY_ROBUSTNESS_ALTERNATIVE

[MauriceHubain]

Depending on if we set BOOST_GEOMETRY_ROBUSTNESS_ALTERNATIVE, either difference or intersection fail for some of our test scenarios.

The test case for difference:

MULTIPOLYGON(((-2.0 -1.5, 2.0 -1.5, 2.0 1.5, -2.0 1.5)))
POLYGON((-0.5 -1.49999999, -2.0 -0.1, -1.99999999 -1.5))

difference(MULTIPOLYGON, POLYGON, result)

result should be like MULTIPOLYGON(((-2.0 -0.1, -0.5 -1.49999999, 2.0 -1.5, 2.0 1.5, -2.0 1.5)))

If we do not define BOOST_GEOMETRY_ROBUSTNESS_ALTERNATIVE, result is empty.
If we define BOOST_GEOMETRY_ROBUSTNESS_ALTERNATIVE, we get the expected result.

The test case for intersection:

MULTIPOLYGON(((1.0 2.0, 3.0 1.0, 4.0 3.0, 2.0 4.0)))
POLYGON((2.0 1.5, 3.0 1.5, 3.0 2.5, 2.0 2.5))

intersection(MULTIPOLYGON, POLYGON, result)

result should be like MULTIPOLYGON(((2.0 1.5, 3.0 1.5, 3.0 2.5, 2.0 2.5))), which is similar to the input POLYGON.

If we do not define BOOST_GEOMETRY_ROBUSTNESS_ALTERNATIVE, we get the expected result.
If we define BOOST_GEOMETRY_ROBUSTNESS_ALTERNATIVE, result contains the bigger input MULTIPOLYGON.

Attached BoostGeometryTests.cpp contains these 2 test cases and some more, which can be enabled by defining ALL_TEST_CASES 1 (in the second line).
BoostGeometryTests.cpp.txt

With Boost 1.71.0 all test cases succeed (in fact all but Intersection_04c1, where Intersection_04c2 is the same test with slightly adapted expectation - 5 points instead of 4).
With Boost 1.80.0 to 1.85.0 multiple tests fail (different, depending on the BOOST_GEOMETRY_ROBUSTNESS_ALTERNATIVE flag).

It is not important, that we get exactly the declared expected results, but reasonable ones.

[vissarion]

Thanks for the detailed examples and code. I can reproduce it. Note that BOOST_GEOMETRY_ROBUSTNESS_ALTERNATIVE is deprecated, and will be removed (probably in the version so we focus on issues without setting this flag. And there are quite a few. At a first look this does not seems as a numerical issue since your instances are failing even with multiprecision arithmetic. Need to look at it in more detail.

[MauriceHubain]

Note that BOOST_GEOMETRY_ROBUSTNESS_ALTERNATIVE is deprecated, and will be removed

This is definitely fine. We just tried that flag, because it seemed to fix our problems with the difference function.

[barendgehrels]

First

Can you next time please provide a minimal scenario? 537 lines is too much to analyze

Second

Herewith an example of such a scenario for the case you list in the description. I cannot reproduce it with Boost 1.85 nor with Boost 1.84

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/geometries.hpp>
#include <iostream>
#include <iomanip>

namespace bg = boost::geometry;

template <typename T>
void test()
{
    using point_t = bg::model::point<T, 2, bg::cs::cartesian>;
    using polygon_t = bg::model::polygon<point_t, false, false>;
    using mp_t = bg::model::multi_polygon<polygon_t>;

    mp_t poly1;
    bg::read_wkt("MULTIPOLYGON(((-2.0 -1.5, 2.0 -1.5, 2.0 1.5, -2.0 1.5)))", poly1);
    bg::correct(poly1);

    mp_t poly2;
    bg::read_wkt("MULTIPOLYGON(((-0.5 -1.49999999, -2.0 -0.1, -1.99999999 -1.5)))", poly2);   
    bg::correct(poly2);

    mp_t diff, inters;
    bg::difference(poly1, poly2, diff);
    bg::intersection(poly1, poly2, inters);

    std::cout << std::setprecision(12) << std::boolalpha;
    std::cout 
        << "Type: " << typeid(T).name()
        << " a: " << bg::area(poly1)
        << " b: " << bg::area(poly2)
        << " diff: " << bg::area(diff)
        << " ints: " << bg::area(inters)
        << " a-b: " << bg::area(poly1) - bg::area(poly2)
        << " d+i " << bg::area(diff) + bg::area(inters)
        << " validity diff: " << bg::is_valid(diff)
        << " ints: " << bg::is_valid(inters)
        << std::endl;
}

int main() {
    test<float>();
    test<double>();
    test<long double>();
}

My output in all cases

Type: f a: 12 b: 1.04999999888 diff: 10.9500000179 ints: 1.04999999888 a-b: 10.9500000011 d+i 12.0000000168 validity diff: true ints: true
Type: d a: 12 b: 1.049999993 diff: 10.95 ints: 1.05 a-b: 10.950000007 d+i 12 validity diff: true ints: true
Type: e a: 12 b: 1.049999993 diff: 10.95 ints: 1.05 a-b: 10.950000007 d+i 12 validity diff: true ints: true

Did you call "correct"?
Can we close it or replace it with a minimal example which can reproduce?

[barendgehrels]

BTW the reason that I checked it again is that I have a concept fix for another recent issue and wanted to check if it would fix this as well.

[MauriceHubain]

If you change your example in a way, that poly2 is a polygon, not a multipolygon (as in our test cases), the test fails for double and long double precision:

[...]
polygon_t poly2;
bg::read_wkt("POLYGON((-0.5 -1.49999999, -2.0 -0.1, -1.99999999 -1.5))", poly2);
bg::correct(poly2);
[...]

Output:

Type: float a: 12 b: 1.04999999888 diff: 10.9500000179 ints: 1.04999999888 a-b: 10.9500000011 d+i 12.0000000168 validity diff: true ints: true
Type: double a: 12 b: 1.049999993 diff: 0 ints: 1.049999993 a-b: 10.950000007 d+i 1.049999993 validity diff: true ints: true
Type: long double a: 12 b: 1.049999993 diff: 0 ints: 1.049999993 a-b: 10.950000007 d+i 1.049999993 validity diff: true ints: true

So at least there is something wrong with this combination of input geometries.
In accordance to the documentation "Behavior: areal (e.g. polygon) - All combinations of: box, ring, polygon, multi_polygon", I would expect, that this combination should be supported.

[MauriceHubain]

Nevertheless, I will check, if changing both input geometries to multipolygon, would solve all our issues.

[barendgehrels]

That should indeed not be necessary, poly should give the same result as multi poly

[MauriceHubain]

Unfortunately I found a first example (from our "real life" issues), where also using mutlipolygon does not give the expected results:

mp_t poly1;
bg::read_wkt("MULTIPOLYGON(((1.2549999979079400 0.85000000411847698, -1.2550000020920500 0.84999999897038103, -1.2549999999999999 -0.85000000102961903, 1.2549999999999999 -0.84999999999999998)))", poly1);
bg::correct(poly1);

mp_t poly2;
bg::read_wkt("MULTIPOLYGON(((-0.87500000000000000 -0.84999999999999998, -0.87500000000000000 -0.070000000000000201, -1.2549999999999999 -0.070000000000000201, -1.2549999999999999 -0.84999999999999998)))", poly2);
bg::correct(poly2);

Output:
Type: float a: 4.26700010347 b: 0.296400005227 diff: 3.97060009825 ints: 0.296400005227 a-b: 3.97060009825 d+i 4.26700010347 validity diff: true ints: true
Type: double a: 4.26700000517 b: 0.2964 diff: 0 ints: 0.2964 a-b: 3.97060000517 d+i 0.2964 validity diff: true ints: true
Type: long double a: 4.26700000517 b: 0.2964 diff: 0 ints: 0.2964 a-b: 3.97060000517 d+i 0.2964 validity diff: true ints: true

[barendgehrels]

That should indeed not be necessary, poly should give the same result as multi poly

and I can reproduce it so it's labeled as a bug now

[barendgehrels]

Unfortunately I found a first example (from our "real life" issues), where also using mutlipolygon does not give the expected results:

mp_t poly1; bg::read_wkt("MULTIPOLYGON(((1.2549999979079400 0.85000000411847698, -1.2550000020920500 0.84999999897038103, -1.2549999999999999 -0.85000000102961903, 1.2549999999999999 -0.84999999999999998)))", poly1); bg::correct(poly1);

mp_t poly2; bg::read_wkt("MULTIPOLYGON(((-0.87500000000000000 -0.84999999999999998, -0.87500000000000000 -0.070000000000000201, -1.2549999999999999 -0.070000000000000201, -1.2549999999999999 -0.84999999999999998)))", poly2); bg::correct(poly2);

Output: Type: float a: 4.26700010347 b: 0.296400005227 diff: 3.97060009825 ints: 0.296400005227 a-b: 3.97060009825 d+i 4.26700010347 validity diff: true ints: true Type: double a: 4.26700000517 b: 0.2964 diff: 0 ints: 0.2964 a-b: 3.97060000517 d+i 0.2964 validity diff: true ints: true Type: long double a: 4.26700000517 b: 0.2964 diff: 0 ints: 0.2964 a-b: 3.97060000517 d+i 0.2964 validity diff: true ints: true

Can you please file other failing cases under new tickets?
Thanks!

[MauriceHubain]

Sure, I can file a new ticket.
But this scenario was already part of the provided test bench.

[barendgehrels]

Sure, I can file a new ticket. But this scenario was already part of the provided test bench.

We can really not work with such a bench, sorry. Because things are fixed one by one, there can (sometimes) be months in between and handled by different persons.

So if one is fixed - the ticket should be closed.

[MauriceHubain]

We can really not work with such a bench, sorry. Because things are fixed one by one, there can (sometimes) be months in between and handled by different persons.

So if one is fixed - the ticket should be closed.

Ok, that's fine.

I filed a new ticket #1299