triangle rect overlap test

tri_rect.cpp

#define _SDL_main_h
#include <SDL/SDL.h>

#include <assert.h>
#include <stdint.h>

#include <array>

struct coord_t {
  float x, y;
};

// dot product between two 2d vectors
constexpr float dot(const coord_t &a, const coord_t &b) {
  return a.x * b.x + a.y * b.y;
}

// squared length of two 2d vectors
constexpr float len_sqr(const coord_t &a, const coord_t &b) {
  return (b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y);
}

constexpr int32_t minv(int32_t a, int32_t b) {
  return a < b ? a : b;
}

constexpr int32_t maxv(int32_t a, int32_t b) {
  return a > b ? a : b;
}

constexpr int32_t clampv(int32_t lo, int32_t v, int32_t hi) {
  return minv(hi, maxv(lo, v));
}

// return true if a triangle is backfacing
bool is_backface(const coord_t &a, const coord_t &b, const coord_t &c) {
  const float v1 = a.x - b.x;
  const float v2 = a.y - b.y;
  const float w1 = c.x - b.x;
  const float w2 = c.y - b.y;
  return 0.f > (v1 * w2 - v2 * w1);
}

// return true if a triangle intersect unit square [0,0] -> [1,1]
bool tri_vis(const coord_t &a, const coord_t &b, const coord_t &c) {

  // cohen-sutherland style trivial box clipping
  {
    const auto classify = [](const coord_t &p) -> int {
      return (p.x < 0.f ? 1 : 0) |
             (p.x > 1.f ? 2 : 0) |
             (p.y < 0.f ? 4 : 0) |
             (p.y > 1.f ? 8 : 0);
    };

    const int ca = classify(a), cb = classify(b), cc = classify(c);

    if (0 == (ca | cb | cc)) {
      // all in center, no clipping
      return true;
    }

    const int code = ca & cb & cc;
    if ((code & 1) || (code & 2) || (code & 4) || (code & 8)) {
      // all outside one plane
      return false;
    }
  }

  // reject backfaces
  // note: they mess up the plane equation test below with inverted normals
  if (is_backface(a, b, c)) {
    return false;
  }

  // triangle edge distance rejection
  {
    struct plane_t {

      plane_t(const coord_t &a, const coord_t &b)
          : nx(b.y - a.y), ny(a.x - b.x), d(a.x * nx + a.y * ny) {
          // fixme: invert to correct normals
      }

      bool test(const coord_t &p) const { return d > (p.x * nx + p.y * ny); }

      // trivial rejection based on edge normal and closest box vertex
      bool trivial() const {
        switch ((nx > 0.f) | ((ny > 0.f) << 1)) {
        case 3: return d > (1.f * nx + 1.f * ny);  // - nx  - ny
        case 2: return d > (0.f * nx + 1.f * ny);  // + nx  - ny
        case 1: return d > (1.f * nx + 0.f * ny);  // - nx  + ny
        case 0: return d > (0.f * nx + 0.f * ny);  // + nx  + ny
        }
        return false;  // keep compiler happy
      }

      float nx, ny, d;
    };

    const plane_t pab{a, b}, pbc{b, c}, pca{c, a};

    // perform trivial reject tests
    if (pab.trivial() || pbc.trivial() || pca.trivial()) {
      return false;
    }

#if 0
    const coord_t ndc[4] = {
      {0.f, 0.f}, {1.f, 0.f},
      {0.f, 1.f}, {1.f, 1.f}};

    if (pab.test(ndc[0]) &&
        pab.test(ndc[1]) &&
        pab.test(ndc[2]) &&
        pab.test(ndc[3])) {
      return false;
    }
    if (pbc.test(ndc[0]) &&
        pbc.test(ndc[1]) &&
        pbc.test(ndc[2]) &&
        pbc.test(ndc[3])) {
      return false;
    }
    if (pca.test(ndc[0]) &&
        pca.test(ndc[1]) &&
        pca.test(ndc[2]) &&
        pca.test(ndc[3])) {
      return false;
    }
#else
    return true;
#endif
  }

  // in, but needs clipping
  return true;
}

// return true if 'p' falls inside triangle 'a, b, c'
bool point_in_tri(const coord_t &a, const coord_t &b, const coord_t &c,
                  const coord_t &p) {
  // Compute vectors
  const auto v0 = coord_t{c.x - a.x, c.y - a.y}; // C - A;
  const auto v1 = coord_t{b.x - a.x, b.y - a.y}; // B - A;
  const auto v2 = coord_t{p.x - a.x, p.y - a.y}; // P - A;

  // Compute dot products
  const float dot00 = dot(v0, v0);
  const float dot01 = dot(v0, v1);
  const float dot02 = dot(v0, v2);
  const float dot11 = dot(v1, v1);
  const float dot12 = dot(v1, v2);

  // Compute barycentric coordinates
  const float denom = (dot00 * dot11 - dot01 * dot01);
  const float u = (dot11 * dot02 - dot01 * dot12) / denom;
  const float v = (dot00 * dot12 - dot01 * dot02) / denom;

  // Check if point is in triangle
  return (u >= 0) && (v >= 0) && (u + v < 1);
}

// plot a pixel to the screen
void plot(SDL_Surface *surf, int32_t x, int32_t y, uint32_t rgb = 0xdadada) {
  assert(surf);
  if (x < 0 || y < 0 || x >= surf->w || y >= surf->h) {
    return;
  }
  uint32_t *pix = (uint32_t *)surf->pixels;
  pix[x + y * surf->w] = rgb;
}

enum clip_span_t { CLIP_SPAN_MIN_X, CLIP_SPAN_MAX_X };

template <clip_span_t CLIP, size_t SIZE>
void scan_convert(coord_t a, coord_t b, std::array<int32_t, SIZE> &span) {

  // XXX: move into global attribute
  static const int32_t screen_w = 511;
  static const int32_t screen_h = 511;

  // assume our vertices are pre-sorted
  __assume(a.y < b.y);

  // scanline rejection
  if (b.y < 0.f || a.y > screen_h) {
    return;
  }

  // reject if vertices are co-linear
  if (b.y <= a.y) {
    return;
  }
  const float dx = (b.x - a.x) / (b.y - a.y);

  // clip to top of window edge
  if (a.y < 0.f) {
    a.x += dx * (0.f - a.y);
    a.y = 0.f;
  } else {
    // align to start of scanline
    const float ceily = ceilf(a.y);
    a.x += dx * (ceily - a.y);
    a.y = ceily;
  }

  const int32_t iay = maxv(int32_t(a.y), 0);
  const int32_t iby = minv(int32_t(b.y), screen_h);

  int32_t x = int32_t(a.x * float(0x10000));
  const int32_t idx = int32_t(dx * float(0x10000));

  switch (CLIP) {
  case CLIP_SPAN_MAX_X:
    for (int32_t y = iay; y < iby; ++y, x += idx) {
      span[y] = std::max<int32_t>(x >> 16, 0);
    }
    break;
  case CLIP_SPAN_MIN_X:
    for (int32_t y = iay; y < iby; ++y, x += idx) {
      span[y] = std::min<int32_t>(x >> 16, screen_w);
    }
    break;
  }
}

// scan convert a triangle
bool scan_triangle(SDL_Surface *surf, std::array<coord_t, 3> v) {

  // sort vertices: top (0), mid (1), bottom (2)
  if (v[1].y < v[0].y)
    std::swap(v[1], v[0]);
  if (v[2].y < v[0].y)
    std::swap(v[2], v[0]);
  if (v[2].y < v[1].y)
    std::swap(v[2], v[1]);

  // check mid vertex side
  const float nx = v[2].y - v[0].y;
  const float ny = v[0].x - v[2].x;
  const float d1 = v[0].x * nx + v[0].y * ny;
  const float d2 = v[1].x * nx + v[1].y * ny;

  if (d1 == d2) {
    // coplanar, so reject triangle
    return false;
  }

  // our y axis span buffers
  std::array<int32_t, 512> lo, hi;

  // scan convert edges
  if (d1 > d2) {
    scan_convert<CLIP_SPAN_MIN_X>(v[0], v[2], hi);
    scan_convert<CLIP_SPAN_MAX_X>(v[0], v[1], lo);
    scan_convert<CLIP_SPAN_MAX_X>(v[1], v[2], lo);
  } else {
    scan_convert<CLIP_SPAN_MAX_X>(v[0], v[2], lo);
    scan_convert<CLIP_SPAN_MIN_X>(v[0], v[1], hi);
    scan_convert<CLIP_SPAN_MIN_X>(v[1], v[2], hi);
  }

  // fill triangle
  {
    const int32_t y0 = std::max(int32_t(ceilf(v[0].y)), 0);
    const int32_t y1 = std::min(int32_t(v[2].y), 511);

    uint32_t *py = (uint32_t *)surf->pixels;
    py += (y0 * surf->pitch) / 4;
    for (int32_t y = y0; y < y1; ++y) {
      // step to edge
      uint32_t *px = py + lo[y];
      // raster scanline
      for (int32_t x = lo[y]; x < hi[y]; ++x, ++px) {
        *px = 0x335577;
      }
      // step scanline
      py += surf->pitch / 4;
    }
  }

  return true;
}

// fixed point line drawing
void draw_line(SDL_Surface *surf, coord_t a, coord_t b, uint32_t rgb) {
  const float dx = b.x - a.x, dy = b.y - a.y;
  const float adx = fabsf(dx), ady = fabs(dy);
  if (fabsf(dx) > fabsf(dy)) {
    if (b.x < a.x)
      std::swap(a, b);
    const int32_t iy = int32_t(float(0x10000) * ((b.y - a.y) / adx));
    int32_t y = int32_t(a.y * float(0x10000));
    for (int32_t x = int32_t(a.x); x < int32_t(b.x); ++x, y += iy) {
      plot(surf, x, y >> 16, rgb);
    }
  } else {
    if (b.y < a.y)
      std::swap(a, b);
    const int32_t ix = int32_t(float(0x10000) * ((b.x - a.x) / ady));
    int32_t x = int32_t(a.x * float(0x10000));
    for (int32_t y = int32_t(a.y); y < int32_t(b.y); ++y, x += ix) {
      plot(surf, x >> 16, y, rgb);
    }
  }
}

// draw a wireframe triangle
void draw_tri(SDL_Surface *surf, const std::array<coord_t, 3> &t,
              uint32_t rgb) {

  for (uint32_t j = 0; j < 3; ++j) {
    const coord_t &a = t[j];
    const coord_t &b = t[j == 2 ? 0 : j + 1];
    draw_line(surf, a, b, rgb);
  }
}

// find the nearest vertex to point p
coord_t *nearest(std::array<coord_t, 3> &tri, const coord_t &p) {
  coord_t *out = nullptr;
  float best = 100000.f;
  for (auto &c : tri) {
    const float dx = p.x - c.x;
    const float dy = p.y - c.y;
    const float ds = dx * dx + dy * dy;
    if (ds < best) {
      best = ds;
      out = &c;
    }
  }
  return best < 256.f ? out : nullptr;
}

// program entry
int main(const int argc, const char **args) {

  if (SDL_Init(SDL_INIT_VIDEO)) {
    return 1;
  }

  SDL_Surface *surf = SDL_SetVideoMode(512, 512, 32, 0);
  if (!surf) {
    return 2;
  }

  std::array<coord_t, 3> tri = {coord_t{256, 64}, coord_t{64, 512 - 64},
                                coord_t{512 - 64, 512 - 64}};
  coord_t *drag = nullptr;
  bool dragall = false;

  bool active = true;
  while (active) {

    SDL_FillRect(surf, nullptr, 0x101010);

    SDL_Event event;
    while (SDL_PollEvent(&event)) {
      switch (event.type) {
      case SDL_QUIT:
        active = false;
        break;
      case SDL_MOUSEBUTTONUP:
        drag = nullptr;
        dragall = false;
        break;
      case SDL_MOUSEBUTTONDOWN:
        if (!drag) {
          const coord_t p =
              coord_t{float(event.button.x), float(event.button.y)};
          // check for drags
          drag = nearest(tri, p);
          if (!drag) {
            dragall = point_in_tri(tri[0], tri[1], tri[2], p);
          }
          // clear relative mouse state
          SDL_GetRelativeMouseState(nullptr, nullptr);
        }
        break;
      }
    }

    if (drag || dragall) {
      int rx, ry;
      SDL_GetRelativeMouseState(&rx, &ry);
      if (drag) {
        drag->x += rx;
        drag->y += ry;
      }
      if (dragall) {
        for (auto &c : tri) {
          c.x += rx;
          c.y += ry;
        }
      }
    }

    uint32_t tri_rgb = 0xff6040;

    // draw screen
    {
      SDL_Rect rect = {(512 - 256) / 2, (512 - 256) / 2, 256, 256};
      SDL_FillRect(surf, &rect, 0x203040);

      std::array<coord_t, 3> xtri = tri;
      for (coord_t &c : xtri) {
        c.x = (c.x - rect.x) / rect.w;
        c.y = (c.y - rect.y) / rect.h;
      }

      if (tri_vis(xtri[0], xtri[1], xtri[2])) {
        tri_rgb = 0x30ff60;
      }

      if (is_backface(tri[0], tri[1], tri[2])) {
        tri_rgb = 0xdadada;
      }
    }

    // draw triangle
    scan_triangle(surf, tri);
    draw_tri(surf, tri, tri_rgb);

    // draw vertex nodes
    for (const auto &p : tri) {
      SDL_Rect r = {int16_t(p.x - 3), int16_t(p.y - 3), 7, 7};
      SDL_FillRect(surf, &r, 0xdadada);
    }

    SDL_Flip(surf);
    SDL_Delay(1);
  }

  return 0;
}