scipy voronoi 3d - показаны не все точки гребня

Question

У меня были проблемы с использованием функции Вороного Сципи. Я следовал 2-му примеру, однако, когда я выполнил аналогичный пример в 3d, вычисляются не все ridge_points. Мои данные - это коробка из 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])

не имеют никаких точек гребня. Я бы предположил (как и в случае 2d), что углы будут иметь ребристые точки. Например, я бы предположил, что точки [6]=[0,0,0] будут иметь точки ребра с [1,0,0], [0,1,0] и [0,0,1]. Разве это невозможно вычислить со Сципи или я думал об этом неправильно?

Answer 1

Сципи использует Qhull для вычислений Делоне, Вороного и Конвексхалла. Данные, содержащиеся в ridge_points это то, что сообщает qvoronoi Fv хотя гребни не обязательно перечислены в том же порядке. (Как проверка: https://gist.github.com/pv/2f756ec83cdf242ce691)

Документация Qhull для Fv ( http://www.qhull.org/html/qh-optf.htm) упоминает предостережение, которое кажется здесь уместным:

Опция 'Fv' не отображает гребни, которые требуют более одной средней точки. Например, на диаграмме воронойых точек Вороного перечислены нулевые гребни (например, "rbox 10 s | qvoronoi Fv Qz"). Другими примерами являются диаграммы Вороного для прямоугольной сетки (например, '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"). Куб должен быть достаточно большим, чтобы ограничить все области Вороного исходного набора точек. Пожалуйста, сообщайте о любых других пропущенных случаях. Если вы можете официально описать эти случаи или написать код для их решения, отправьте электронное письмо по адресу [e-mail].

rbox 27 M1,0 Упомянутый в тексте набор точно такой же, как и в вашем примере (в другом порядке).

Как правило, Qhull имеет проблемы с геометрическими вырождениями, которые встречаются, например, в прямоугольных сетках. Обходной путь должен установить qhull_options="QJ" который говорит ему, чтобы добавить случайные возмущения к точкам данных, пока вырождения не разрешены. Обычно это генерирует диаграммы тесселяций / вороной с несколькими дополнительными симплексами / гребнями, но может решить проблемы такого типа.

Answer 2

У меня тоже была такая же проблема. Затем я использовал 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))