src/math/hit_test.cpp

/*
 * Copyright 2022 Rive
 */

#include "rive/math/hit_test.hpp"

#include "rive/math/mat2d.hpp"
#include <algorithm>
#include <assert.h>
#include <cmath>

// Should we make this an option at runtime?
#define CULL_BOUNDS true

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 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;
}

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

bool HitTester::testMesh(Vec2D pt, Span<Vec2D> verts, Span<uint16_t> indices)
{
  if (verts.size() < 3)
  {
    return false;
  }

  // Test against the bounds of the entire mesh
  // Make this optional?
  if (CULL_BOUNDS)
  {
    const auto bounds = AABB(verts);

    if (bounds.bottom() < pt.y || pt.y < bounds.top() ||
        bounds.right() < pt.x || pt.x < bounds.left())
    {
      return false;
    }
  }

  for (size_t i = 0; i < indices.size(); i += 3)
  {
    const auto a = verts[indices[i + 0]];
    const auto b = verts[indices[i + 1]];
    const auto c = verts[indices[i + 2]];
    // A collapsed triangle covers nothing, but every point passes the sign
    // test below against one: its cross products are all zero.
    if (Vec2D::cross(b - a, c - a) == 0.0f)
    {
      continue;
    }

    auto pa = a - pt;
    auto pb = b - pt;
    auto pc = c - pt;

    const float ab = Vec2D::cross(pa, pb);
    const float bc = Vec2D::cross(pb, pc);
    const float ca = Vec2D::cross(pc, pa);

    // Inside or on an edge, whichever way the triangle winds, so triangles
    // sharing an edge leave no seam between them.
    if ((ab >= 0.0f && bc >= 0.0f && ca >= 0.0f) ||
        (ab <= 0.0f && bc <= 0.0f && ca <= 0.0f))
    {
      return true;
    }
  }
  return false;
}

bool HitTester::testMesh(
    const IAABB& area,
    Span<Vec2D> verts,
    Span<uint16_t> indices)
{
  // this version can give slightly different results, so perhaps we should do
  // this automatically, ... its just much faster if we do.
  if (area.width() * area.height() == 1)
  {
    return testMesh(Vec2D((float)area.left, (float)area.top), verts, indices);
  }

  if (verts.size() < 3)
  {
    return false;
  }

  // Test against the bounds of the entire mesh
  // Make this optional?
  if (CULL_BOUNDS)
  {
    const auto bounds = AABB(verts);

    if (bounds.bottom() <= area.top || area.bottom <= bounds.top() ||
        bounds.right() <= area.left || area.right <= bounds.left())
    {
      return false;
    }
  }

  std::vector<int> windings(area.width() * area.height());
  const auto offset = Vec2D((float)area.left, (float)area.top);
  int* deltas = windings.data();

  for (size_t i = 0; i < indices.size(); i += 3)
  {
    const auto a = verts[indices[i + 0]] - offset;
    const auto b = verts[indices[i + 1]] - offset;
    const auto c = verts[indices[i + 2]] - offset;

    clip_line((float)area.height(), a, b, deltas, area.width());
    clip_line((float)area.height(), b, c, deltas, area.width());
    clip_line((float)area.height(), c, a, deltas, area.width());

    int nonzero = 0;
    for (auto w : windings)
    {
      nonzero |= w;
    }
    if (nonzero)
    {
      return true;
    }
  }
  return false;
}

// Barycentric weights of `p` in triangle abc; false when it's degenerate.
// From cross products in doubles: the Gram determinant squares the area,
// which a long, thin triangle then loses in float.
static bool barycentric(Vec2D p, Vec2D a, Vec2D b, Vec2D c, float weights[3])
{
  const double abx = (double)b.x - a.x;
  const double aby = (double)b.y - a.y;
  const double acx = (double)c.x - a.x;
  const double acy = (double)c.y - a.y;
  const double apx = (double)p.x - a.x;
  const double apy = (double)p.y - a.y;
  const double area = abx * acy - aby * acx;
  if (!(std::abs(area) >
        1e-6 * std::sqrt((abx * abx + aby * aby) * (acx * acx + acy * acy))))
  {
    return false;
  }
  const double wb = (apx * acy - apy * acx) / area;
  const double wc = (abx * apy - aby * apx) / area;
  weights[0] = (float)(1.0 - wb - wc);
  weights[1] = (float)wb;
  weights[2] = (float)wc;
  return true;
}

// The affine map taking triangle abc onto ABC, which barycentric found
// isn't degenerate.
static Mat2D triangleMap(Vec2D a, Vec2D b, Vec2D c, Vec2D A, Vec2D B, Vec2D C)
{
  const double abx = (double)b.x - a.x;
  const double aby = (double)b.y - a.y;
  const double acx = (double)c.x - a.x;
  const double acy = (double)c.y - a.y;
  const double area = abx * acy - aby * acx;
  // How the weights of B and C change with x and with y.
  const double bx = acy / area;
  const double by = -acx / area;
  const double cx = -aby / area;
  const double cy = abx / area;
  const double ABx = (double)B.x - A.x;
  const double ABy = (double)B.y - A.y;
  const double ACx = (double)C.x - A.x;
  const double ACy = (double)C.y - A.y;
  const double xx = ABx * bx + ACx * cx;
  const double xy = ABy * bx + ACy * cx;
  const double yx = ABx * by + ACx * cy;
  const double yy = ABy * by + ACy * cy;
  return Mat2D(
      (float)xx,
      (float)xy,
      (float)yx,
      (float)yy,
      (float)(A.x - (xx * a.x + yx * a.y)),
      (float)(A.y - (xy * a.x + yy * a.y)));
}

static float segmentDistanceSquared(Vec2D p, Vec2D a, Vec2D b)
{
  const Vec2D ab = b - a;
  const float lengthSquared = ab.lengthSquared();
  const float t = lengthSquared > 0.0f ? Vec2D::dot(p - a, ab) / lengthSquared
                                       : 0.0f;
  return Vec2D::distanceSquared(p, a + ab * std::min(std::max(t, 0.0f), 1.0f));
}

bool HitTester::mapMesh(
    Vec2D point,
    Span<const Vec2D> from,
    Span<const Vec2D> to,
    Span<const uint16_t> indices,
    Vec2D* out,
    bool* inside,
    Mat2D* map)
{
  auto land = [&](const uint16_t triangle[3], const float weights[3]) {
    *out = to[triangle[0]] * weights[0] + to[triangle[1]] * weights[1] +
           to[triangle[2]] * weights[2];
    if (map != nullptr)
    {
      *map = triangleMap(
          from[triangle[0]],
          from[triangle[1]],
          from[triangle[2]],
          to[triangle[0]],
          to[triangle[1]],
          to[triangle[2]]);
    }
  };
  // Lets a point on the seam between two triangles land in one of them.
  constexpr float edgeTolerance = 1e-5f;
  bool hasNearest = false;
  float nearestDistance = 0.0f;
  uint16_t nearest[3] = {0, 0, 0};
  float nearestWeights[3] = {0.0f, 0.0f, 0.0f};
  // Later triangles draw over earlier ones, so walk back to find the top one.
  for (size_t i = indices.size() - indices.size() % 3; i >= 3; i -= 3)
  {
    const uint16_t triangle[3] = {
        indices[i - 3],
        indices[i - 2],
        indices[i - 1]};
    bool inRange = true;
    for (auto index : triangle)
    {
      inRange = inRange && index < from.size() && index < to.size();
    }
    float weights[3];
    if (!inRange || !barycentric(
                        point,
                        from[triangle[0]],
                        from[triangle[1]],
                        from[triangle[2]],
                        weights))
    {
      continue;
    }
    if (weights[0] >= -edgeTolerance && weights[1] >= -edgeTolerance &&
        weights[2] >= -edgeTolerance)
    {
      land(triangle, weights);
      *inside = true;
      return true;
    }
    const float distance = std::min(
        {segmentDistanceSquared(point, from[triangle[0]], from[triangle[1]]),
         segmentDistanceSquared(point, from[triangle[1]], from[triangle[2]]),
         segmentDistanceSquared(point, from[triangle[2]], from[triangle[0]])});
    if (!hasNearest || distance < nearestDistance)
    {
      hasNearest = true;
      nearestDistance = distance;
      std::copy(triangle, triangle + 3, nearest);
      std::copy(weights, weights + 3, nearestWeights);
    }
  }
  if (!hasNearest)
  {
    return false;
  }
  // Off the mesh: extend the nearest triangle past its edges, which keeps
  // the mapping continuous as the point crosses them.
  land(nearest, nearestWeights);
  *inside = false;
  return true;
}