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 40 forks source link

DelaunayTriangulation.DivideAndConquer nontermination #173

Open aavogt opened 2 years ago

aavogt commented 2 years ago

If I load this file here with cabal repl, bd and br get stuck after printing out _neighbours =. DelaunayTriangulation.Naive works, but it is too slow to help me find the two closest points between two sets of points.

I am using hgeometry-0.13 plus this patch https://github.com/noinia/hgeometry/commit/a1d6f20a98669a1dc1493bb47f00da90e32c29e7

noinia commented 2 years ago

Thanks for reporting this! I'll try to have a look at this this weekend.

noinia commented 2 years ago

Seems this issue is related to having multiple colinear points. I've managed to reduce it to a point set with only 5 points. Let's see if I can fix it :).

useronym commented 2 years ago

I think I stumbled upon the same issue. I am trying to compute the minimum spanning tree (that's using triangulation in the process afaik) but it's not terminating. I can provide an example with 10 points, in case that's of any help.

hgeometry-0.13

noinia commented 2 years ago

Ok, I finally had some time to think a bit more about this again. The core of the issue is that the input point set has colinear triples, for example (54,178), (56,176), and (58,174). (See also the minimal example below). For such a set of input points there does not exist a valid Delaunay triangulation. I assume you are hoping that the function constructs the delaunay graph (i.e. the dual graph of the Voronoi Diagram) in this situation?

buggyM :: NonEmpty.NonEmpty (Point 2 R :+ ())
buggyM = fmap (\(a,b) -> Point2 (fromIntegral a) (fromIntegral b) :+ ()) $ NonEmpty.fromList
  [ (33,261)
  , (54,178)
  , (56,176)
  , (58,174)
  , (61,174)
  ]
useronym commented 2 years ago

Well, here is for example what scipy does. You can notice that there is an extra edge between the middle colinear point and bottom right that one could do without. I think as long as it does something that is a valid triangulation instead of crashing, it's acceptable.

import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

x = [33, 54, 56, 58, 61]
y = [261, 178, 176, 174, 174]

tri = Delaunay(list(zip(x, y)))
plt.triplot(x, y, tri.simplices)

plt.plot(x, y, 'o')

plt.savefig("delaunay.png", dpi=200)

delaunay