I am trying to render the results of a GA approach to a covering problem. Earlier, I found this question extremely helpful in generating boundaries. Some edge cases do exist:
I seem to need (often/always? I'm not positive) 4 or more points to generate the zoning (seems like 2 should be the minimum, its a minor detail though). I also, every once in a while, fail to cover the surface with my polygons. Here's an example:
with this output
...
>79, new best f(8) = 1532561.736
>84, new best f(8) = 1532586.256
>97, new best f(8) = 1532910.361
Score: 1532910.361458032
points: [(7.385603823101929, 3.8405732260590417), (7.385603823101929, 1.0891597118471825), (6.97653903235321, 6.537437011319787), (6.911445327302822, 7.404222303328822)]
polys: [[(10.0, 5.6165848996388315), (10.0, 3.7685057514544313), (3.394616991360788, 2.46486646895311), (0.0, 2.46486646895311), (0.0, 4.0997682439946175), (10.0, 5.6165848996388315)], [(10.0, 2.46486646895311), (10.0, 0.0), (0.0, 0.0), (0.0, 2.46486646895311), (10.0, 2.46486646895311)], [(0.0, 4.0997682439946175), (0.0, 6.4493508286073835), (10.0, 7.200329250259113), (10.0, 5.6165848996388315), (0.0, 4.0997682439946175)], [(0.0, 6.4493508286073835), (0.0, 10.0), (10.0, 10.0), (10.0, 7.200329250259113), (0.0, 6.4493508286073835)]]
The entire square should be covered. How can I correct this issue?
The code I am using to generate it is:
from scipy.spatial import Voronoi #, voronoi_plot_2d
from shapely.geometry import Polygon
from voronoi_polygons import voronoi_polygons
def render_score(score):
print('Score:', score)
def render_zones(points):
print("points:", points)
vor = Voronoi(points)
#fig = voronoi_plot_2d(vor)
boundary = np.array([[0,0], [0,10], [10,10], [10,0]])
boundary_polygon = Polygon(boundary)
x, y = boundary.T
plt.xlim(round(x.min() - 1), round(x.max() + 1))
plt.ylim(round(y.min() - 1), round(y.max() + 1))
plt.plot(*np.array(points).T, 'b.')
diameter = np.linalg.norm(boundary.ptp(axis=0))
polys = []
for p in voronoi_polygons(vor, diameter):
polys.append(list(p.intersection(boundary_polygon).exterior.coords))
x, y = zip(*p.intersection(boundary_polygon).exterior.coords)
plt.plot(x, y, 'r-')
plt.show()
print("polys:", polys)
def render_plan(centers):
pass
def render_results(agent):
render_score(agent.evaluate())
zone_centers = agent.predict()
render_zones(zone_centers)
render_plan(zone_centers)
with voronoi_polygons.py:
# src https://stackoverflow.com/questions/23901943/voronoi-compute-exact-boundaries-of-every-region#answer-52727406
from collections import defaultdict
from shapely.geometry import Polygon
import numpy as np
def voronoi_polygons(voronoi, diameter):
"""
Generate shapely.geometry.Polygon objects corresponding to the
regions of a scipy.spatial.Voronoi object, in the order of the
input points. The polygons for the infinite regions are large
enough that all points within a distance 'diameter' of a Voronoi
vertex are contained in one of the infinite polygons.
"""
centroid = voronoi.points.mean(axis=0)
# Mapping from (input point index, Voronoi point index) to list of
# unit vectors in the directions of the infinite ridges starting
# at the Voronoi point and neighbouring the input point.
ridge_direction = defaultdict(list)
for (p, q), rv in zip(voronoi.ridge_points, voronoi.ridge_vertices):
u, v = sorted(rv)
if u == -1:
# Infinite ridge starting at ridge point with index v,
# equidistant from input points with indexes p and q.
t = voronoi.points[q] - voronoi.points[p] # tangent
n = np.array([-t[1], t[0]]) / np.linalg.norm(t) # normal
midpoint = voronoi.points[[p, q]].mean(axis=0)
direction = np.sign(np.dot(midpoint - centroid, n)) * n
ridge_direction[p, v].append(direction)
ridge_direction[q, v].append(direction)
for i, r in enumerate(voronoi.point_region):
region = voronoi.regions[r]
if -1 not in region:
# Finite region.
yield Polygon(voronoi.vertices[region])
continue
# Infinite region.
inf = region.index(-1) # Index of vertex at infinity.
j = region[(inf - 1) % len(region)] # Index of previous vertex.
k = region[(inf + 1) % len(region)] # Index of next vertex.
if j == k:
# Region has one Voronoi vertex with two ridges.
dir_j, dir_k = ridge_direction[i, j]
else:
# Region has two Voronoi vertices, each with one ridge.
dir_j, = ridge_direction[i, j]
dir_k, = ridge_direction[i, k]
# Length of ridges needed for the extra edge to lie at least
# 'diameter' away from all Voronoi vertices.
length = 2 * diameter / np.linalg.norm(dir_j + dir_k)
# Polygon consists of finite part plus an extra edge.
finite_part = voronoi.vertices[region[inf + 1:] + region[:inf]]
extra_edge = [voronoi.vertices[j] + dir_j * length,
voronoi.vertices[k] + dir_k * length]
yield Polygon(np.concatenate((finite_part, extra_edge)))
In fact, this may be a problem with scipy.spatial.Voronoi itself, as the built-in 2d plotter produces a strange looking result (it is missing a solid-or-dotted line):
edit:
I notice that I can render it correctly by using different bounds for determining the dotted lines ("diameter" property) than are used for the bounds of the shape and the render. Ie, this (here) works:
points = agent.predict()
vor = Voronoi(points)
determine_boundary = np.array([[-1000,-1000], [-1000,1010], [1010,1010], [1010,-1000]])
render_boundary = np.array([[0,0], [0,10], [10,10], [10,0]])
boundary_polygon = Polygon(render_boundary)
x, y = render_boundary.T
plt.xlim(round(x.min() - 1), round(x.max() + 1))
plt.ylim(round(y.min() - 1), round(y.max() + 1))
plt.plot(*np.array(points).T, 'b.')
diameter = np.linalg.norm(determine_boundary.ptp(axis=0))
polys = []
for p in voronoi_polygons(vor, diameter):
polys.append(list(p.intersection(boundary_polygon).exterior.coords))
x, y = zip(*p.intersection(boundary_polygon).exterior.coords)
plt.plot(x, y, 'r-')
plt.show()
print(polys)
Because the final shapes (the polys) maintain the original boundary, the sizes and positioning of the polys remain normal, this is just what I want. Now, I just want a way to know how large the diameter polygon really has to be to produce an error free coordinate set within the granularity of the floats that compose it.
My instinct was the the angles would be different the wider you go out, but this is not correct. In fact, this produces identical results with only slightly larger bounds: -2.5,-2.5,12.5,12.5.
Two details that could help me here are that the number of coordinates always shrinks, and the sum of the areas of the polys equals the sum of the manifold. But you cant just take measurements of the area and if less, increase the diameter a bit, measure again, then use the delta to calculate the limit; that would work if there was only one surface that was not covered, but there could be multiple faults.
