Python: Calcule la Tesselación de Voronoi de la Triangulación de Delaunay de Scipy en 3D

Question

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

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

... Me preguntaba si hay una manera más fácil de llegar a la "grafo dual" de esta triangulación, que es el de Voronoi Tesselación

¿Alguna pista? Mi búsqueda en este parece no mostrar pre-construido en funciones scipy, que me parece casi extraño!

Gracias, Edward

Answer 1

La información de adyacencia se puede encontrar en el atributo neighbors del objeto Delaunay. Desafortunadamente, el código no expone los circuncentros al usuario en este momento, por lo que deberá recalcularlos usted mismo.

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

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

Answer 2

no sé de una función para hacer esto, pero no parece ser una tarea excesivamente complicada.

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

Así que podría comenzar con una función que encuentre el centro de las circunferencias de un triángulo, que es matemáticas básicas (http://en.wikipedia.org/wiki/Circumscribed_circle).

Luego, solo une los centros de triángulos adyacentes.

Answer 3

me encontré con el mismo problema y construyeron una solución de pv. De respuesta y otros fragmentos de código que he encontrado en la web. La solución devuelve un diagrama de Voronoi completo, incluidas las líneas externas donde no están presentes los vecinos de triángulo.

#!/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 = diagrama de Voronoi, las líneas rojas = Delauny triángulos

Answer 4

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

El código está en https://gist.github.com/letmaik/8803860 y se extiende en 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 en índices en lugar de coordenadas de puntos.

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

El último paso es crear los polígonos reales (consulte el método voronoi_polygons).Primero, las celdas externas que tienen bordes colgantes necesitan estar cerradas. Esto es tan simple como mirar a través de todos los bordes y verificar cuáles tienen solo un borde contiguo. Puede haber cero o dos de esos bordes. En el caso de dos, los conecto introduciendo un borde adicional.

Finalmente, los bordes desordenados en cada celda deben ponerse 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 emite:

El código no es probablemente adecuada para un gran número de puntos de entrada y se puede mejorar en algunas áreas. Sin embargo, puede ser útil para otras personas que tienen problemas similares.