Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Get point associated with Voronoi region (scipy.spatial.Voronoi)

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.

like image 938
Gabriel Avatar asked Aug 14 '15 23:08

Gabriel


People also ask

How do you find Voronoi vertices?

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.

How do you find the center of a Voronoi diagram?

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.

What is the Voronoi diagram for a set of three points?

 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.


2 Answers

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')
like image 199
aqueiros Avatar answered Oct 19 '22 23:10

aqueiros


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]
like image 30
Gabriel Avatar answered Oct 19 '22 23:10

Gabriel