boost.geometry 多边形差异返回自相交结果或错误地检测到它

如何解决boost.geometry 多边形差异返回自相交结果或错误地检测到它

我正在尝试将 boost.geometry 用于 polygon subtraction。多边形可能是凹面，但没有 interior rings。

在大多数情况下，这很有效，但我发现至少有一种情况，boost 几何本身检测到的结果为 self intersecting。

Polygon Concept 定义了一些我认为我完全满足的规则，但至少对于峰值我不完全确定。 Boost.geometry 指的是 OGC Simple Feature Specification，它注意以下关于切割线、尖峰和穿孔的内容：

d) 多边形不得有切割线、尖刺或穿孔，例如：

∀ P ∈ 多边形，P = P.Interior.Closure;

以下示例失败：

我正在尝试从红色多边形中减去白色：

结果是看起来不错的绿色虚线多边形

近距离观察有趣的部分

在这里：

放大到标记的角落时

（虽然诚然受到查看器浮点精度的限制）

似乎有两个非常接近的点可能重叠也可能不重叠

（我没有做计算；重点是 boost.geometry 认为它们重叠）

这是a godbolt showing the behavior。

在结果多边形中形成 Z 的四条线段的坐标是

-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

虽然小数点后有一些差异，但感觉它们应该仍然足够大，不会导致浮点精度问题。

1. 是否有更详细的解释说明文档中提到的尖峰的定义 - “不应有切割线、尖峰或穿孔”
2. boost.geometry 中有什么东西吗？我不太了解的 strategies 或其他任何可以完全避免这种情况的方法？
3. 如果没有其他帮助，切换到整数会完全解决问题还是使用 boost.polygon 是唯一稳定的 boost 选项？

编辑 1：

我没有问一个可能再次归结为同一问题的类似问题，而是在这里添加另一个复制品，它没有在对交叉点的调用中显示问题，但在结果中表示应该有一个洞的地方。

以下示例失败：

我正在尝试从红色多边形中减去白色。 这应该会产生一个与红色几乎相同的多边形，并且没有任何孔。相反，结果是原始红色多边形作为外环，白色多边形作为内环。

添加和删除看似无关的点，例如红色多边形中的第 7、8 和 9 点会改变行为并使其正常工作。

据说增加更高的精度可以解决这个例子，但我怀疑这是算法固有的问题。

这是a godbolt showing the behavior。

将poly2的点向右旋转一位后，问题消失，如this godbolt所示。 covered_by 似乎与此行为相符。

解决方法

准备好：首先我确认假设（这是一个浮点精度限制/缺陷）。接下来我想出了最简单的解决方法......显然有效:)

检查前提

Firat 我简化了测试仪，添加了很多诊断：

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},-29.970448623772313},{-38.166710648232517,-29.970440761860232},{-46.094160894726727,-23.318520183850637},{-46.499997761818364,}},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,}};

    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";
}

印刷品

---
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,-38.1667 -29.9704,-46.5 -22.978)))
diff is self-intersecting
diff area: 3468.78
---
result#1: POLYGON((-46.5 -22.978,-3
8.1667 -29.9704,-46.5 -22.978))
result#1 is self-intersecting
result#1 area: 3468.78
Diff in areas: 0.0690986

从而确认前提。

在黑暗中拍摄

作为盲注，我尝试用 double 替换 long double。这导致没有明显的输出差异。

用多精度类型替换它，例如：

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

确实显示出不同：

如您所见，该区域没有显着差异，但区域 delta 发生了变化（0.0690986 变为 0.0690985），更有趣的是，“误诊”自相交已经消失。

问题

上面的分析有一个问题：不改变代码中的几个地方就不能使用

• constexpr 被错误地假定为计算类型（阻止编译）
• std::abs 是限定的，而不是从 Boost Multiprecision 调用支持 ADL 的 abs

如果您愿意，可以在此处查看相关补丁（相对于 1.76.0），但这意味着您可能会遇到更多类似的问题（就像我之前遇到的那样）：https://github.com/sehe/geometry/commit/31a7ccf1730b09b827ba6cc4dabcb845c3582a9b

commit 31a7ccf1730b09b827ba6cc4dabcb845c3582a9b
Author: sehe <[e-mail]>
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>

总结/解决方法

我不推荐使用补丁。我建议得出结论，这确实是一个精度问题。如果您认为这是库的缺陷，请考虑向库维护者报告问题。

作为一种解决方法，您可以尝试缩放输入：

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

这也消除了症状：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 -22978)))
diff area: 3.46878e+09
---
result#1: POLYGON((-46500 -22978,-46500 -22978))
result#1 area: 3.46878e+09
Diff in areas: 69098.1