7

ポイントリスト=[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。なぜなら、船体だけを計算することは、原則として三角測量を計算するよりも速く実行できるからです。残念ながら、Scipyには、Qhullを使用して船体を直接計算するためのインターフェースはまだありません。

于 2011-06-30T22:19:29.150 に答える