search

noinia / hgeometry

HGeometry is a library for computing with geometric objects in Haskell. It defines basic geometric types and primitives, and it implements some geometric data structures and algorithms.
122 stars 41 forks source link

DelaunayTriangulation.DivideAndConquer sometimes loops forever/gets stuck #28

Open patrickt opened 5 years ago

patrickt commented 5 years ago

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:

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 commented 5 years ago

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

glittershark commented 4 years ago

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 commented 4 years ago

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

noinia commented 4 years ago

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 commented 4 years ago

they're Doubles.

glittershark commented 4 years ago

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 :smile:

patrickt commented 4 years ago

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

lemmih commented 4 years ago

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

noinia commented 4 years ago

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 :).