I'm trying to compute the exact boundaries of every region of a Voronoi Diagram using scipy.spatial.Voronoi, in the case when all the points are inside a pre-defined polygon.
For example, using the example in the documentation,
http://docs.scipy.org/doc/scipy-dev/reference/generated/scipy.spatial.Voronoi.html what if I need to compute Voroni with the same points but inside a rectangle with the following boundariesglobal_boundaries = np.array([[-2, -2], [4, -2], [4, 4], [-2, 4], [-2, -2]])
and I need to compute the precise boundaries of every voronoi region, like that?
voronoi_region_1_boundaries = [[-2, -2], [0.5, -2], [0.5, 0.5], [-2, 0-5], [-2, -2]]
voronoi_region_2_boundaries = [[-2, 1.5], [0.5, 1.5], [0.5, 4], [-2, 4], [-2, 1.5]]
voronoi_region_3_boundaries = [[-2, 0.5], [0.5, 0.5], [0.5, 1.5], [-2, 1.5], [-2, 0.5]]
and so on for all the 9 regions, instead of
vor.regions
[[], [-1, 0], [-1, 1], [1, -1, 0], [3, -1, 2], [-1, 3], [-1, 2], [3, 2, 0, 1], [2, -1, 0], [3, -1, 1]]
How do I compute the missing endpoint of an infinite ridge?
I've tried to adapt this code http://nbviewer.ipython.org/gist/pv/8037100
related to this problem Colorize Voronoi Diagram
but it's working only for rounded boundaries. I've modified it considering a radius such that my area is completely inside the circle, and then computing the intersection between the line connecting the points and the circumference and the boundaries. It works, but only for the first point and after that I have "GEOMETRYCOLLECTION EMPTY" as a result.
direction = np.sign(np.dot(midpoint - center, n)) * n
super_far_point = vor.vertices[v2] + direction * radius
line_0 = LineString([midpoint, super_far_point])
for i in range(0, len(map_boundaries)-1):
i += 1
line_i = LineString([(map_boundaries[i-1]), (map_boundaries[i])])
if line_0.intersection(line_i) != 0:
far_point = line_0.intersection(line_i)
new_region.append(len(new_vertices))
new_vertices.append(far_point.tolist())
Has anyone ever solved a similar problem?
Can anyone help?
I took the
voronoi_plot_2d
and modified it. see below.Note: There are two simple constraints:
colors:
1. Algorithm
I suggest the following two-step approach:
First, make a convex polygon for each of the Voronoi regions. In the case of the infinite regions, do this by splitting the point at infinity into two points that are far enough away, joined by an edge. ("Far enough away" means that the extra edge passes entirely outside the bounding polygon.)
Intersect each of the polygons from step (1) with the bounding polygon, using shapely's
intersection
method.The benefit of this approach over Ophir Cami's answer is that it works with non-convex bounding polygons, and the code is a little simpler.
2. Example
Let's start with the Voronoi diagram for the points from Ophir Cami's answer. The infinite ridges are shown as dashed lines by
scipy.spatial.voronoi_plot_2d
:Then for each Voronoi region we construct a convex polygon. This is easy for the finite regions, but we have to zoom out a long way to see what happens to the infinite Voronoi regions. The polygons corresponding to these regions have an extra edge that is far enough away that it lies entirely outside the bounding polygon:
Now we can intersect the polygons for each Voronoi region with the bounding polygon:
In this case, all the Voronoi polygons have non-empty intersection with the bounding polygon, but in the general case some of them might vanish.
3. Code
The first step is to generate polygons corresponding to the Voronoi regions. Like Ophir Cami, I derived this from the implementation of
scipy.spatial.voronoi_plot_2d
.The second step is to intersect the Voronoi polygons with the bounding polygon. We also need to choose an appropriate diameter to pass to
voronoi_polygons
.This plots the last of the figures in §2 above.