问题描述
我正在尝试将 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
虽然小数点后有一些差异,但感觉它们应该仍然足够大,不会导致浮点精度问题。
- 是否有更详细的解释说明文档中提到的尖峰的定义 - “不应有切割线、尖峰或穿孔”
- boost.geometry 中有什么东西吗?我不太了解的 strategies 或其他任何可以完全避免这种情况的方法?
- 如果没有其他帮助,切换到整数会完全解决问题还是使用 boost.polygon 是唯一稳定的 boost 选项?
编辑 1:
我没有问一个可能再次归结为同一问题的类似问题,而是在这里添加另一个复制品,它没有在对交叉点的调用中显示问题,但在结果中表示应该有一个洞的地方。
以下示例失败:
我正在尝试从红色多边形中减去白色。 这应该会产生一个与红色几乎相同的多边形,并且没有任何孔。相反,结果是原始红色多边形作为外环,白色多边形作为内环。
添加和删除看似无关的点,例如红色多边形中的第 7、8 和 9 点会改变行为并使其正常工作。
据说增加更高的精度可以解决这个例子,但我怀疑这是算法固有的问题。
这是a godbolt showing the behavior。
将poly2的点向右旋转一位后,问题消失,如this godbolt所示。
covered_by
似乎与此行为相符。
解决方法
准备好:首先我确认假设(这是一个浮点精度限制/缺陷)。接下来我想出了最简单的解决方法......显然有效:)
检查前提
Firat 我简化了测试仪,添加了很多诊断:
#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 <sgheeren@gmail.com>
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