7

我有一个 pointlist=[p1,p2,p3...] 其中 p1 = [x1,y1],p2=[x2,y2] ...

我想使用 scipy.spatial.Delaunay 对这些点云进行三角剖分，然后绘制它

我怎样才能做到这一点 ？

Delaunay 的文档真的很稀缺

到目前为止我有这个代码

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()

谢谢

4

1 回答 1

18

编辑：也绘制凸包

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()

请注意，scipy.spatial.Delaunay仅用于计算复杂的 hull 可能是多余的，因为原则上仅计算 hull 可以比计算三角剖分更快。不幸的是，Scipy 中还没有直接使用 Qhull 计算外壳的接口。

于 2011-06-30T22:19:29.150 回答