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?