How to best create a circulair voronoi with python and matplotlib | R8

How to best create a voronoi pattern inside a circle with Python? I now have this

import numpy as np
from scipy.spatial import Voronoi, voronoi_plot_2d
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
import Rhino.Geometry as rg
import scriptcontext as sc

def poisson_disc_samples(radius, k=30, width=1, height=1):
    """
    Generate Poisson disc samples with the given radius and k (limit of samples to choose before rejection).
    """
    def in_circle(point):
        return np.linalg.norm(point - np.array([0.5, 0.5])) <= 0.5
    
    cell_size = radius / np.sqrt(2)
    grid_width = int(np.ceil(width / cell_size))
    grid_height = int(np.ceil(height / cell_size))
    grid = [None] * (grid_width * grid_height)
    points = []
    spawns = []

    def grid_coords(p):
        return int(p[0] / cell_size), int(p[1] / cell_size)
    
    def fits(p):
        return 0 <= p[0] < width and 0 <= p[1] < height and in_circle(p)
    
    p0 = np.array([np.random.uniform(), np.random.uniform()]) * np.array([width, height])
    if not in_circle(p0 + np.array([0.5, 0.5])):
        return points

    spawns.append(p0)
    grid_x, grid_y = grid_coords(p0)
    grid[grid_x + grid_y * grid_width] = p0
    points.append(p0)

    while spawns:
        idx = np.random.randint(0, len(spawns))
        base = spawns[idx]
        for _ in range(k):
            angle = np.random.uniform(0, 2 * np.pi)
            dist = np.random.uniform(radius, 2 * radius)
            new_point = base + dist * np.array([np.cos(angle), np.sin(angle)])
            new_grid_x, new_grid_y = grid_coords(new_point)
            if fits(new_point) and new_grid_x + new_grid_y * grid_width < len(grid) and grid[new_grid_x + new_grid_y * grid_width] is None:
                grid[new_grid_x + new_grid_y * grid_width] = new_point
                points.append(new_point)
                spawns.append(new_point)
                break
        else:
            spawns.pop(idx)

    # Normalize the points to [-1, 1] range to fit the circle of radius 1 centered at (0, 0)
    scaler = MinMaxScaler(feature_range=(-1, 1))
    points = scaler.fit_transform(points)

    return points

def main():
    points = poisson_disc_samples(0.1, width=2, height=2)
    vor = Voronoi(points)

    # Plotting
    fig, ax = plt.subplots()
    voronoi_plot_2d(vor, ax=ax, show_vertices=True, line_colors='orange', line_width=2, line_alpha=0.6, point_size=2)
    circle = plt.Circle((0, 0), 1, edgecolor='b', facecolor='none')
    ax.add_artist(circle)
    ax.set_xlim(-1, 1)
    ax.set_ylim(-1, 1)
    ax.set_aspect('equal')
    plt.show()

if __name__ == "__main__":
    main()