scipy voronoi 3d - not all ridge points are shown

Question

I have been having trouble using scipy's Voronoi function. I have followed the 2d example, however when I performed a similar example in 3d not all the ridge_points are computed. My data is a box of 27 points in a [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)

I noticed the points on the corners:

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

do not have any ridge points. I would assume (like the 2d case) that corners would have ridge points. For example I would assume points[6]=[0,0,0] would have ridge points with [1,0,0], [0,1,0], and [0,0,1]. Is this not possible to compute with scipy or have I been thinking about this wrong?

Answer 1

I also had the same problem. Then I used Delaunay to get all rigde points in 3D. As follow:

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 uses Qhull for the Delaunay/Voronoi/Convexhull computations. The data contained in ridge_points is what is reported by qvoronoi Fv although the ridges are not necessarily listed in the same order. (As a check: https://gist.github.com/pv/2f756ec83cdf242ce691)

The Qhull documentation for Fv (http://www.qhull.org/html/qh-optf.htm#Fv2) mentions a caveat that seems relevant here:

Option 'Fv' does not list ridges that require more than one midpoint. For example, the Voronoi diagram of cospherical points lists zero ridges (e.g., 'rbox 10 s | qvoronoi Fv Qz'). Other examples are the Voronoi diagrams of a rectangular mesh (e.g., 'rbox 27 M1,0 | qvoronoi Fv') or a point set with a rectangular corner (e.g., 'rbox P4,4,4 P4,2,4 P2,4,4 P4,4,2 10 | qvoronoi Fv'). Both cases miss unbounded rays at the corners. To determine these ridges, surround the points with a large cube (e.g., 'rbox 10 s c G2.0 | qvoronoi Fv Qz'). The cube needs to be large enough to bound all Voronoi regions of the original point set. Please report any other cases that are missed. If you can formally describe these cases or write code to handle them, please send email to qhull@qhull.org.

The rbox 27 M1,0 mentioned in the text is exactly the same set of points as in your example (in different order).

Typically, Qhull has problems in dealing with geometrical degeneracies, which occur for example in rectangular meshes. A generic workaround is to set qhull_options="QJ" which tells it to add random perturbations to the data points until degeneracies are resolved. This typically generates tesselations/voronoi diagrams with several additional simplices/ridges, but may resolve issues of this type.