Search code examples
pythonscipytriangulationdelaunay

python scipy Delaunay plotting point cloud


I have a pointlist=[p1,p2,p3...] where p1 = [x1,y1],p2=[x2,y2] ...

I want to use scipy.spatial.Delaunay to do trianglation on these point clouds and then plot it

How can i do this ?

The documentation for the Delaunay is really scarce

so far i have this code

from subprocess import Popen, PIPE
import os


os.environ['point_num'] = "2000"

cmd = 'rbox $point_num D2 | tail -n $point_num'
sub_process = Popen(cmd, shell=True,stdout=PIPE,stderr=PIPE)
output = sub_process.communicate()
points = [line.split() for line in output[0].split('\n') if line]
x = [p[0] for p in points if p]
y = [p[1] for p in points if p]

import matplotlib.pyplot as plt
plt.plot(x,y,'bo')

from scipy.spatial import Delaunay

dl = Delaunay(points)
convex = dl.convex_hull

from numpy.core.numeric import reshape,shape
convex = reshape(convex,(shape(convex)[0]*shape(convex)[1],1))
convex_x = [x[i] for i in convex]
convex_y = [y[i] for i in convex]

plt.plot(convex_x,convex_y,'r')
plt.show()

Thanks


Solution

  • EDIT: plot also the convex hull

    import numpy as np
    from scipy.spatial import Delaunay
    
    points = np.random.rand(30, 2) # 30 points in 2-d
    tri = Delaunay(points)
    
    # Make a list of line segments: 
    # edge_points = [ ((x1_1, y1_1), (x2_1, y2_1)),
    #                 ((x1_2, y1_2), (x2_2, y2_2)),
    #                 ... ]
    edge_points = []
    edges = set()
    
    def add_edge(i, j):
        """Add a line between the i-th and j-th points, if not in the list already"""
        if (i, j) in edges or (j, i) in edges:
            # already added
            return
        edges.add( (i, j) )
        edge_points.append(points[ [i, j] ])
    
    # loop over triangles: 
    # ia, ib, ic = indices of corner points of the triangle
    for ia, ib, ic in tri.vertices:
        add_edge(ia, ib)
        add_edge(ib, ic)
        add_edge(ic, ia)
    
    # plot it: the LineCollection is just a (maybe) faster way to plot lots of
    # lines at once
    import matplotlib.pyplot as plt
    from matplotlib.collections import LineCollection
    
    lines = LineCollection(edge_points)
    plt.figure()
    plt.title('Delaunay triangulation')
    plt.gca().add_collection(lines)
    plt.plot(points[:,0], points[:,1], 'o', hold=1)
    plt.xlim(-1, 2)
    plt.ylim(-1, 2)
    
    # -- the same stuff for the convex hull
    
    edges = set()
    edge_points = []
    
    for ia, ib in tri.convex_hull:
        add_edge(ia, ib)
    
    lines = LineCollection(edge_points)
    plt.figure()
    plt.title('Convex hull')
    plt.gca().add_collection(lines)
    plt.plot(points[:,0], points[:,1], 'o', hold=1)
    plt.xlim(-1, 2)
    plt.ylim(-1, 2)
    plt.show()
    

    Note that using scipy.spatial.Delaunay just for computing the complex hull is probably overkill, because computing just the hull can in principle done faster than computing the triangulation. Unfortunately, there's no interface in Scipy yet for computing hulls directly with Qhull.