scipy voronoi 3d - 并未显示所有山脊点答案

【问题标题】：scipy voronoi 3d - not all ridge points are shownscipy voronoi 3d - 并未显示所有山脊点
【发布时间】：2014-09-09 21:53:57
【问题描述】：

我在使用 scipy 的 Voronoi 函数时遇到了问题。我遵循了 2d 示例，但是当我在 3d 中执行类似示例时，并没有计算所有 ridge_points。我的数据是 [0,2]x[0,2]x[0,2] 中的 27 个点的框：

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

没有任何脊点。我会假设（如 2d 案例）角落会有脊点。例如，我假设 points[6]=[0,0,0] 的脊点为 [1,0,0]、[0,1,0] 和 [0,0,1]。这不可能用 scipy 计算还是我一直在想这个错误？

【问题讨论】：

FWIW，Scipy 输出与qvoronoi Fv < pts 一致，请参阅gist.github.com/pv/2f756ec83cdf242ce691 - 但是，我不完全理解为什么 Qhull 报告的 Voronoi 山脊是这样的。

标签： python numpy scipy voronoi

【解决方案1】：

Scipy 使用 Qhull 进行 Delaunay/Voronoi/Convexhull 计算。 ridge_points 中包含的数据是 qvoronoi Fv 报告的数据，尽管这些脊不一定以相同的顺序列出。（作为支票：https://gist.github.com/pv/2f756ec83cdf242ce691）

Fv (http://www.qhull.org/html/qh-optf.htm#Fv2) 的 Qhull 文档提到了一个似乎与此处相关的警告：

选项“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 图，但可能会解决此类问题。

【讨论】：

【解决方案2】：

我也遇到了同样的问题。然后我使用 Delaunay 来获得 3D 中的所有 rigde 点。如下：

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

【讨论】：