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 != 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, t]) # 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, vs[:,0] - c) 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.1 max_x = vor.max_bound + 0.1 min_y = vor.min_bound - 0.1 max_y = vor.max_bound + 0.1 mins = np.tile((min_x, min_y), (vertices.shape, 1)) bounded_vertices = np.max((vertices, mins), axis=0) maxs = np.tile((max_x, max_y), (vertices.shape, 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.1, vor.max_bound + 0.1) plt.ylim(vor.min_bound - 0.1, vor.max_bound + 0.1) plt.savefig('voro.png') plt.show()

### 13cailla commented Oct 19, 2018

I might be wrong, but the line 53 can throw an error if p1 not in all_ridges. Adding a line `if p1 in all_ridges:` simply solves the problem I think

### Javier1314 commented Nov 28, 2020

Hi, In which list are the generated vertices stored? I managed to graph such points by modifying some attributes, but I would like to know the position (x, y) of each generated vertices.

### Sklavit commented Nov 29, 2020

@Javier1314 It is hard to recall and test now for me, but I guess that generated vertices (which are for infinity) is in the end of `new_vertices`.
So you should take `new_vertices[len(vor.vertices):]`. In this list tuples (or lists) are stored, so I suppose `point` will be for `x` and `point` will be for `y`.

### Vosatorp commented Mar 2, 2021

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