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

Question

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.

Explicable

Voronoi instances have a point_region attribute that does exactly what you are after. Read the docs! – Funicular

@Funicular I did try that attribute, please see what I added to the question. – Explicable

Answer 1

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')

Zante

Answer 2

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]

Explicable

Did you figure out why certain regions are empty? This documentation says "Note here that due to similar numerical precision issues as in Delaunay triangulation above, there may be fewer Voronoi regions than input points." So I'm guessing that means qhull calculates Voronoi using Delaunay? And "such degeneracies can occur not only because of duplicated points, but also for more complicated geometrical reasons, even in point sets that at first sight seem well-behaved" – Appointor

But although I've seen how to hack around this problem, I'm still confused whether those solutions give you real answers or just keep "regions" from being empty. I will keep looking. Perhaps for my case I will just throw out the points that give empty regions – Appointor

Unfortunately it's been a while since I used Voronoi regions in my research and I ended up not actually using them, so I can not give you any insight into what's going on. Perhaps you could ask over at math.stackexchange.com since it would appear to be an issue related to the properties of Voronoi regions rather than its implementation in Scipy? – Explicable

Thanks for your quick answer. I believe it's actually the opposite, haha. Voronoi should always have a region associated with a given point so long as that point is not a duplicate. It seems to have something to do with qhull's implementation – Appointor

I hope @Zante ' answer is correct, and these empty regions just belong to the points added at infinity – Appointor

Answer 3

Here is a solution:

import numpy as np
from scipy.spatial import Voronoi
import matplotlib.pyplot as plt
from plotutils_moje import voronoi_plot_2d


class VoronoiRegion:
    def __init__(self, region_id):
        self.id = region_id
        self.vertices = []
        self.is_inf = False
        self.point_inside = None

    def __str__(self):
        text = f'region id={self.id}'
        if self.point_inside:
            point_idx, point = self.point_inside
            text = f'{text}[point:{point}(point_id:{point_idx})]'
        text += ', vertices: '
        if self.is_inf:
            text += '(inf)'
        for v in self.vertices:
            text += f'{v}'
        return text

    def __repr__(self):
        return str(self)

    def add_vertex(self, vertex, vertices):
        if vertex == -1:
            self.is_inf = True
        else:
            point = vertices[vertex]
            self.vertices.append(point)

def voronoi_to_voronoi_regions(voronoi):
    voronoi_regions = []

    for i, point_region in enumerate(voronoi.point_region):
        region = voronoi.regions[point_region]
        vr = VoronoiRegion(point_region)
        for r in region:
            vr.add_vertex(r, voronoi.vertices)
        vr.point_inside = (i, voronoi.points[i])
        voronoi_regions.append(vr)
    return voronoi_regions


points = np.array([[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1], [2, 2]])
vor = Voronoi(points)
regions = voronoi_to_voronoi_regions(vor)
for r in regions:
    print(r)

and result for sample:

region id=1[point:[0. 0.](point_id:0)], vertices: (inf)[0.5 0.5]
region id=3[point:[0. 1.](point_id:1)], vertices: (inf)[0.5 1.5][0.5 0.5]
region id=2[point:[0. 2.](point_id:2)], vertices: (inf)[0.5 1.5]
region id=8[point:[1. 0.](point_id:3)], vertices: (inf)[1.5 0.5][0.5 0.5]
region id=7[point:[1. 1.](point_id:4)], vertices: [0.5 0.5][0.5 1.5][1.5 1.5][1.5 0.5]
region id=9[point:[1. 2.](point_id:5)], vertices: (inf)[1.5 1.5][0.5 1.5]
region id=6[point:[2. 0.](point_id:6)], vertices: (inf)[1.5 0.5]
region id=4[point:[2. 1.](point_id:7)], vertices: (inf)[1.5 1.5][1.5 0.5]
region id=5[point:[2. 2.](point_id:8)], vertices: (inf)[1.5 1.5]

Streamliner