Python的:从SciPy的的Delaunay三角在3D的Voronoi计算Tesselation的SciPy、Python、Delaunay、Tesselation

2023-09-08 00:14:28 作者:月稀掩朦胧

我对我所来自新SciPy的(我使用的是0.10),这给了我一个非常有用的三角运行scipy.spatial.Delaunay在3D约50,000个数据点。

I have about 50,000 data points in 3D on which I have run scipy.spatial.Delaunay from the new scipy (I'm using 0.10) which gives me a very useful triangulation.

根据http://en.wikipedia.org/wiki/Delaunay_triangulation (一节的Voronoi图的关系)

Based on: http://en.wikipedia.org/wiki/Delaunay_triangulation (section "Relationship with the Voronoi diagram")

......我在想,如果有一个简单的方法去双图这个三角的,这是沃罗诺伊Tesselation的。

...I was wondering if there is an easy way to get to the "dual graph" of this triangulation, which is the Voronoi Tesselation.

任何线索?我摸索这个似乎没有显示pre-内置SciPy的功能,我觉得这几乎是陌生的!

Any clues? My searching around on this seems to show no pre-built in scipy functions, which I find almost strange!

谢谢, 爱德华

推荐答案

邻接信息可以在邻居中找到的德劳内对象的属性。不幸的是,code没有的circumcenters暴露此刻的用户,所以你必须要重新计算的自己。

The adjacency information can be found in the neighbors attribute of the Delaunay object. Unfortunately, the code does not expose the circumcenters to the user at the moment, so you'll have to recompute those yourself.

此外,该沃罗诺伊边缘延伸到无穷不直接以这种方式获得。它仍然可能是可行的,但需要更多的思考。

Also, the Voronoi edges that extend to infinity are not directly obtained in this way. It's still probably possible, but needs some more thinking.

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