Вопрос по python, scipy – Python: вычисление тесселяции Вороного по триангуляции Делоне Сципи в 3D

15

У меня есть около 50 000 точек данных в 3D, для которых я запустил scipy.spatial. Delaunay из нового scipy (я использую 0.10), который дает мне очень полезную триангуляцию.

На основании:http://en.wikipedia.org/wiki/Delaunay_triangulation (раздел "Связь с диаграммой Вороного")

... Мне было интересно, есть ли простой способ добраться до "двойного графа"? этой триангуляции, которая является тесселяцией Вороного.

Есть какие-нибудь подсказки? Мой поиск по этому вопросу, кажется, не показывает никаких встроенных функций scipy, что я нахожу почти странным!

Спасибо, Эдвард

Ваш Ответ

4   ответа
17

neighbors атрибут объекта Делоне. К сожалению, в данный момент код не предоставляет пользователю окружности, поэтому вам придется пересчитать их самостоятельно.

Кроме того, ребра Вороного, которые простираются до бесконечности, не получаются напрямую таким образом. Это, вероятно, все еще возможно, но требует дополнительного размышления.

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()
Error: User Rate Limit Exceeded
Error: User Rate Limit Exceededncross2Error: User Rate Limit ExceededuError: User Rate Limit ExceededvError: User Rate Limit ExceededaError: User Rate Limit ExceededbError: User Rate Limit ExceededaError: User Rate Limit ExceededbError: User Rate Limit ExceededuError: User Rate Limit Exceededv?
Error: User Rate Limit Exceeded EdwardAndo
Error: User Rate Limit Exceeded
Error: User Rate Limit Exceeded
7

зователя pv. И других фрагментов кода, которые я нашел в Интернете. Решение возвращает полную диаграмму Вороного, включая внешние линии, где нет соседей треугольника.

#!/usr/bin/env python
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

def voronoi(P):
    delauny = Delaunay(P)
    triangles = delauny.points[delauny.vertices]

    lines = []

    # Triangle vertices
    A = triangles[:, 0]
    B = triangles[:, 1]
    C = triangles[:, 2]
    lines.extend(zip(A, B))
    lines.extend(zip(B, C))
    lines.extend(zip(C, A))
    lines = matplotlib.collections.LineCollection(lines, color='r')
    plt.gca().add_collection(lines)

    circum_centers = np.array([triangle_csc(tri) for tri in triangles])

    segments = []
    for i, triangle in enumerate(triangles):
        circum_center = circum_centers[i]
        for j, neighbor in enumerate(delauny.neighbors[i]):
            if neighbor != -1:
                segments.append((circum_center, circum_centers[neighbor]))
            else:
                ps = triangle[(j+1)%3] - triangle[(j-1)%3]
                ps = np.array((ps[1], -ps[0]))

                middle = (triangle[(j+1)%3] + triangle[(j-1)%3]) * 0.5
                di = middle - triangle[j]

                ps /= np.linalg.norm(ps)
                di /= np.linalg.norm(di)

                if np.dot(di, ps) < 0.0:
                    ps *= -1000.0
                else:
                    ps *= 1000.0
                segments.append((circum_center, circum_center + ps))
    return segments

def triangle_csc(pts):
    rows, cols = pts.shape

    A = np.bmat([[2 * np.dot(pts, pts.T), np.ones((rows, 1))],
                 [np.ones((1, rows)), np.zeros((1, 1))]])

    b = np.hstack((np.sum(pts * pts, axis=1), np.ones((1))))
    x = np.linalg.solve(A,b)
    bary_coords = x[:-1]
    return np.sum(pts * np.tile(bary_coords.reshape((pts.shape[0], 1)), (1, pts.shape[1])), axis=0)

if __name__ == '__main__':
    P = np.random.random((300,2))

    X,Y = P[:,0],P[:,1]

    fig = plt.figure(figsize=(4.5,4.5))
    axes = plt.subplot(1,1,1)

    plt.scatter(X, Y, marker='.')
    plt.axis([-0.05,1.05,-0.05,1.05])

    segments = voronoi(P)
    lines = matplotlib.collections.LineCollection(segments, color='k')
    axes.add_collection(lines)
    plt.axis([-0.05,1.05,-0.05,1.05])
    plt.show()

Черные линии = диаграмма Вороного, красные линии = треугольники Делоне Black lines = Voronoi diagram, Red lines = Delauny triangles

Error: User Rate Limit ExceededdiError: User Rate Limit ExceededpsError: User Rate Limit Exceeded
7

я хотел бы поделиться своим решением о том, как получить Вороной.polygons а не только края.

Код находится наhttps://gist.github.com/letmaik/8803860 и распространяется на решениеТоран.

Во-первых, я изменил код для предоставления мне вершин и (пар) индексов (= ребер) по отдельности, так как многие вычисления могут быть упрощены при работе с индексами вместо координат точек.

Затем вvoronoi_cell_lines Метод, который я определяю, какие ребра принадлежат к каким ячейкам. Для этого я использую предложенное решениеСсылка из смежного вопроса. То есть для каждого ребра найдите две ближайшие точки ввода (= ячейки) и создайте отображение из этого.

Последний шаг - создать фактические полигоны (см.voronoi_polygons метод). Во-первых, внешние ячейки, у которых есть свисающие края, должны быть закрыты. Это так же просто, как просмотреть все ребра и проверить, какие из них имеют только одно соседнее ребро. Таких ребер может быть либо ноль, либо два. В случае двух я тогда соединяю их, вводя дополнительное ребро.

Наконец, неупорядоченные ребра в каждой ячейке должны быть расположены в правильном порядке, чтобы вывести из них многоугольник.

Использование:

P = np.random.random((100,2))

fig = plt.figure(figsize=(4.5,4.5))
axes = plt.subplot(1,1,1)

plt.axis([-0.05,1.05,-0.05,1.05])

vertices, lineIndices = voronoi(P)        
cells = voronoi_cell_lines(P, vertices, lineIndices)
polys = voronoi_polygons(cells)

for pIdx, polyIndices in polys.items():
    poly = vertices[np.asarray(polyIndices)]
    p = matplotlib.patches.Polygon(poly, facecolor=np.random.rand(3,1))
    axes.add_patch(p)

X,Y = P[:,0],P[:,1]
plt.scatter(X, Y, marker='.', zorder=2)

plt.axis([-0.05,1.05,-0.05,1.05])
plt.show()

какие выводы:

Voronoi polygons

Код, вероятно, не подходит для большого количества точек ввода и может быть улучшен в некоторых областях. Тем не менее, это может быть полезно для других, у которых есть подобные проблемы.

Error: User Rate Limit Exceeded
Error: User Rate Limit Exceeded
0

но это не кажется слишком сложной задачей.

Граф Вороного - это соединение окружностей, как описано в статье в Википедии.

Таким образом, вы можете начать с функции, которая находит центр окружностей треугольника, который является основной математикой (http://en.wikipedia.org/wiki/Circumscribeed_circle).

Тогда просто соединяйте центры соседних треугольников.

Error: User Rate Limit Exceeded

Похожие вопросы