import bpy
import bmesh
import numpy as np
import triangle as tr
from scipy.spatial import ConvexHull

def get_points():
    """ Extracts 3D points from the 'Point' object in Blender. """
    
    obj_points = bpy.data.objects.get("Point")  # Ensure point object exists
    
    if not obj_points:
        print("Error: Object 'Point' not found in the scene!")
        return None

    # Extract 3D points
    points = []
    for vert in obj_points.data.vertices:
        world_co = obj_points.matrix_world @ vert.co  # Convert to global coordinates
        points.append((world_co.x, world_co.y, world_co.z))

    return points

def triangulate_points():
    """ Performs triangulation of the given points, preserving the convex hull as a boundary. """

    points = get_points()
    if points is None:
        return

    # Convert points into numpy array for Triangle
    vertices = np.array([(x, y) for x, y, _ in points])

    # Compute convex hull using SciPy
    hull = ConvexHull(vertices)
    segments = [(hull.vertices[i], hull.vertices[(i + 1) % len(hull.vertices)]) for i in range(len(hull.vertices))]

    print(f"Convex hull computed with {len(hull.vertices)} boundary points.")

    # Prepare input for Triangle
    triangulation_data = {'vertices': vertices, 'segments': np.array(segments)}

    # Perform constrained triangulation
    tri = tr.triangulate(triangulation_data, 'p')  # 'p' = preserve input points

    if 'triangles' not in tri:
        print("ERROR: Triangle did not return any triangles!")
        print("Returned data:", tri)
        return

    # Create a new triangulated mesh in Blender
    mesh = bpy.data.meshes.new("TriangulatedMesh")
    obj = bpy.data.objects.new("TriangulatedTopo", mesh)
    bpy.context.collection.objects.link(obj)

    bm = bmesh.new()
    vert_list = {}

    # Create triangulated vertices
    for i, v in enumerate(tri['vertices']):
        x, y = v
        # Find the corresponding Z value from the original points
        z = next((z for x2, y2, z in points if round(x2, 5) == round(x, 5) and round(y2, 5) == round(y, 5)), 0)
        vert_list[i] = bm.verts.new((x, y, z))

    # Create triangular faces
    for tri in tri['triangles']:
        try:
            bm.faces.new([vert_list[i] for i in tri])
        except ValueError:
            print(f" Warning: Skipping invalid face {tri}")

    bm.to_mesh(mesh)
    bm.free()

    print("Triangulation of points completed WITH convex hull boundary!")

# Run the triangulation
triangulate_points()