Obtención de coordenadas poligonales acotadas de celdas Voronoi

Tengo puntos (por ejemplo, lat, lon pares de ubicaciones de torres celulares) y necesito obtener el polígono de las celdas Voronoi que forman.

from scipy.spatial import Voronoi tower = [[ 24.686 , 46.7081], [ 24.686 , 46.7081], [ 24.686 , 46.7081]] c = Voronoi(towers) 

Ahora, necesito obtener los límites de polígono en lat, lon coordenadas para cada celda (y cuál fue el centroide que rodea este polígono). Necesito que este Voronoi también esté limitado. Lo que significa que los límites no van al infinito, sino más bien dentro de un cuadro delimitador.

Dado un cuadro de límite rectangular, mi primera idea fue definir un tipo de operación de intersección entre este cuadro de límite y el diagtwig de Voronoï producido por scipy.spatial.Voronoi . Una idea no necesariamente buena, ya que esto requiere codificar una gran cantidad de funciones básicas de geometría computacional.

Sin embargo, aquí está la segunda idea (¿hack?) Que se me ocurrió: los algoritmos para calcular el diagtwig Voronoï de un conjunto de n puntos en el plano tienen una complejidad temporal de O(n ln(n)) . ¿Qué hay de agregar puntos para restringir las celdas Voronoï de los puntos iniciales para que se encuentren en el cuadro delimitador?

Solución para un diagtwig de Voronoï acotado.

Una imagen vale un gran discurso:

introduzca la descripción de la imagen aquí

¿Qué hice aquí? ¡Eso es bastante simple! Los puntos iniciales (en azul) se encuentran en [0.0, 1.0] x [0.0, 1.0] . Luego obtengo los puntos (en azul) a la izquierda (es decir, [-1.0, 0.0] x [0.0, 1.0] ) mediante una simetría de reflexión de acuerdo con x = 0.0 (borde izquierdo del cuadro delimitador). Con simetrías de reflexión de acuerdo con x = 1.0 , y = 0.0 y y = 1.0 (otros bordes del cuadro delimitador), obtengo todos los puntos (en azul) que necesito para hacer el trabajo.

Entonces corro scipy.spatial.Voronoi . La imagen anterior muestra el diagtwig de Voronoï resultante (uso scipy.spatial.voronoi_plot_2d ).

¿Qué hacer a continuación? Solo filtre los puntos, bordes o caras según el cuadro delimitador. Y obtenga el centroide de cada cara de acuerdo con la conocida fórmula para calcular el centroide del polígono . Aquí hay una imagen del resultado (los centroides están en rojo):

introduzca la descripción de la imagen aquí

Un poco de diversión antes de mostrar el código

¡Genial! Parece funcionar. ¿Qué sucede si después de una iteración trato de volver a ejecutar el algoritmo en los centroides (en rojo) en lugar de los puntos iniciales (en azul)? ¿Y si lo bash una y otra vez?

Paso 2

introduzca la descripción de la imagen aquí

Paso 10

introduzca la descripción de la imagen aquí

Paso 25

introduzca la descripción de la imagen aquí

¡Guay! Las células Voronoï tienden a minimizar su energía

Aquí está el código

 import matplotlib.pyplot as pl import numpy as np import scipy as sp import scipy.spatial import sys eps = sys.float_info.epsilon n_towers = 100 towers = np.random.rand(n_towers, 2) bounding_box = np.array([0., 1., 0., 1.]) # [x_min, x_max, y_min, y_max] def in_box(towers, bounding_box): return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0], towers[:, 0] <= bounding_box[1]), np.logical_and(bounding_box[2] <= towers[:, 1], towers[:, 1] <= bounding_box[3])) def voronoi(towers, bounding_box): # Select towers inside the bounding box i = in_box(towers, bounding_box) # Mirror points points_center = towers[i, :] points_left = np.copy(points_center) points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0]) points_right = np.copy(points_center) points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0]) points_down = np.copy(points_center) points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2]) points_up = np.copy(points_center) points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1]) points = np.append(points_center, np.append(np.append(points_left, points_right, axis=0), np.append(points_down, points_up, axis=0), axis=0), axis=0) # Compute Voronoi vor = sp.spatial.Voronoi(points) # Filter regions regions = [] for region in vor.regions: flag = True for index in region: if index == -1: flag = False break else: x = vor.vertices[index, 0] y = vor.vertices[index, 1] if not(bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and bounding_box[2] - eps <= y and y <= bounding_box[3] + eps): flag = False break if region != [] and flag: regions.append(region) vor.filtered_points = points_center vor.filtered_regions = regions return vor def centroid_region(vertices): # Polygon's signed area A = 0 # Centroid's x C_x = 0 # Centroid's y C_y = 0 for i in range(0, len(vertices) - 1): s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1]) A = A + s C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s A = 0.5 * A C_x = (1.0 / (6.0 * A)) * C_x C_y = (1.0 / (6.0 * A)) * C_y return np.array([[C_x, C_y]]) vor = voronoi(towers, bounding_box) fig = pl.figure() ax = fig.gca() # Plot initial points ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.') # Plot ridges points for region in vor.filtered_regions: vertices = vor.vertices[region, :] ax.plot(vertices[:, 0], vertices[:, 1], 'go') # Plot ridges for region in vor.filtered_regions: vertices = vor.vertices[region + [region[0]], :] ax.plot(vertices[:, 0], vertices[:, 1], 'k-') # Compute and plot centroids centroids = [] for region in vor.filtered_regions: vertices = vor.vertices[region + [region[0]], :] centroid = centroid_region(vertices) centroids.append(list(centroid[0, :])) ax.plot(centroid[:, 0], centroid[:, 1], 'r.') ax.set_xlim([-0.1, 1.1]) ax.set_ylim([-0.1, 1.1]) pl.savefig("bounded_voronoi.png") sp.spatial.voronoi_plot_2d(vor) pl.savefig("voronoi.png")