boost.geometry polygon difference returns self intersecting result (or is wrongly detecting it as such)

Question

I'm trying to use boost.geometry for polygon subtraction. The polygons might be concave but don't have interior rings.

Most of the time this works well but i found at least one occasion where the result is detected by boost geometry itself as self intersecting.

The Polygon Concept defines some rules which i think i'm all fulfilling, but at least about the spikes i'm not entirely sure. Boost.geometry refers to the OGC Simple Feature Specification which notes the following about cut lines, spikes and punctures:

d) A Polygon may not have cut lines, spikes or punctures e.g.:

∀ P ∈ Polygon, P = P.Interior.Closure;

The following example fails:

I'm trying to subtract the white from the red polygon:

The result is the green dashed polygon which looks fine

A closer view of the interesting part

is here:

When zooming in to the marked corner

(while admittedly being somewhat limited by the floating point accuracy of the viewer)

it appears that there are two very close points which might or might not overlap

(i didn't do the calculations; the main point is that boost.geometry thinks they overlap)

Here's a godbolt showing the behavior.

The coordinates of the four line segments forming the Z in the result polygon are

-38.166710648232516689, -29.97044076186023176
-46.093710179049615761, -23.318898379209066718 // these two points are notably
-46.093707539928615802, -23.318900593694593226 // but not awfully close
-46.499997777166534263, -22.977982153068655435

While the differences are somewhat after the decimal point, it feels like they should be still more than big enough to not cause floating point accuracy problems.

1. Is there a more detailed explanation of what defines a spike as mentioned in the docs - "There should be no cut lines, spikes or punctures"
2. Is there something in boost.geometry, f.e. the strategies which i don't know much about, or anything else which can be used to completely avoid such situations?
3. If nothing else helps, would switching to integers completely solve the problem or is using boost.polygon the only stable boost option?

EDIT 1:

Instead of asking a similar question which likely boils down to the same problem again, i add another reproduction here which doesn't show the problem in the call to intersection but in a result which expresses a hole where there should be one.

The following example fails:

I'm trying to subtract the white from the red polygon. This should result in a polygon almost identical with the red one and without any holes. Instead the result is the original red polygon as outer ring and the white polygon as inner ring.

Adding and removing seemingly unrelated points, f.e. points 7, 8 and 9 from the red polygon, changes the behavior and makes it work correctly.

Adding more precision supposedly could solve this example but i suspect it to be a problem inherent to the algorithm.

Here's a godbolt showing the behavior.

After rotating the points of poly2 by one to the right, the problem disappears as shown in this godbolt. covered_by seems to match with this behavior.

Answer 1

Brace for a ride: first I confirm the hypotheses (it's a floating point accuracy limitation/defect). Next I come up with the simplest of workarounds that... apparently works :)

Checking The Premise

Firat I simplified the tester adding a lot of diagnostics:

Live On Wandbox

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/geometries.hpp>
#include <iostream>
#include <iomanip>
namespace bg = boost::geometry;

void report(auto& geo, std::string_view caption) {
    std::cout << "---\n";
    std::cout << caption << ": " << bg::wkt(geo) << "\n";

    std::string reason;
    if (!bg::is_valid(geo, reason)) {
        std::cout << caption << ": " << reason << "\n";

        bg::correct(geo);
        std::cout << caption << " corrected: " << bg::wkt(geo) << "\n";
    }

    if (bg::intersects(geo)) {
        std::cout << caption << " is self-intersecting\n";
    }
    std::cout << caption << " area: " << bg::area(geo) << "\n";
}

int main() {
    using point_t         = bg::model::d2::point_xy<double>;
    using polygon_t       = bg::model::polygon<point_t /*, true, true*/>;
    using multi_polygon_t = bg::model::multi_polygon<polygon_t>;

    polygon_t poly1{{
        {-46.499997761818364, -23.318506263117456},
        {-46.499999998470159, 26.305250946791375},
        {-5.3405104310993323, 15.276598956337693},
        {37.500000001521741, -9.4573812741570009},
        {37.500000001521741, -29.970448623772313},
        {-38.166710648232517, -29.970440761860232},
        {-46.094160894726727, -23.318520183850637},
        {-46.499997761818364, -23.318506263117456},
    }},
        poly2{{
            {-67.554314795325794, -23.318900735763236},
            {-62.596713294359084, -17.333596950467950},
            {-60.775620215083222, -15.852879652420938},
            {-58.530163386780792, -15.186307709861694},
            {-56.202193256444019, -15.435360658555282},
            {-54.146122173314907, -16.562122444733632},
            {-46.093707539928616, -23.318900593694593},
            {-67.554314795325794, -23.318900735763236},
        }};

    report(poly1, "poly1");
    report(poly2, "poly2");

    multi_polygon_t diff;
    bg::difference(poly1, poly2, diff);

    if (diff.empty()) {
        std::cout << "difference is empty\n";
    } else {
        report(diff, "diff");

        for (size_t i = 0; i < diff.size(); ++i) {
            report(diff[i], "result#" + std::to_string(i+1));
        }
    }

    std::cout << "Diff in areas: " << (bg::area(poly1) - bg::area(diff)) << "\n";
}

Prints

---
poly1: POLYGON((-46.5 -23.3185,-46.5 26.3053,-5.34051 15.2766,37.5 -9.45738,37.5 -29.9704,-38.
1667 -29.9704,-46.0942 -23.3185,-46.5 -23.3185))
poly1 area: 3468.84
---
poly2: POLYGON((-67.5543 -23.3189,-62.5967 -17.3336,-60.7756 -15.8529,-58.5302 -15.1863,-56.20
22 -15.4354,-54.1461 -16.5621,-46.0937 -23.3189,-67.5543 -23.3189))
poly2 area: 105.495
---
diff: MULTIPOLYGON(((-46.5 -22.978,-46.5 26.3053,-5.34051 15.2766,37.5 -9.45738,37.5 -29.9704,
-38.1667 -29.9704,-46.0937 -23.3189,-46.0937 -23.3189,-46.5 -22.978)))
diff is self-intersecting
diff area: 3468.78
---
result#1: POLYGON((-46.5 -22.978,-46.5 26.3053,-5.34051 15.2766,37.5 -9.45738,37.5 -29.9704,-3
8.1667 -29.9704,-46.0937 -23.3189,-46.0937 -23.3189,-46.5 -22.978))
result#1 is self-intersecting
result#1 area: 3468.78
Diff in areas: 0.0690986

Thus confirming the premise.

Shot In The Dark

As a blind shot I tried replacing double with long double. This leads to no discernible difference in output.

Replacing it with a multi-precision type like:

using Coord = boost::multiprecision::number<
    boost::multiprecision::backends::cpp_dec_float<100>,
    boost::multiprecision::et_off>;

does show a difference:

As you can see there is no significant difference in the area, but the area delta changed (0.0690986 became 0.0690985) and more interestingly, the "mis-diagnosed" self-intersection has vanished.

THE PROBLEM

There's a problem with the above analysis: you can't use it without changing a few spots in the code where

• constexpr is wrongly assumed for the calculation type (preventing compilation)
• std::abs is qualified, instead of invoking ADL-enabled abs from Boost Multiprecision

If you want you can see the relevant patch (relative to 1.76.0) here, but it will by implication mean you may run into more similar issues (like I have before): https://github.com/sehe/geometry/commit/31a7ccf1730b09b827ba6cc4dabcb845c3582a9b

commit 31a7ccf1730b09b827ba6cc4dabcb845c3582a9b
Author: sehe <[email protected]>
Date:   Wed Jun 9 16:35:53 2021 +0200

    Minimal patch for https://stackoverflow.com/q/67904576/85371
    
    Allows compilation of bg::difference (specifically, sort_by_side) using
    Multiprecision number type. Expressions templates hav already been
    disabled to sidestep the bulk of the issues.

diff --git a/include/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp b/include/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
index f65c8ebae..72f4aa724 100644
--- a/include/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
+++ b/include/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
@@ -299,7 +299,7 @@ public :
         // Including distance would introduce cyclic dependencies.
         using coor_t = typename select_coordinate_type<Point1, Point2>::type;
         using calc_t = typename geometry::select_most_precise <coor_t, T>::type;
-        constexpr calc_t machine_giga_epsilon = 1.0e9 * std::numeric_limits<calc_t>::epsilon();
+        calc_t machine_giga_epsilon = 1.0e9 * std::numeric_limits<calc_t>::epsilon();
 
         calc_t const& a0 = geometry::get<0>(a);
         calc_t const& b0 = geometry::get<0>(b);
@@ -310,9 +310,10 @@ public :
 
         // The maximum limit is avoid, for floating point, large limits like 400
         // (which are be calculated using eps)
-        constexpr calc_t maxlimit = 1.0e-3;
+        calc_t maxlimit = 1.0e-3;
         auto const limit = (std::min)(maxlimit, limit_giga_epsilon * machine_giga_epsilon * c);
-        return std::abs(a0 - b0) <= limit && std::abs(a1 - b1) <= limit;
+        using std::abs;
+        return abs(a0 - b0) <= limit && abs(a1 - b1) <= limit;
     }
 
     template <typename Operation, typename Geometry1, typename Geometry2>

SUMMARY / WORKAROUND

I donot recommend using the patch. I recommend concluding that indeed it's a precision problem. If you think that is a defect with the library, consider reporting the problem to the library maintainers.

As a workaround, you can try scaling your inputs:

for (auto& pt: poly1.outer()) bg::multiply_value(pt, 1'000);
for (auto& pt: poly2.outer()) bg::multiply_value(pt, 1'000);

This also removes the symptom: Live On Wandbox

---
poly1: POLYGON((-46500 -23318.5,-46500 26305.3,-5340.51 15276.6,37500 -9457.38,37500 -29970.4,-38166.7 -29970.4,-46094.2 -23318.5,-46500 -23318.5))
poly1 area: 3.46884e+09
---
poly2: POLYGON((-67554.3 -23318.9,-62596.7 -17333.6,-60775.6 -15852.9,-58530.2 -15186.3,-56202.2 -15435.4,-54146.1 -16562.1,-46093.7 -23318.9,-67554.3 -23318.9))
poly2 area: 1.05495e+08
---
diff: MULTIPOLYGON(((-46500 -22978,-46500 26305.3,-5340.51 15276.6,37500 -9457.38,37500 -29970.4,-38166.7 -29970.4,-46093.7 -23318.9,-46500 -22978)))
diff area: 3.46878e+09
---
result#1: POLYGON((-46500 -22978,-46500 26305.3,-5340.51 15276.6,37500 -9457.38,37500 -29970.4,-38166.7 -29970.4,-46093.7 -23318.9,-46500 -22978))
result#1 area: 3.46878e+09
Diff in areas: 69098.1