Wrong result of sym_difference.

[awulkiew]

To reproduce (polygon's point order doesn't matter):

namespace bg = boost::geometry;
namespace bgm = boost::geometry::model;

using cpt_t = bgm::point<double, 2, bg::cs::cartesian>;
using cpo_t = bgm::polygon<cpt_t, false, true>;
using cmpo_t = bgm::multi_polygon<cpo_t>;

int main()
{
    cpo_t a, b, c, d;
    bg::read_wkt("POLYGON((25 0,0 15,30 15,22 10,25 0))", a);
    bg::read_wkt("POLYGON((5 0,15 25,25 0,15 5,5 0))", b);
    bg::read_wkt("POLYGON((15 0,17 10,10 15,20 15,25 25,30 15,40 15,32 10,35 0,25 5,15 0))", c);
    bg::read_wkt("POLYGON((5 0,7 10,0 15,10 15,15 25,20 15,30 15,22 10,25 0,15 5,5 0))", d);
    bg::correct(a);
    bg::correct(b);
    bg::correct(c);
    bg::correct(d);

    cmpo_t ab, abc, dc, abcdc;
    bg::sym_difference(a, b, ab);       // ok
    bg::sym_difference(ab, c, abc);     // ok
    bg::sym_difference(d, c, dc);       // ok
    bg::sym_difference(abc, dc, abcdc); // wrong!
}

Wrong result for geometries:

abc = MULTIPOLYGON(((20 2.5,25 0,25 -8.8817841970012523e-16,20.454545454545453 2.7272727272727271,20 2.5)),((20 2.5,15.909090909090908 4.5454545454545459,15 0,20 2.5)),((20.454545454545453 2.7272727272727271,23.333333333333332 4.166666666666667,19 15,15 25,11 15,10.777777777777779 14.444444444444443,17 10,16.071428571428573 5.3571428571428577,20.454545454545453 2.7272727272727271)),((16.071428571428573 5.3571428571428577,8.870967741935484 9.67741935483871,5 0,15 5,15.90909090909091 4.5454545454545459,16.071428571428573 5.3571428571428577)),((23.333333333333332 4.166666666666667,25 0,23.695652173913043 4.3478260869565215,23.333333333333332 4.166666666666667)),((10.777777777777779 14.444444444444443,10 15,0 15,8.870967741935484 9.67741935483871,10.777777777777779 14.444444444444443)),((23.695652173913043 4.3478260869565215,25 5,35 0,32 10,40 15,30 15,22 10,23.695652173913043 4.3478260869565215)),((30 15,25 25,20 15,30 15)))
dc = MULTIPOLYGON(((20 15,15 25,10 15,20 15)),((20 15,30 15,25 25,20 15)),((15.90909090909091 4.5454545454545459,17 10,10 15,0 15,7 10,5 0,15 5,15.90909090909091 4.5454545454545459)),((15.90909090909091 4.5454545454545459,15 0,20 2.5,15.90909090909091 4.5454545454545459)),((20 2.5,25 0,23.695652173913043 4.3478260869565215,20 2.5)),((23.695652173913043 4.3478260869565215,25 5,35 0,32 10,40 15,30 15,22 10,23.695652173913043 4.3478260869565215)))

Visual help:

[awulkiew]

@barendgehrels I'm debugging the case with robustness enabled. sym_difference() calculates difference two times and then union of the results here: https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/algorithms/sym_difference.hpp#L134-L152 In all of these calculations the same rescale policy is used. The one that was calculated for the initial input geometries. (with diff12 as blue and diff21 as orange)

Inside union turns and self-turns of geometry2 are caluclated here: https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/algorithms/detail/overlay/overlay.hpp#L296 and here: https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/algorithms/detail/overlay/overlay.hpp#L316 (with geometry1 as blue and geometry2 as orange)

Note that there several wrong turns here. Turn 2 (regular turn) and 8 (self turn) are caused by a small segment generated by difference. These are the points of the nearly rectangular polygon in the middle:

{10.777778284173866, 14.444445710434662}
{8.8709677419354840, 9.6774193548387100}
{16.071428650510203, 5.3571432525510234}
{17.000000000000000, 10.000000000000000}
{10.777777750841748, 14.444444377104368}
{10.777778284173866, 14.444445710434662}

So when the rescale policy of the initial geometries is used this polygon is detected as self-intersecting. So this is one issue. Difference should not generate an invalid polygon. Note that if you check validity of diff21 or this single polygon the geometry is considered to be valid because of different rescale policy. This however shouldn't prevent the geometry to be generated because the algorithm follows u.

Turns 0 and 1 would probably cause a spike but probably shouldn't prevent the geometry to be generated. Or am I wrong?

Anyway AFAIU the problem are turns 5, 6, and 7. It looks like if the operations were reversed for them. Furthermore if the operations were reversed for turns 0 and 1 there wouldn't be a spike. Do you have an idea what may be wrong?

Btw here are the points of geometry1:

{20.454545435878362, 2.7272729639275339}
{22.056447156464479, 3.5282235782322404}
{23.333333333333332, 4.1666666666666670}
{19.000000000000000, 15.000000000000000}
{11.000000000000000, 15.000000000000000}
{10.777778284173866, 14.444445710434662}
{17.000000000000000, 10.000000000000000}
{16.071429972789407, 5.3571424591834358}
{20.454545435878362, 2.7272729639275339}

and the lower-right triangle of geometry2:

{20.454545435878362, 2.7272729639275339}
{25.000000000000000, 0.0000000000000000}
{23.333333398692986, 4.1666665032675381}
{22.056447156464479, 3.5282235782322404}
{20.454545435878362, 2.7272729639275339}

[vissarion]

I have tried running the case with different number types than double and got different results for all the combinations below. Note that the result in the first row (i.e. geometry mp1) seems correct and it is the only valid geometry among the four.

number type	BOOST_GEOMETRY_NO_ROBUSTNESS	result
float	no	mp1
float	yes	mp2
long double	yes	mp3
long double	no	mp4

mp1 = MULTIPOLYGON(((8.87097 9.67742,16.0714 5.35714,20.4545 2.72727,25 0,23.3333 4.16667,23.3333 4.16667,19 15,11 15,8.87097 9.67742)),((11 15,15 25,10 15,11 15)),((19 15,20 15,15 25,19 15)),((8.87097 9.67742,0 15,7 10,5 0,8.87097 9.67742)))
mp2 = MULTIPOLYGON(((23.3333 4.1666700000000105,19 15,11 15,10.7778 14.444400000000002,8.87097 9.677419999999998,16.0714 5.357140000000001,20.4545 2.7272700000000043,20.7244 2.8622199999999935,23.3333 4.1666700000000105)),((11 15,15 24.999999999999986,10 15,11 15)),((19 15,20 15,15 24.999999999999986,19 15)),((20.4545 2.7272700000000043,20 2.5,25 0,25 5.117280039712568e-7,20.4545 2.7272700000000043)),((25 9.536740037674463e-7,23.3333 4.1666700000000105,20.4545 2.7272700000000043,25 0.000001254829996355511,25 9.536740037674463e-7)),((8.87097 9.677419999999998,0 15,7.000000000000001 10,5 0,8.87097 9.677419999999998)))
mp3 = MULTIPOLYGON(((20.4545 2.72727,25 0,25 4.33681e-19,23.3333 4.16667,20.4545 2.72727)),((20.4545 2.72727,23.3333 4.16667,19 15,11 15,10.7778 14.4444,8.87097 9.67742,16.0714 5.35714,15.9091 4.54545,15 0,20 2.5,20.4545 2.72727)),((11 15,15 25,10 15,11 15)),((19 15,20 15,15 25,19 15)),((8.87097 9.67742,0 15,7 10,5 0,8.87097 9.67742)))
mp4 = MULTIPOLYGON(((16.0714 5.357140000000001,20.4545 2.7272700000000043,25 0,23.3333 4.1666700000000105,19 15,11 15,10.7778 14.444400000000002,17 10,10.7778 14.444400000000002,8.87097 9.677419999999998,16.0714 5.357140000000001)),((11 15,15 24.999999999999986,10 15,11 15)),((19 15,20 15,15 24.999999999999986,19 15)),((8.87097 9.677419999999998,0 15,7.000000000000001 10,5 0,8.87097 9.677419999999998)))

EDIT: the above results are with 1.77 (develop produces different results that also depend on number types).

[barendgehrels]

Thanks. I can (edit: partly) reproduce it with current branch (fix/other_types).

@vissarion for me long double seems always correct, but I get the difference in double (wrong with rescaling) and float (wrong without rescaling). Note that I'm not on 1.77. Because rescaling will be removed (hopefully in next version), and we don't trust float in all cases, it's becoming a non-issue but still there are some weird things. I'll add a testcase anyway.

Areas are of abc dc abcdc respectively

Without rescaling (BOOST_GEOMETRY_NO_ROBUSTNESS) Areas 363.96 363.933 23.2258 type f Areas 363.96 363.933 136.452 type d Areas 363.96 363.933 136.452 type e Areas 243.451 363.933 101.015 type m

With rescaling Areas 363.96 363.933 136.452 type f Areas 363.96 363.933 23.2258 type d Areas 363.96 363.933 136.452 type e Areas 243.451 363.933 101.015 type m

where m = boost::multiprecision bm::number<bm::cpp_bin_float<100>>

m is supposed to give better results, but there is something wrong, already in the steps before. I'll continue next week.

[barendgehrels]

My code here https://godbolt.org/z/9T6PoePoz

[vissarion]

Areas may give results that seems correct but the generated geometries does not look correct. In some cases there are missing parts (this is reflected in area computations) but is some other cases there are "spikes" that slightly affect the area but return weird geometries. See bellow illustrations for all combinations of number types with and without rescaling.

EDIT: figures are ploted with WKT-playground for simplicity, the geometries even if plotted on map are cartesian

	float	double	long double	m
With rescaling
Without rescaling

[barendgehrels]

Areas may give results that seems correct but the generated geometries does not look correct. In some cases there are missing parts (this is reflected in area computations) but is some other cases there are "spikes" that slightly affect the area but return weird geometries. See bellow illustrations for all combinations of number types with and without rescaling.

Indeed! Thanks. There are some spikes. We should also record the perimeter (to detect spikes - and still have a relatively simple check). But that's a lot of work for all cases.

Anyway, apart from the spikes, using area still seems to filter out the "most" wrong cases - but to be sure I need more figures.

[awulkiew]

@barendgehrels is it fixed now after rescaling has been removed?

[vissarion]

I am testing again this issue with current develop https://github.com/boostorg/geometry/commit/2ee0967344121c71a7e95bae84bf821e98a32cf3 (i.e. after the removal of rescaling and the merging of https://github.com/boostorg/geometry/pull/1010). I am getting the following:

For floats the issue is solved, for (long) double there is some garbage line in the bottom right corner and the multiprecision case is still wrong.

float	double	long double	m

[tinko92]

@vissarion Thanks for the test results with long double and multi-precision. I tested it with double before mentioning it in the PR and saw the line in the bottom right corner. I attributed this to precision limitations in the construction of intersection points in the result and I believed it was no longer an algorithmic issue with the turn traversal after #1010 but only a precision issue in the construction of intersection points.

Does multi-precision refer to something like boost::multiprecision::cpp_bin_float_100 ? This will still lose precision in the division, I think. In a rescaled version such that no imprecise divisions occur, e.g. https://gist.github.com/tinko92/3ee03341fa1d3ef85c11321bdd3a44be the correct result is computed as far as I can see.

With CORE::Expr from CGAL (https://doc.cgal.org/latest/Number_types/classCORE_1_1Expr.html), which provides exactness for +,-,/,,sqrt() the result is also correct, as far as I can see (MULTIPOLYGON(((25. 0,19.0000 15.0000,11.0000 15.0000,8.87097 9.67742,25. 0)),((11.0000 15.0000,15. 25.,10. 15.,11.0000 15.0000)),((19.0000 15.0000,20. 15.,15. 25.,19.0000 15.0000)),((8.87097 9.67742,0 15.,7. 10.,5.00000 0,8.87097 9.67742)))).

I thought that #1010 fixes an algorithmic issue here and just leaves an issue that requires exact constructions but maybe I was wrong and the improvement is just incidental. I will look into it again, soon when I look into exact constructions for overlay.