Siéntete en la libertad de compartir nuestra página y códigos con tus amigos, ayúdanos a hacer crecer nuestra comunidad.
Solución:
Dado un cuadro delimitador rectangular, mi primera idea fue definir una especie de operación de intersección entre este cuadro delimitador y el diagrama de Voronoï producido por scipy.spatial.Voronoi. Una idea no necesariamente genial, ya que requiere codificar una gran cantidad de funciones básicas de geometría computacional.
Sin embargo, aquí está la segunda idea (¿truco?) que me vino a la mente: los algoritmos para calcular el diagrama de Voronoï de un conjunto de n puntos en el plano tienen una complejidad temporal de O(n ln(n)). ¿Qué pasa con la adición de puntos para restringir las celdas de Voronoï de los puntos iniciales para que se encuentren en el cuadro delimitador?
Solución para un diagrama de Voronoï acotado
Una imagen vale un gran discurso:
¿Qué hice aquí? ¡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]) por una simetría de reflexión según x = 0.0 (borde izquierdo del cuadro delimitador). Con simetrías de reflexión según 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 diagrama de Voronoï resultante (yo uso scipy.spatial.voronoi_plot_2d).
¿Qué hacer a continuación? Simplemente filtre puntos, bordes o caras de acuerdo con el cuadro delimitador. Y obtenga el centroide de cada cara de acuerdo con la conocida fórmula para calcular el centroide de un polígono. Aquí hay una imagen del resultado (los centroides están en rojo):
Un poco de diversión antes de mostrar el código
¡Genial! Parece funcionar. ¿Qué pasa 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)? ¿Qué pasa si lo intento una y otra vez?
Paso 2
Paso 10
Paso 25
¡Frio! Las células de 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")
Si conservas alguna duda y capacidad de ascender nuestro artículo puedes escribir una anotación y con deseo lo estudiaremos.