Keeping SciPy Voronoi Iteration Bounded

I’m looking at various means of balancing points in a given space, and voronoi iteration / relaxation (aka Lloyds algorithm) has been beckoning me.

However, when using SciPy Voronoi, the points seem to be leaking out of the boundary, and spreading into the known universe, which doesn’t work for me at all!

The following generates ten points in the region [-0.5,0.5] but after 100 generations they have spread to a region -800k …600k (depending on start conditions.
What I am hoping for is an evenly spaced set of points within [-0.5,0.5].

import numpy as np
from scipy.spatial import Voronoi, voronoi_plot_2d
from random import random


if __name__ == '__main__':

    points = [[random()-0.5, random()-0.5] for pt in range(100)]
    voronoi = None

    for i in range(60):
        voronoi = Voronoi(np.array(points))
        points = []

        for reg in voronoi.point_region:
            if len(voronoi.regions[reg]) > 0:
                points.append(np.mean([voronoi.vertices[i] for i in voronoi.regions[reg] if i >= 0], axis=0))

    plt = voronoi_plot_2d(voronoi, show_vertices=False, line_colors="orange", line_width=2, line_alpha=0.6, point_size=2)
    plt.show()

You need to sign in to view this answers