One Voronoi cell created from a set of 10 points is non-convex
jaakkor2 opened this issue · 2 comments
Voronoi diagram of this WKB produces one non-convex cell. It is a ten-point multipoint.
UInt8[0x01, 0x04, 0x00, 0x00, 0x00, 0x0a, 0x00, 0x00, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x58, 0x69, 0xc8, 0xb9, 0xcf, 0xae, 0x32, 0x40, 0x1f, 0x85, 0xeb, 0x51, 0xb8, 0x06, 0x59, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0xc0, 0xe2, 0x2b, 0xbe, 0xe2, 0x0b, 0x30, 0x40, 0x1f, 0x85, 0xeb, 0x51, 0xb8, 0x06, 0x59, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x52, 0xb8, 0x1e, 0x85, 0xeb, 0xd1, 0x2a, 0x40, 0x1f, 0x85, 0xeb, 0x51, 0xb8, 0x06, 0x59, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x52, 0xb8, 0x1e, 0x85, 0xeb, 0xd1, 0x2a, 0x40, 0x13, 0x83, 0xc0, 0xca, 0xa1, 0x9d, 0x59, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x52, 0xb8, 0x1e, 0x85, 0xeb, 0xd1, 0x2a, 0x40, 0x07, 0x81, 0x95, 0x43, 0x8b, 0x34, 0x5a, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x52, 0xb8, 0x1e, 0x85, 0xeb, 0xd1, 0x2a, 0x40, 0xfa, 0x7e, 0x6a, 0xbc, 0x74, 0xcb, 0x5a, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x52, 0xb8, 0x1e, 0x85, 0xeb, 0xd1, 0x2a, 0x40, 0xee, 0x7c, 0x3f, 0x35, 0x5e, 0x62, 0x5b, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x52, 0xb8, 0x1e, 0x85, 0xeb, 0xd1, 0x2a, 0x40, 0xe2, 0x7a, 0x14, 0xae, 0x47, 0xf9, 0x5b, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0xbe, 0xe2, 0x2b, 0xbe, 0xe2, 0x0b, 0x30, 0x40, 0xe2, 0x7a, 0x14, 0xae, 0x47, 0xf9, 0x5b, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x54, 0x69, 0xc8, 0xb9, 0xcf, 0xae, 0x32, 0x40, 0xe2, 0x7a, 0x14, 0xae, 0x47, 0xf9, 0x5b, 0x40]
Here as WKT
MULTIPOINT ((18.68285714285716 100.105), (16.046428571428578 100.105), (13.41 100.105), (13.41 102.46300000000001), (13.41 104.82100000000001), (13.41 107.179), (13.41 109.537), (13.41 111.89500000000001), (16.04642857142857 111.89500000000001), (18.682857142857145 111.89500000000001))
Hopefully #953 fixes this, but I cannot easily test. Maybe this could be a test case.
The test code was made with Julia, LibGEOS.jl v0.8.5 and GEOS_jll v3.12.0.
import LibGEOS
wkb = UInt8[0x01, 0x04, 0x00, 0x00, 0x00, 0x0a, 0x00, 0x00, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x58, 0x69, 0xc8, 0xb9, 0xcf, 0xae, 0x32, 0x40, 0x1f, 0x85, 0xeb, 0x51, 0xb8, 0x06, 0x59, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0xc0, 0xe2, 0x2b, 0xbe, 0xe2, 0x0b, 0x30, 0x40, 0x1f, 0x85, 0xeb, 0x51, 0xb8, 0x06, 0x59, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x52, 0xb8, 0x1e, 0x85, 0xeb, 0xd1, 0x2a, 0x40, 0x1f, 0x85, 0xeb, 0x51, 0xb8, 0x06, 0x59, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x52, 0xb8, 0x1e, 0x85, 0xeb, 0xd1, 0x2a, 0x40, 0x13, 0x83, 0xc0, 0xca, 0xa1, 0x9d, 0x59, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x52, 0xb8, 0x1e, 0x85, 0xeb, 0xd1, 0x2a, 0x40, 0x07, 0x81, 0x95, 0x43, 0x8b, 0x34, 0x5a, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x52, 0xb8, 0x1e, 0x85, 0xeb, 0xd1, 0x2a, 0x40, 0xfa, 0x7e, 0x6a, 0xbc, 0x74, 0xcb, 0x5a, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x52, 0xb8, 0x1e, 0x85, 0xeb, 0xd1, 0x2a, 0x40, 0xee, 0x7c, 0x3f, 0x35, 0x5e, 0x62, 0x5b, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x52, 0xb8, 0x1e, 0x85, 0xeb, 0xd1, 0x2a, 0x40, 0xe2, 0x7a, 0x14, 0xae, 0x47, 0xf9, 0x5b, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0xbe, 0xe2, 0x2b, 0xbe, 0xe2, 0x0b, 0x30, 0x40, 0xe2, 0x7a, 0x14, 0xae, 0x47, 0xf9, 0x5b, 0x40, 0x01, 0x01, 0x00, 0x00, 0x00, 0x54, 0x69, 0xc8, 0xb9, 0xcf, 0xae, 0x32, 0x40, 0xe2, 0x7a, 0x14, 0xae, 0x47, 0xf9, 0x5b, 0x40]
mp = LibGEOS.readgeom(wkb)
# setting precision fixes the problem
#mp = LibGEOS.setPrecision(mp, 1e-7)
# GEOSVoronoiDiagram_r is not yet wrapped by LibGEOS.jl-interface
ctx = LibGEOS.get_context(mp)
env = LibGEOS.envelope(mp)
tol = 0.0
flags = 0x2
vrn = LibGEOS.geomFromGEOS(LibGEOS.GEOSVoronoiDiagram_r(ctx.ptr, mp.ptr, env.ptr, tol, flags))
# plot
using GLMakie, GeoInterfaceMakie
GeoInterfaceMakie.@enable(LibGEOS.AbstractGeometry)
v = LibGEOS.getGeometries(vrn)
fig = Figure()
ax = Axis(fig[1,1], aspect = DataAspect())
plot!(ax, mp, color = :red)
for v in v
plot!(ax, LibGEOS.boundary(v))
end
plot!(ax, LibGEOS.boundary(v[5]), color = :green, linewidth = 5)
#save("nonconvex.png", fig)
The non-convex cell is plotted with thick green line
Setting the precision to 1e-7
makes the problem go away, the line is commented out in the code.
Inputs which lie on a regular grid are challenging for the Voronoi computation. Due to numerical roundoff the circumcentre computation can produce slightly different values for circumcentres which in fact should be identical. This is what is happening here, and unfortunately the slightly different values end up creating an invalid polygon for one of the cells. The cell-clipping then runs into a known bug in OverlayNG, and results in a malformed cell polygon.
Unfortunately #953 does not fix this. At this point I'm not sure what a suitable fix is.
Perturb the inputs ever so slightly? Kind of like with buffer, there is a "close enough" aspect to the outputs of this operation.