VoronoiLines and VoronoiPolygons hangs indefinitely on multipoint

[robe]

I tested on 3.8.0 and works fine there (still have to test on 3.8.1) but definitely an issue on GEOS head.

SELECT ST_VoronoiLines('MULTIPOINT(-10 40,5 40,20 40,35 40,50 40,-10 55,5 55,20 55,35 55,50 55,-10 70,5 70,20 70,35 70,50 70)'::geometry, 20.1, 'POINT(-11.1111111 40)'::geometry)

See related PostGIS ticket

​https://trac.osgeo.org/postgis/ticket/4726

Raul's trace

#0  0x00007fc3eacf126c in geos::triangulate::quadedge::Vertex::isCCW (this=<optimized out>, b=..., c=...) at ../../../include/geos/triangulate/quadedge/Vertex.h:226
#1  geos::triangulate::quadedge::Vertex::rightOf (this=0x55bba48b9720, e=...) at Vertex.cpp:83
#2  0x00007fc3eace807c in geos::triangulate::IncrementalDelaunayTriangulator::insertSite (this=<optimized out>, v=...) at IncrementalDelaunayTriangulator.cpp:90
#3  0x00007fc3eace7e4b in geos::triangulate::IncrementalDelaunayTriangulator::insertSites (this=0x7ffe3cf96408, vertices=...) at IncrementalDelaunayTriangulator.cpp:40
#4  0x00007fc3eace9f5d in geos::triangulate::VoronoiDiagramBuilder::create (this=0x7ffe3cf964b8) at VoronoiDiagramBuilder.cpp:93
#5  0x00007fc3eacea45b in geos::triangulate::VoronoiDiagramBuilder::getDiagramEdges (this=0x7ffe3cf964b8, geomFact=...) at VoronoiDiagramBuilder.cpp:122
#6  0x00007fc3eba9c83e in GEOSVoronoiDiagram_r::$_188::operator() (this=<optimized out>) at geos_ts_c.cpp:3261
#7  execute<GEOSVoronoiDiagram_r::$_188, (decltype(nullptr))0>(GEOSContextHandle_HS*, GEOSVoronoiDiagram_r::$_188&&) (extHandle=0x55bba48b88b0, f=...) at geos_ts_c.cpp:379
#8  GEOSVoronoiDiagram_r (extHandle=0x55bba48b88b0, g1=0x55bba48b6840, env=0x55bba48e6f40, tolerance=20.100000000000001, onlyEdges=1) at geos_ts_c.cpp:3253
#9  0x00007fc3ebc02ecb in lwgeom_voronoi_diagram (g=<optimized out>, env=0x7ffe3cf965c0, tolerance=20.100000000000001, output_edges=-1534355264) at lwgeom_geos.c:1926
#10 0x00007fc3ebb3267e in ST_Voronoi (fcinfo=0x55bba48fbb20) at lwgeom_geos.c:3465
#11 0x000055bba3c126e1 in ExecInterpExpr (state=<optimized out>, econtext=<optimized out>, isnull=0x7ffe3cf9675f) at execExprInterp.c:625
#12 0x000055bba3c45874 in ExecEvalExprSwitchContext (state=0x55bba48fba48, econtext=0x55bba48fb770, isNull=0x7ffe3cf9675f) at ../../../src/include/executor/executor.h:307
[...]

(gdb) c
Continuing.
^C
Program received signal SIGINT, Interrupt.
geos::triangulate::quadedge::Vertex::isCCW (this=<optimized out>, b=..., c=...) at ../../../include/geos/triangulate/quadedge/Vertex.h:227
227                    > (b.p.y - p.y) * (c.p.x - p.x);
(gdb) bt
#0  geos::triangulate::quadedge::Vertex::isCCW (this=<optimized out>, b=..., c=...) at ../../../include/geos/triangulate/quadedge/Vertex.h:227
#1  geos::triangulate::quadedge::Vertex::rightOf (this=0x55bba48b9ee0, e=...) at Vertex.cpp:83
#2  0x00007fc3eace807c in geos::triangulate::IncrementalDelaunayTriangulator::insertSite (this=<optimized out>, v=...) at IncrementalDelaunayTriangulator.cpp:90
#3  0x00007fc3eace7e4b in geos::triangulate::IncrementalDelaunayTriangulator::insertSites (this=0x7ffe3cf96408, vertices=...) at IncrementalDelaunayTriangulator.cpp:40
#4  0x00007fc3eace9f5d in geos::triangulate::VoronoiDiagramBuilder::create (this=0x7ffe3cf964b8) at VoronoiDiagramBuilder.cpp:93
#5  0x00007fc3eacea45b in geos::triangulate::VoronoiDiagramBuilder::getDiagramEdges (this=0x7ffe3cf964b8, geomFact=...) at VoronoiDiagramBuilder.cpp:122
#6  0x00007fc3eba9c83e in GEOSVoronoiDiagram_r::$_188::operator() (this=<optimized out>) at geos_ts_c.cpp:3261
#7  execute<GEOSVoronoiDiagram_r::$_188, (decltype(nullptr))0>(GEOSContextHandle_HS*, GEOSVoronoiDiagram_r::$_188&&) (extHandle=0x55bba48b88b0, f=...) at geos_ts_c.cpp:379
#8  GEOSVoronoiDiagram_r (extHandle=0x55bba48b88b0, g1=0x55bba48b6840, env=0x55bba48e6f40, tolerance=20.100000000000001, onlyEdges=1) at geos_ts_c.cpp:3253
#9  0x00007fc3ebc02ecb in lwgeom_voronoi_diagram (g=<optimized out>, env=0x7ffe3cf965c0, tolerance=20.100000000000001, output_edges=-1534355264) at lwgeom_geos.c:1926
#10 0x00007fc3ebb3267e in ST_Voronoi (fcinfo=0x55bba48fbb20) at lwgeom_geos.c:3465
[...]

Same with VoronoiPolygons

SELECT ST_VoronoiPolygons('MULTIPOINT(-10 40,5 40,20 40,35 40,50 40,-10 55,5 55,20 55,35 55,50 55,-10 70,5 70,20 70,35 70,50 70)'::geometry, 20.1, 'POINT(-11.1111111 40)'::geometry)

[mdavis]

See also #976 and #859.. Those issues are a different bug.

This case works fine in JTS. It seems that the fix that worked in JTS (computing triangle circumcentres) does not work in GEOS?

[mdavis]

The tolerance value of 20.1 is pretty massive compared to the scale of the input.

It looks like a non-zero tolerance value is the cause of the infinite loop (even 0.1 causes it). Not sure why that is at the moment.

[mdavis]

This case also creates an infinite loop when evaluating a Delaunay Triangulation.

[mdavis]

A slightly smaller reproducer is:

MULTIPOINT ((-10 40), (5 40), (20 40), (-10 55), (5 55), (20 55), (-10 70), (5 70), (20 70))

[mdavis]

After bisecting it appears (not surprisingly) that the cause of this bug is the commit:

​https://github.com/libgeos/geos/commit/8b433f0887633e320053e169e5bf71c5bbaecb1f

Not sure why, however. I suspect that there's an issue in the QuadEdgeSubdivision::remove method, since that is only called when there is a non-zero tolerance being used. Further evidence: if the removal code is disabled, then the case works.

[mdavis]

Found the fix for this. See ​https://github.com/libgeos/geos/pull/342

The problem is caused by two bugs in the QuadEdgeSubdivision::remove method:

• Since the method argument can be an edge in either base or sym orientation, the code to remove the corresponding QuadEdgeQuartet needs to check if the quartet base matches the arg edge OR its sym
• The lambda capture clause must specify the edge parameter as e, not &e, since it is already a reference

[mdavis]

Fixed by 1198e06/git