提问人：jesses 提问时间：6/9/2021 最后编辑：jesses 更新时间：7/1/2021 访问量：810
boost.geometry 多边形差异返回自相交结果（或错误地检测到它）
boost.geometry polygon difference returns self intersecting result (or is wrongly detecting it as such)
问：
我正在尝试使用 boost.geometry 进行多边形减法。多边形可能是凹形的，但没有内环。
大多数情况下，这效果很好，但我发现至少有一种情况，结果被提升几何体本身检测为自相交。
多边形概念定义了一些规则，我认为我都满足了，但至少关于尖峰，我不完全确定。 Boost.geometry 指的是 OGC 简单特征规范，该规范记录了以下有关切割线、尖峰和穿刺的信息：
d） 多边形可能没有切割线、尖刺或穿孔，例如：
∀ P ∈ 多边形，P = P.Interior.Closure;
以下示例失败：
我正在尝试从红色多边形中减去白色：
结果是绿色虚线多边形，看起来很好
有趣部分的近距离观察
在这里：
放大到标记的角落时
（诚然，这在一定程度上受到查看器浮点精度的限制）
似乎有两个非常接近的点，它们可能会重叠，也可能不会重叠
（我没有做计算;重点是boost.geometry认为它们重叠）
这是一个显示该行为的 godbolt。
在结果多边形中形成 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
虽然差异在小数点后一点，但感觉它们应该仍然足够大，不会导致浮点精度问题。
- 是否有更详细的解释，说明文档中提到的尖峰的定义 - “不应该有切割线、尖峰或穿孔”
- boost.geometry 中是否有我不太了解的策略，或者其他任何可以用来完全避免这种情况的东西？
- 如果没有其他帮助，切换到整数会完全解决问题，还是使用 boost.polygon 是唯一稳定的提升选项？
编辑1：
我没有提出一个可能再次归结为相同问题的类似问题，而是在这里添加了另一个复制品，该复制品没有在对交集的调用中显示问题，而是在结果中表示应该有一个洞。
以下示例失败：
我正在尝试从红色多边形中减去白色。 这应该会产生一个与红色多边形几乎相同且没有任何孔的多边形。相反，结果是原始红色多边形作为外环，白色多边形作为内环。
添加和删除看似不相关的点，例如红色多边形中的点 7、8 和 9，会改变行为并使其正常工作。
增加更多的精度应该可以解决这个例子，但我怀疑这是算法固有的问题。
这是一个显示该行为的 godbolt。
将 poly2 的点向右旋转 1 后，问题消失了，如此 godbolt 所示。 似乎与这种行为相匹配。covered_by
答：
2赞
sehe
6/9/2021
#1
准备搭便车：首先我确认假设（这是一个浮点精度限制/缺陷）。接下来，我想出最简单的解决方法......显然:)
检查前提
Firat I 简化了测试器，添加了大量诊断功能：
#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";
}
指纹
---
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
从而证实了前提。
在黑暗中拍摄
作为盲目拍摄，我尝试用 .这导致输出没有明显的差异。doublelong double
将其替换为多精度类型，例如：
using Coord = boost::multiprecision::number<
boost::multiprecision::backends::cpp_dec_float<100>,
boost::multiprecision::et_off>;
确实显示了差异：
正如你所看到的，面积没有显着差异，但面积增量发生了变化（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 <[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>
摘要/解决方法
我不建议使用补丁。我建议得出结论，这确实是一个精度问题。如果您认为这是库的缺陷，请考虑将问题报告给库维护人员。
作为解决方法，您可以尝试缩放输入：
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 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
评论
0赞
jesses
6/30/2021
感谢您的详细回答以及您所做的工作。我真的很感激。很抱歉耽搁了。我还有很多其他事情要做，但现在可以测试更多。我尝试使用不同的比例，切换到整数，甚至尝试了 boost：:p olygon，但最后当我将此示例反馈给 intersect 方法时，它们似乎都存在问题。有些秤甚至根本无法提供结果。我还发现了一个类似的问题，我在上面更新了这个问题，其中差异导致一个显然不应该有的洞。
1赞
jesses
6/30/2021
虽然提高准确性可以解决这个例子，但我认为它并不能解决所有问题，而只是让它们不太可能出现。所以我听从了你的建议，在 github.com/boostorg/geometry/issues/875 和 github.com/boostorg/geometry/issues/876 创建了两个错误报告
评论