¿Cómo encontrar todos los vecinos de un punto dado en una triangulación de delaunay usando scipy.spatial.Delaunay?

Question

He estado buscando una respuesta a esta pregunta pero no encuentro nada útil.

estoy trabajando con la pila de la computación científica pitón (scipy, numpy, matplotlib) y que tiene un conjunto de 2 puntos dimensionales, para lo cual calculamos el traingulation Delaunay (wiki) usando scipy.spatial.Delaunay.

I necesidad de escribir una función que, dado cualquier punto a, devolverá todos los otros puntos que son vértices de cualquier simplex (triángulo es decir) que a es también un vértice del (los vecinos de a en la triangulación). Sin embargo, la documentación para scipy.spatial.Delaunay (here) es bastante mala, y no puedo entender cómo se especifican las simplicidades o si voy a hacer esto. Incluso solo una explicación de cómo se organizan las matrices neighbors, vertices y vertex_to_simplex en la salida de Delaunay sería suficiente para ponerme en marcha.

Muchas gracias por cualquier ayuda.

Answer 1

Lo descubrí por mi cuenta, así que aquí hay una explicación para cualquier persona futura que esté confundida por esto.

A modo de ejemplo, vamos a usar la sencilla entramado de puntos que los que estaba trabajando en mi código, que genero la siguiente manera

import numpy as np 
import itertools as it 
from matplotlib import pyplot as plt 
import scipy as sp 

inputs = list(it.product([0,1,2],[0,1,2])) 
i = 0 
lattice = range(0,len(inputs)) 
for pair in inputs: 
    lattice[i] = mksite(pair[0], pair[1]) 
    i = i +1 

detalles aquí realmente no es importante, es suficiente decir que genera un triangular periódica celosía en el que la distancia entre un punto y cualquiera de sus seis vecinos más cercanos es 1.

Para representar que

plt.plot(*np.transpose(lattice), marker = 'o', ls = '') 
axes().set_aspect('equal') 

Ahora calcular la triangulación:

dela = sp.spatial.Delaunay 
triang = dela(lattice) 

Veamos lo que esto nos da.

triang.points 

salida:

array([[ 0.  , 0.  ], 
     [ 0.5  , 0.8660254 ], 
     [ 1.  , 1.73205081], 
     [ 1.  , 0.  ], 
     [ 1.5  , 0.8660254 ], 
     [ 2.  , 1.73205081], 
     [ 2.  , 0.  ], 
     [ 2.5  , 0.8660254 ], 
     [ 3.  , 1.73205081]]) 

simple, sólo un conjunto de los nueve puntos de la retícula se ilustra arriba. Como vamos a ver:

triang.vertices 

de salida:

array([[4, 3, 6], 
     [5, 4, 2], 
     [1, 3, 0], 
     [1, 4, 2], 
     [1, 4, 3], 
     [7, 4, 6], 
     [7, 5, 8], 
     [7, 5, 4]], dtype=int32) 

En esta matriz, cada fila representa una simple (triángulo) en la triangulación. Las tres entradas en cada fila son los índices de los vértices de ese símplex en la matriz de puntos que acabamos de ver. Así, por ejemplo la primera simplex en esta matriz, [4, 3, 6] se compone de los puntos:

[ 1.5  , 0.8660254 ] 
[ 1.  , 0.  ] 
[ 2.  , 0.  ] 

Es fácil ver esto dibujando la red en un trozo de papel, el etiquetado de cada punto de acuerdo a su índice, y luego trazando a través de cada fila en triang.vertices.

Esta es toda la información que necesitamos para escribir la función que especifiqué en mi pregunta. Parece que :

def find_neighbors(pindex, triang): 
    neighbors = list() 
    for simplex in triang.vertices: 
     if pindex in simplex: 
      neighbors.extend([simplex[i] for i in range(len(simplex)) if simplex[i] != pindex]) 
      ''' 
      this is a one liner for if a simplex contains the point we`re interested in, 
      extend the neighbors list by appending all the *other* point indices in the simplex 
      ''' 
    #now we just have to strip out all the dulicate indices and return the neighbors list: 
    return list(set(neighbors)) 

Y eso es todo! Estoy seguro de que la función anterior podría funcionar con alguna optimización, es solo lo que se me ocurrió en unos minutos. Si alguien tiene alguna sugerencia, no dude en publicarla. Espero que esto ayude a alguien en el futuro que está tan confundido acerca de esto como yo.

Answer 2

Aquí también es una versión simple de una línea propia respuesta de James Porter usando listas por comprensión:

find_neighbors = lambda x,triang: list(set(indx for simplex in triang.simplices if x in simplex for indx in simplex if indx !=x)) 

Answer 3

yo también necesitaba esto y me encontré con la siguiente respuesta. Resulta que si necesita los vecinos para todos puntos iniciales, es mucho más eficiente para producir un diccionario de vecinos de una sola vez (el ejemplo siguiente es para 2D):

def find_neighbors(tess, points): 

    neighbors = {} 
    for point in range(points.shape[0]): 
     neighbors[point] = [] 

    for simplex in tess.simplices: 
     neighbors[simplex[0]] += [simplex[1],simplex[2]] 
     neighbors[simplex[1]] += [simplex[2],simplex[0]] 
     neighbors[simplex[2]] += [simplex[0],simplex[1]] 

    return neighbors 

Los vecinos de punto v son entonces neighbors[v]. Para 10,000 puntos en esto se necesitan 370ms para correr en mi computadora portátil. ¿Tal vez otros tienen ideas sobre cómo optimizar esto aún más?

Answer 4

Los métodos descritos anteriormente recorren todas las simplicidades, lo que podría llevar mucho tiempo, en caso de que haya una gran cantidad de puntos. Una mejor manera podría ser usar Delaunay.vertex_neighbor_vertices, que ya contiene toda la información sobre los vecinos. Por desgracia, la extracción de la información

def find_neighbors(pindex, triang): 

    return triang.vertex_neighbor_vertices[1][triang.vertex_neighbor_vertices[0][pindex]:triang.vertex_neighbor_vertices[0][pindex+1]] 

El código siguiente muestra cómo obtener los índices de algunos vértices (número 17, en este ejemplo):

import scipy.spatial 
import numpy 
import pylab 

x_list = numpy.random.random(200) 
y_list = numpy.random.random(200) 

tri = scipy.spatial.Delaunay(numpy.array([[x,y] for x,y in zip(x_list, y_list)])) 

pindex = 17 

neighbor_indices = find_neighbors(pindex,tri) 

pylab.plot(x_list, y_list, 'b.') 
pylab.plot(x_list[pindex], y_list[pindex], 'dg') 
pylab.plot([x_list[i] for i in neighbor_indices], 
      [y_list[i] for i in neighbor_indices], 'ro')  

pylab.show() 

Answer 5

Aquí es una elaboracion de respuesta @astrofrog. Esto también funciona en más de 2D.

Tardó unos 300 ms en un conjunto de 2430 puntos en 3D (alrededor de 16000 simplices).

from collections import defaultdict 

def find_neighbors(tess): 
    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