Overlapping Voronoi output from four input points

[mappingvermont]

👋 hi geos folks!

Thanks for writing an amazing library. I don't know where we'd be without it 🙏

I'm using this to build voronoi polygons in PostGIS, and found an issue when using four very specific points as an input:

create table _tmp1 as
select
(ST_Dump(ST_VoronoiPolygons(ST_Collect(ST_SetSRID(geom, 4326))))).geom  geom
from (
  values
    (ST_GeomFromWKB(E'\\x0101000000e0be0e9c33a21a40b5c876be9fca4a40')),
    (ST_GeomFromWKB(E'\\x0101000000fe65f7e461a11a400af9a067b3ca4a40')),
    (ST_GeomFromWKB(E'\\x010100000054e3a59bc4a01a40b459f5b9daca4a40')),
    (ST_GeomFromWKB(E'\\x01010000008db96b09f9a01a405f29cb10c7ca4a40'))
) s(geom);

Here's the output in QGIS, showing the overlapping polygons along with the gap over one of the input points:

I asked about this on the spatial community slack and @dbaston graciously took my example and was able to reproduce the issue using a geos CLI tool:

echo '0101000000e0be0e9c33a21a40b5c876be9fca4a40\n0101000000fe65f7e461a11a400af9a067b3ca4a40\n010100000054e3a59bc4a01a40b459f5b9daca4a40\n01010000008db96b09f9a01a405f29cb10c7ca4a40' > /tmp/example.wkb
geosop -c -a /tmp/example.wkb -f wkt voronoi
POLYGON ((6.655200000000004 53.5866, 6.660199999999996 53.5866, 6.660199999999996 53.58584444444445, 6.659499999999997 53.585300000000004, 6.655200000000004 53.583866666666665, 6.655200000000004 53.5866))
POLYGON ((6.655200000000004 53.5866, 6.660199999999996 53.5866, 6.660199999999996 53.58120000000001, 6.655200000000004 53.58120000000001, 6.655200000000004 53.582433333333334, 6.659499999999997 53.585300000000004, 6.655200000000004 53.583866666666665, 6.655200000000004 53.5866))
POLYGON ((6.660199999999996 53.58120000000001, 6.656425000000008 53.58120000000001, 6.659499999999996 53.585300000000004, 6.659499999999997 53.585300000000004, 6.660199999999996 53.58584444444445, 6.660199999999996 53.58120000000001))
POLYGON ((6.655200000000004 53.58120000000001, 6.655200000000004 53.582433333333334, 6.659499999999996 53.585300000000004, 6.656425000000008 53.58120000000001, 6.655200000000004 53.58120000000001))

I'm not sure where to go from here, but happy to help test however I can. Thank you!

[mwtoews]

Seems to reproduce with current main branch.

This is the unexpected result in the issue:

$ ./bin/geosop -a "MULTIPOINT ((6.6584 53.583000000000006), (6.6576 53.583600000000004), (6.657 53.5848), (6.6572000000000005 53.5842))" -f wkt voronoi
POLYGON ((6.655200000000004 53.5866, 6.660199999999996 53.5866, 6.660199999999996 53.58584444444445, 6.659499999999997 53.585300000000004, 6.655200000000004 53.583866666666665, 6.655200000000004 53.5866))
POLYGON ((6.655200000000004 53.5866, 6.660199999999996 53.5866, 6.660199999999996 53.58120000000001, 6.655200000000004 53.58120000000001, 6.655200000000004 53.582433333333334, 6.659499999999997 53.585300000000004, 6.655200000000004 53.583866666666665, 6.655200000000004 53.5866))
POLYGON ((6.660199999999996 53.58120000000001, 6.656425000000008 53.58120000000001, 6.659499999999996 53.585300000000004, 6.659499999999997 53.585300000000004, 6.660199999999996 53.58584444444445, 6.660199999999996 53.58120000000001))
POLYGON ((6.655200000000004 53.58120000000001, 6.655200000000004 53.582433333333334, 6.659499999999996 53.585300000000004, 6.656425000000008 53.58120000000001, 6.655200000000004 53.58120000000001))

Chopping a few decimal places from the input, here is an expected result from a very similar example:

$ ./bin/geosop -a "MULTIPOINT ((6.6584 53.583), (6.6576 53.5836), (6.657 53.5848), (6.6572 53.5842))" -f wkt voronoi
POLYGON ((6.655199999999997 53.586600000000004, 6.660200000000003 53.586600000000004, 6.660200000000003 53.58584444444444, 6.65950000000001 53.585300000000004, 6.655199999999997 53.583866666666665, 6.655199999999997 53.586600000000004))
POLYGON ((6.655199999999997 53.583866666666665, 6.65950000000001 53.585300000000004, 6.659499999999982 53.585299999999975, 6.655199999999997 53.58243333333334, 6.655199999999997 53.583866666666665))
POLYGON ((6.660200000000003 53.581199999999995, 6.656425000000003 53.581199999999995, 6.659499999999982 53.585299999999975, 6.65950000000001 53.585300000000004, 6.660200000000003 53.58584444444444, 6.660200000000003 53.581199999999995))
POLYGON ((6.655199999999997 53.581199999999995, 6.655199999999997 53.58243333333334, 6.659499999999982 53.585299999999975, 6.656425000000003 53.581199999999995, 6.655199999999997 53.581199999999995))

Given the sensitivity to a few decimal places, it's probably a robustness issue.