// A modified version of the bezier code in ImGui with extra features for bitmap and path interaction.

// Function to convert a ratio back into a point for the bezier handles.
v2 Line_RatioToPoint(v2 a, v2 b, real32 ratio)
{
    v2 ab_dir = b - a;
    real32 ab_len_sqr = ab_dir.x * ab_dir.x + ab_dir.y * ab_dir.y;
    real32 dot = ratio*ab_len_sqr;
    return a + ab_dir * dot / V2(ab_len_sqr);
}

// The ratio here is just the dot product divided by the squared length.
v2 Bezier_LineClosestPoint2(v2 a, v2 b, v2 p, real32 *ratio)
{
    v2 ap = p - a;
    v2 ab_dir = b - a;
    real32 dot = ap.x * ab_dir.x + ap.y * ab_dir.y;
    if (dot < 0.0f) {
        *ratio = 0.0f;
        return a;
    }
    real32 ab_len_sqr = ab_dir.x * ab_dir.x + ab_dir.y * ab_dir.y;
    if (dot > ab_len_sqr) {
        *ratio = 1.0f;
        return b;
    }
    *ratio = dot / ab_len_sqr;
    return a + ab_dir * dot / ab_len_sqr;
}

// Following the algorithm, we take the ratio from the _leftmost_ point in each
// subdivision of the cubic spline until we're within tess_tol, and then we
// interpolate it with the subdivision's rightmost point in the ClosestPoint call.
// The pow(0.5, level) represents the ratio of the next subdivision's leftmost
// point (AKA the rightmost point of the current subdivision).

// finds the point closest to p and also fills out its ratio along the curve
static void Bezier_CubicClosestPointCasteljauStep(v2 p, v2 *p_closest, real32 ratio, real32 *r_closest, v2 *p_last, real32 *p_closest_dist2,
                                                    real32 x1, real32 y1, real32 x2, real32 y2, real32 x3, real32 y3, real32 x4, real32 y4, real32 tess_tol, int level)
{
    real32 dx = x4 - x1;
    real32 dy = y4 - y1;
    real32 d2 = ((x2 - x4) * dy - (y2 - y4) * dx);
    real32 d3 = ((x3 - x4) * dy - (y3 - y4) * dx);
    d2 = (d2 >= 0) ? d2 : -d2;
    d3 = (d3 >= 0) ? d3 : -d3;
    if ((d2 + d3) * (d2 + d3) < tess_tol * (dx * dx + dy * dy))
    {
        v2 p_current = V2(x4, y4);
        real32 added_ratio;
        v2 p_line = Bezier_LineClosestPoint2(*p_last, p_current, p, &added_ratio);
        real32 dist2 = LengthSq(p - p_line);
        if (dist2 < *p_closest_dist2)
        {
            *p_closest = p_line;
            *p_closest_dist2 = dist2;
            *r_closest = ratio + pow(0.5, level)*added_ratio;
        }
        *p_last = p_current;
    }
    else if (level < 10)
    {
        real32 x12 = (x1 + x2)*0.5f,       y12 = (y1 + y2)*0.5f;
        real32 x23 = (x2 + x3)*0.5f,       y23 = (y2 + y3)*0.5f;
        real32 x34 = (x3 + x4)*0.5f,       y34 = (y3 + y4)*0.5f;
        real32 x123 = (x12 + x23)*0.5f,    y123 = (y12 + y23)*0.5f;
        real32 x234 = (x23 + x34)*0.5f,    y234 = (y23 + y34)*0.5f;
        real32 x1234 = (x123 + x234)*0.5f, y1234 = (y123 + y234)*0.5f;
        Bezier_CubicClosestPointCasteljauStep(p, p_closest, ratio, r_closest, p_last, p_closest_dist2, x1, y1, x12, y12, x123, y123, x1234, y1234, tess_tol, level + 1);
        Bezier_CubicClosestPointCasteljauStep(p, p_closest, ratio + pow(0.5, level+1), r_closest, p_last, p_closest_dist2, x1234, y1234, x234, y234, x34, y34, x4, y4, tess_tol, level + 1);
    }
}

// finds the min/max bounds of the curve
static void Bezier_CubicMinMaxCasteljauStep(v2 *p_min, v2 *p_max, real32 x1, real32 y1, real32 x2, real32 y2, real32 x3, real32 y3, real32 x4, real32 y4, real32 tess_tol, int level)
{
    real32 dx = x4 - x1;
    real32 dy = y4 - y1;
    real32 d2 = ((x2 - x4) * dy - (y2 - y4) * dx);
    real32 d3 = ((x3 - x4) * dy - (y3 - y4) * dx);
    d2 = (d2 >= 0) ? d2 : -d2;
    d3 = (d3 >= 0) ? d3 : -d3;
    if ((d2 + d3) * (d2 + d3) < tess_tol * (dx * dx + dy * dy))
    {
        v2 p_current = V2(x4, y4);
        if (p_current.x < p_min->x) {
            p_min->x = p_current.x;
        }
        if (p_current.y < p_min->y) {
            p_min->y = p_current.y;
        }
        if (p_current.x > p_max->x) {
            p_max->x = p_current.x;
        }
        if (p_current.y > p_max->y) {
            p_max->y = p_current.y;
        }
    }
    else if (level < 10)
    {
        real32 x12 = (x1 + x2)*0.5f,       y12 = (y1 + y2)*0.5f;
        real32 x23 = (x2 + x3)*0.5f,       y23 = (y2 + y3)*0.5f;
        real32 x34 = (x3 + x4)*0.5f,       y34 = (y3 + y4)*0.5f;
        real32 x123 = (x12 + x23)*0.5f,    y123 = (y12 + y23)*0.5f;
        real32 x234 = (x23 + x34)*0.5f,    y234 = (y23 + y34)*0.5f;
        real32 x1234 = (x123 + x234)*0.5f, y1234 = (y123 + y234)*0.5f;
        Bezier_CubicMinMaxCasteljauStep(p_min, p_max, x1, y1, x12, y12, x123, y123, x1234, y1234, tess_tol, level + 1);
        Bezier_CubicMinMaxCasteljauStep(p_min, p_max, x1234, y1234, x234, y234, x34, y34, x4, y4, tess_tol, level + 1);
    }
}

// return all points
static void Bezier_CubicCalcPointsCasteljauStep(void *Data, uint32 *Increment, real32 x1, real32 y1, real32 x2, real32 y2, real32 x3, real32 y3, real32 x4, real32 y4, real32 tess_tol, int level)
{
    real32 dx = x4 - x1;
    real32 dy = y4 - y1;
    real32 d2 = ((x2 - x4) * dy - (y2 - y4) * dx);
    real32 d3 = ((x3 - x4) * dy - (y3 - y4) * dx);
    d2 = (d2 >= 0) ? d2 : -d2;
    d3 = (d3 >= 0) ? d3 : -d3;
    if ((d2 + d3) * (d2 + d3) < tess_tol * (dx * dx + dy * dy))
    {
        *((real32 *)Data + *Increment*3) = x4;
        *((real32 *)Data + *Increment*3 + 1) = y4;
        *((real32 *)Data + *Increment*3 + 2) = 0;
        *Increment += 1;
    }
    else if (level < 10)
    {
        real32 x12 = (x1 + x2)*0.5f,       y12 = (y1 + y2)*0.5f;
        real32 x23 = (x2 + x3)*0.5f,       y23 = (y2 + y3)*0.5f;
        real32 x34 = (x3 + x4)*0.5f,       y34 = (y3 + y4)*0.5f;
        real32 x123 = (x12 + x23)*0.5f,    y123 = (y12 + y23)*0.5f;
        real32 x234 = (x23 + x34)*0.5f,    y234 = (y23 + y34)*0.5f;
        real32 x1234 = (x123 + x234)*0.5f, y1234 = (y123 + y234)*0.5f;
        Bezier_CubicCalcPointsCasteljauStep(Data, Increment, x1, y1, x12, y12, x123, y123, x1234, y1234, tess_tol, level + 1);
        Bezier_CubicCalcPointsCasteljauStep(Data, Increment, x1234, y1234, x234, y234, x34, y34, x4, y4, tess_tol, level + 1);
    }
}

real32 Bezier_CubicRatioOfPoint(v2 p1, v2 p2, v2 p3, v2 p4, v2 p)
{
    real32 tess_tol = TESS_TOL;
    v2 p_last = p1;
    v2 p_closest;
    real32 p_closest_dist2 = FLT_MAX;
    real32 ratio = 0;
    Bezier_CubicClosestPointCasteljauStep(p, &p_closest, 0, &ratio, &p_last, &p_closest_dist2, p1.x, p1.y, p2.x, p2.y, p3.x, p3.y, p4.x, p4.y, tess_tol, 0);
    return ratio;
}

void Bezier_CubicCalcPoints(v2 p1, v2 p2, v2 p3, v2 p4, void *Data, uint32 *Increment)
{
    real32 tess_tol = TESS_TOL;
    void *Pointer = Data;
    Bezier_CubicCalcPointsCasteljauStep(Pointer, Increment, p1.x, p1.y, p2.x, p2.y, p3.x, p3.y, p4.x, p4.y, tess_tol, 0);
}

// These functions will become more generalized as shapes are added.

static void
Mask_AddPointToLine(mask *Mask, uint16 Index, v2 Pos)
{
    for (int i = Mask->NumberOfPoints - 1; i > Index; i--) {
        Mask->Point[i+1] = Mask->Point[i];
    }
    mask_point *PointToAdd = &Mask->Point[Index+1];
    PointToAdd->Pos = Pos;
    Mask->NumberOfPoints++;
}

static void
Mask_PushPoint(mask *Mask, v2 Pos)
{
    mask_point *PointToAdd = &Mask->Point[Mask->NumberOfPoints];
    PointToAdd->Pos = Pos;
    Mask->NumberOfPoints++;
}

static void
Mask_AddPointToCurve(mask *Mask, uint16 Index, real32 ratio)
{
    mask_point *Point0 = &Mask->Point[Index];
    mask_point *Point1 = &Mask->Point[Index+1];

    v2 Point0_Pos_Right = Point0->Pos + Point0->TangentRight;
    v2 Point1_Pos_Left = Point1->Pos + Point1->TangentLeft;
    v2 Handle0_Half   = Line_RatioToPoint(Point0->Pos,      Point0_Pos_Right, ratio);
    v2 Handle1_Half   = Line_RatioToPoint(Point1_Pos_Left,  Point1->Pos,      ratio);
    v2 Top_Half       = Line_RatioToPoint(Point0_Pos_Right, Point1_Pos_Left,  ratio);
    v2 NewHandleLeft  = Line_RatioToPoint(Handle0_Half,     Top_Half,         ratio);
    v2 NewHandleRight = Line_RatioToPoint(Top_Half,         Handle1_Half,     ratio);
    v2 NewPos         = Line_RatioToPoint(NewHandleLeft,    NewHandleRight,   ratio);

    Point0->TangentRight = -(Point0->Pos - Handle0_Half);
    Point1->TangentLeft  = -(Point1->Pos - Handle1_Half);

    for (int i = Mask->NumberOfPoints - 1; i > Index; i--) {
        Mask->Point[i+1] = Mask->Point[i];
    }
    mask_point *PointToAdd = &Mask->Point[Index+1];

    PointToAdd->Pos = NewPos;
    PointToAdd->TangentLeft  = -(NewPos - NewHandleLeft);
    PointToAdd->TangentRight = -(NewPos - NewHandleRight);
    Mask->NumberOfPoints++;
}