I'm generating a simple 2D Voronoi tessellation, using the scipy.spatial.Voronoi function. I use a random 2D distribution of points (see MCVE below).
I need a way to go through each defined region (defined by scipy.spatial.Voronoi
) and get the coordinates of the point associated to it (ie: the point that said region encloses).
The issue is that there are N+1
regions (polygons) defined for the N
points, and I'm not sure what this means.
Here's a MCVE that will fail when it gets to the last region:
from scipy.spatial import Voronoi
import numpy as np
# Generate random data.
N = 10
x = [np.random.random() for i in xrange(N)]
y = [np.random.random() for i in xrange(N)]
points = zip(x, y)
# Obtain Voronoi regions.
vor = Voronoi(points)
# Loop through each defined region/polygon
for i, reg in enumerate(vor.regions):
print 'Region:', i
print 'Indices of vertices of Voronoi region:', reg
print 'Associated point:', points[i], '\n'
Another thing I don't understand is why are there empty vor.regions
stored? According to the docs:
regions: Indices of the Voronoi vertices forming each Voronoi region. -1 indicates vertex outside the Voronoi diagram.
What does an empty region mean?
Add
I tried the point_region
attribute but apparently I don't understand how it works. It returns indexes outside of the range of the points
list. For example: in the MCVE above it will always show an index 10
for a list of 10 points, which is clearly out of range.
The Voronoi diagram is just the dual graph of the Delaunay triangulation. So, the edges of the Voronoi diagram are along the perpendicular bisectors of the edges of the Delaunay triangulation, so compute those lines. Then, compute the vertices of the Voronoi diagram by finding the intersections of adjacent edges.
First determine all the edges for a particular cell then then for that cell take the x component average separate then the y component average separate then that linear combination is the 'center of mass' of a cell. Then just do the same for each cell in the Voronoi diagram. This method works for triangles.
The set with three or more nearest neighbors make up the vertices of the diagram. The points �� are called the sites of the Voronoi diagram. The three bisectors intersect at a point The intersection can be outside the triangle. The point of intersection is center of the circle passing through the three points.
For your first question:
The issue is that there are N+1 regions (polygons) defined for the N points, and I'm not sure what this means.
This is because your vor.regions will always have an empty array. Something like
[[],[0, 0],[0, 1],[1, 1]]
This is related to your second question:
Another thing I don't understand is why are there empty vor.regions stored? According to the docs: regions: Indices of the Voronoi vertices forming each Voronoi region. -1 indicates vertex outside the Voronoi diagram. What does an empty region mean?
By default Voronoi() uses QHull with options 'Qbb Qc Qz Qx' enabled (qhull.org/html/qvoronoi.htm). This inserts a "point-at-infinity" which is used to improve precision on circular inputs. Therefore, being a "fake" point, it has no region. If you want to get rid of this, try removing the Qz option:
vor = Voronoi(points, qhull_options='Qbb Qc Qx')
I was misreading the docs. It says:
point_region: Index of the Voronoi region for each input point.
and I was using point_region
it as if it were the: "Index of the input point for each Voronoi region".
Instead of using:
points[i]
the correct point coordinates for each region can be obtained with:
np.where(vor.point_region == i)[0][0]
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With