Opened 4 months ago
Last modified 8 weeks ago
#976 reopened defect
Voronoi polygons robustness issue
Reported by: | komzpa | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | 3.9.0 |
Component: | Default | Version: | master |
Severity: | Unassigned | Keywords: | |
Cc: |
Description
https://trac.osgeo.org/postgis/ticket/4335
SELECT ST_VoronoiPolygons(ST_GeomFromWKB( E'\\001\\354\\003\\000\\000\\027\\000\\000\\000\\001\\351\\003\\000\\000\\354Q\\270\\036\\021\\242\\007A\\354Q\\270\\036\\205\\245\\034A\\\\\\217\\302\\365(\\3345@\\001\\351\\003\\000\\000\\037\\205\\353Q\\024\\242\\007A\\\\\\217\\302\\365\\205\\245\\034A{\\024\\256G\\341\\2725@\\001\\351\\003\\000\\000\\205\\353Q\\270\\030\\242\\007A\\250\\306K7\\206\\245\\034A>\\012\\327\\243p\\2355@\\001\\351\\003\\000\\000\\000\\000\\000\\000\\033\\242\\007A\\376\\324x\\351\\204\\245\\034A>\\012\\327\\243p\\2355@\\001\\351\\003\\000\\000\\205\\353Q\\270\\030\\242\\007A\\376\\324x\\351\\204\\245\\034A>\\012\\327\\243p\\2355@\\001\\351\\003\\000\\000\\012\\327\\243p\\026\\242\\007A\\376\\324x\\351\\204\\245\\034A>\\012\\327\\243p\\2355@\\001\\351\\003\\000\\000\\000\\000\\000\\000\\033\\242\\007AT\\343\\245\\233\\203\\245\\034A>\\012\\327\\243p\\2355@\\001\\351\\003\\000\\000\\205\\353Q\\270\\030\\242\\007AT\\343\\245\\233\\203\\245\\034A>\\012\\327\\243p\\2355@\\001\\351\\003\\000\\000\\012\\327\\243p\\026\\242\\007AT\\343\\245\\233\\203\\245\\034A>\\012\\327\\243p\\2355@\\001\\351\\003\\000\\000\\000\\000\\000\\000\\033\\242\\007A\\252\\361\\322M\\202\\245\\034A>\\012\\327\\243p\\2355@\\001\\351\\003\\000\\000\\205\\353Q\\270\\030\\242\\007A\\252\\361\\322M\\202\\245\\034A>\\012\\327\\243p\\2355@\\001\\351\\003\\000\\000\\366(\\\\\\217\\022\\242\\007A\\354Q\\270\\036\\210\\245\\034AA`\\345\\320"\\3335@\\001\\351\\003\\000\\000""""\\020\\242\\007AR\\270\\036\\305\\206\\245\\034AA`\\345\\320"\\3335@\\001\\351\\003\\000\\000O\\033\\350\\264\\015\\242\\007AR\\270\\036\\305\\206\\245\\034AA`\\345\\320"\\3335@\\001\\351\\003\\000\\000{\\024\\256G\\013\\242\\007AR\\270\\036\\305\\206\\245\\034AA`\\345\\320"\\3335@\\001\\351\\003\\000\\000""""\\020\\242\\007A\\270\\036\\205k\\205\\245\\034AA`\\345\\320"\\3335@\\001\\351\\003\\000\\000O\\033\\350\\264\\015\\242\\007A\\270\\036\\205k\\205\\245\\034AA`\\345\\320"\\3335@\\001\\351\\003\\000\\000{\\024\\256G\\013\\242\\007A\\270\\036\\205k\\205\\245\\034AA`\\345\\320"\\3335@\\001\\351\\003\\000\\000\\247\\015t\\332\\010\\242\\007A\\270\\036\\205k\\205\\245\\034AA`\\345\\320"\\3335@\\001\\351\\003\\000\\000\\324\\006:m\\006\\242\\007A\\270\\036\\205k\\205\\245\\034AA`\\345\\320"\\3335@\\001\\351\\003\\000\\000{\\024\\256G\\013\\242\\007A\\037\\205\\353\\021\\204\\245\\034AA`\\345\\320"\\3335@\\001\\351\\003\\000\\000\\247\\015t\\332\\010\\242\\007A\\037\\205\\353\\021\\204\\245\\034AA`\\345\\320"\\3335@\\001\\351\\003\\000\\000\\324\\006:m\\006\\242\\007A\\037\\205\\353\\021\\204\\245\\034AA`\\345\\320"\\3335@' ))
Attachments (3)
Change History (26)
comment:1 Changed 4 months ago by
comment:2 Changed 4 months ago by
0104000080170000000101000080EC51B81E11A20741EC51B81E85A51C415C8FC2F528DC354001010000801F85EB5114A207415C8FC2F585A51C417B14AE47E1BA3540010100008085EB51B818A20741A8C64B3786A51C413E0AD7A3709D35400101000080000000001BA20741FED478E984A51C413E0AD7A3709D3540010100008085EB51B818A20741FED478E984A51C413E0AD7A3709D354001010000800AD7A37016A20741FED478E984A51C413E0AD7A3709D35400101000080000000001BA2074154E3A59B83A51C413E0AD7A3709D3540010100008085EB51B818A2074154E3A59B83A51C413E0AD7A3709D354001010000800AD7A37016A2074154E3A59B83A51C413E0AD7A3709D35400101000080000000001BA20741AAF1D24D82A51C413E0AD7A3709D3540010100008085EB51B818A20741AAF1D24D82A51C413E0AD7A3709D35400101000080F6285C8F12A20741EC51B81E88A51C414160E5D022DB354001010000802222222210A2074152B81EC586A51C414160E5D022DB354001010000804F1BE8B40DA2074152B81EC586A51C414160E5D022DB354001010000807B14AE470BA2074152B81EC586A51C414160E5D022DB354001010000802222222210A20741B81E856B85A51C414160E5D022DB354001010000804F1BE8B40DA20741B81E856B85A51C414160E5D022DB354001010000807B14AE470BA20741B81E856B85A51C414160E5D022DB35400101000080A70D74DA08A20741B81E856B85A51C414160E5D022DB35400101000080D4063A6D06A20741B81E856B85A51C414160E5D022DB354001010000807B14AE470BA207411F85EB1184A51C414160E5D022DB35400101000080A70D74DA08A207411F85EB1184A51C414160E5D022DB35400101000080D4063A6D06A207411F85EB1184A51C414160E5D022DB3540
comment:3 Changed 4 months ago by
Thanks!
This is a robustness issue in the current JTS overlay algorithm. Might be able to be fixed with a smarter implementation of Voronoi polygons... but more likely will be fixed with the coming smarter implementation of overlay. I will log this as a JTS issue, to surface it at source.
comment:4 Changed 4 months ago by
This is due to a couple of the generated Voronoi polygons being topologically invalid, and because of that failing during the intersection computation (when clipping the raw Voronoi polygons to a surrounding rectangle).
The reason the polygons are invalid is that they contain two points which are almost identical, one of which happens to lie on another edge of the polygon, thus creating a self-intersection.
The reason the very close points are present is because the input points (Voronoi sites) generating them lie on a regular grid pattern. This produces two Delaunay triangles with *almost* identical circumcentres. They are not 100% identical because of round-off error. The DT circumcentres form the vertices of the Voronoi polygons, and so polygons are created with almost-but-not-identical points in them.
See attached diagrams for images.
Changed 4 months ago by
Attachment: | sites-and-delaunay.png added |
---|
Shows some input site vertices lying on a regular grid, and the Delaunay triangles containing them which have identical circumcentres
Changed 4 months ago by
Attachment: | invalid-voronoi-polygons.png added |
---|
Shows one of the invalid Voronoi polygons, containing a self-intersection
Changed 4 months ago by
Attachment: | invalid-polygon-detail.png added |
---|
Detail of an invalid polygon, showing very close vertices causing self-intersection
comment:5 Changed 4 months ago by
At the moment I'm not sure what an appropriate fix would be. Some sort of slight snapping, probably, but not sure where or how to determine an appropriate tolerance.
comment:6 Changed 4 months ago by
Perhaps another way of fixing this is to detect pairs of triangles which should have identical circumcentres and force this to be the case. I'm not sure it's possible to do this in full generality, but for points falling exactly on a square (or diamond) it is possible to robustly detect that the circumcentres should be the same.
comment:7 Changed 4 months ago by
Can longer arithmetics solve that? There's a patch of mine in master that promotes double to long double - is it reproducing on master?
comment:8 Changed 4 months ago by
Reported upstream in JTS: https://github.com/locationtech/jts/issues/447
comment:10 Changed 4 months ago by
This is the place to see if extra precision can provide a more stable circumcentre result:
https://github.com/libgeos/geos/blob/master/src/triangulate/quadedge/QuadEdgeSubdivision.cpp#L424
The circumcentre call could be replaced with one that uses higher precision arithmetic.
comment:11 Changed 4 months ago by
Here's a simple test case to reproduce the problem. It's two of the Delaunay triangles.
GEOMETRYCOLLECTION (POLYGON ((193600.80333333334 469345.355, 193600.80333333334 469345.0175, 193601.10666666666 469345.0175, 193600.80333333334 469345.355)), POLYGON ((193600.80333333334 469345.355, 193601.10666666666 469345.0175, 193601.10666666666 469345.355, 193600.80333333334 469345.355)))
Computing the circumcentres produces slightly different points:
GEOMETRYCOLLECTION (POINT (193600.95500000002 469345.18625), POINT (193600.95500000002 469345.18624999997) )
This can be tried with extended-precision to see if it produces identical points.
comment:12 Changed 4 months ago by
Good news... I implemented Triangle.circumcentreDD in JTS, using the JTS extended precision floating point provided by DD. The test case above now works, and so does the original Voronoi calculation.
I'll be committing this to JTS very shortly.
Based on the thread from Dec 2018 about DoubleDouble? vs long double [1], it seems like perhaps DoubleDouble? is the way to go for this in GEOS?
[1] http://osgeo-org.1560.x6.nabble.com/double-double-td5388315.html
comment:13 Changed 4 months ago by
It also looks like just about any set of points on a grid will trigger this issue, and is handled by the fix in the works. So this is definitely worthwhile to implement.
DD arithmetic is slower than double precision, of course, but this code is only used when generating Voronoi diagrams, and the slowdown is likely a small percentage of the runtime.
comment:15 follow-up: 16 Changed 4 months ago by
Based on the thread from Dec 2018 about DoubleDouble?? vs long double [1], it seems like perhaps DoubleDouble?? is the way to go for this in GEOS?
going to long double is just a matter of adding "long" all over the function, and has almost no performance penalty. Worth trying that first?
comment:16 Changed 4 months ago by
Replying to komzpa:
going to long double is just a matter of adding "long" all over the function, and has almost no performance penalty. Worth trying that first?
My reading suggests that compilers can choose how to implement long double, and some may choose to implement it as plain double (including clang and MSVC?). So that would not make it useful for GEOS in general- although we can try it to see if it works in some cases.
comment:17 Changed 4 months ago by
long double
is also not the same as DoubleDouble
even if it available. Given that GEOS is already set up to use the ttmath
library for extended-precision, I can't see a reason not to use it.
comment:21 Changed 2 months ago by
Milestone: | → 3.8.0 |
---|---|
Resolution: | fixed |
Status: | closed → reopened |
Version: | 3.6.2 → master |
Actually, I just ran this particular SQL test against master in PostGIS and it still fails. So the fix needs more unit testing in GEOS land and confirmation it actually does something useful.
comment:23 Changed 8 weeks ago by
Milestone: | 3.8.0 → 3.9.0 |
---|
Can you post plain 'ol WKB?