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()

The simplest solution is to not update points, where the new position is outside the desired bounding box.

#!/usr/bin/env python
"""
Implementation of a constrained Lloyds algorithm.
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi


def lloyds(positions, origin=(0,0), scale=(1,1), total_iterations=3):
    positions = positions.copy()
    for _ in range(total_iterations):
        centroids = _get_voronoi_centroids(positions)
        is_valid = _is_within_bbox(centroids, origin, scale)
        positions[is_valid] = centroids[is_valid]
    return positions


def _get_voronoi_centroids(positions):
    voronoi = Voronoi(positions)
    centroids = np.zeros_like(positions)
    for ii, idx in enumerate(voronoi.point_region):
        region = [jj for jj in voronoi.regions[idx] if jj != -1]
        centroids[ii] = np.mean(voronoi.vertices[region], axis=0)
    return centroids


def _is_within_bbox(points, origin, scale):
    minima = np.array(origin)
    maxima = minima + np.array(scale)
    return np.all(np.logical_and(points >= minima, points <= maxima), axis=1)


if __name__ == '__main__':

    positions = np.random.rand(200, 2) - 0.5
    adjusted = lloyds(positions, origin=(-0.5, -0.5), scale=(1, 1))

    fig, axes = plt.subplots(1, 2, sharex=True, sharey=True)
    axes[0].plot(*positions.T, 'ko')
    axes[1].plot(*adjusted.T, 'ro')

    for ax in axes:
        ax.set_aspect('equal')

    fig.tight_layout()
    plt.show()