Voronoi tesselation python

Question

I want to do a voronoi tessellation in a closed domain but I can't identify where I'm going wrong. The problem is that the domain where I want to do the tessellation is a rectangle

here is the code if anyone can help i appreciate it.

I can't identify which are some parts of the code that are causing problems, like the generation of ghost points and the filtering of valid rules. I've already looked at expanding the points to make sure they're inside the rectangular domain.

import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys
from scipy.spatial import Voronoi, voronoi_plot_2d

eps = sys.float_info.epsilon # variável que armazena o valor da menor diferença,
                              # entre números de ponto flutuante representáveis

n_geradores = 5 # quantidade de geradores da tesselação
a = 0  # Limite inferior
b = 6  # Limite superior

c = 0
d = 10 

def gerar_geradores(n, a, b, c, d):
    geradores = np.empty((n, 2))
    for i in range(n):
        geradores[i, 0] = np.random.uniform(a, b + np.finfo(float).eps)
        geradores[i, 1] = np.random.uniform(c, d + np.finfo(float).eps)
    return geradores

# Gera uma matriz com n_geradores linhas e 2 colunas, com valores aleatórios no intervalo [a, b]
geradores = gerar_geradores(n_geradores, a, b, c, d)
#print(geradores)
#print(geradores[0, :])
dominio = np.array([0., 6., 0., 10.]) 


def verificar_limites_do_dominio(geradores, dominio):
    return np.logical_and(np.logical_and(dominio[0] <= geradores[:, 0],
                                         geradores[:, 0] <= dominio[1]),
                          np.logical_and(dominio[2] <= geradores[:, 1],
                                         geradores[:, 1] <= dominio[3]))


def voronoi(geradores, dominio):
    i = verificar_limites_do_dominio(geradores, dominio)
    #print(i)
    points_center = geradores[i, :]
    points_left = np.copy(points_center)
    #rint(points_left)
    points_left[:, 0] = dominio[0] - (points_left[:, 0] - dominio[0])
    #print(points_left[:, 0])
    #print(points_left)
    points_right = np.copy(points_center)
    points_right[:, 0] = dominio[1] + (dominio[1] - points_right[:, 0])
    #print(points_right[:, 0])
    points_down = np.copy(points_center)
    points_down[:, 1] = dominio[2] - (points_down[:, 1] - dominio[2])
    #print(points_down)
    points_up = np.copy(points_center)
    points_up[:, 1] = dominio[3] + (dominio[3] - points_up[:, 1])
    points = np.append(points_center,
                       np.append(np.append(points_left,
                                           points_right,
                                           axis=0),
                                 np.append(points_down,
                                           points_up,
                                           axis=0),
                                 axis=0),
                       axis=0)
    #print(points)
    # Compute Voronoi
    vor = sp.spatial.Voronoi(points)
    fig, ax = plt.subplots()
    voronoi_plot_2d(vor, ax = ax)
    ax.grid(True)
    print(vor.vertices)
    # Filter regions
    regions = []
    points_to_filter = [] # nós também precisaremos coletar pontos
    ind = np.arange(points.shape[0])
    ind = np.expand_dims(ind,axis= 1)


    for i,region in enumerate(vor.regions): # Enumerar as regiõe
        if not region: # Melhor pular completamente a região vazia
            continue

        flag = True
        for index in region:
            if index == -1:
                flag = False
                break
            else:
                x = vor.vertices[index, 0]
                y = vor.vertices[index, 1]
                if not(dominio[0] - eps <= x and x <= dominio[1] + eps and
                      dominio[2] - eps <= y and y <= dominio[3] + eps):
                    flag = False
                    print(i)
                    break
        if flag:
            regions.append(region)

            # Encontrar o ponto que está dentro
            points_to_filter.append(vor.points[vor.point_region == i][0,:])

    vor.filtered_points = np.array(points_to_filter)
    vor.filtered_regions = regions
    return vor

def centroid_region(vertices):

    A = 0

    C_x = 0

    C_y = 0
    for i in range(0, len(vertices) - 1):
        s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
        A = A + s
        C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
        C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
    A = 0.5 * A
    C_x = (1.0 / (6.0 * A)) * C_x
    C_y = (1.0 / (6.0 * A)) * C_y
    return np.array([[C_x, C_y]])

def setor_cor(y, dominio):
    terco = (dominio[3] - dominio[2]) / 3
    if y <= terco:
        return 'red'
    elif terco < y <= 2 * terco:
        return 'green'
    else:
        return 'orange'

def plot(r,index):
    vor = voronoi(r, dominio)

    fig = pl.figure()
    ax = fig.gca()

    # Plot initial points
    ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')

    # Plot ridges
    for i, region in enumerate(vor.filtered_regions):
        vertices = vor.vertices[region + [region[0]], :]
        print(vertices)
        color = setor_cor(centroid_region(vertices)[0, 1], dominio)
        polygon = pl.Polygon(vertices, closed=True, facecolor=color, edgecolor='none')
        ax.add_patch(polygon)
        ax.plot(vertices[:, 0], vertices[:, 1], 'k-')

        # Compute and plot centroids
        centroid = centroid_region(vertices)
        #ax.plot(centroid[:, 0], centroid[:, 1], 'w.')

    ger = np.copy(vor.filtered_points)

    # Ajustar a relação de aspecto e os limites dos eixos
    ax.set_aspect('equal', adjustable='box')
    ax.set_xlim([-0.5, 6.5])
    ax.set_ylim([-0.5, 10.5])

    pl.savefig("voronoi" + str(index) + ".png")
    return vor.filtered_regions

if __name__ == '__main__':
    for i in range(1):
        regions = plot(geradores, i)