Calculate and draw bounded (closed) Voronoi diagrams in Python

Introduction

The libraries for calculating and drawing Voronoi diagrams in Python include scipy.spatial.Voronoi and scipy.spatial.voronoi_plot_2d. For details on these, refer to the following sites.

• Official documentation for scipy.spatial.Voronoi
• Official documentation for scipy.spatial.voronoi_plot_2d
• scipy Voronoi

However, since the Voronoi diagrams created by these libraries are calculated in an unbounded (unbounded) space, it is not possible to obtain information on the outer Voronoi side or Voronoi region (polygon). Therefore, this time, I will introduce a method of dividing a closed area into Voronoi diagrams.

Implementation

For the part that makes the Voronoi region bounded, I referred to this article. The general flow is as follows.

1. Add three dummy mother points to make all Voronoi regions bounded.
2. Calculate the Voronoi diagram with scipy.spatial.Voronoi.
3. Calculate the intersection of the area to be divided and each Voronoi area with shapely. 4.3 Draw the polygon obtained in 3.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import PolyCollection
from scipy.spatial import Voronoi, voronoi_plot_2d
from shapely.geometry import Polygon

def bounded_voronoi(bnd, pnts):
    """
A function that calculates and draws bounded Voronoi diagrams.
    """
    
    #Added 3 dummy mother points to make the Voronoi region of all mother points bounded
    gn_pnts = np.concatenate([pnts, np.array([[100, 100], [100, -100], [-100, 0]])])
    
    #Voronoi diagram calculation
    vor = Voronoi(gn_pnts)
    
    #Polygon for the area to be divided
    bnd_poly = Polygon(bnd)
    
    #List to store each Voronoi area
    vor_polys = []
    
    #Repeat for non-dummy mother points
    for i in range(len(gn_pnts) - 3):
        
        #Voronoi region that does not consider closed space
        vor_poly = [vor.vertices[v] for v in vor.regions[vor.point_region[i]]]
        #Calculate the intersection of the Voronoi area for the area to be divided
        i_cell = bnd_poly.intersection(Polygon(vor_poly))
        
        #Stores the vertex coordinates of the Voronoi region considering the closed space
        vor_polys.append(list(i_cell.exterior.coords[:-1]))
    
    
    #Drawing a Voronoi diagram
    fig = plt.figure(figsize=(7, 6))
    ax = fig.add_subplot(111)
    
    #Mother point
    ax.scatter(pnts[:,0], pnts[:,1])
    
    #Voronoi area
    poly_vor = PolyCollection(vor_polys, edgecolor="black",
                              facecolors="None", linewidth = 1.0)
    ax.add_collection(poly_vor)
    
    xmin = np.min(bnd[:,0])
    xmax = np.max(bnd[:,0])
    ymin = np.min(bnd[:,1])
    ymax = np.max(bnd[:,1])
    
    ax.set_xlim(xmin-0.1, xmax+0.1)
    ax.set_ylim(ymin-0.1, ymax+0.1)
    ax.set_aspect('equal')
    
    plt.show()
    
    return vor_polys

Divide the unit square

Dividing the unit square using the above function gives the following.

#Voronoi division area
bnd = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])

#Number of mother points
n = 30
#Mother point coordinates
pnts = np.random.rand(n, 2)

#Calculation and drawing of Voronoi diagram
vor_polys = bounded_voronoi(bnd, pnts)

Divide a general convex polygon

A general convex polygon that is not a square can be divided in the same way.

from scipy.spatial import ConvexHull

def points_in_convex_polygon(bnd, n):
    """
A function that randomly generates n points inside a convex polygon.
    """
    
    #Creating a matrix that represents the boundaries of a region
    bndhull = ConvexHull(bnd)
    bndTmp = bndhull.equations
    bndMat = np.matrix(bndTmp)
    Abnd = np.array(bndMat[:,0:2])
    bbnd = np.array(bndMat[:,2])
    
    #Rectangle surrounding the area
    xmin = np.min(bnd[:,0])
    xmax = np.max(bnd[:,0])
    ymin = np.min(bnd[:,1])
    ymax = np.max(bnd[:,1])
    
    #For repetition
    i = 0
    pnts = []
    
    while i < n:
        
        #Generate points
        pnt = np.random.rand(2)
        pnt[0] = xmin + (xmax - xmin) * pnt[0]
        pnt[1] = ymin + (ymax - ymin) * pnt[1]
        
        #If the point is inside a convex polygon
        if (np.round(np.dot(Abnd,pnt.transpose()),7) <= np.round(-bbnd.transpose(),7)).all():

            pnts.append(pnt.tolist())
            i += 1

    return np.array(pnts)


#Voronoi division area
bnd = np.array([[0.1,0.4],[0.3,0.2],[0.8,0.3],[0.9,0.6],[0.7,0.7],[0.4,0.7],[0.2,0.6]])

#Number of mother points
n = 10
#Mother point coordinates
pnts = points_in_convex_polygon(bnd, n)

#Calculation and drawing of Voronoi diagram
vor_polys = bounded_voronoi(bnd, pnts)