/*
 * Copyright 2022 Rive
 */

#include "rive/math/hit_test.hpp"

#include <algorithm>
#include <assert.h>
#include <cmath>

using namespace rive;

static inline float graphics_roundf(float x) {
    return std::floor(x + 0.5f);
}

static inline int graphics_round(float x) {
    return (int)graphics_roundf(x);
}

struct Point {
    float x, y;
    
    Point() {}
    Point(float xx, float yy) : x(xx), y(yy) {}
    Point(const Vec2D& src) : x(src.x()), y(src.y()) {}

    Point operator+(Point v) const {
        return { x + v.x, y + v.y };
    }
    Point operator-(Point v) const {
        return { x - v.x, y - v.y };
    }

    Point& operator+=(Point v) {
        *this = *this + v;
        return *this;
    }
    Point& operator-=(Point v) {
        *this = *this - v;
        return *this;
    }

    friend Point operator*(Point v, float s) {
        return { v.x * s, v.y * s };
    }
    friend Point operator*(float s, Point v) {
        return { v.x * s, v.y * s };
    }
};

template <typename T> T lerp(T a, T b, float t) {
    return a + (b - a) * t;
}

template <typename T> T ave(T a, T b) {
    return lerp(a, b, 0.5f);
}

static void append_line(const float height, Point p0, Point p1, float m,
                        int winding, int delta[], int iwidth) {
    assert(winding == 1 || winding == -1);

    int top = graphics_round(p0.y);
    int bottom = graphics_round(p1.y);
    if (top == bottom) {
        return;
    }

    assert(top < bottom);
    assert(top >= 0);
    assert((float)bottom <= height);

    // we add 0.5 at the end to pre-round the values
    float x = p0.x + m * (top - p0.y + 0.5f) + 0.5f;

    int* row = delta + top * iwidth;
    for (int y = top; y < bottom; ++y) {
        int ix = (int)std::max(x, 0.0f);
        if (ix < iwidth) {
            row[ix] += winding;
        }
        x += m;
        row += iwidth;
    }
}

static void clip_line(const float height, Point p0, Point p1, int delta[], const int iwidth) {
    if (p0.y == p1.y) {
        return;
    }

    int winding = 1;
    if (p0.y > p1.y) {
        winding = -1;
        std::swap(p0, p1);
    }
    // now we're monotonic in Y: p0 <= p1
    if (p1.y <= 0 || p0.y >= height) {
        return;
    }

    const float m = (float)(p1.x - p0.x) / (p1.y - p0.y);
    if (p0.y < 0) {
        p0.x += m * (0 - p0.y);
        p0.y = 0;
    }
    if (p1.y > height) {
        p1.x += m * (height - p1.y);
        p1.y = height;
    }
    
    assert(p0.y <= p1.y);
    assert(p0.y >= 0);
    assert(p1.y <= height);

    append_line(height, p0, p1, m, winding, delta, iwidth);
}

#define MAX_CURVE_SEGMENTS  (1 << 8)

static int compute_cubic_segments(Point a, Point b, Point c, Point d) {
    Point abc = a - b - b + c;
    Point bcd = b - c - c + d;
    float dx = std::max(std::abs(abc.x), std::abs(bcd.x));
    float dy = std::max(std::abs(abc.y), std::abs(bcd.y));
    float dist = sqrtf(dx*dx + dy*dy);
    // count = sqrt(6*dist / 8*tol)
    // tol = 0.25
    // count = sqrt(3*dist)
    float count = sqrtf(3*dist);
    return std::max(1, std::min((int)ceilf(count), MAX_CURVE_SEGMENTS));
}

// cubic a(1-t)^3 + 3bt(1-t)^2 + 3c(1-t)t^2 + dt^3
// becomes
// At^3 + Bt^2 + Ct + D
//
struct CubicCoeff {
    Point A, B, C, D;
    
    // a(1-t)^3 + 3bt(1-t)^2 + 3ct^2(1-t) + dt^3
    // a - 3at + 3at^2 -  at^3      a(1 - 3t + 3t^2 - t^3)
    //     3bt - 6bt^2 + 3bt^3     3b(     t - 2t^2 + t^3)
    //           3ct^2 - 3ct^3     3c(          t^2 - t^3)
    //                    dt^3      d(                t^3)
    // ...
    // D  + Ct  + Bt^2  + At^3
    //
    CubicCoeff(Point a, Point b, Point c, Point d) {
        A = (d - a) + 3.0f*(b - c);
        B = 3.0f*((c - b) + (a - b));
        C = 3.0f*(b - a);
        D = a;
    }
    
    Point eval(float t) const {
        return ((A*t + B)*t + C)*t + D;
    }
};

////////////////////////////////////////////

void HitTester::reset() {
    m_DW.clear();
}

void HitTester::reset(const IAABB& clip) {
    m_offset = Vec2D{ (float)clip.left, (float)clip.top };
    m_height = (float)clip.height();

    m_IWidth = clip.width();
    m_IHeight = clip.height();
    m_DW.resize(m_IWidth * m_IHeight);
    for (size_t i = 0; i < m_DW.size(); ++i) {
        m_DW[i] = 0;
    }

    m_ExpectsMove = true;
}

void HitTester::move(Vec2D v) {
    if (!m_ExpectsMove) {
        this->close();
    }
    m_First = m_Prev = v - m_offset;
    m_ExpectsMove = false;
}

void HitTester::line(Vec2D v) {
    assert(!m_ExpectsMove);

    v = v - m_offset;
    clip_line(m_height, m_Prev, v, m_DW.data(), m_IWidth);
    m_Prev = v;
}

void HitTester::quad(Vec2D b, Vec2D c) {
    assert(!m_ExpectsMove);


    m_Prev = c;
}

static bool quickRejectCubic(float height, Point a, Point b, Point c, Point d) {
    const float h = height;
    return (a.y <= 0 && b.y <= 0 && c.y <= 0 && d.y <= 0) ||
           (a.y >= h && b.y >= h && c.y >= h && d.y >= h);
}

struct CubicChop {
    Vec2D storage[7];
    
    CubicChop(Vec2D a, Vec2D b, Vec2D c, Vec2D d) {
        auto ab = ave(a, b);
        auto bc = ave(b, c);
        auto cd = ave(c, d);
        auto abc = ave(ab, bc);
        auto bcd = ave(bc, cd);
        
        storage[0] = a;
        storage[1] = ab;
        storage[2] = abc;
        storage[3] = ave(abc, bcd);
        storage[4] = bcd;
        storage[5] = cd;
        storage[6] = d;
    }
    
    Vec2D operator[](unsigned index) const {
        assert(index < 7);
        return storage[index];
    }
};

// Trial and error to pick a good value for this.
//
// Subdivision and recursion have their own cost, so at some point
// just evaluating the cubic (count) times is cheaper than continuing
// to chop.
//
// The key win is quickRejectCubic. This is how we save time over just
// evaluating the cubic up front.
//
#define MAX_LOCAL_SEGMENTS  16

void HitTester::recurse_cubic(Vec2D b, Vec2D c, Vec2D d, int count) {
    if (quickRejectCubic(m_height, m_Prev, b, c, d)) {
        m_Prev = d;
        return;
    }
    
    if (count > MAX_LOCAL_SEGMENTS) {
        CubicChop chop(m_Prev, b, c, d);
        const int newCount = (count + 1) >> 1;
        assert(newCount < count);
        this->recurse_cubic(chop[1], chop[2], chop[3], newCount);
        this->recurse_cubic(chop[4], chop[5], chop[6], newCount);
    } else {
        const float dt = 1.0f / (float)count;
        float t = dt;

        CubicCoeff cube(m_Prev, b, c, d);
        // we don't need the first point eval(0) or the last eval(1)
        Point prev = m_Prev;
        for (int i = 1; i < count - 1; ++i) {
            auto next = cube.eval(t);
            clip_line(m_height, prev, next, m_DW.data(), m_IWidth);
            prev = next;
            t += dt;
        }
        clip_line(m_height, prev, d, m_DW.data(), m_IWidth);
        m_Prev = d;
    }
}
void HitTester::cubic(Vec2D b, Vec2D c, Vec2D d) {
    assert(!m_ExpectsMove);

    b = b - m_offset;
    c = c - m_offset;
    d = d - m_offset;
    
    if (quickRejectCubic(m_height, m_Prev, b, c, d)) {
        m_Prev = d;
        return;
    }
    
    const int count = compute_cubic_segments(m_Prev, b, c, d);

    this->recurse_cubic(b, c, d, count);
}

void HitTester::close() {
    assert(!m_ExpectsMove);

    clip_line(m_height, m_Prev, m_First, m_DW.data(), m_IWidth);
    m_ExpectsMove = true;
}

void HitTester::addRect(const AABB& rect, const Mat2D& xform, PathDirection dir) {
    const Vec2D pts[] = {
        xform * Vec2D{rect.left(),  rect.top()},
        xform * Vec2D{rect.right(), rect.top()},
        xform * Vec2D{rect.right(), rect.bottom()},
        xform * Vec2D{rect.left(),  rect.bottom()},
    };

    move(pts[0]);
    if (dir == PathDirection::clockwise) {
        line(pts[1]);
        line(pts[2]);
        line(pts[3]);
    } else {
        line(pts[3]);
        line(pts[2]);
        line(pts[1]);
    }
    close();
}

bool HitTester::test(FillRule rule) {
    if (!m_ExpectsMove) {
        this->close();
    }

    const int mask = (rule == rive::FillRule::nonZero) ? -1 : 1;

    int nonzero = 0;
    for (auto m : m_DW) {
        nonzero |= (m & mask);
    }
    return nonzero != 0;
}