Handling Infinite Voronoi Cells in Point Cloud

[tnakaicode]

I tried to generate a Voronoi diagram from a point cloud and convert it into a TopoDS_Shape.
The code works well for finite cells, but I'm struggling with handling infinite Voronoi cells.
Specifically, I want to close these infinite cells to make them finite.

How can I modify this code to handle infinite Voronoi cells and close them to make them finite?

Here's a snippet of my code:

import numpy as np
from scipy.spatial import Voronoi, ConvexHull
from OCC.Core.BRep import BRep_Builder
from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakeVertex, BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeWire, BRepBuilderAPI_MakeFace
from OCC.Core.gp import gp_Pnt
from OCC.Core.TopoDS import TopoDS_Compound
from OCC.Display.SimpleGui import init_display


def generate_points_on_plane(num_points):
    """
    Function to generate random points on a plane

    Args:
        num_points: Number of points to generate. Example: 30

    Returns:
        Generated points. Shape is (num_points, 3).
    """
    points = np.random.rand(num_points, 2) * 100
    # Place points within a 100x100 range
    return points


def compute_voronoi_region(points, index, hull):
    """
    Function to compute the polyhedron defining the Voronoi region around a specific point

    Args:
        points: Point cloud. Shape is (num_points, 2).
        index: Index of the point to compute the Voronoi region for. Example: 0
        hull: ConvexHull object of the point cloud.

    Returns:
        Vertices of the region and the region index.
    """
    vor = Voronoi(points)
    region_index = vor.point_region[index]
    region = vor.regions[region_index]

    if -1 in region:
        # If the region is infinite, connect the vertices on both sides of -1 to make it finite
        finite_vertices = []
        for i in range(len(region)):
            if region[i] == -1:
                print(vor.vertices[-1])
                start_vertex = vor.vertices[region[i - 1]]
                end_vertex = vor.vertices[region[(i + 1) % len(region)]]
                if not np.array_equal(start_vertex, end_vertex):
                    # finite_vertices.append(start_vertex)
                    # finite_vertices.append(end_vertex)
                    pass
            else:
                finite_vertices.append(vor.vertices[region[i]])
        return finite_vertices, 1

    vertices = [vor.vertices[i] for i in region]
    return vertices, 1


def create_shape_from_points(points):
    """
    Function to create a Voronoi diagram from a point cloud and generate a TopoDS_Shape

    Args:
        points: Point cloud. Shape is (num_points, 2).

    Returns:
        Generated shape.
    """
    vor = Voronoi(points)
    hull = ConvexHull(points)
    builder = BRep_Builder()
    compound = TopoDS_Compound()
    builder.MakeCompound(compound)

    # Compute the Voronoi region for each point and add vertices and edges
    for i, point in enumerate(points):
        region_vertices, region = compute_voronoi_region(points, i, hull)

        if region is None:
            continue  # Skip if the region is infinite

        # Create shape from the vertices of the region and add it
        wire_builder = BRepBuilderAPI_MakeWire()
        for j in range(len(region_vertices)):
            start_vertex = region_vertices[j]
            end_vertex = region_vertices[(j + 1) % len(region_vertices)]
            # Set z-coordinate to 0
            start_pnt = gp_Pnt(start_vertex[0], start_vertex[1], 0)
            # Set z-coordinate to 0
            end_pnt = gp_Pnt(end_vertex[0], end_vertex[1], 0)
            edge_shape = BRepBuilderAPI_MakeEdge(start_pnt, end_pnt).Edge()
            wire_builder.Add(edge_shape)

        wire = wire_builder.Wire()
        face_shape = BRepBuilderAPI_MakeFace(wire).Face()
        builder.Add(compound, face_shape)

    return compound


if __name__ == "__main__":
    # Initialize display
    display, start_display, add_menu, add_function_to_menu = init_display()

    # Generate points (restricted to a plane)
    num_points = 30
    points = generate_points_on_plane(num_points)

    # Create shape from points
    shape = create_shape_from_points(points)

    # Draw points
    for point in points:
        pnt = gp_Pnt(point[0], point[1], 0)
        display.DisplayShape(
            BRepBuilderAPI_MakeVertex(pnt).Shape(), update=True)

    # Display shape
    display.DisplayShape(shape, update=True)
    start_display()

By the way, most of this code was generated using Copilot.