RGeo 交叉函数的问题

Problems with RGeo intersection functions

我从 RGeo 多边形相交函数(Ruby 2.3.0,RGeo 0.5.3)strange/incorrect 得到结果

示例 1:

我有两个多边形,我认为它们共享一个边界但不共享任何内部 space(即它们 touch but do not overlap):

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)

当我们检查它们之间的交点时,它 returns 一条线,正如仅共享边界的几何图形所预期的那样:

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 段线,两者都不是。

在第一个示例中,您为两个多边形提供了一些浮点数。是什么让您认为它们接触但不重叠?你怎么知道?在一般情况下,如果两条线段完全重叠或不重叠,则不能绝对地说(不使用公差)。您可以提供具有圆角坐标的线,在这种情况下可以为相同的参数提出论点。

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

我没有使用过 RGeo,但是有没有办法在计算交叉点几何图形时调整公差?尝试降低公差值(使该阈值更小)。这将允许函数在不合并彼此太近的点的情况下计算几何。

简答

不幸的是, 在使用 Float 坐标时,您无法期望 touches?overlaps? 得到可靠和正确的输出。

它不依赖于 RGeo 或 GEOS 版本(或者就此而言,JTS,GEOS 所基于的项目)。

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

Geometry#intersectiontouches?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?都在[=30之间切换=] 和 false.

交叉点要么是空的,要么是一条线,要么是一个多边形,但至少 intersectionoverlaps?touches? 之间的结果看起来是一致的。

在您的情况下,具有倾斜边的更复杂的多边形会使问题变得更难。在计算两个线段的交点时,很难保持 1 埃的精度!

是不是 RGeo 的问题?

RGeo uses GEOS, which is a C++ port of JTS (Java Topology Suite).

为了检查问题是否特定于 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 相同的错误结果。

为什么结果不一致?

为什么 intersectiontouches? return 结果不一致?

DE-9IM

touches?intersects?overlaps?等谓词都是从the Dimensionally Extended nine-Intersection Model (DE-9IM)派生出来的。它是一个描述几何内部、边界和外部交集维度的矩阵。

这个矩阵是在 GEOS src/operation/relate/RelateComputer.cpp 的第 72 行计算的:

IntersectionMatrix* RelateComputer::computeIM()

该算法似乎 ,但在您的任何示例中都不是这种情况。

我能找到的所有 JTS 测试都使用整数坐标,甚至是名为 "complex polygons touching" 的测试:

  # 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_3wkt_4 在您的示例中具有圆角坐标,但计算它们的差异会创建坐标不精确的多边形:x1 of diff_a-8243077.837796713.

几何#intersection

Geometry#intersection 是在 GEOS 中 src/operation/overlay/OverlayOp.cpp 的第 670 行计算的:

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

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

这个方法似乎比RelateComputer::computeIM更稳健。

使用 GeometryPrecisionReducer?

使用 GeometryPrecisionReducerPrecisionModel,可以为所有几何体指定允许点的网格。

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

它也改变了边界和它们的交叉点。根据 PrecisionModel,第一个示例的交点可以是空的、线或多边形。