Python: Calcule el Voronoi Tesselation a partir de la Triangulación Delaunay de Scipy en 3D

Tengo unos 50,000 puntos de datos en 3D en los que he ejecutado scipy.spatial.Delaunay desde el nuevo scipy (estoy usando 0.10), lo que me da una triangulación muy útil.

Basado en: http://en.wikipedia.org/wiki/Delaunay_triangulation (sección “Relación con el diagtwig de Voronoi”)

… Me preguntaba si hay una manera fácil de llegar al “gráfico doble” de esta triangulación, que es la Voronoi Tesselation.

¿Alguna pista? Mi búsqueda alrededor de esto parece no mostrar funciones de scipy pre-construidas, ¡lo cual me parece casi extraño!

Gracias edward

La información de adyacencia se puede encontrar en el atributo de neighbors del objeto Delaunay. Desafortunadamente, el código no expone los circuncéntricos al usuario en este momento, por lo que tendrá que volver a calcularlos usted mismo.

Además, los bordes Voronoi que se extienden hasta el infinito no se obtienen directamente de esta manera. Todavía es posible, pero necesita un poco más de reflexión.

 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): """ux (vxw)""" return dot2(u, w)*v - dot2(u, v)*w def ncross2(u, v): """|| uxv ||^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() 

Encontré el mismo problema y construí una solución a partir de la respuesta de pv. Y otros fragmentos de código que encontré en la web. La solución devuelve un diagtwig completo de Voronoi, incluidas las líneas externas donde no están presentes los triangularjs vecinos.

 #!/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() 

Líneas negras = diagtwig de Voronoi, líneas rojas = triangularjs de Delauny Líneas negras = diagrama de Voronoi, líneas rojas = triángulos de Delauny

Como pasé una cantidad considerable de tiempo en esto, me gustaría compartir mi solución sobre cómo obtener los polígonos Voronoi en lugar de solo los bordes.

El código está en https://gist.github.com/letmaik/8803860 y se extiende sobre la solución de tauran .

Primero, cambié el código para darme vértices y (pares de) índices (= bordes) por separado, ya que muchos cálculos se pueden simplificar cuando se trabaja con índices en lugar de coordenadas de puntos.

Luego, en el método voronoi_cell_lines , determino qué bordes pertenecen a qué celdas. Para eso uso la solución propuesta de Alink a partir de una pregunta relacionada. Es decir, para cada borde, encuentre los dos puntos de entrada más cercanos (= celdas) y cree una asignación a partir de eso.

El último paso es crear los polígonos reales (consulte el método voronoi_polygons ). En primer lugar, las celdas exteriores que tienen bordes colgantes deben cerrarse. Esto es tan simple como mirar a través de todos los bordes y verificar cuáles tienen un solo borde adyacente. Puede haber cero o dos de estos bordes. En el caso de dos, los conecto introduciendo una ventaja adicional.

Finalmente, los bordes desordenados en cada celda deben colocarse en el orden correcto para derivar un polígono de ellos.

El uso es:

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

que produce:

Polígonos voronoi

El código probablemente no sea adecuado para grandes cantidades de puntos de entrada y puede mejorarse en algunas áreas. Sin embargo, puede ser útil para otros que tienen problemas similares.

No conozco una función para hacer esto, pero no parece una tarea demasiado complicada.

El gráfico de Voronoi es la unión de los circuncírculos, como se describe en el artículo de wikipedia.

Por lo tanto, podría comenzar con una función que encuentre el centro de los círculos circulares de un triángulo, que es matemática básica (http://en.wikipedia.org/wiki/Circumscribed_circle).

Entonces, simplemente une los centros de los triangularjs adyacentes.