DelaunayTriangulation.DivideAndConquer sometimes loops forever/gets stuck

[patrickt]

I’ve been working on a little project that uses hgeometry, and so far I’m loving the library. In using it, however, I’ve found some circumstances where the divide-and-conquer algorithm for Delaunay triangulation loops infinitely, such as when given these points:

Point2 [217.44781269876754,249.24741543498274] :+ 'a'
Point2 [237.91428506927295,105.8082929316906] :+ 'b'
Point2 [51.46936876163245,193.21960885915342] :+ 'c'
Point2 [172.55365082143922,2.8346743864823387] :+ 'd'
Point2 [250.55083565080437,93.13205719006257] :+ ‘e’

To reproduce:

• git clone [email protected]:patrickt/voronoid.git
• git checkout a2d95e00f4774b19b03c640442d7c52d3f871c49
• cabal new-build
• cabal new-run

Approximately 75% of the time, this will loop infinitely (or at least for a very long time). If it doesn’t, run it a few times.

I ran some profiling tests, and it appears to be stuck in Data.Geometry.Ball.disk. I’m new to the codebase so I wasn’t able to debug any further. For the time being I’m using the Naive module, which is good enough for my purposes, but I figured you’d want a bug report.

voronoid.prof.txt

[noinia]

Thanks for the kind words and the report. I'll try to figure out what is happening exactly!

[glittershark]

Not sure if we necessarily need other test cases, but just in case it helps this also reproduces the issue:

[ Point2 [38.5,3.5] :+ 0
, Point2 [67.0,33.0] :+ 1
, Point2 [46.0,45.5] :+ 2
, Point2 [55.5,42.0] :+ 3
, Point2 [36.0,25.0] :+ 4
, Point2 [76.5,12.0] :+ 5
, Point2 [29.0,26.5] :+ 6
, Point2 [55.0,10.5] :+ 7
]

[glittershark]

Actually maybe it would be useful to look at my example, since I've actually seen that hang 100% of the time

[noinia]

Thanks for reporting the specific testcase! What "real number type" (i.e. type 'r') are you using? If you are using Doubles or Floats you may want to try Rational numbers instead; some of the geometric predicates used don't go that well with floating point numbers (see also the README).

[glittershark]

they're Doubles.

[glittershark]

I don't necessarily need the exact right answers in this case - an approximation is fine - though I would like them within a reasonable amount of time 😄

[patrickt]

Aha, interesting! I’ll try Rational and see if that helps.

[lemmih]

@patrickt Did the problem go away when you used Rationals?

[noinia]

Just as a heads-up; I've been working on an implementation of an O(n log n) time algorithm to compute the convex hull in R^3 (see hgeometry-devel). That then will also give easy O(n log n) time algorithms for delaunay triangulation and Voronoi diagrams (in R^2). This could then be an alternative to the current implementation.

In trying to deal with all sorts of degeneracies in the 3D CH properly this all has turned into some sort of yak-shaving project.... causing me to implement some stuff on symbolic computation, and reimmplement some geometric primitives. I'm hoping I'm close to the end of the tunnel for that one now though :).