search

boostorg / geometry

Boost.Geometry - Generic Geometry Library | Requires C++14 since Boost 1.75
http://boost.org/libs/geometry
Boost Software License 1.0
460 stars 215 forks source link

Wrong result of sym_difference. #900

Open awulkiew opened 3 years ago

awulkiew commented 3 years ago

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:

image871

awulkiew commented 3 years ago

@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. obraz (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 obraz (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 commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

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

vissarion commented 3 years ago

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](https://clydedacruz.github.io/openstreetmap-wkt-playground/) for simplicity, the geometries even if plotted on map are cartesian
barendgehrels commented 3 years ago

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 commented 2 years ago

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

vissarion commented 2 years ago

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:

float double long double m
With rescaling
Without rescaling
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.
tinko92 commented 2 years ago

@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.

  • © Githubissues.
  • Githubissues is a development platform for aggregating issues.
float double long double m