Opened 4 months ago

Last modified 8 weeks ago

#976 reopened defect

Voronoi polygons robustness issue

Reported by: komzpa Owned by: geos-devel@…
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)

sites-and-delaunay.png (12.7 KB) - added by mdavis 4 months ago.
Shows some input site vertices lying on a regular grid, and the Delaunay triangles containing them which have identical circumcentres
invalid-voronoi-polygons.png (17.1 KB) - added by mdavis 4 months ago.
Shows one of the invalid Voronoi polygons, containing a self-intersection
invalid-polygon-detail.png (46.1 KB) - added by mdavis 4 months ago.
Detail of an invalid polygon, showing very close vertices causing self-intersection

Download all attachments as: .zip

Change History (26)

comment:1 Changed 4 months ago by mdavis

Can you post plain 'ol WKB?

comment:2 Changed 4 months ago by komzpa

0104000080170000000101000080EC51B81E11A20741EC51B81E85A51C415C8FC2F528DC354001010000801F85EB5114A207415C8FC2F585A51C417B14AE47E1BA3540010100008085EB51B818A20741A8C64B3786A51C413E0AD7A3709D35400101000080000000001BA20741FED478E984A51C413E0AD7A3709D3540010100008085EB51B818A20741FED478E984A51C413E0AD7A3709D354001010000800AD7A37016A20741FED478E984A51C413E0AD7A3709D35400101000080000000001BA2074154E3A59B83A51C413E0AD7A3709D3540010100008085EB51B818A2074154E3A59B83A51C413E0AD7A3709D354001010000800AD7A37016A2074154E3A59B83A51C413E0AD7A3709D35400101000080000000001BA20741AAF1D24D82A51C413E0AD7A3709D3540010100008085EB51B818A20741AAF1D24D82A51C413E0AD7A3709D35400101000080F6285C8F12A20741EC51B81E88A51C414160E5D022DB354001010000802222222210A2074152B81EC586A51C414160E5D022DB354001010000804F1BE8B40DA2074152B81EC586A51C414160E5D022DB354001010000807B14AE470BA2074152B81EC586A51C414160E5D022DB354001010000802222222210A20741B81E856B85A51C414160E5D022DB354001010000804F1BE8B40DA20741B81E856B85A51C414160E5D022DB354001010000807B14AE470BA20741B81E856B85A51C414160E5D022DB35400101000080A70D74DA08A20741B81E856B85A51C414160E5D022DB35400101000080D4063A6D06A20741B81E856B85A51C414160E5D022DB354001010000807B14AE470BA207411F85EB1184A51C414160E5D022DB35400101000080A70D74DA08A207411F85EB1184A51C414160E5D022DB35400101000080D4063A6D06A207411F85EB1184A51C414160E5D022DB3540

comment:3 Changed 4 months ago by mdavis

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 mdavis

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 mdavis

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 mdavis

Shows one of the invalid Voronoi polygons, containing a self-intersection

Changed 4 months ago by mdavis

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 mdavis

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 mdavis

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 komzpa

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 mdavis

comment:9 Changed 4 months ago by mdavis

Extra precision might fix the problem. Have to try it and see.

comment:10 Changed 4 months ago by mdavis

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 mdavis

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 mdavis

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 mdavis

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 Changed 4 months ago by komzpa

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 in reply to:  15 Changed 4 months ago by mdavis

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 dbaston

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:18 Changed 4 months ago by mdavis

Agreed, ttmath seems like the best solution.

comment:20 Changed 4 months ago by pramsey

Resolution: fixed
Status: newclosed

Pushed fix to master

comment:21 Changed 2 months ago by pramsey

Milestone: 3.8.0
Resolution: fixed
Status: closedreopened
Version: 3.6.2master

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.

Last edited 2 months ago by pramsey (previous) (diff)

comment:22 Changed 2 months ago by pramsey

The test case from #859 should be checked when this case passes.

comment:23 Changed 8 weeks ago by pramsey

Milestone: 3.8.03.9.0
Note: See TracTickets for help on using tickets.