Opened 5 years ago
Closed 4 years ago
#976 closed defect (fixed)
Voronoi polygons robustness issue
Reported by: | komzpa | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | 3.9.0 |
Component: | Default | Version: | main |
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 (29)
comment:1 by , 5 years ago
comment:2 by , 5 years ago
0104000080170000000101000080EC51B81E11A20741EC51B81E85A51C415C8FC2F528DC354001010000801F85EB5114A207415C8FC2F585A51C417B14AE47E1BA3540010100008085EB51B818A20741A8C64B3786A51C413E0AD7A3709D35400101000080000000001BA20741FED478E984A51C413E0AD7A3709D3540010100008085EB51B818A20741FED478E984A51C413E0AD7A3709D354001010000800AD7A37016A20741FED478E984A51C413E0AD7A3709D35400101000080000000001BA2074154E3A59B83A51C413E0AD7A3709D3540010100008085EB51B818A2074154E3A59B83A51C413E0AD7A3709D354001010000800AD7A37016A2074154E3A59B83A51C413E0AD7A3709D35400101000080000000001BA20741AAF1D24D82A51C413E0AD7A3709D3540010100008085EB51B818A20741AAF1D24D82A51C413E0AD7A3709D35400101000080F6285C8F12A20741EC51B81E88A51C414160E5D022DB354001010000802222222210A2074152B81EC586A51C414160E5D022DB354001010000804F1BE8B40DA2074152B81EC586A51C414160E5D022DB354001010000807B14AE470BA2074152B81EC586A51C414160E5D022DB354001010000802222222210A20741B81E856B85A51C414160E5D022DB354001010000804F1BE8B40DA20741B81E856B85A51C414160E5D022DB354001010000807B14AE470BA20741B81E856B85A51C414160E5D022DB35400101000080A70D74DA08A20741B81E856B85A51C414160E5D022DB35400101000080D4063A6D06A20741B81E856B85A51C414160E5D022DB354001010000807B14AE470BA207411F85EB1184A51C414160E5D022DB35400101000080A70D74DA08A207411F85EB1184A51C414160E5D022DB35400101000080D4063A6D06A207411F85EB1184A51C414160E5D022DB3540
comment:3 by , 5 years ago
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 by , 5 years ago
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.
by , 5 years ago
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
by , 5 years ago
Attachment: | invalid-voronoi-polygons.png added |
---|
Shows one of the invalid Voronoi polygons, containing a self-intersection
by , 5 years ago
Attachment: | invalid-polygon-detail.png added |
---|
Detail of an invalid polygon, showing very close vertices causing self-intersection
comment:5 by , 5 years ago
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 by , 5 years ago
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 by , 5 years ago
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 by , 5 years ago
Reported upstream in JTS: https://github.com/locationtech/jts/issues/447
comment:10 by , 5 years ago
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 by , 5 years ago
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 by , 5 years ago
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 by , 5 years ago
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.
follow-up: 16 comment:15 by , 5 years ago
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 by , 5 years ago
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 by , 5 years ago
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 by , 5 years ago
Milestone: | → 3.8.0 |
---|---|
Resolution: | fixed |
Status: | closed → reopened |
Version: | 3.6.2 → master |
Actually, I just ran this particular test against master and it still fails. So the fix needs more testing in GEOS land.
comment:23 by , 5 years ago
Milestone: | 3.8.0 → 3.9.0 |
---|
Can you post plain 'ol WKB?