import numpy as np
from scipy.spatial import Delaunay, ConvexHull
from shapely.geometry import Polygon, Point


def alpha_shape(points, alpha):
    """
    Compute the alpha shape (concave hull) of a set of points.
    """

    def add_edge(edges, edge_points, points, i, j):
        """Add a line between the i-th and j-th points if not in the list already"""
        if (i, j) in edges or (j, i) in edges:
            return
        edges.add((i, j))
        edge_points.append(points[[i, j]])

    tri = Delaunay(points)
    edges = set()
    edge_points = []

    # Loop over triangles:
    for ia, ib, ic in tri.simplices:
        pa = points[ia]
        pb = points[ib]
        pc = points[ic]
        # Lengths of sides of triangle
        a = np.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
        b = np.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
        c = np.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)
        # Semiperimeter of triangle
        s = (a + b + c) / 2.0
        # Area of triangle by Heron's formula
        area = np.sqrt(s * (s - a) * (s - b) * (s - c))
        circum_r = a * b * c / (4.0 * area)
        # Here's the radius filter.
        if circum_r < 1.0 / alpha:
            add_edge(edges, edge_points, points, ia, ib)
            add_edge(edges, edge_points, points, ib, ic)
            add_edge(edges, edge_points, points, ic, ia)

    m = np.array(edge_points)
    return m


def generate_mesh(points, depth, alpha=0.1):
    """
    Generate a mesh by extruding a 2D shape along the Z-axis, automatically detecting holes.

    :param points: List of (x, y) tuples representing all points of the 2D shape, including potential holes.
    :param depth: Extrusion depth along the Z-axis.
    :param alpha: Alpha value for the alpha shape algorithm (controls the "tightness" of the boundary).
    :return: Tuple of vertices and faces.
    """
    # Convert points to a numpy array
    points_2d = np.array(points)

    # Compute the alpha shape (outer boundary)
    boundary_edges = alpha_shape(points_2d, alpha)

    # Create a Polygon from the boundary
    boundary_polygon = Polygon(boundary_edges)

    # Separate points into boundary and interior
    boundary_points = []
    interior_points = []
    for point in points:
        if Point(point).touches(boundary_polygon) or Point(point).within(boundary_polygon):
            if Point(point).touches(boundary_polygon):
                boundary_points.append(point)
            else:
                interior_points.append(point)

    # Perform Delaunay triangulation on all points
    tri = Delaunay(points_2d)

    # Generate the top and bottom faces
    bottom_face = np.hstack((tri.points, np.zeros((tri.points.shape[0], 1))))
    top_face = np.hstack((tri.points, np.ones((tri.points.shape[0], 1)) * depth))

    # Combine top and bottom vertices
    vertices_array = np.vstack((bottom_face, top_face))

    # Create faces
    faces = []

    # Bottom face triangulation
    for simplex in tri.simplices:
        faces.append(simplex.tolist())

    # Top face triangulation (with an offset)
    top_offset = len(tri.points)
    for simplex in tri.simplices:
        faces.append([i + top_offset for i in simplex])

    # Side faces for the outer boundary
    for i in range(len(boundary_points)):
        next_i = (i + 1) % len(boundary_points)
        current = points.index(boundary_points[i])
        next_point = points.index(boundary_points[next_i])
        faces.append([current, top_offset + current, top_offset + next_point])
        faces.append([current, top_offset + next_point, next_point])

    # Convert vertices to the desired format: list of tuples
    vertices = [tuple(vertex) for vertex in vertices_array]

    return vertices, faces

def generate_mesh_wholes(points, holes, depth):
    """
    Generate a mesh by extruding a 2D shape along the Z-axis, including holes.

    :param points: List of (x, y) tuples representing the outer boundary of the 2D shape.
    :param holes: List of lists, where each inner list contains (x, y) tuples representing a hole.
    :param depth: Extrusion depth along the Z-axis.
    :return: Tuple of vertices and faces.
    """
    # Convert points to a numpy array
    points_2d = np.array(points)

    # Prepare points for triangulation
    triangulation_points = points_2d.tolist()
    for hole in holes:
        triangulation_points.extend(hole)

    # Perform Delaunay triangulation
    tri = Delaunay(np.array(triangulation_points))

    # Generate the top and bottom faces
    bottom_face = np.hstack((tri.points, np.zeros((tri.points.shape[0], 1))))
    top_face = np.hstack((tri.points, np.ones((tri.points.shape[0], 1)) * depth))

    # Combine top and bottom vertices
    vertices_array = np.vstack((bottom_face, top_face))

    # Create faces
    faces = []

    # Bottom face triangulation
    for simplex in tri.simplices:
        faces.append(simplex.tolist())

    # Top face triangulation (with an offset)
    top_offset = len(tri.points)
    for simplex in tri.simplices:
        faces.append([i + top_offset for i in simplex])

    # Side faces
    for i in range(len(points)):
        next_i = (i + 1) % len(points)
        faces.append([i, top_offset + i, top_offset + next_i])
        faces.append([i, top_offset + next_i, next_i])

    # Side faces for holes
    start_index = len(points)
    for hole in holes:
        for i in range(len(hole)):
            current = start_index + i
            next_i = start_index + (i + 1) % len(hole)
            faces.append([current, top_offset + next_i, top_offset + current])
            faces.append([current, next_i, top_offset + next_i])
        start_index += len(hole)

    # Convert vertices to the desired format: list of tuples
    vertices = [tuple(vertex) for vertex in vertices_array]

    return vertices, faces

def generate_mesh_simple(points, depth):
    """
    Generate a mesh by extruding a 2D shape along the Z-axis.

    :param points: List of (x, y) tuples representing the 2D shape.
    :param depth: Extrusion depth along the Z-axis.
    :return: Tuple of vertices and faces.
    """
    # Convert points to a numpy array
    points_2d = np.array(points)

    # Get the convex hull of the points to ensure they form a proper polygon
    hull = ConvexHull(points_2d)
    hull_points = points_2d[hull.vertices]

    # Generate the top and bottom faces
    bottom_face = np.hstack((hull_points, np.zeros((hull_points.shape[0], 1))))
    top_face = np.hstack((hull_points, np.ones((hull_points.shape[0], 1)) * depth))

    # Combine top and bottom vertices
    vertices_array = np.vstack((bottom_face, top_face))

    # Create faces
    faces = []

    # Bottom face triangulation (counter-clockwise)
    for i in range(len(hull_points) - 2):
        faces.append([0, i + 2, i + 1])

    # Top face triangulation (counter-clockwise, with an offset)
    top_offset = len(hull_points)
    for i in range(len(hull_points) - 2):
        faces.append([top_offset, top_offset + i + 1, top_offset + i + 2])

    # Side faces (ensure counter-clockwise order)
    for i in range(len(hull_points)):
        next_i = (i + 1) % len(hull_points)
        faces.append([i, top_offset + i, top_offset + next_i])
        faces.append([i, top_offset + next_i, next_i])

    # Convert vertices to the desired format: list of tuples
    vertices = [tuple(vertex) for vertex in vertices_array]

    return vertices, faces