pythoncomputational-geometrypolygonsvoronoi

Getting a bounded polygon coordinates from Voronoi cells


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.


Solution

  • 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")