recortar un diagtwig voronoi python

Estoy calculando un diagtwig voronoi a partir de un conjunto de puntos como sigue:

from scipy.spatial import Voronoi import numpy as np np.random.seed(0) points = np.random.uniform(-0.5, 0.5, (100, 2)) // Compute Voronoi v = Voronoi(points) voronoi_plot_2d(v) plt.show() 

Esto crea una imagen de la siguiente manera:

Voronoi

Como se puede ver, esto es crear vértices que van hacia el infinito (líneas discontinuas) y también más allá del cuadro delimitador original para los puntos, que es:

  bbox = np.array([[-0.5, -0.5], [0.5, -0.5], [0.5, 0.5], [-0.5, 0.5]]) 

Lo que me gustaría hacer es recortar el diagtwig voronoi a este cuadro delimitador, es decir, proyectar los vértices fuera de límites e infinitos en las ubicaciones apropiadas en este cuadro delimitador. Por lo tanto, los vértices tendrían que reorganizarse y proyectarse de nuevo a los puntos de intersección adecuados desde el infinito o los vértices finitos, pero que están fuera de los límites de mi región de recorte.

Se puede hacer fácilmente con Shapely . Puede instalarlo desde Conda Forge: conda install shapely -c conda-forge

Código que necesita en github.gist , basado en la respuesta de @Gabriel y @pv. :

 # coding=utf-8 import numpy as np import matplotlib.pyplot as plt from scipy.spatial import Voronoi from shapely.geometry import Polygon def voronoi_finite_polygons_2d(vor, radius=None): """ Reconstruct infinite voronoi regions in a 2D diagram to finite regions. Parameters ---------- vor : Voronoi Input diagram radius : float, optional Distance to 'points at infinity'. Returns ------- regions : list of tuples Indices of vertices in each revised Voronoi regions. vertices : list of tuples Coordinates for revised Voronoi vertices. Same as coordinates of input vertices, with 'points at infinity' appended to the end. """ if vor.points.shape[1] != 2: raise ValueError("Requires 2D input") new_regions = [] new_vertices = vor.vertices.tolist() center = vor.points.mean(axis=0) if radius is None: radius = vor.points.ptp().max()*2 # Construct a map containing all ridges for a given point all_ridges = {} for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices): all_ridges.setdefault(p1, []).append((p2, v1, v2)) all_ridges.setdefault(p2, []).append((p1, v1, v2)) # Reconstruct infinite regions for p1, region in enumerate(vor.point_region): vertices = vor.regions[region] if all(v >= 0 for v in vertices): # finite region new_regions.append(vertices) continue # reconstruct a non-finite region ridges = all_ridges[p1] new_region = [v for v in vertices if v >= 0] for p2, v1, v2 in ridges: if v2 < 0: v1, v2 = v2, v1 if v1 >= 0: # finite ridge: already in the region continue # Compute the missing endpoint of an infinite ridge t = vor.points[p2] - vor.points[p1] # tangent t /= np.linalg.norm(t) n = np.array([-t[1], t[0]]) # normal midpoint = vor.points[[p1, p2]].mean(axis=0) direction = np.sign(np.dot(midpoint - center, n)) * n far_point = vor.vertices[v2] + direction * radius new_region.append(len(new_vertices)) new_vertices.append(far_point.tolist()) # sort region counterclockwise vs = np.asarray([new_vertices[v] for v in new_region]) c = vs.mean(axis=0) angles = np.arctan2(vs[:,1] - c[1], vs[:,0] - c[0]) new_region = np.array(new_region)[np.argsort(angles)] # finish new_regions.append(new_region.tolist()) return new_regions, np.asarray(new_vertices) # make up data points np.random.seed(1234) points = np.random.rand(15, 2) # compute Voronoi tesselation vor = Voronoi(points) # plot regions, vertices = voronoi_finite_polygons_2d(vor) min_x = vor.min_bound[0] - 0.1 max_x = vor.max_bound[0] + 0.1 min_y = vor.min_bound[1] - 0.1 max_y = vor.max_bound[1] + 0.1 mins = np.tile((min_x, min_y), (vertices.shape[0], 1)) bounded_vertices = np.max((vertices, mins), axis=0) maxs = np.tile((max_x, max_y), (vertices.shape[0], 1)) bounded_vertices = np.min((bounded_vertices, maxs), axis=0) box = Polygon([[min_x, min_y], [min_x, max_y], [max_x, max_y], [max_x, min_y]]) # colorize for region in regions: polygon = vertices[region] # Clipping polygon poly = Polygon(polygon) poly = poly.intersection(box) polygon = [p for p in poly.exterior.coords] plt.fill(*zip(*polygon), alpha=0.4) plt.plot(points[:, 0], points[:, 1], 'ko') plt.axis('equal') plt.xlim(vor.min_bound[0] - 0.1, vor.max_bound[0] + 0.1) plt.ylim(vor.min_bound[1] - 0.1, vor.max_bound[1] + 0.1) plt.savefig('voro.png') plt.show()