scipy voronoi 3d - 并非所有脊点都显示

Question

我在使用scipy的Voronoi函数时遇到了一些问题。我按照二维的例子进行了操作，但在三维的类似例子中，并不是所有的边界点都被计算出来。我的数据是一组27个点，分布在一个[0,2]x[0,2]x[0,2]的立方体里：

points = np.array([
    # bottom plane
    [0,2,0], [1,2,0], [2,2,0],
    [0,1,0], [1,1,0], [2,1,0],
    [0,0,0], [1,0,0], [2,0,0],
    # middle plane
    [0,2,1], [1,2,1], [2,2,1],
    [0,1,1], [1,1,1], [2,1,1],
    [0,0,1], [1,0,1], [2,0,1],
    # top plane
    [0,2,2], [1,2,2], [2,2,2],
    [0,1,2], [1,1,2], [2,1,2],
    [0,0,2], [1,0,2], [2,0,2]
    ])

vor = Voronoi(points)

print vor.ridge_points
# outputed
array([[ 4,  7],
       [ 4,  5],
       [ 4,  3],
       [ 4,  1],
       [ 4, 13],
       [ 3, 12],
       [ 7, 16],
       [15, 12],
       [15, 16],
       [ 9, 12],
       [ 9, 10],
       [ 1, 10],
       [12, 21],
       [12, 13],
       [23, 14],
       [23, 22], 
       [14, 17],
       [14, 11],
       [14,  5],
       [14, 13],
       [22, 19],
       [22, 21],
       [22, 13],
       [22, 25],
       [17, 16],
       [11, 10],
       [25, 16],
       [16, 13],
       [13, 10],
       [19, 10], dtype=int32)

我注意到角落上的点：

points[0] = array([0, 2, 0])
points[2] = array([2, 2, 0])
points[6] = array([0, 0, 0])
points[8] = array([2, 0, 0])
points[18] = array([0, 2, 2])
points[20] = array([2, 2, 2])
points[24] = array([0, 0, 2])
points[26] = array([2, 0, 2])

没有任何边界点。我本以为（就像二维的情况一样）角落的点应该会有边界点。例如，我认为points[6]=[0,0,0]应该会和[1,0,0]、[0,1,0]以及[0,0,1]有边界点。难道这在scipy中无法计算，还是我想错了？

Answer 1

我也遇到过同样的问题。然后我使用了Delaunay算法来获取3D中的所有边界点。具体代码如下：

def find_neighbors(tess):
"""
Parameters
----------
tess : Delaunay

Returns
-------
neighbors : neighbors in defaultdict type

"""

neighbors = defaultdict(set)

for simplex in tess.simplices:
    for idx in simplex:
        other = set(simplex)
        other.remove(idx)
        neighbors[idx] = neighbors[idx].union(other)
return neighbors

import scipy.spatial
from collections import defaultdict

x_list = np.random.random(8)
y_list = np.random.random(8)
z_list = np.random.random(8)

tri = scipy.spatial.Delaunay(np.array([[x,y,z] for x,y,z in zip(x_list, y_list, z_list)])) # create the Delaunay triangles
print(find_neighbors(tri))

Answer 2

Scipy使用了Qhull这个工具来进行Delaunay、Voronoi和凸包的计算。ridge_points里面的数据是由qvoronoi Fv提供的，不过这些山脊的顺序可能不一样。（你可以通过这个链接检查一下：https://gist.github.com/pv/2f756ec83cdf242ce691）

Qhull的文档中提到的Fv选项（http://www.qhull.org/html/qh-optf.htm#Fv2）有一个需要注意的地方：

选项'Fv'不会列出需要多个中点的山脊。例如，针对共球点的Voronoi图会列出零个山脊（比如说，'rbox 10 s | qvoronoi Fv Qz'）。其他例子包括矩形网格的Voronoi图（比如说，'rbox 27 M1,0 | qvoronoi Fv'）或者有矩形角的点集（比如说，'rbox P4,4,4 P4,2,4 P2,4,4 P4,4,2 10 | qvoronoi Fv'）。在这两种情况下，角落的无界射线会被遗漏。为了找出这些山脊，可以用一个大立方体包围这些点（比如说，'rbox 10 s c G2.0 | qvoronoi Fv Qz'）。这个立方体需要足够大，以包住原始点集的所有Voronoi区域。如果你发现了其他遗漏的情况，请反馈。如果你能正式描述这些情况或者写代码来处理它们，请发邮件到qhull@qhull.org。

文中提到的rbox 27 M1,0实际上和你的例子是同一组点，只是顺序不同。

通常情况下，Qhull在处理几何退化问题时会遇到困难，比如在矩形网格中就会出现这种情况。一个通用的解决办法是设置qhull_options="QJ"，这会告诉它对数据点添加随机扰动，直到解决退化问题。这样通常会生成带有多个额外单纯形/山脊的剖分/Voronoi图，但可能会解决这类问题。