Instantly share code, notes, and snippets.

# Sklavit/colorized_voronoi_with_clipping.py

Forked from pv/colorized_voronoi.py
Created March 26, 2017 00:36
Star You must be signed in to star a gist
Colorized Voronoi diagram with Scipy, in 2D, including infinite regions which are clipped to given box.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 # 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()

### Vosatorp commented Mar 2, 2021

What is poly.exterior. coords? My code don't work because <<AttributeError: 'GeometryCollection' object has no attribute 'exterior'>>