README.md

Lloyd Relaxation on a voronoi diagram in Python

Lloyd_relaxation.py

# MIT License

# Copyright (c) 2023 Dylan Missuwe

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

import numpy as np
from scipy.spatial import Voronoi
from matplotlib.animation import FuncAnimation
from matplotlib import pyplot as plt, animation
from matplotlib.collections import PathCollection
from matplotlib.path import Path
from shapely.geometry import Polygon
from shapely.ops import unary_union

# This method clips the voronoi polygons that cross the bounding box. 
# It then returns a new clipped polygon
def clip_polygon_to_box(polygon_vertices, clip_box):
    # Convert the bounding box to a polygon
    clip_polygon = Polygon(clip_box)

    # Create a Polygon object from the input polygon_vertices
    polygon = Polygon(polygon_vertices)

    # Perform the intersection/clipping using Shapely
    clipped_polygon = polygon.intersection(clip_polygon)

    # If the intersection resulted in multiple polygons, we take their union
    if clipped_polygon.geom_type == 'MultiPolygon':
        clipped_polygon = unary_union(clipped_polygon)

    # return coords of the newly created clipped polygon
    return list(clipped_polygon.exterior.coords)

# https://gist.github.com/pv/8036995
def voronoi_finite_polygons_2d(vor, radius=None):
    if vor.points.shape[1] != 2:
        raise ValueError("Requires 2D input")

    new_regions = []
    new_vertices = vor.vertices.tolist()

    center = vor.points.mean(axis=0)
    if radius is None:
        radius = vor.points.ptp().max() * 2

    # Construct a map containing all ridges for a given point
    all_ridges = {}
    for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
        all_ridges.setdefault(p1, []).append((p2, v1, v2))
        all_ridges.setdefault(p2, []).append((p1, v1, v2))

    # Reconstruct infinite regions
    for p1, region in enumerate(vor.point_region):
        vertices = vor.regions[region]

        if all(v >= 0 for v in vertices):
            # finite region
            new_regions.append(vertices)
            continue

        # reconstruct a non-finite region
        ridges = all_ridges[p1]
        new_region = [v for v in vertices if v >= 0]

        for p2, v1, v2 in ridges:
            if v2 < 0:
                v1, v2 = v2, v1
            if v1 >= 0:
                # finite ridge: already in the region
                continue

            # Compute the missing endpoint of an infinite ridge
            t = vor.points[p2] - vor.points[p1]  # tangent
            t /= np.linalg.norm(t)
            n = np.array([-t[1], t[0]])  # normal

            midpoint = vor.points[[p1, p2]].mean(axis=0)
            direction = np.sign(np.dot(midpoint - center, n)) * n
            far_point = vor.vertices[v2] + direction * radius

            new_region.append(len(new_vertices))
            new_vertices.append(far_point.tolist())

        # sort region counterclockwise
        vs = np.asarray([new_vertices[v] for v in new_region])
        c = vs.mean(axis=0)
        angles = np.arctan2(vs[:, 1] - c[1], vs[:, 0] - c[0])
        new_region = np.array(new_region)[np.argsort(angles)]

        # finish
        new_regions.append(new_region.tolist())

    return new_regions, np.asarray(new_vertices)

# compute the polygon centroid (thx ChatGPT)
def calculate_polygon_centroid(polygon_vertices):
    # calculate the centroid of a polygon given its vertices
    x = np.array([v[0] for v in polygon_vertices])
    y = np.array([v[1] for v in polygon_vertices])
    area = 0.5 * np.sum(x[:-1] * y[1:] - x[1:] * y[:-1])
    cx = np.sum((x[:-1] + x[1:]) * (x[:-1] * y[1:] - x[1:] * y[:-1])) / (6 * area)
    cy = np.sum((y[:-1] + y[1:]) * (x[:-1] * y[1:] - x[1:] * y[:-1])) / (6 * area)
    return [cx, cy]

# perform Lloyd iterations on the points
def animate_voronoi(i):
    global points

    # Generate Voronoi diagram
    vor = Voronoi(points)

    # convert the voronoi diagram to a list of polygons
    regions, vertices = voronoi_finite_polygons_2d(vor, None)

    distance_sum = 0
    polygons = []
    centroids = []

    # for every polygon in the voronoid diagram:
    for j, region in enumerate(regions):
        polygon_indices = region

        # clip the voronoi cell with the bounderies of the voronoi diagram 
        clipped_polygon = clip_polygon_to_box(vertices[polygon_indices], bounding_box)
        polygons.append(Path(clipped_polygon, closed=True))

        # find the geometric centroids of the voronoi cells (not the seed points)
        centroid = calculate_polygon_centroid(clipped_polygon)
        centroids.append(centroid)

        # calculate the distance between the seed points and the geometric centroids of the voronoi cells
        # if this distance is zero for all points, the lloyd relaxation has been fully converged
        distance = np.sqrt((points[j, 0] - centroid[0])**2 + (points[j, 1] - centroid[1])**2)
        distance_sum += distance

    # convert centroids to np array
    centroids = np.array(centroids, dtype=float)

    # calculate and print the distance sum
    print("Sum of the moving distance on iteration " + str(i) + " is: " + str(distance_sum))

    plt.cla()  # Clear the previous frame

    # Use PathCollection to draw all voronoi polygons at once
    collection = PathCollection(polygons, alpha=1, facecolors='none', edgecolors='black')
    ax.add_collection(collection)

    # draw the seed points in red
    plt.scatter(points[:, 0], points[:, 1], c='red', s=1, zorder=2)

    # draw the centroids in green
    plt.scatter(centroids[:, 0], centroids[:, 1], c='green', s=1, zorder=2)

    # assign the old geometric centroids to the new seed points to perform a Lloid iteration
    points = centroids

    plt.xlim(-1.5, 1.5)
    plt.ylim(-1.5, 1.5)
    plt.axis('equal')

# Number of Lloyd iterations
n_it = 10000

# Number of points
count = 200

# generate random points within a circular distribution
angle = (2 * np.pi) * np.random.random(count)
points = np.stack([np.cos(angle), np.sin(angle)], axis=1)
points *= np.sqrt(np.random.random(count)*0.1)[:, None]

# Define the bounding box vertices (clockwise order)
bounding_box = [(-1.5, -1.5), (-1.5, 1.5), (1.5, 1.5), (1.5, -1.5)]

# Create the figure and axis
fig, ax = plt.subplots()

# Create the animation
anim = animation.FuncAnimation(fig, animate_voronoi,
                               frames=300, interval=20, blit=True)

#anim.save('animation.gif', writer='imagemagick', fps=60)

# Display the animation
plt.show()