How do you create outlines to the polygons of a Voronoi tessellation in Python using matplotlib?

Question

I'm trying to create a voronoi tessellation from a trajectory.xyz file with 305 particles in it and multiple frames of data. I have used this script I found on GitHub here: https://gist.github.com/pv/8036995 and retooled some stuff to more accurately fit my goal of plotting a single frame's tessellation (code below). My problem is, however, that whenever I run the code, there's no border/outline for any of the polygons generated.

This makes the whole thing look kind of messy and muddled. I am fairly confident it has to do with the graphing of the tessellation, but I'm unsure of how to fix it. Even when I run the GitHub script exactly as it is with the random points, my polygons don't have outlines. Any ideas?

Code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi
from curved_analysis import read_xyz, read_nfo, project_frame

def voronoi(traj):
    voronoi_frames = []
    coords = np.array(read_xyz("traj0.xyz"))
    for k in range(coords.shape[0]):
        points = project_frame(coords[k])
        vor = Voronoi(points[:, :2], qhull_options=('Qz'))
        voronoi_frames.append(vor) 
    return points, voronoi_frames


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)

# make up data points


def voronoi_plot(vor, points):
    regions, vertices = voronoi_finite_polygons_2d(vor)


    for region in regions:
        polygon = vertices[region]
        plt.fill(*zip(*polygon), alpha=0.4)

    plt.plot(points[:,0], points[:,1], 'ko', marker=".", markersize=2)
    plt.axis('equal')
    plt.xlim(vor.min_bound[0] - 0.1, vor.max_bound[0] + 0.1)
    plt.ylim(vor.min_bound[1] - 0.1, vor.max_bound[1] + 0.1)

    plt.savefig('voro.png')
    plt.show()

if __name__ == "__main__":
    frame_num = 166
    points, vor = voronoi("traj0.xyz")
    voronoi_plot(vor[frame_num], points)

Answer 1

Replace the line

plt.fill(*zip(*polygon), alpha=0.4)

with:

plt.fill(*zip(*polygon), alpha=0.4, lw=1, ec='k')

lw sets the width of polygon edges, and ec their color.