Garden hang on ST_VoronoiLines and ST_VoronoiPolygons

[robe]

This hangs indefinitely on my PostgreSQL 13 Geos 3.9

POSTGIS="3.1.0alpha2 3.1.0alpha1-168-gceca2a80d" [EXTENSION] PGSQL="130" GEOS="3.9.0-CAPI-1.14.0" SFCGAL="1.3.8" PROJ="6.2.1" GDAL="GDAL 3.0.2, released 2019/10/28" LIBXML="2.9.9" LIBJSON="0.12" LIBPROTOBUF="1.2.1" WAGYU="0.5.0 (Internal)" RASTER PostgreSQL 13beta2 on x86_64-w64-mingw32, compiled by gcc.exe (x86_64-posix-seh-rev0, Built by MinGW-W64 project) 8.1.0, 64-bit

SELECT ST_VoronoiLines(foo1.the_geom, 20.1, foo2.the_geom) As result, ST_AsText(foo1.the_geom) As ref1_geom, ST_AsText(foo2.the_geom) As ref2_geom FROM ((SELECT ST_Collect(s.the_geom) As the_geom FROM (SELECT ST_SetSRID(ST_Point(i,j),4326) As the_geom FROM generate_series(-10,50,15) As i CROSS JOIN generate_series(40,70, 15) j ) As s)) As foo1 CROSS JOIN ((SELECT ST_SetSRID(ST_Point(i,j),4326) As the_geom FROM (SELECT a*1.11111111 FROM generate_series(-10,50,2) As a) As i(i) CROSS JOIN generate_series(40,70, 5) j ORDER BY i,j )) As foo2 LIMIT 10;

So does this reduced one:

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)

and it's none-cancellable.

I tried it on

POSTGIS="2.5.3 r17699" [EXTENSION] PGSQL="110" GEOS="3.8.0-CAPI-1.13.1 " PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.4.4, released 2020/01/08" LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER

and doesn't hang.

and as text rep is

MULTILINESTRING((-2.5 47.5,-2.5 62.5),(-2.5 62.5,5 70),(5 70,20 55),(20 55,12.5 47.5),(12.5 47.5,-2.5 47.5),(-2.5 47.5,-70 47.5),(-70 62.5,-2.5 62.5),(-2.5 47.5,-2.5 -20),(12.5 47.5,12.5 -20),(5 130,5 70),(20 55,27.5 55),(27.5 55,27.5 -20),(27.5 130,27.5 55),(42.5 130,42.5 55),(42.5 55,27.5 55),(42.5 55,42.5 -20),(110 55,42.5 55))

I'll try the latest GEOS 3.9 head next

[Algunenano]

Can confirm it hangs in current GEOS master / HEAD.

I've attached gdb to the process and it's doing stuff:

#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
[...]

[robe]

related ticket on GEOS

​https://trac.osgeo.org/geos/ticket/1035

[robe]

No surprise I guess same issue with ST_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]

This is fixed in GEOS.