Python：從3D中的Scipy's Delaunay三角測量計算Voronoi Tesselation

Question

我有大約50,000個3D數據點，我已經從新的scipy（我使用0.10）運行scipy.spatial.Delaunay，它給了我一個非常有用的三角剖分。

基於：

...我想知道如果有一個簡單的方法來得到這個三角的「雙圖」，這是維諾http://en.wikipedia.org/wiki/Delaunay_triangulation（節「與Voronoi圖的關係」）鑲嵌。

任何線索？我在這方面的搜索似乎沒有預先建立在scipy函數中，我發現它幾乎很奇怪！

感謝， 愛德華

Answer 1

可以在Delaunay對象的neighbors屬性中找到鄰接信息。不幸的是，此時代碼並未向用戶公開外核，因此您必須自己重新計算這些內容。

另外，延伸到無窮大的Voronoi邊緣不是直接以這種方式獲得的。這仍然可能，但需要更多的思考。

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

Answer 2

我不知道一個函數來做到這一點，但它似乎並不像一個過於複雜的任務。

Voronoi圖是外接圓的交點，如維基百科文章中所述。

所以，你可以從一個函數開始，找到三角形的外接圓的中心，這是一個基本的數學（http://en.wikipedia.org/wiki/Circumscribed_circle）。

然後，只需加入相鄰三角形的中心。

Answer 3

我遇到了同樣的問題，並在pv。的答案和我在網上找到的其他代碼片段中構建了一個解決方案。該解決方案返回一個完整的Voronoi圖，包括沒有三角形鄰居的外線。

#!/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() 

黑線= Voronoi圖，紅線= Delauny三角形

Answer 4

正如我花在這個相當長的時間，我想與大家分享我對如何讓維諾解決方案多邊形而不只是邊緣。

該代碼爲https://gist.github.com/letmaik/8803860，並擴展到tauran的解決方案。

首先，我改變了代碼以分別給出頂點和（對）索引（=邊），因爲在處理索引而不是點座標時可以簡化許多計算。

然後，在voronoi_cell_lines方法中，我確定哪些邊緣屬於哪些單元。爲此，我使用相關問題中提出的Alink解決方案。也就是說，對於每個邊找到兩個最近的輸入點（=單元格）並從中創建一個映射。

最後一步是創建實際的多邊形（參見voronoi_polygons方法）。首先，需要關閉具有懸掛邊緣的外部單元。這與查看所有邊緣並檢查哪些邊緣只有一個相鄰邊緣一樣簡單。可以有零個或兩個這樣的邊緣。如果是兩個，我通過引入一個額外的邊緣來連接它們。

最後，每個單元格中的無序邊緣需要按照正確的順序放置，以便從中導出多邊形。

的用法是：

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

，其輸出：

的代碼可能是不適合於大量的輸入點，並且可以在某些方面得到改善。不過，對於有類似問題的其他人可能會有所幫助。