18

新しい scipy (0.10 を使用) から scipy.spatial.Delaunay を実行した 3D には約 50,000 のデータ ポイントがあり、非常に便利な三角形分割が得られます。

基: http://en.wikipedia.org/wiki/Delaunay_triangulation (セクション「ボロノイ図との関係」)

...ボロノイ分割であるこの三角形分割の「デュアルグラフ」に到達する簡単な方法があるかどうか疑問に思っていました。

手がかりはありますか？これを調べてみると、組み込みの scipy 関数が表示されないようです。これはほとんど奇妙です!

ありがとう、エドワード

4

4 に答える 4

19

隣接情報はneighborsDelaunay オブジェクトの属性にあります。残念ながら、現時点ではコードは外心をユーザーに公開していないため、自分で再計算する必要があります。

また、無限に広がるボロノイ エッジは、この方法では直接得られません。おそらくまだ可能ですが、もう少し考える必要があります。

import numpy as np
from scipy.spatial import Delaunay

points = np.random.rand(30, 2)
tri = Delaunay(points)

p = tri.points[tri.vertices]

# Triangle vertices
A = p[:,0,:].T
B = p[:,1,:].T
C = p[:,2,:].T

# See http://en.wikipedia.org/wiki/Circumscribed_circle#Circumscribed_circles_of_triangles
# The following is just a direct transcription of the formula there
a = A - C
b = B - C

def dot2(u, v):
    return u[0]*v[0] + u[1]*v[1]

def cross2(u, v, w):
    """u x (v x w)"""
    return dot2(u, w)*v - dot2(u, v)*w

def ncross2(u, v):
    """|| u x v ||^2"""
    return sq2(u)*sq2(v) - dot2(u, v)**2

def sq2(u):
    return dot2(u, u)

cc = cross2(sq2(a) * b - sq2(b) * a, a, b) / (2*ncross2(a, b)) + C

# Grab the Voronoi edges
vc = cc[:,tri.neighbors]
vc[:,tri.neighbors == -1] = np.nan # edges at infinity, plotting those would need more work...

lines = []
lines.extend(zip(cc.T, vc[:,:,0].T))
lines.extend(zip(cc.T, vc[:,:,1].T))
lines.extend(zip(cc.T, vc[:,:,2].T))

# Plot it
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection

lines = LineCollection(lines, edgecolor='k')

plt.hold(1)
plt.plot(points[:,0], points[:,1], '.')
plt.plot(cc[0], cc[1], '*')
plt.gca().add_collection(lines)
plt.axis('equal')
plt.xlim(-0.1, 1.1)
plt.ylim(-0.1, 1.1)
plt.show()
于 2012-05-18T17:18:09.233 に答える
10

これにかなりの時間を費やしたので、エッジだけでなくボロノイポリゴンを取得する方法についての解決策を共有したいと思います。

コードはhttps://gist.github.com/letmaik/8803860にあり、 tauranのソリューションを拡張しています。

まず、ポイント座標ではなくインデックスを操作するときに多くの計算を簡略化できるため、コードを変更して、頂点と (ペアの) インデックス (= エッジ) を別々に提供するようにしました。

次に、voronoi_cell_linesメソッドで、どのエッジがどのセルに属しているかを判断します。そのために、関連する質問からAlinkの提案されたソリューションを使用します。つまり、エッジごとに最も近い 2 つの入力ポイント (=セル) を見つけ、そこからマッピングを作成します。

最後のステップは、実際のポリゴンを作成することです (voronoi_polygons方法を参照)。まず、ダングリング エッジを持つ外側のセルを閉じる必要があります。これは、すべてのエッジを調べて、隣接するエッジが 1 つしかないエッジを確認するのと同じくらい簡単です。そのようなエッジはゼロまたは 2 つ存在する可能性があります。2 つの場合は、追加のエッジを導入してこれらを接続します。

最後に、各セルの順序付けられていないエッジを正しい順序に並べて、そこからポリゴンを派生させる必要があります。

使用法は次のとおりです。

P = np.random.random((100,2))

fig = plt.figure(figsize=(4.5,4.5))
axes = plt.subplot(1,1,1)

plt.axis([-0.05,1.05,-0.05,1.05])

vertices, lineIndices = voronoi(P)        
cells = voronoi_cell_lines(P, vertices, lineIndices)
polys = voronoi_polygons(cells)

for pIdx, polyIndices in polys.items():
    poly = vertices[np.asarray(polyIndices)]
    p = matplotlib.patches.Polygon(poly, facecolor=np.random.rand(3,1))
    axes.add_patch(p)

X,Y = P[:,0],P[:,1]
plt.scatter(X, Y, marker='.', zorder=2)

plt.axis([-0.05,1.05,-0.05,1.05])
plt.show()

出力：

ボロノイ多角形

このコードは、多数の入力ポイントにはおそらく適しておらず、いくつかの領域で改善することができます。それでも、同様の問題を抱えている他の人にとっては役立つかもしれません。

于 2014-02-04T14:17:20.050 に答える
8

私は同じ問題に遭遇し、ウェブ上で見つけた pv. の回答と他のコード スニペットからソリューションを構築しました。解は、隣接する三角形が存在しない外側の線を含む完全なボロノイ図を返します。

#!/usr/bin/env python
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

def voronoi(P):
    delauny = Delaunay(P)
    triangles = delauny.points[delauny.vertices]

    lines = []

    # Triangle vertices
    A = triangles[:, 0]
    B = triangles[:, 1]
    C = triangles[:, 2]
    lines.extend(zip(A, B))
    lines.extend(zip(B, C))
    lines.extend(zip(C, A))
    lines = matplotlib.collections.LineCollection(lines, color='r')
    plt.gca().add_collection(lines)

    circum_centers = np.array([triangle_csc(tri) for tri in triangles])

    segments = []
    for i, triangle in enumerate(triangles):
        circum_center = circum_centers[i]
        for j, neighbor in enumerate(delauny.neighbors[i]):
            if neighbor != -1:
                segments.append((circum_center, circum_centers[neighbor]))
            else:
                ps = triangle[(j+1)%3] - triangle[(j-1)%3]
                ps = np.array((ps[1], -ps[0]))

                middle = (triangle[(j+1)%3] + triangle[(j-1)%3]) * 0.5
                di = middle - triangle[j]

                ps /= np.linalg.norm(ps)
                di /= np.linalg.norm(di)

                if np.dot(di, ps) < 0.0:
                    ps *= -1000.0
                else:
                    ps *= 1000.0
                segments.append((circum_center, circum_center + ps))
    return segments

def triangle_csc(pts):
    rows, cols = pts.shape

    A = np.bmat([[2 * np.dot(pts, pts.T), np.ones((rows, 1))],
                 [np.ones((1, rows)), np.zeros((1, 1))]])

    b = np.hstack((np.sum(pts * pts, axis=1), np.ones((1))))
    x = np.linalg.solve(A,b)
    bary_coords = x[:-1]
    return np.sum(pts * np.tile(bary_coords.reshape((pts.shape[0], 1)), (1, pts.shape[1])), axis=0)

if __name__ == '__main__':
    P = np.random.random((300,2))

    X,Y = P[:,0],P[:,1]

    fig = plt.figure(figsize=(4.5,4.5))
    axes = plt.subplot(1,1,1)

    plt.scatter(X, Y, marker='.')
    plt.axis([-0.05,1.05,-0.05,1.05])

    segments = voronoi(P)
    lines = matplotlib.collections.LineCollection(segments, color='k')
    axes.add_collection(lines)
    plt.axis([-0.05,1.05,-0.05,1.05])
    plt.show()

黒線＝ボロノイ図、赤線＝ドロニー三角形 黒線＝ボロノイ図、赤線＝ドロニー三角形

于 2013-04-03T09:27:58.450 に答える
0

これを行う関数はわかりませんが、それほど複雑な作業ではないようです。

ウィキペディアの記事で説明されているように、ボロノイグラフは外接円の接合部です。

したがって、基本的な数学である三角形の外接円の中心を見つける関数から始めることができます（http://en.wikipedia.org/wiki/Circumscribed_circle）。

次に、隣接する三角形の中心を結合します。

于 2012-05-18T10:12:47.157 に答える