Python: Calcola la Tesselazione Voronoi dalla triangolazione di Delaunay di Scipy in 3D

Question

Ho circa 50.000 punti dati in 3D su cui ho eseguito scipy.spatial.Delaunay dal nuovo scipy (sto usando 0.10) che mi dà una triangolazione molto utile .

Sulla base: http://en.wikipedia.org/wiki/Delaunay_triangulation (sezione "Rapporti con il diagramma di Voronoi")

... Mi chiedevo se c'è un modo semplice per raggiungere il "grafo duale" di questo triangolazioni, che è il Voronoi tesselation.

Eventuali indizi? La mia ricerca su questo sembra non mostrare funzioni pre-costruite in scipy, che trovo quasi strano!

Grazie, Edward

Answer 1

Le informazioni di adiacenza sono disponibili nell'attributo neighbors dell'oggetto Delaunay. Sfortunatamente, il codice non espone i circumcenter all'utente, in questo momento, quindi dovrai ricalcolare te stesso.

Inoltre, i bordi di Voronoi che si estendono all'infinito non vengono ottenuti direttamente in questo modo. Probabilmente è ancora possibile, ma ha bisogno di pensare ancora.

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

Non so di una funzione per fare questo, ma non sembra un compito troppo complicato.

Il grafico di Voronoi è la giunzione dei circumcircles, come descritto nell'articolo di Wikipedia.

Quindi è possibile iniziare con una funzione che trova il centro dei circoncelli di un triangolo, che è la matematica di base (http://en.wikipedia.org/wiki/Circumscribed_circle).

Quindi, basta unire i centri dei triangoli adiacenti.

Answer 3

mi sono imbattuto lo stesso problema e costruirono una soluzione di pv. Di risposta e di altri frammenti di codice che ho trovato sul web. La soluzione restituisce un diagramma di Voronoi completo, incluse le linee esterne dove non sono presenti i vicini del triangolo.

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

linee nere = diagramma di Voronoi, Linee rosse = Delauny triangoli

Answer 4

Come ho trascorso una notevole quantità di tempo su questo, mi piacerebbe condividere la mia soluzione su come ottenere il Voronoi poligoni anziché solo i bordi.

Il codice è a https://gist.github.com/letmaik/8803860 e si estende sulla soluzione di tauran.

Innanzitutto, ho modificato il codice per fornirmi vertici e (coppie di) indici (= spigoli) separatamente, poiché molti calcoli possono essere semplificati quando si lavora su indici anziché coordinate punto.

Quindi, nel metodo voronoi_cell_lines, determino quali bordi appartengono a quali celle. Per questo io uso la soluzione proposta di Alink da una domanda correlata. Cioè, per ogni spigolo trova i due punti di input più vicini (= celle) e crea una mappatura da quella.

L'ultimo passaggio consiste nel creare i poligoni effettivi (vedere metodo voronoi_polygons).Innanzitutto, le celle esterne che hanno i bordi penzolanti devono essere chiuse. Questo è semplice come guardare attraverso tutti i bordi e verificare quali hanno solo un bordo vicino. Possono esserci zero o due di questi bordi. In caso di due, li collego introducendo un ulteriore margine.

Infine, i bordi non ordinati di ogni cella devono essere inseriti nell'ordine giusto per ricavarne un poligono.

L'utilizzo è:

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

cui uscite:

Il codice non è probabilmente adatto per un gran numero di punti di ingresso e può essere migliorata in alcune aree. Tuttavia, può essere utile per gli altri che hanno problemi simili.