ruby - RGeo 交集函数的问题

Question

我从 RGeo 多边形相交函数（Ruby 2.3.0，RGeo 0.5.3）得到奇怪/不正确的结果

示例 1：

我有两个多边形，我相信它们共享一个边界但不共享任何内部空间（即它们接触但不重叠）：

wkt_1 = "POLYGON ((-8226874.27782158 4962626.76394919, -8223358.174520462 4961756.817075645, -8223358.174520462 4960289.557693501, -8224471.369428394 4960289.557693501, -8226874.27782158 4962253.674727506, -8226874.27782158 4962626.76394919))"
wkt_2 = "POLYGON ((-8224757.546680832 4960523.476563589, -8225269.1002275925 4959296.105368667, -8226993.791361805 4959219.668340384, -8226420.900079966 4961883.087589158, -8224757.546680832 4960523.476563589))"
poly_1 = RGeo::Geos.factory(:srid => 3857).parse_wkt(wkt_1)
poly_2 = RGeo::Geos.factory(:srid => 3857).parse_wkt(wkt_2)

当我们检查它们之间的交点时，它会返回一条线，正如预期的仅共享边界的几何图形：

poly_1.intersection poly_2
=> #<RGeo::Geos::CAPILineStringImpl:0x3fc0249af168 "LINESTRING (-8224757.546680832 4960523.476563589, -8225598.074380083 4961210.51680879)">

但是，在运行以下检查时，我们得到的结果与预期相反：

poly_1.overlaps? poly_2
=> true
poly_1.touches? poly_2
=> false

示例 2：

我们采用两个合法重叠的多边形：

wkt_3 = "POLYGON ((-8243237.0 4970203.0, -8243237.0 4968735.0, -8242123.0 4968735.0, -8242123.0 4970203.0, -8243237.0 4970203.0))"
wkt_4 = "POLYGON ((-8244765.0 4966076.0, -8244765.0 4964608.0, -8243652.0 4964608.0, -8243652.0 4966076.0, -8242680.0 4969362.0, -8244765.0 4966076.0))"
poly_3 = RGeo::Geos.factory(:srid => 3857).parse_wkt(wkt_3)
poly_4 = RGeo::Geos.factory(:srid => 3857).parse_wkt(wkt_4)

并计算两个方向的差异：

diff_a = poly_3.difference poly_4
=> #<RGeo::Geos::CAPIPolygonImpl:0x3fe3fca26028 "POLYGON ((-8243077.837796713 4968735.0, -8243237.0 4968735.0, -8243237.0 4970203.0, -8242123.0 4970203.0, -8242123.0 4968735.0, -8242865.466828971 4968735.0, -8242680.0 4969362.0, -8243077.837796713 4968735.0))">
diff_b = poly_4.difference poly_3
=> #<RGeo::Geos::CAPIPolygonImpl:0x3fe3fd1dda28 "POLYGON ((-8242865.466828971 4968735.0, -8243652.0 4966076.0, -8243652.0 4964608.0, -8244765.0 4964608.0, -8244765.0 4966076.0, -8243077.837796713 4968735.0, -8242865.466828971 4968735.0))">

现在我们根据减法器检查剩余的多边形：

diff_b.touches? poly_3
=> true
diff_b.overlaps? poly_3
=> false

这很好。

diff_a.touches? poly_4
=> false
diff_a.overlaps? poly_4
=> true

这是不可能的-从另一个多边形中减去的多边形的其余部分不可能与该多边形重叠！

为什么它只发生在一个方向？

示例 3：

怪事还在继续。现在让我们得到 poly_3 和 poly_4 的交集

intersection_a = poly_3.intersection poly_4
=> #<RGeo::Geos::CAPIPolygonImpl:0x3fd1084ece88 "POLYGON ((-8242865.724520766 4968734.582337743, -8243078.32501591 4968734.582337743, -8242680.062418439 4969362.301390027, -8242865.724520766 4968734.582337743))">

现在，既然这是应该从 poly_3 中减去得到 diff_a 的值，因此 diff_a 应该与 intersection_a 相交的方式与它与 poly_4（减法器）相交的方式完全相同。

除非它没有：

diff_a.touches? poly_4
=> false
diff_a.touches? intersection_a
=> true
diff_a.intersection poly_4
=> #<RGeo::Geos::CAPILineStringImpl:0x3fe3f98fb854 "LINESTRING (-8242680.0 4969362.0, -8243077.837796713 4968735.0)">
diff_a.intersection intersection_a
=> #<RGeo::Geos::CAPIMultiLineStringImpl:0x3fe3fca157b4 "MULTILINESTRING ((-8242865.466828971 4968735.0, -8242680.0 4969362.0), (-8242680.0 4969362.0, -8243077.837796713 4968735.0))">

更糟糕的是，这两个交集结果都没有意义。它应该是一条单独的 2 段线，但两者都不是。

Answer 1

简短的回答

不幸的是，在使用浮点坐标时，您似乎不能期望得到可靠和正确的输出touches?。overlaps?

它不依赖于 RGeo 或 GEOS 版本（或者就此而言，JTS是 GEOS 所基于的项目）。

如果您需要有关两个多边形位置的信息，您可以使用Geometry#distance并检查它是否小于给定的 epsilon。

Geometry#intersection比touches?and更可靠一点overlaps?，但也不能保证适用于每个示例。

---

是精度问题吗？

理论

touches?根据定义非常敏感：多边形必须在其边界上共享一个点或一条线，但不应有任何共同的内部点。

• 如果它们彼此相距太远，则它们没有任何共同点，并且touches?是错误的。
• 如果它们太靠近，它们的交点是一个多边形，并且touches?是错误的。

敏感性分析

它到底有多敏感？

让我们考虑两个多边形，一个在另一个之上：

POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))
POLYGON((0 1, 0 2, 1 2, 1 1, 0 1))

我们可以将上多边形的下边界移动一个 epsilon，看看它有什么影响：

require 'rgeo'

epsilon = 1E-15

deltas = [-epsilon, 0, epsilon]

deltas.each do |delta|
  puts "--------------------------------"
  puts "Delta : #{delta}\n\n"
  simple1 = 'POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))'
  simple2 = "POLYGON((0 #{1+delta}, 0 2, 1 2, 1 #{1+delta}, 0 #{1+delta}))"

  puts "Polygon #1 : #{simple1}\n"
  puts "Polygon #2 : #{simple2}\n\n"

  poly_1 = RGeo::Geos.factory(:srid => 4326).parse_wkt(simple1)
  poly_2 = RGeo::Geos.factory(:srid => 4326).parse_wkt(simple2)

  puts "Intersection : #{poly_1.intersection poly_2}"
  puts

  puts "Distance between polygons :"
  puts format('%.30f°', poly_1.distance(poly_2))
  puts

  puts "Overlaps? : #{poly_1.overlaps? poly_2}"
  puts "Touches? : #{poly_1.touches? poly_2}"
end

它输出：

--------------------------------
Delta : -1.0e-15

Polygon #1 : POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))
Polygon #2 : POLYGON((0 0.999999999999999, 0 2, 1 2, 1 0.999999999999999, 0 0.999999999999999))

Intersection : POLYGON ((0.0 0.999999999999999, 0.0 1.0, 1.0 1.0, 1.0 0.999999999999999, 0.0 0.999999999999999))

Distance between polygons :
0.000000000000000000000000000000°

Overlaps? : true
Touches? : false
--------------------------------
Delta : 0

Polygon #1 : POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))
Polygon #2 : POLYGON((0 1, 0 2, 1 2, 1 1, 0 1))

Intersection : LINESTRING (0.0 1.0, 1.0 1.0)

Distance between polygons :
0.000000000000000000000000000000°

Overlaps? : false
Touches? : true
--------------------------------
Delta : 1.0e-15

Polygon #1 : POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))
Polygon #2 : POLYGON((0 1.000000000000001, 0 2, 1 2, 1 1.000000000000001, 0 1.000000000000001))

Intersection : GEOMETRYCOLLECTION EMPTY

Distance between polygons :
0.000000000000001110223024625157°

Overlaps? : false
Touches? : false

这些是表现良好的多边形，具有

• 平行于轴线的边
• 相同大小
• 整个长度的共享边

一个度数的差异1E-15（大约1 Ångström）足以完全改变结果，两者overlaps?和之间touches?切换。truefalse

交点要么是空的，要么是一条线，要么是一个多边形，但至少结果在 和 之间看起来是一致intersection的。overlaps?touches?

在您的情况下，具有倾斜边的更复杂的多边形会使问题变得更加困难。在计算两条线段的交点时，很难保持 1-Ångström 的精度！

RGeo 有问题吗？

RGeo使用GEOS，它是JTS (Java Topology Suite)的 C++ 端口。

为了检查问题不是特定于 RGeo 或 GEOS，我使用 JTS 1.14 计算了示例 1：

WKTReader wktReader = new WKTReader();
String wkt1 = "POLYGON ((-8226874.27782158 4962626.76394919, -8223358.174520462 4961756.817075645, -8223358.174520462 4960289.557693501, -8224471.369428394 4960289.557693501, -8226874.27782158 4962253.674727506, -8226874.27782158 4962626.76394919))";
String wkt2 = "POLYGON ((-8224757.546680832 4960523.476563589, -8225269.1002275925 4959296.105368667, -8226993.791361805 4959219.668340384, -8226420.900079966 4961883.087589158, -8224757.546680832 4960523.476563589))";
Polygon poly1 = (Polygon) wktReader.read(wkt1);
Polygon poly2 = (Polygon) wktReader.read(wkt2);
System.out.println("Intersection : " + poly1.intersection(poly2));
System.out.println("Overlaps?    : " + poly1.overlaps(poly2));
System.out.println("Intersects?  : " + poly1.intersects(poly2));
System.out.println("Touches?     : " + poly1.touches(poly2));
showMatrixWith(poly1, poly2);

它输出：

Intersection : LINESTRING (-8224757.546680832 4960523.476563589, -8225598.074380083 4961210.51680879)
Overlaps?    : true
Intersects?  : true
Touches?     : false

212
101
212

交点与您的示例完全相同，intersects?并touches? 输出与 RGeo 相同的错误结果。

为什么结果不一致？

为什么intersection和touches?返回不一致的结果？

DE-9IM

touches?, intersects?,overlaps?和其他谓词派生自维度扩展的九交集模型 (DE-9IM)。它是一个矩阵，描述几何内部、边界和外部之间的交集维度。

该矩阵src/operation/relate/RelateComputer.cpp在 GEOS 的第 72 行计算：

IntersectionMatrix* RelateComputer::computeIM()

该算法似乎需要精确的 noding，在您的任何示例中都不是这种情况。

我可以为 JTS 找到的所有测试都使用整数坐标，甚至是一个名为“复杂多边形接触”的测试：

  # line 477 in jts-1.14/testxml/general/TestFunctionAA.xml
  <desc>mAmA - complex polygons touching</desc>
  <a>
    MULTIPOLYGON(
      (
        (100 200, 100 180, 120 180, 120 200, 100 200)),
      (
        (60 240, 60 140, 220 140, 220 160, 160 160, 160 180, 200 180, 200 200, 160 200,
        160 220, 220 220, 220 240, 60 240),
        (80 220, 80 160, 140 160, 140 220, 80 220)),
      (
        (280 220, 240 180, 260 160, 300 200, 280 220)))
  </a>
  <b>
    MULTIPOLYGON(
      (
        (80 220, 80 160, 140 160, 140 220, 80 220),
        (100 200, 100 180, 120 180, 120 200, 100 200)),
      (
        (220 240, 220 220, 160 220, 160 200, 220 200, 220 180, 160 180, 160 160, 220 160,
        220 140, 320 140, 320 240, 220 240),
        (240 220, 240 160, 300 160, 300 220, 240 220)))
  </b>

没有一个测试用例为浮点坐标计算谓词。

请注意，即使wkt_3并且wkt_4在您的示例中具有舍入坐标，计算它们的差异也会创建具有不精确坐标的多边形：x1of diff_ais -8243077.837796713。

几何#intersection

Geometry#intersectionsrc/operation/overlay/OverlayOp.cpp在 GEOS的第 670 行计算：

void OverlayOp::computeOverlay(OverlayOp::OpCode opCode)

代码中的注释似乎表明不需要精确的点头，并且有多个 if 语句来检查结果是否正确。

这种方法似乎比RelateComputer::computeIM.

使用 GeometryPrecisionReducer？

使用GeometryPrecisionReducer和 a PrecisionModel，可以为所有几何体指定允许点的网格。

GeometryPrecisionReducer在 GEOS 中实现，但在 RGeo 中不可用。在 JTS 中使用您的第一个示例进行的测试表明，它无论如何都不能解决您的问题：不精确的坐标被捕捉到网格上最近的点。每个点都向北/南/东或西移动一点，这会改变每一侧的坡度。

它还改变了边界及其交叉点。根据PrecisionModel，您的第一个示例的交点可能是空的、线或多边形。

Answer 2

在第一个示例中，您为两个多边形提供了一些浮点数。是什么让你认为它们接触而不重叠？你怎么知道？在一般情况下，您不能用绝对术语（不使用公差）说两条线段是否完全重叠。您可以提供具有圆角坐标的线，在这种情况下，可以为其设置一个参数。

浮点数的使用是实数空间离散化的一种形式。当使用不同的算法或方法来计算相同的结果时，这会导致结果不一致。例如，考虑 2D 空间中三条线的交点，它们应该在同一点相交。现在，每次使用不同的线对独立计算交点 3 次。您每次都可能得到相似但不平等的答案。

我没有使用 RGeo，但是在计算相交几何形状时有没有办法调整公差？尝试降低容差值（使该阈值更小）。这将允许函数在不合并彼此太近的点的情况下计算几何。