Opened 5 years ago

Closed 3 years ago

#976 closed defect (fixed)

Voronoi polygons robustness issue

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

sites-and-delaunay.png (12.7 KB ) - added by mdavis 5 years 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 5 years ago.
Shows one of the invalid Voronoi polygons, containing a self-intersection
invalid-polygon-detail.png (46.1 KB ) - added by mdavis 5 years ago.
Detail of an invalid polygon, showing very close vertices causing self-intersection

Download all attachments as: .zip

Change History (29)

comment:1 by mdavis, 5 years ago

Can you post plain 'ol WKB?

comment:2 by komzpa, 5 years ago

0104000080170000000101000080EC51B81E11A20741EC51B81E85A51C415C8FC2F528DC354001010000801F85EB5114A207415C8FC2F585A51C417B14AE47E1BA3540010100008085EB51B818A20741A8C64B3786A51C413E0AD7A3709D35400101000080000000001BA20741FED478E984A51C413E0AD7A3709D3540010100008085EB51B818A20741FED478E984A51C413E0AD7A3709D354001010000800AD7A37016A20741FED478E984A51C413E0AD7A3709D35400101000080000000001BA2074154E3A59B83A51C413E0AD7A3709D3540010100008085EB51B818A2074154E3A59B83A51C413E0AD7A3709D354001010000800AD7A37016A2074154E3A59B83A51C413E0AD7A3709D35400101000080000000001BA20741AAF1D24D82A51C413E0AD7A3709D3540010100008085EB51B818A20741AAF1D24D82A51C413E0AD7A3709D35400101000080F6285C8F12A20741EC51B81E88A51C414160E5D022DB354001010000802222222210A2074152B81EC586A51C414160E5D022DB354001010000804F1BE8B40DA2074152B81EC586A51C414160E5D022DB354001010000807B14AE470BA2074152B81EC586A51C414160E5D022DB354001010000802222222210A20741B81E856B85A51C414160E5D022DB354001010000804F1BE8B40DA20741B81E856B85A51C414160E5D022DB354001010000807B14AE470BA20741B81E856B85A51C414160E5D022DB35400101000080A70D74DA08A20741B81E856B85A51C414160E5D022DB35400101000080D4063A6D06A20741B81E856B85A51C414160E5D022DB354001010000807B14AE470BA207411F85EB1184A51C414160E5D022DB35400101000080A70D74DA08A207411F85EB1184A51C414160E5D022DB35400101000080D4063A6D06A207411F85EB1184A51C414160E5D022DB3540

comment:3 by mdavis, 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 mdavis, 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 mdavis, 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 mdavis, 5 years ago

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

by mdavis, 5 years ago

Attachment: invalid-polygon-detail.png added

Detail of an invalid polygon, showing very close vertices causing self-intersection

comment:5 by mdavis, 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 mdavis, 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 komzpa, 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:9 by mdavis, 5 years ago

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

comment:10 by mdavis, 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 mdavis, 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 mdavis, 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 mdavis, 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.

comment:15 by komzpa, 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?

in reply to:  15 comment:16 by mdavis, 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 dbaston, 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:18 by mdavis, 5 years ago

Agreed, ttmath seems like the best solution.

comment:20 by pramsey, 5 years ago

Resolution: fixed
Status: newclosed

Pushed fix to master

comment:21 by pramsey, 5 years ago

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 5 years ago by pramsey (previous) (diff)

comment:22 by pramsey, 5 years ago

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

comment:23 by pramsey, 5 years ago

Milestone: 3.8.03.9.0

comment:25 by Sandro Santilli <strk@…>, 3 years ago

In 26cd478/git:

Remove incorrect isoceles heuristic in QuadEdgeSubdivision::TriangleCircumcentreVisitor, fix unit test

References #976 in master (3.9.0dev)

comment:26 by Sandro Santilli <strk@…>, 3 years ago

Resolution: fixed
Status: reopenedclosed

In e4a0126/git:

Remove incorrect isoceles heuristic in QuadEdgeSubdivision::TriangleCircumcentreVisitor

Fix related unit test

Closes #976 in 3.8 branch (3.8.2dev)

The heuristic was added with 56417ce55

Note: See TracTickets for help on using tickets.