Search code examples
pythonscipydelaunayscipy-spatialqhull

Python Scipy Delaunay 2D doesn't connect some points

enter image description here

I generated 20 random points and applied Delaunay triangulation and created a visualization. I can't understand why there are always some points on the outside that aren't connected. I've looked in the qhull documentation, but I can't find the setting that would solve my issue. I also tried changing the seed, but even with different numbers, the problem persists. I tried the opitions in scipy.spatial.Delaunay like furthest_site, incremental and qhull_options="QJ" but they didn't solved the issue.

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

#The points that I used

np.random.seed(10)  
points = 10 * np.random.rand(20, 2)


#qhull_options = "QJ"

tri = Delaunay(points)



plt.plot(points[:, 0], points[:, 1], 'o', markersize=8, label='Points')

for simplex in tri.simplices:
    plt.plot(points[simplex, 0], points[simplex, 1], 'k-')
    #point1 = points[simplex, 0]
    #point2 = points[simplex, 1]
    #distance = np.linalg.norm(point1 - point2)
    #point_pairs.append((point1, point2))
    #distances.append(distance)
    

plt.xlabel('X')
plt.ylabel('Y')
plt.title('Delaunay Triangulation ')
plt.legend()
plt.axis('equal')  

plt.show()


"""
[[7.71320643 0.20751949]
 [6.33648235 7.48803883]
 [4.98507012 2.24796646]
 [1.98062865 7.60530712]
 [1.69110837 0.88339814]
 [6.85359818 9.53393346]
 [0.03948266 5.12192263]
 [8.12620962 6.12526067]
 [7.21755317 2.91876068]
 [9.17774123 7.14575783]
 [5.42544368 1.42170048]
 [3.7334076  6.74133615]
 [4.41833174 4.34013993]
 [6.17766978 5.13138243]
 [6.50397182 6.01038953]
 [8.05223197 5.21647152]
 [9.08648881 3.19236089]
 [0.90459349 3.00700057]
 [1.13984362 8.28681326]
 [0.46896319 6.26287148]]
"""

'''
Solution

  • A Minimal example

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

np.random.seed(1234)
points=np.random.rand(10,2)

tri=Delaunay(points)

plt.scatter(points[:,0], points[:,1])
for s in tri.simplices:
    plt.plot(points[s,0], points[s,1], 'k-')

plt.show()

    Its shows:

    enter image description here

    I was lucky with my seed, because (unlike yours — even tho, without bragging, I knew anyway what was the problem from yours too) it make it obvious what the problem is: that single line at the bottom cannot be part of a triangle (on your example, you seem to have only closed triangles

    Plot just one triangle

    The reason why, your, otherwise quite clever, way to draw triangles, just draw 2 lines of each triangle. Each simplex is a triplet of indexes. You use fancy indexing to turn those into a triplet of x and a triplet of y. And plot those. But in matplotlib, when you plot 3 x and 3 y you get 2 lines, not 3.

    Just draw the 1st triangle only to be sure (we are lucky again, that 1st triangle happen to be the one we already noticed)

    plt.scatter(points[:,0], points[:,1])
for s in tri.simplices:
    plt.plot(points[s,0], points[s,1], 'k-')
    break

    enter image description here

    So, you see, 2 lines (tho I am sure you didn't really need convincing at this point. That is the kind of bug sometimes hard to find, but that seem obvious once seen).

    Close the loops

    So, one, solution could be to use plt.fill instead (with option fill=False, which may seem oxymoronic, since that means fill, but don't fill. But fill role is to close the loop, when drawing, and to fill the inside, unless asked not to)

    plt.scatter(points[:,0], points[:,1])
for s in tri.simplices:
    plt.fill(points[s,0], points[s,1], fill=False)

    enter image description here

    Another option could have been to alter all s to add the first point at the end

    plt.scatter(points[:,0], points[:,1])
for s in tri.simplices:
    t=s.tolist()+[s[0]]
    plt.plot(points[t,0], points[t,1], 'k-')

    showing the exact same result.

    But fill, without fill=False can also make it quite clear where the triangles are:

    for s in tri.simplices:
    plt.fill(points[s,0], points[s,1])
plt.scatter(points[:,0], points[:,1]) # I need to draw the points after, or otherwise they are covered by filled triangles

    enter image description here

    Your seed

    I missed the fact that, while I was playing with my own example, you edited your post to create a MRE. So now I can show what happens with your seed enter image description here