import numpy as np
import matplotlib.pyplot as plt

def line_sdf(P, x1, y1, x2, y2):
    """
    Calculate the signed distance from a point P to a line segment
    defined by two endpoints (x1, y1) and (x2, y2).
    """

    a = np.array([x2 - x1, y2 - y1])
    a = a / np.linalg.norm(a)
    b = P - np.array([x1, y1])
    d = np.dot(b, np.array([a[1], -a[0]]))
    return d

def circle_sdf(P, xc, yc, r):
    """
    Calculate the signed distance from a point P to a circle defined by its center (xc, yc) and radius r.
    """
    return np.sqrt((P[0] - xc) ** 2 + (P[1] - yc) ** 2) - r

def intersect_sdf(d1, d2):
    """
    Calculate the signed distance field resulting from the intersection of two distance fields (d1 and d2).
    """
    return max(d1, d2)

def union_sdf(d1, d2):
    """
    Calculate the signed distance field resulting from the union of two distance fields (d1 and d2).
    """
    return min(d1, d2)

def find_tangent_point(x1, y1, d1, x2, y2, d2):
    """
    Find the tangent point that has radius(d1) distance from center of circle
    and calculated distance (d2) from the bottom tip of the heart. The result is
    limited to desired tangent point.
    """

    centerdx = x1 - x2
    centerdy = y1 - y2
    R = np.sqrt(centerdx**2 + centerdy**2)

    if not (abs(d1 - d2) <= R and R <= d1 + d2):
        # No intersections
        return []

    d1d2 = d1**2 - d2**2
    a = d1d2 / (2 * R**2)

    c = np.sqrt(2 * (d1**2 + d2**2) / R**2 - (d1d2**2) / R**4 - 1)

    fx = (x1 + x2) / 2 + a * (x2 - x1)
    gx = c * (y2 - y1) / 2
    # ix1 = fx + gx
    ix2 = fx - gx

    fy = (y1 + y2) / 2 + a * (y2 - y1)
    gy = c * (x1 - x2) / 2
    # iy1 = fy + gy
    iy2 = fy - gy

    return [ix2, iy2]

def heart_sdf(p, r=4):
    """
    Calculate the signed distance from a point P to a heart shape defined by two circles and three line segments.

    The heart shape consists of two circles representing the left and right sides, and three line segments representing the bottom tip and the tangent lines connecting the circles.

    Parameters:
        p (tuple): A tuple containing the coordinates (x, y) of the point P.
        r (float): The radius of the heart shape. Default is 4.

    Returns:
        float: The signed distance from the point P to the heart shape. Negative values indicate that the point is inside the heart shape, zero indicates that the point is on the boundary, and positive values indicate that the point is outside the heart shape.

    Note:
        The heart shape is constructed based on mathematical equations and geometric calculations. The function utilizes signed distance functions (SDFs) for circles and lines to determine the distance from the point P to various components of the heart shape.
    """
    # Some experimental ratio for the center of circles based on their radius
    a = r * 3 / 4

    circle1 = circle_sdf(p, -a, 0, r)  # Left Circle
    circle2 = circle_sdf(p, a, 0, r)  # Right Circle

    # Distance from bottom tip to center of circles
    d2c = np.sqrt((a - 0) ** 2 + (0 - 2 * a - r) ** 2)  
    
    # Distance to tangent point of circle using Pythagorean theorem
    d2t = np.sqrt(d2c**2 - r**2)  

    tpr = find_tangent_point(
        a, 0, r, 0, -2 * a - r, d2t
    )  # Tangent point on right circle

    line1 = line_sdf(p, -tpr[0], tpr[1], 0, -2 * a - r)
    line2 = line_sdf(p, 0, -2 * a - r, tpr[0], tpr[1])
    line3 = line_sdf(p, tpr[0], tpr[1], -tpr[0], tpr[1])

    #    Create a triangle which base is the line that
    #    connects tangent points and the vertex point is heart bottom tip

    dl = intersect_sdf(intersect_sdf(line1, line2), line3)

    d = union_sdf(union_sdf(circle1, circle2), dl)

    return d

def plot_sdf(SDF, BdBox=[-10, 10, -12, 8], n=300):
    """
    Plots the signed distance function.
    """
    x, y = np.meshgrid(
        np.linspace(BdBox[0], BdBox[1], n),
        np.linspace(BdBox[2], BdBox[3], n),
    )
    points = np.hstack([x.reshape((-1, 1)), y.reshape((-1, 1))])

    sdf = np.fromiter(map(SDF, points), dtype=float)

    inner = np.where(sdf <= 0, 1, 0)

    _, ax = plt.subplots(figsize=(8, 6))
    _ = ax.imshow(
        inner.reshape((n, n)),
        extent=(BdBox[0], BdBox[1], BdBox[2], BdBox[3]),
        origin="lower",
        cmap="Reds",
        alpha=0.8,
    )
    ax.contour(x, y, sdf.reshape((n, n)), levels=[0], colors="gold", linewidths=2)
    ax.set_xlabel("X", fontweight="bold")
    ax.set_ylabel("Y", fontweight="bold")
    ax.set_title("SDF Visualization", fontweight="bold", fontsize=16)
    ax.set_aspect("equal")

    plt.show()

plot_sdf(heart_sdf)