Last active
October 9, 2017 12:19
-
-
Save manuelep/ead2fded9ebbd812ef8ed99230ba502d to your computer and use it in GitHub Desktop.
How to divide poolygon into voronoi sub-polygons
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
from shapely.geometry import MultiPoint, Point, Polygon | |
from scipy.spatial import Voronoi | |
import numpy as np | |
import matplotlib.pyplot as plt | |
""" | |
1. buffer di max 20 m intorno al poligono building | |
2. selezione dei punti contenuti nel buffer | |
3. suddivizione del poligono originario in aree di voronoi come qui di seguito | |
4. assegnazione per ogni poligono di voronoi delle caratteristiche del punto | |
contenuto | |
""" | |
def get_bigger(polygon): | |
return Polygon(polygon.exterior.coords) | |
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() | |
# 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) | |
def get_max_dim(poly): | |
cc = [i for i in poly.boundary.coords] | |
return max([max([i[0] for i in cc])-min([i[0] for i in cc]), max([i[1] for i in cc])-min([i[1] for i in cc])]) | |
points = [[1,1], [3,1.1], [2,.5]] | |
points = np.array(points) | |
vor = Voronoi(points) | |
rect = np.array([[0,0], [4,0], [4,1], [0,1]]) | |
pts = MultiPoint([Point(i) for i in rect]) | |
mask = pts.convex_hull | |
regions, vertices = voronoi_finite_polygons_2d(vor, 2*get_max_dim(mask)) | |
new_vertices = [] | |
for region in regions: | |
polygon = vertices[region] | |
shape = list(polygon.shape) | |
shape[0] += 1 | |
p = Polygon(np.append(polygon, polygon[0]).reshape(*shape)).intersection(mask) | |
poly = np.array(list(zip(p.boundary.coords.xy[0][:-1], p.boundary.coords.xy[1][:-1]))) | |
new_vertices.append(poly) | |
plt.fill(*zip(*poly), alpha=0.4) | |
plt.plot(points[:,0], points[:,1], 'ko') | |
plt.title("Clipped Voronois") | |
plt.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment