17

I have points (e.g., lat, lon pairs of cell tower locations) and I need to get the polygon of the Voronoi cells they form.

from scipy.spatial import Voronoi

tower = [[ 24.686 ,  46.7081],
       [ 24.686 ,  46.7081],
       [ 24.686 ,  46.7081]]

c = Voronoi(towers)

Now, I need to get the polygon boundaries in lat,lon coordinates for each cell (and what was the centroid this polygon surrounds). I need this Voronoi to be bounded as well. Meaning that the boundaries don't go to infinity, but rather within a bounding box.

amaatouq
  • 2,297
  • 5
  • 29
  • 50

2 Answers2

33

Given a rectangular bounding box, my first idea was to define a kind of intersection operation between this bounding box and the Voronoï diagram produce by scipy.spatial.Voronoi. An idea not necessarily great, since this requires to code a large number of basic functions of computational geometry.

However, here is the second idea (hack?) that came to my mind: the algorithms to compute the Voronoï diagram of a set of n points in the plane have a time complexity of O(n ln(n)). What about adding points to constraint the Voronoï cells of the initial points to lie in the bounding box?

Solution for a bounded Voronoï diagram

A picture is worth a great speech:

enter image description here

What I did here? That's pretty simple! The initial points (in blue) lie in [0.0, 1.0] x [0.0, 1.0]. Then I get the points (in blue) on the left (i.e. [-1.0, 0.0] x [0.0, 1.0]) by a reflection symmetry according to x = 0.0 (left edge of the bounding box). With reflection symmetries according to x = 1.0, y = 0.0 and y = 1.0 (other edges of the bounding box), I get all the points (in blue) I need to do the job.

Then I run scipy.spatial.Voronoi. The previous image depicts the resulting Voronoï diagram (I use scipy.spatial.voronoi_plot_2d).

What to do next? Just filter points, edges or faces according to the bounding box. And get the centroid of each face according to the well-known formula to compute centroid of polygon. Here is an image of the result (centroids are in red):

enter image description here

Some fun before showing code

Great! It seems to work. What if after one iteration I try to re-run the algorithm on the centroids (in red) rather than the initial points (in blue)? What if I try it again and again?

Step 2

enter image description here

Step 10

enter image description here

Step 25

enter image description here

Cool! Voronoï cells tend to minimize their energy...

Here is the code

import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys

eps = sys.float_info.epsilon

n_towers = 100
towers = np.random.rand(n_towers, 2)
bounding_box = np.array([0., 1., 0., 1.]) # [x_min, x_max, y_min, y_max]

def in_box(towers, bounding_box):
    return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0],
                                         towers[:, 0] <= bounding_box[1]),
                          np.logical_and(bounding_box[2] <= towers[:, 1],
                                         towers[:, 1] <= bounding_box[3]))


def voronoi(towers, bounding_box):
    # Select towers inside the bounding box
    i = in_box(towers, bounding_box)
    # Mirror points
    points_center = towers[i, :]
    points_left = np.copy(points_center)
    points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])
    points_right = np.copy(points_center)
    points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])
    points_down = np.copy(points_center)
    points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])
    points_up = np.copy(points_center)
    points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])
    points = np.append(points_center,
                       np.append(np.append(points_left,
                                           points_right,
                                           axis=0),
                                 np.append(points_down,
                                           points_up,
                                           axis=0),
                                 axis=0),
                       axis=0)
    # Compute Voronoi
    vor = sp.spatial.Voronoi(points)
    # Filter regions
    regions = []
    for region in vor.regions:
        flag = True
        for index in region:
            if index == -1:
                flag = False
                break
            else:
                x = vor.vertices[index, 0]
                y = vor.vertices[index, 1]
                if not(bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and
                       bounding_box[2] - eps <= y and y <= bounding_box[3] + eps):
                    flag = False
                    break
        if region != [] and flag:
            regions.append(region)
    vor.filtered_points = points_center
    vor.filtered_regions = regions
    return vor

def centroid_region(vertices):
    # Polygon's signed area
    A = 0
    # Centroid's x
    C_x = 0
    # Centroid's y
    C_y = 0
    for i in range(0, len(vertices) - 1):
        s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
        A = A + s
        C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
        C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
    A = 0.5 * A
    C_x = (1.0 / (6.0 * A)) * C_x
    C_y = (1.0 / (6.0 * A)) * C_y
    return np.array([[C_x, C_y]])

vor = voronoi(towers, bounding_box)

fig = pl.figure()
ax = fig.gca()
# Plot initial points
ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')
# Plot ridges points
for region in vor.filtered_regions:
    vertices = vor.vertices[region, :]
    ax.plot(vertices[:, 0], vertices[:, 1], 'go')
# Plot ridges
for region in vor.filtered_regions:
    vertices = vor.vertices[region + [region[0]], :]
    ax.plot(vertices[:, 0], vertices[:, 1], 'k-')
# Compute and plot centroids
centroids = []
for region in vor.filtered_regions:
    vertices = vor.vertices[region + [region[0]], :]
    centroid = centroid_region(vertices)
    centroids.append(list(centroid[0, :]))
    ax.plot(centroid[:, 0], centroid[:, 1], 'r.')

ax.set_xlim([-0.1, 1.1])
ax.set_ylim([-0.1, 1.1])
pl.savefig("bounded_voronoi.png")

sp.spatial.voronoi_plot_2d(vor)
pl.savefig("voronoi.png")
Flabetvibes
  • 3,066
  • 20
  • 32
  • why do you subtract epsilon when checking if points are within the bounding box? To make sure floating point rounding errors don't make the point appear outside the bounding box when it's in fact inside? – duhaime Sep 25 '18 at 14:32
  • Also, why do you multiple A by 6.0 when getting the centroid coordinates? Any help you can offer on these questions would be greatly appreciated! – duhaime Sep 25 '18 at 16:02
  • Ah, this answers my second question: https://en.wikipedia.org/wiki/Centroid#Of_a_polygon – duhaime Sep 25 '18 at 16:04
  • 1
    To answer your first question: yes, this substraction is present only to handle floating point rounding errors. If I remember well, I've started without but was not able to get the expected result. If you try the code without the substraction, tell me if it works properly. – Flabetvibes Sep 25 '18 at 16:17
  • 1
    Thanks very much for this code & insight! I have an improvement for your procedure though: by construction, the regions you want to keep belong to the first 1/5th of the points you gave to `scipy.spatial.Voronoi()`. This means you can use the `Voronoi.point_region` attribute that gives for each point, the index of the region it belongs to: `vor.filtered_regions = np.array(vor.regions)[vor.point_region[:vor.npoints//5]]`. This a) saves a bunch of manual checking, and b) guarantees the regions are returned in the same order as their matching points (which was crucial for my use-case :) ) – Energya May 28 '19 at 00:22
  • Thanks for this great answer @Flabetvibes! I've used it here https://stackoverflow.com/a/61893015/139144 to address an artifact with KDE plotting that was bothering me. – Gabriel May 19 '20 at 14:04
1

I was having a lot of trouble using scipy's voronoi function and creating CVD's, so these wonderful posts and comments helped immensely. As a programming novice, I tried to understand the code from Flabetvvibes answer and I'm going to share my interpretation of how it works along with Energya's and my own modifications. I have also posted my version of the code in its entirety at the bottom of this answer

import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys
import copy

eps = sys.float_info.epsilon

# Returns a new np.array of towers that within the bounding_box
def in_box(towers, bounding_box):
    return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0],
                                         towers[:, 0] <= bounding_box[1]),
                          np.logical_and(bounding_box[2] <= towers[:, 1],
                                         towers[:, 1] <= bounding_box[3]))

the in_box function uses numpy's logical_and method to return a boolean array that represents which coordinates from towers are in the bounding box.

# Generates a bounded vornoi diagram with finite regions in the bounding box
def bounded_voronoi(towers, bounding_box):
    # Select towers inside the bounding box
    i = in_box(towers, bounding_box)

    # Mirror points left, right, above, and under to provide finite regions for the
    # edge regions of the bounding box
    points_center = towers[i, :]

    points_left = np.copy(points_center)
    points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])

    points_right = np.copy(points_center)
    points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])

    points_down = np.copy(points_center)
    points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])

    points_up = np.copy(points_center)
    points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])

    points = np.append(points_center,
                       np.append(np.append(points_left,
                                           points_right,
                                           axis=0),
                                 np.append(points_down,
                                           points_up,
                                           axis=0),
                                 axis=0),
                       axis=0)

Flabetvvibes mirrors the points to allow the regions along the inside edges of the bounding box to be finite. Scipy's voronoi method returns -1 for vertices that are undefined, so mirroring the points allows all the regions inside the bounding box to be finite and all the infinite regions are in the mirrored areas outside of the bounding box that will be discarded later.

# Compute Voronoi
vor = sp.spatial.Voronoi(points)

# creates a new attibute for points that form the diagram within the region
vor.filtered_points = points_center 
# grabs the first fifth of the regions, which are the original regions
vor.filtered_regions = np.array(vor.regions)[vor.point_region[:vor.npoints//5]]

return vor

this last bit of the bounded_voronoi method calls scipy's voronoi function and adds new attributes for filtered points and and regions that are within the bounding box. Energya suggested removing Flabetvvibe's code that manually found all the finite regions within the bounding box with a one liner that gets the first fifth of the regions, which are the original inputs along with points making up the bounding box.

def generate_CVD(points, iterations, bounding_box):
    p = copy.copy(points)

    for i in range(iterations):
        vor = bounded_voronoi(p, bounding_box)
        centroids = []

        for region in vor.filtered_regions:
            # grabs vertices for the region and adds a duplicate
            # of the first one to the end
            vertices = vor.vertices[region + [region[0]], :] 
            centroid = centroid_region(vertices)
            centroids.append(list(centroid[0, :]))

        p = np.array(centroids)

    return bounded_voronoi(p, bounding_box)

I took Flabetvvibe's code that perform an iteration of loyd's algorithm and formed it into a method for ease of use. For each iteration the previous bounded_voronoi function is called and then the centroids are found for each cell and they become the new set of points for the next iteration. vertices = vor.vertices[region + [region[0]], :] simply grabs all the vertices for the current region and duplicates the first vertices to the end, so that the first and last vertices are the same for calculating centroids.

Thanks to Flabetvvibes and Energya. Your posts/answers taught me how to use scipy's voronoi method better than the its documentation. I have also posted the code as a single body underneath to any other looking for a copy/paste.

import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys
import copy

eps = sys.float_info.epsilon

# Returns a new np.array of towers that within the bounding_box
def in_box(towers, bounding_box):
    return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0],
                                         towers[:, 0] <= bounding_box[1]),
                          np.logical_and(bounding_box[2] <= towers[:, 1],
                                         towers[:, 1] <= bounding_box[3]))


# Generates a bounded vornoi diagram with finite regions
def bounded_voronoi(towers, bounding_box):
    # Select towers inside the bounding box
    i = in_box(towers, bounding_box)

    # Mirror points left, right, above, and under to provide finite regions for the edge regions of the bounding box
    points_center = towers[i, :]

    points_left = np.copy(points_center)
    points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])

    points_right = np.copy(points_center)
    points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])

    points_down = np.copy(points_center)
    points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])

    points_up = np.copy(points_center)
    points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])

    points = np.append(points_center,
                       np.append(np.append(points_left,
                                           points_right,
                                           axis=0),
                                 np.append(points_down,
                                           points_up,
                                           axis=0),
                                 axis=0),
                       axis=0)

    # Compute Voronoi
    vor = sp.spatial.Voronoi(points)

    vor.filtered_points = points_center # creates a new attibute for points that form the diagram within the region
    vor.filtered_regions = np.array(vor.regions)[vor.point_region[:vor.npoints//5]] # grabs the first fifth of the regions, which are the original regions

    return vor


# Finds the centroid of a region. First and last point should be the same.
def centroid_region(vertices):
    # Polygon's signed area
    A = 0
    # Centroid's x
    C_x = 0
    # Centroid's y
    C_y = 0
    for i in range(0, len(vertices) - 1):
        s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
        A = A + s
        C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
        C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
    A = 0.5 * A
    C_x = (1.0 / (6.0 * A)) * C_x
    C_y = (1.0 / (6.0 * A)) * C_y
    return np.array([[C_x, C_y]])


# Performs x iterations of loyd's algorithm to calculate a centroidal vornoi diagram
def generate_CVD(points, iterations, bounding_box):
    p = copy.copy(points)

    for i in range(iterations):
        vor = bounded_voronoi(p, bounding_box)
        centroids = []

        for region in vor.filtered_regions:
            vertices = vor.vertices[region + [region[0]], :] # grabs vertices for the region and adds a duplicate of the first one to the end
            centroid = centroid_region(vertices)
            centroids.append(list(centroid[0, :]))

        p = np.array(centroids)

    return bounded_voronoi(p, bounding_box)


# returns a pyplot of given voronoi data
def plot_vornoi_diagram(vor, bounding_box, show_figure):
    # Initializes pyplot stuff
    fig = pl.figure()
    ax = fig.gca()

    # Plot initial points
    ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')

    # Plot ridges points
    for region in vor.filtered_regions:
        vertices = vor.vertices[region, :]
        ax.plot(vertices[:, 0], vertices[:, 1], 'go')

    # Plot ridges
    for region in vor.filtered_regions:
        vertices = vor.vertices[region + [region[0]], :]
        ax.plot(vertices[:, 0], vertices[:, 1], 'k-')

    # stores references to numbers for setting axes limits
    margin_percent = .1
    width = bounding_box[1]-bounding_box[0]
    height = bounding_box[3]-bounding_box[2]

    ax.set_xlim([bounding_box[0]-width*margin_percent, bounding_box[1]+width*margin_percent])
    ax.set_ylim([bounding_box[2]-height*margin_percent, bounding_box[3]+height*margin_percent])

    if show_figure:
        pl.show()
    return fig
Nick Maclean
  • 11
  • 2