Instantly share code, notes, and snippets.

# pv/colorized_voronoi.py

Last active August 7, 2023 16:39
Star You must be signed in to star a gist
Colorized Voronoi diagram with Scipy, in 2D, including infinite regions
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
 import numpy as np import matplotlib.pyplot as plt from scipy.spatial import Voronoi 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) print "--" print regions print "--" print vertices # colorize for region in regions: polygon = vertices[region] 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()

### liubenyuan commented Sep 19, 2015

and `edge_project` is a function projecting points back on the boundary of shape using numerical gradient

```def edge_project(pts, fd, h0=1.0):
"""
project points back on the boundary (where fd=0) using numerical gradient

note : you should specify h0 according to your actual mesh size
"""
deps = sqrt(np.finfo(float).eps)*h0
d = fd(pts)
dgradx = (fd(pts + [deps, 0]) - d) / deps
dgrady = (fd(pts + [0, deps]) - d) / deps

### geertn444 commented Dec 17, 2019 • edited

The procedure voronoi_finite_polygons_2d creates remote points, but it creates duplicate remote points.
For example:

35 [-0.60794158 -0.24345993]
36 [-3.41872139 -0.21721891]
37 [-3.41872139 -0.21721891]
38 [-0.99550093 0.83040914]

p36 and p37 are the same.

p37 belongs to polygon 8:
8 [37, 11, 1, 5, 7, 38]

p36 belongs to polygon 6:
6 [11, 36, 35, 9, 10]

Because they also share p11, in reality polygon 6 & 8 are adjacent. However, this point is missed and this influences the coloring !
Would it be possible to run the procedure with a parameter to "merge" common points together ?

### Magnuti commented Apr 15, 2022

I agree with @geertn444, the code generates duplicate points.