'Centroidal Voronoi tesselations around obstacles
I am trying to implement a Centroidal Voronoi tessellation algorithm in a bounded rectangular space such that, the bounding rectangle has a number of obstacles(polygons) in it.
The below code gives the centroidal voronoi tessellations in a bounding box without the presence of obstacles(polygons). The blue points are the generators, red points the centroids and yellow points are the in between points between the blue and red points. 
import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys
np.random.seed(1)
eps = sys.float_info.epsilon
n_robots = 10
robots = np.random.rand(n_robots, 2)
#print(robots)
bounding_box = np.array([0., 1., 0., 1.])
def in_box(robots, bounding_box):
return np.logical_and(np.logical_and(bounding_box[0] <= robots[:, 0],
robots[:, 0] <= bounding_box[1]),
np.logical_and(bounding_box[2] <= robots[:, 1],
robots[:, 1] <= bounding_box[3]))
def voronoi(robots, bounding_box):
i = in_box(robots, bounding_box)
points_center = robots[i, :]
points_left = np.copy(points_center)
points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])
points_right = np.copy(points_center)
points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])
points_down = np.copy(points_center)
points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])
points_up = np.copy(points_center)
points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])
points = np.append(points_center,
np.append(np.append(points_left,
points_right,
axis=0),
np.append(points_down,
points_up,
axis=0),
axis=0),
axis=0)
# Compute Voronoi
vor = sp.spatial.Voronoi(points)
# Filter regions and select corresponding points
regions = []
points_to_filter = [] # we'll need to gather points too
ind = np.arange(points.shape[0])
ind = np.expand_dims(ind,axis= 1)
for i,region in enumerate(vor.regions): # enumerate the regions
if not region: # nicer to skip the empty region altogether
continue
flag = True
for index in region:
if index == -1:
flag = False
break
else:
x = vor.vertices[index, 0]
y = vor.vertices[index, 1]
if not(bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and
bounding_box[2] - eps <= y and y <= bounding_box[3] + eps):
flag = False
break
if flag:
regions.append(region)
# find the point which lies inside
points_to_filter.append(vor.points[vor.point_region == i][0,:])
vor.filtered_points = np.array(points_to_filter)
vor.filtered_regions = regions
return vor
def centroid_region(vertices):
A = 0
C_x = 0
C_y = 0
for i in range(0, len(vertices) - 1):
s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
A = A + s
C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
A = 0.5 * A
C_x = (1.0 / (6.0 * A)) * C_x
C_y = (1.0 / (6.0 * A)) * C_y
return np.array([[C_x, C_y]])
def plot(r,index):
vor = voronoi(r, bounding_box)
fig = pl.figure()
ax = fig.gca()
#ax.plot(pol2[:,0],pol2[:,1],'k-')
# Plot initial points
ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')
print("initial",vor.filtered_points)
# Plot ridges points
for region in vor.filtered_regions:
vertices = vor.vertices[region, :]
ax.plot(vertices[:, 0], vertices[:, 1], 'go')
# Plot ridges
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :]
ax.plot(vertices[:, 0], vertices[:, 1], 'k-')
# Compute and plot centroids
centroids = []
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :]
centroid = centroid_region(vertices)
centroids.append(list(centroid[0, :]))
ax.plot(centroid[:, 0], centroid[:, 1], 'r.')
centroids = np.asarray(centroids)
rob = np.copy(vor.filtered_points)
# the below code is for the plotting purpose the update happens in the update function
interim_x = np.asarray(centroids[:,0] - rob[:,0])
interim_y = np.asarray(centroids[:,1] - rob[:,1])
magn = [np.linalg.norm(centroids[i,:] - rob[i,:]) for i in range(rob.shape[0])]
x = np.copy(interim_x)
x = np.asarray([interim_x[i]/magn[i] for i in range(interim_x.shape[0])])
y = np.copy(interim_y)
y = np.asarray([interim_y[i]/magn[i] for i in range(interim_y.shape[0])])
nor = np.copy(rob)
for i in range(x.shape[0]):
nor[i,0] = x[i]
nor[i,1] = y[i]
temp = np.copy(rob)
temp[:,0] = [rob[i,0] + 0.5*interim_x[i] for i in range(rob.shape[0])]
temp[:,1] = [rob[i,1] + 0.5*interim_y[i] for i in range(rob.shape[0])]
pol = [[]]
ax.plot(temp[:,0] ,temp[:,1], 'y.' )
ax.set_xlim([-0.1, 1.1])
ax.set_ylim([-0.1, 1.1])
pl.savefig("voronoi" + str(index) + ".png")
return centroids
def update(rob,centroids):
interim_x = np.asarray(centroids[:,0] - rob[:,0])
interim_y = np.asarray(centroids[:,1] - rob[:,1])
magn = [np.linalg.norm(centroids[i,:] - rob[i,:]) for i in range(rob.shape[0])]
x = np.copy(interim_x)
x = np.asarray([interim_x[i]/magn[i] for i in range(interim_x.shape[0])])
y = np.copy(interim_y)
y = np.asarray([interim_y[i]/magn[i] for i in range(interim_y.shape[0])])
nor = [np.linalg.norm([x[i],y[i]]) for i in range(x.shape[0])]
temp = np.copy(rob)
temp[:,0] = [rob[i,0] + 0.5*interim_x[i] for i in range(rob.shape[0])]
temp[:,1] = [rob[i,1] + 0.5*interim_y[i] for i in range(rob.shape[0])]
return np.asarray(temp)
if __name__ == '__main__':
for i in range(1):
centroids = plot(robots,i)
robots = update(robots,centroids)
Now I want to extend this particular code to the case with obstacles i.e I want to do something like 
except I don't want anything in the red polygons.
One approach I tried was to partition the space using the free area which did not go well. 
.
Code for approach one :
import random
from shapely.geometry import Polygon, Point
import numpy as np
import matplotlib.pyplot as pl
def get_random_point_in_polygon(poly,polygons,num):
(minx, miny, maxx, maxy) = poly.bounds
points =[]
while num != 0:
p = Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
if any(poly.contains(p) for poly in polygons):
continue
else:
num = num-1
#print(num)
points.append([p.x,p.y])
return np.asarray(points)
def polysplit(poly,polygons):
(minx, miny, maxx, maxy) = poly.bounds
pols =[]
return pols
def randomRects(p,poly):
(minx, miny, maxx, maxy) = poly.bounds
rect = []
while True:
w = round(random.uniform(0, 1),3)
h = round(random.uniform(0, 1),3)
if (((p[:,0]+w) < maxx) and ((p[:,1]+h) < maxy)):
rect.append(np.squeeze([np.squeeze(p[:,0]),np.squeeze(p[:,1])]))
rect.append(np.squeeze([np.squeeze(p[:,0]+w),np.squeeze(p[:,1])]))
rect.append(np.squeeze([np.squeeze(p[:,0]+w),np.squeeze(p[:,1]+h)]))
rect.append(np.squeeze([np.squeeze(p[:,0]),np.squeeze(p[:,1]+h)]))
rect.append(np.squeeze([np.squeeze(p[:,0]),np.squeeze(p[:,1])]))
break
else:
continue
return np.asarray(rect)
def rect(poly,polygons):
rec =[]
area = poly.area
areas = 0
for i in polygons:
areas = areas+i.area
#print(area - areas)
flag = False
while (area - areas) > 0.4:
p = get_random_point_in_polygon(poly,polygons,1)
#print(p)
rect = randomRects(p,poly)
if any(poly.intersects(Polygon(rect)) for poly in polygons):
continue
#elif any(poly.intersects(Polygon(rect)) for poly in rec):
#continue
else:
if rec == []:
rec.append(Polygon(rect))
print("hi")
elif any(pol.intersects(Polygon(rect)) for pol in rec):
continue
else:
areas = areas+Polygon(rect).area
print(area-areas)
rec.append(Polygon(rect))
return rec
p = Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
p2 = Polygon([(0, 0), (.2,0), (.2,.2), (0, 0.2), (0,0)])
p3 = Polygon([(0.4, 0.4), (0.8,0.4), (.8,.8), (0.4, 0.8), (0.4,0.4)])
p4 = Polygon([(0.1,0.6),(0.3,.6),(0.3,0.9),(0.1,0.9),(0.1,0.6)])
p5 = Polygon([(0.25,0.25),(0.85,.25),(0.85,0.35),(0.25,0.35),(0.25,0.25)])
polygons = []
polygons.append(p2)
polygons.append(p3)
polygons.append(p4)
polygons.append(p5)
point_in_poly = get_random_point_in_polygon(p,polygons,10000)
fig = pl.figure()
ax = fig.gca()
#ax.plot(point_in_poly[:,0],point_in_poly[:,1],'b.')
area = 0
for po in polygons:
#area = area +po.area
x,y = po.exterior.xy
#print [x,y]
ax.plot(x,y,'r-')
#print(p.area - area)
r = rect(p,polygons)
for rr in r:
#area = area +po.area
x,y = rr.exterior.xy
#print [x,y]
ax.plot(x,y,'b-')
ax.set_xlim([-0.1, 1.1])
ax.set_ylim([-0.1, 1.1])
pl.savefig("test1.png")
Second approach i thought to use binary space partitioning to divide the free area into rectangles and apply the above code to each of these free area rectangles. But I am not sure how to do this in a python.
Third approach : I used the Python triangle library to calculate the conforming constrained delaunay triangulation of the free space and tried porting it back to a voronoi diagram. The results were not as expected.
and it's 
The code below is a compilation of all the methods i tried so it might be messy. I tried using the Voronoi function in Scipy,Triangle libraries and also a tried converting the triangulation to voronoi using a custom approach. The code doesn't work well and also has some errors .
from numpy import array
import numpy as np
def read_poly(file_name):
"""
Simple poly-file reader, that creates a python dictionary
with information about vertices, edges and holes.
It assumes that vertices have no attributes or boundary markers.
It assumes that edges have no boundary markers.
No regional attributes or area constraints are parsed.
"""
output = {'vertices': None, 'holes': None, 'segments': None}
# open file and store lines in a list
file = open(file_name, 'r')
lines = file.readlines()
file.close()
lines = [x.strip('\n').split() for x in lines]
# Store vertices
vertices= []
N_vertices, dimension, attr, bdry_markers = [int(x) for x in lines[0]]
# We assume attr = bdrt_markers = 0
for k in range(N_vertices):
label, x, y = [items for items in lines[k+1]]
vertices.append([float(x), float(y)])
output['vertices']=array(vertices)
# Store segments
segments = []
N_segments, bdry_markers = [int(x) for x in lines[N_vertices+1]]
for k in range(N_segments):
label, pointer_1, pointer_2 = [items for items in lines[N_vertices+k+2]]
segments.append([int(pointer_1)-1, int(pointer_2)-1])
output['segments'] = array(segments)
# Store holes
N_holes = int(lines[N_segments+N_vertices+2][0])
holes = []
for k in range(N_holes):
label, x, y = [items for items in lines[N_segments + N_vertices + 3 + k]]
holes.append([float(x), float(y)])
output['holes'] = array(holes)
print(holes)
return output
from triangle import triangulate,voronoi, plot as tplot
import matplotlib.pyplot as plt
image = read_poly("/home/pranav/catkin_ws/src/beginner_tutorials/scripts/test.poly")
cncfq20adt = triangulate(image, 'pq20a.01D')
#print(cncfq20adt['vertices'])
#print(cncfq20adt['triangles'])
plt.figure(figsize=(10, 10))
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax, **cncfq20adt)
plt.savefig("image.png")
import triangle
from scipy.spatial import Delaunay
pts = cncfq20adt['vertices']
tri = Delaunay(pts)
p = tri.points[tri.vertices]
#print(pts)
# Triangle vertices
A = p[:,0,:].T
B = p[:,1,:].T
C = p[:,2,:].T
print(C)
# See http://en.wikipedia.org/wiki/Circumscribed_circle#Circumscribed_circles_of_triangles
# The following is just a direct transcription of the formula there
a = A - C
b = B - C
def dot2(u, v):
return u[0]*v[0] + u[1]*v[1]
def cross2(u, v, w):
"""u x (v x w)"""
return dot2(u, w)*v - dot2(u, v)*w
def ncross2(u, v):
"""|| u x v ||^2"""
return sq2(u)*sq2(v) - dot2(u, v)**2
def sq2(u):
return dot2(u, u)
cc = cross2(sq2(a) * b - sq2(b) * a, a, b) / (2*ncross2(a, b)) + C
# Grab the Voronoi edges
vc = cc[:,tri.neighbors]
vc[:,tri.neighbors == -1] = np.nan # edges at infinity, plotting those would need more work...
lines = []
lines.extend(zip(cc.T, vc[:,:,0].T))
lines.extend(zip(cc.T, vc[:,:,1].T))
lines.extend(zip(cc.T, vc[:,:,2].T))
# Plot it
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
lines = LineCollection(lines, edgecolor='b')
#plt.hold(1)
plt.plot(pts[:,0], pts[:,1], '.')
plt.plot(cc[0], cc[1], '*')
plt.gca().add_collection(lines)
plt.axis('equal')
plt.xlim(-0.1, 1.1)
plt.ylim(-0.1, 1.1)
plt.savefig("vor2.png")
ax1 = plt.subplot(121, aspect='equal')
triangle.plot.plot(ax1, vertices=pts)
lim = ax1.axis()
points, edges, ray_origin, ray_direct = triangle.voronoi(pts)
d = dict(vertices=points, edges=edges, ray_origins=ray_origin, ray_directions=ray_direct)
ax2 = plt.subplot(111, aspect='equal')
triangle.plot.plot(ax2, **d)
ax2.axis(lim)
plt.savefig("vor.png")
import matplotlib.pyplot as pl
import scipy as sp
import scipy.spatial
import sys
from shapely.geometry import Polygon,Point
import random
np.random.seed(1)
eps = sys.float_info.epsilon
"""
n_robots = 50
#robots = np.random.rand(n_robots, 2)
def get_random_point_in_polygon(poly,polygons,num):
(minx, miny, maxx, maxy) = poly.bounds
points =[]
while num != 0:
p = Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
if any(poly.contains(p) for poly in polygons):
continue
else:
num = num-1
print(num)
points.append([p.x,p.y])
return np.asarray(points)
def polysplit(poly,polygons):
rectangles = []
return rectangles
p = Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
p2 = Polygon([(0, 0), (.2,0), (.2,.2), (0, 0.2), (0,0)])
p3 = Polygon([(0.4, 0.4), (0.7,0.4), (.7,.7), (0.4, 0.7), (0.4,0.4)])
polygons = []
polygons.append(p2)
polygons.append(p3)
#point_in_poly = get_random_point_in_polygon(p,polygons,10)
robots = get_random_point_in_polygon(p,polygons,n_robots)
#print(sampl)
print(robots)
bounding_box = np.array([0., 1, 0., 1])
box = np.array([0.2, 0.6, 0, 0.6])
box2 = np.array([0, 0.6, 0.2, 0.6])
boxes =[]
boxes.append(box)
boxes.append(box2)
"""
robots = cncfq20adt['vertices']
print("length",len(robots))
bounding_box = np.array([0., 1., 0., 1.])
def in_box(robots, bounding_box):
return np.logical_and(np.logical_and(bounding_box[0] <= robots[:, 0],
robots[:, 0] <= bounding_box[1]),
np.logical_and(bounding_box[2] <= robots[:, 1],
robots[:, 1] <= bounding_box[3]))
def voronoi(robots, bounding_box):
i = in_box(robots, bounding_box)
points_center = robots[i, :]
points_left = np.copy(points_center)
points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])
points_right = np.copy(points_center)
points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])
points_down = np.copy(points_center)
points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])
points_up = np.copy(points_center)
points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])
points = np.append(points_center,
np.append(np.append(points_left,
points_right,
axis=0),
np.append(points_down,
points_up,
axis=0),
axis=0),
axis=0)
# Compute Voronoi
vor = sp.spatial.Voronoi(points)
# Filter regions and select corresponding points
regions = []
points_to_filter = [] # we'll need to gather points too
ind = np.arange(points.shape[0])
ind = np.expand_dims(ind,axis= 1)
for i,region in enumerate(vor.regions): # enumerate the regions
if not region: # nicer to skip the empty region altogether
continue
flag = True
for index in region:
if index == -1:
flag = False
break
else:
x = vor.vertices[index, 0]
y = vor.vertices[index, 1]
if not(bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and
bounding_box[2] - eps <= y and y <= bounding_box[3] + eps):
flag = False
break
if flag:
regions.append(region)
# find the point which lies inside
points_to_filter.append(vor.points[vor.point_region == i][0,:])
vor.filtered_points = np.array(points_to_filter)
vor.filtered_regions = regions
return vor
def centroid_region(vertices):
A = 0
C_x = 0
C_y = 0
for i in range(0, len(vertices) - 1):
s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
A = A + s
C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
A = 0.5 * A
C_x = (1.0 / (6.0 * A)) * C_x
C_y = (1.0 / (6.0 * A)) * C_y
return np.array([[C_x, C_y]])
def plot(r,index):
vor = voronoi(r, bounding_box)
fig = pl.figure()
ax = fig.gca()
#ax.plot(pol2[:,0],pol2[:,1],'k-')
# Plot initial points
ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')
print("initial",vor.filtered_points)
# Plot ridges points
for region in vor.filtered_regions:
vertices = vor.vertices[region, :]
ax.plot(vertices[:, 0], vertices[:, 1], 'go')
# Plot ridges
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :]
ax.plot(vertices[:, 0], vertices[:, 1], 'k-')
# Compute and plot centroids
centroids = []
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :]
centroid = centroid_region(vertices)
centroids.append(list(centroid[0, :]))
ax.plot(centroid[:, 0], centroid[:, 1], 'r.')
centroids = np.asarray(centroids)
rob = np.copy(vor.filtered_points)
# the below code is for the plotting purpose the update happens in the update function
interim_x = np.asarray(centroids[:,0] - rob[:,0])
interim_y = np.asarray(centroids[:,1] - rob[:,1])
magn = [np.linalg.norm(centroids[i,:] - rob[i,:]) for i in range(rob.shape[0])]
x = np.copy(interim_x)
x = np.asarray([interim_x[i]/magn[i] for i in range(interim_x.shape[0])])
y = np.copy(interim_y)
y = np.asarray([interim_y[i]/magn[i] for i in range(interim_y.shape[0])])
nor = np.copy(rob)
for i in range(x.shape[0]):
nor[i,0] = x[i]
nor[i,1] = y[i]
temp = np.copy(rob)
temp[:,0] = [rob[i,0] + 0.5*interim_x[i] for i in range(rob.shape[0])]
temp[:,1] = [rob[i,1] + 0.5*interim_y[i] for i in range(rob.shape[0])]
pol = [[]]
ax.plot(temp[:,0] ,temp[:,1], 'y.' )
ax.set_xlim([-0.1, 1.1])
ax.set_ylim([-0.1, 1.1])
pl.savefig("voronoi" + str(index) + ".png")
return centroids
def update(rob,centroids):
interim_x = np.asarray(centroids[:,0] - rob[:,0])
interim_y = np.asarray(centroids[:,1] - rob[:,1])
magn = [np.linalg.norm(centroids[i,:] - rob[i,:]) for i in range(rob.shape[0])]
x = np.copy(interim_x)
x = np.asarray([interim_x[i]/magn[i] for i in range(interim_x.shape[0])])
y = np.copy(interim_y)
y = np.asarray([interim_y[i]/magn[i] for i in range(interim_y.shape[0])])
nor = [np.linalg.norm([x[i],y[i]]) for i in range(x.shape[0])]
temp = np.copy(rob)
temp[:,0] = [rob[i,0] + 0.5*interim_x[i] for i in range(rob.shape[0])]
temp[:,1] = [rob[i,1] + 0.5*interim_y[i] for i in range(rob.shape[0])]
return np.asarray(temp)
if __name__ == '__main__':
for i in range(1):
centroids = plot(robots,i)
robots = update(robots,centroids)
I will be really really grateful if someone can help me.
Solution 1:[1]
I was looking for a solution to this, and found that several very clever people had done most of the work for me! There is a useful package called geovoronoi, that performs voronoi calculations on constrained spaces via shapely polygons, and then these can be plotted using code from: https://sgillies.net/2010/04/06/painting-punctured-polygons-with-matplotlib.html.
I put together the following code, which should help
from geovoronoi import voronoi_regions_from_coords
import numpy as np
from shapely.geometry import MultiPolygon, Polygon, Point
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon as mPolygon
from shapely.geometry.polygon import LinearRing
from matplotlib.patches import PathPatch
from matplotlib.path import Path
def RingCoding(ob):
# The codes will be all "LINETO" commands, except for "MOVETO"s at the
# beginning of each subpath
n = len(ob.coords)
codes = np.ones(n, dtype=Path.code_type) * Path.LINETO
codes[0] = Path.MOVETO
return codes
def Pathify(polygon):
# Convert coordinates to path vertices. Objects produced by Shapely's
# analytic methods have the proper coordinate order, no need to sort.
vertices = np.concatenate(
[np.asarray(polygon.exterior)]
+ [np.asarray(r) for r in polygon.interiors])
codes = np.concatenate(
[RingCoding(polygon.exterior)]
+ [RingCoding(r) for r in polygon.interiors])
return Path(vertices, codes)
def CreatePatch(poly,area_override=None):
MAX_DENSITY = 0.75
area = poly.area
if area_override is not None:
area=area_override
density = 1 / area
color = (min(1, density / MAX_DENSITY), max(0, (MAX_DENSITY - density) / MAX_DENSITY), 0, 0.5)
region_external_coords = list(poly.exterior.coords)
if len(poly.interiors) > 0:
path = Pathify(poly)
patch = PathPatch(path, facecolor=color, edgecolor=color)
else:
patch = mPolygon(region_external_coords, True)
patch.set_color(color)
return patch
def main():
coords = np.array([[-29, 4], [-6, 3], [-1, -1], [-1.5, 0], [-9, -2],
[0, 0], [-1, 0], [3, 7], [3.2, 6.8],
[3.5, 7.2], [0.1, 2], [-3, 3],
[10, 10], [7, 15]
])
#DEFINE EXTERIOR POLYGONS HERE
a = [(-40, -4), (-40, 6), (2,6), (2, -4), (-40, -4)]
b = [(2, 6), (14, 6), (14, 19), (2, 19), (2, 6)]
#DEFINE INTERNAL HOLES HERE
hole_a_1 = LinearRing([(-20,-2), (-25,-2), (-25, 2),(-20, 2), (-20, -2)])
hole_a_2 = LinearRing([(-30, -2), (-35, -2), (-35, 2), (-30, 2), (-30, -2)])
hole_a_3 = LinearRing([(-20, 4.9), (-10, 4.9), (-10, 2), (-20, 4.9),])
shapely_poly = MultiPolygon([[a, [hole_a_1, hole_a_2, hole_a_3]], [b,[]]])
min_x, min_y = np.inf, np.inf
max_x, max_y = -np.inf, -np.inf
for poly in shapely_poly:
bounds=poly.bounds
min_x=min(bounds[0], min_x)
max_x=max(bounds[2], max_x)
min_y=min(bounds[1], min_y)
max_y=max(bounds[3], max_y)
fig, ax = plt.subplots()
ax.set_xlim(min_x-5, max_x+5)
ax.set_ylim(min_y-5, max_y+5)
# this creates a dictionary of polygons/multipolygons
# and a dictionary of lists, indicating which point is in those polygons
# (if there are identical points, those lists might have 2+ numbers in them)
region_polys, region_pts = voronoi_regions_from_coords(coords, shapely_poly)
for i in region_polys:
if type(region_polys[i]) is MultiPolygon:
# this means that the voronoi cell is technically a multipolygon.
# while you could argue whether this should ever occur, the current implementation
# does this.
# so, we should probably check which polygon actually contains the point.
point=region_pts[i][0]
temp_point=Point(coords[point])
for poly in region_polys[i]:
if poly.contains(temp_point):
#this is the one.
patch=CreatePatch(poly)
ax.add_patch(patch)
temp_area=poly.area
for poly in region_polys[i]:
if not poly.contains(temp_point):
patch=CreatePatch(poly, temp_area)
ax.add_patch(patch)
else:
patch=CreatePatch(region_polys[i])
ax.add_patch(patch)
points=list(zip(*coords))
plt.scatter(points[0], points[1])
plt.show()
if __name__=="__main__":
main()
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
| Solution | Source |
|---|---|
| Solution 1 |