import math


def cotangent_poisson_field(
    planar_points,
    planar_faces,
    boundary_mask,
    tolerance=1e-8,
):
    import numpy as np

    points = np.asarray(planar_points, dtype=np.float64)
    boundary_mask = np.asarray(boundary_mask, dtype=bool)
    loads = np.zeros(len(points), dtype=np.float64)
    edge_weights = {}

    def add_edge(start, end, weight):
        edge = tuple(sorted((int(start), int(end))))
        edge_weights[edge] = edge_weights.get(edge, 0.0) + weight

    for face in planar_faces:
        first, second, third = (int(index) for index in face)
        a, b, c = points[[first, second, third]]
        cross = abs(
            (b[0] - a[0]) * (c[1] - a[1])
            - (b[1] - a[1]) * (c[0] - a[0])
        )
        if cross <= 1e-14:
            continue
        area = cross * 0.5
        loads[[first, second, third]] += area / 3.0
        cotangent_a = float((b - a).dot(c - a) / cross)
        cotangent_b = float((a - b).dot(c - b) / cross)
        cotangent_c = float((a - c).dot(b - c) / cross)
        add_edge(second, third, cotangent_a * 0.5)
        add_edge(first, third, cotangent_b * 0.5)
        add_edge(first, second, cotangent_c * 0.5)

    interior = np.flatnonzero(~boundary_mask)
    field = np.zeros(len(points), dtype=np.float64)
    if not len(interior):
        return field
    compact = np.full(len(points), -1, dtype=np.int64)
    compact[interior] = np.arange(len(interior), dtype=np.int64)
    diagonal = np.zeros(len(interior), dtype=np.float64)
    edge_starts = []
    edge_ends = []
    weights = []
    for (start, end), weight in edge_weights.items():
        start_compact = compact[start]
        end_compact = compact[end]
        if start_compact >= 0:
            diagonal[start_compact] += weight
        if end_compact >= 0:
            diagonal[end_compact] += weight
        if start_compact >= 0 and end_compact >= 0:
            edge_starts.append(start_compact)
            edge_ends.append(end_compact)
            weights.append(weight)
    if np.any(diagonal <= 1e-14):
        raise ValueError("Unable to build a stable Poisson Balloon field")

    edge_starts = np.asarray(edge_starts, dtype=np.int64)
    edge_ends = np.asarray(edge_ends, dtype=np.int64)
    weights = np.asarray(weights, dtype=np.float64)

    def multiply(values):
        result = diagonal * values
        if len(weights):
            np.add.at(result, edge_starts, -weights * values[edge_ends])
            np.add.at(result, edge_ends, -weights * values[edge_starts])
        return result

    right_hand_side = loads[interior]
    solution = np.zeros(len(interior), dtype=np.float64)
    residual = right_hand_side.copy()
    preconditioned = residual / diagonal
    direction = preconditioned.copy()
    residual_product = float(residual.dot(preconditioned))
    target = max(
        float(np.linalg.norm(right_hand_side)) * tolerance,
        1e-12,
    )
    max_iterations = max(64, min(len(interior) * 4, 4000))
    for _iteration in range(max_iterations):
        if float(np.linalg.norm(residual)) <= target:
            break
        multiplied = multiply(direction)
        denominator = float(direction.dot(multiplied))
        if denominator <= 1e-20 or not math.isfinite(denominator):
            raise ValueError("Unable to solve the Poisson Balloon field")
        step = residual_product / denominator
        solution += step * direction
        residual -= step * multiplied
        preconditioned = residual / diagonal
        next_product = float(residual.dot(preconditioned))
        if not math.isfinite(next_product):
            raise ValueError("Unable to solve the Poisson Balloon field")
        if abs(residual_product) <= 1e-30:
            break
        direction = preconditioned + (
            next_product / residual_product
        ) * direction
        residual_product = next_product

    field[interior] = np.maximum(solution, 0.0)
    return field


def poisson_balloon_heights(
    planar_points,
    planar_faces,
    boundary_mask,
    inflation,
):
    import numpy as np

    field = cotangent_poisson_field(
        planar_points,
        planar_faces,
        boundary_mask,
    )
    return 2.0 * inflation * np.sqrt(field)