2 files changed, 201 insertions, 44 deletions

diff --git a/include/swgfx.h b/include/swgfx.h
index 9905382..cdeefe8 100644
--- a/include/swgfx.h
+++ b/include/swgfx.h
@@ -34,11 +34,13 @@ typedef sgVec3 sgNormal;
 typedef struct sgVert2i { sgVec2i pos; sgVec2 uv; } sgVert2i;
 typedef struct sgVert2  { sgVec2  pos; sgVec2 uv; } sgVert2;
 typedef struct sgVert3  { sgVec3  pos; sgVec2 uv; } sgVert3;
+typedef struct sgVert4  { sgVec4  pos; sgVec2 uv; } sgVert4;
 
 typedef struct sgQuadi { sgVert2i p0, p1; }     sgQuadi;
 typedef struct sgQuad  { sgVert2  p0, p1; }     sgQuad;
 typedef struct sgTri2  { sgVert2  p0, p1, p2; } sgTri2;
 typedef struct sgTri3  { sgVert3  p0, p1, p2; } sgTri3;
+typedef struct sgTri4  { sgVert4  p0, p1, p2; } sgTri4;
 
 typedef uint16_t sgIdx;
 typedef struct   sgVertIdx { sgIdx pos, uv, normal; } sgVertIdx;
diff --git a/src/swgfx.c b/src/swgfx.c
index c54fdc1..2262dd8 100644
--- a/src/swgfx.c
+++ b/src/swgfx.c
@@ -44,39 +44,44 @@ typedef struct swgfx {
   // Make it so that changing the model matrix only requires one matrix
   // multiplication (mvp = model * viewProj) and not two (mvp = model * view * projection)
   // before rendering the model's triangles.
-  sgMat4       viewProj; // View-projection matrix.
-  sgMat4       mvp;      // Model-view-projection matrix.
-  const sgImage* texture;// User-specified texture.
+  sgMat4         viewProj;       // View-projection matrix.
+  sgMat4         mvp;            // Model-view-projection matrix.
+  const sgImage* texture;        // User-specified texture.
   sgImage        defaultTexture; // A default for when no texture is provided.
-  sgPixel            defaultPixel;   // The single-pixel of the default texture.
+  sgPixel        defaultPixel;   // The single-pixel of the default texture.
 } swgfx;
 
 static inline int mod(int a, int m) { return (m + (a % m)) % m; }
+static inline R   frac(R a) { return a - (R)((int)a); }
 
-static inline R frac(R a) { return a - (R)((int)a); }
-
-static inline R   rmin(R   a, R   b) { return (a <= b) ? a : b; }
-static inline R   rmax(R   a, R   b) { return (a >= b) ? a : b; }
 static inline int imin(int a, int b) { return (a <= b) ? a : b; }
 static inline int imax(int a, int b) { return (a >= b) ? a : b; }
+static inline R   rmin(R   a, R   b) { return (a <= b) ? a : b; }
+static inline R   rmax(R   a, R   b) { return (a >= b) ? a : b; }
+static inline R   lerp(R   a, R   b, R t) { return a + t*(b-a); }
 
-static inline sgVec2  min2(sgVec2 a, sgVec2 b) { return (sgVec2){.x = rmin(a.x, b.x), .y = rmin(a.y, b.y) }; }
-static inline sgVec2  max2(sgVec2 a, sgVec2 b) { return (sgVec2){.x = rmax(a.x, b.x), .y = rmax(a.y, b.y) }; }
 static inline sgVec2i min2i(sgVec2i a, sgVec2i b) { return (sgVec2i){.x = imin(a.x, b.x), .y = imin(a.y, b.y) }; }
 static inline sgVec2i max2i(sgVec2i a, sgVec2i b) { return (sgVec2i){.x = imax(a.x, b.x), .y = imax(a.y, b.y) }; }
+static inline sgVec2  min2(sgVec2 a, sgVec2 b) { return (sgVec2){.x = rmin(a.x, b.x), .y = rmin(a.y, b.y) }; }
+static inline sgVec2  max2(sgVec2 a, sgVec2 b) { return (sgVec2){.x = rmax(a.x, b.x), .y = rmax(a.y, b.y) }; }
+static inline sgVec2  add2(sgVec2 a, sgVec2 b) { return (sgVec2){a.x + b.x, a.y + b.y}; }
+static inline sgVec2  sub2(sgVec2 a, sgVec2 b) { return (sgVec2){a.x - b.x, a.y - b.y}; }
+static inline sgVec2  scale2(sgVec2 v, R s) { return (sgVec2){v.x * s, v.y * s}; }
+static inline sgVec2  frac2(sgVec2 v) { return (sgVec2){frac(v.x), frac(v.y)}; }
+static inline sgVec2  lerp2(sgVec2 a, sgVec2 b, R t) { return add2(a, scale2(sub2(b,a), t)); }
 
-static inline sgVec2 frac2(sgVec2 v) { return (sgVec2){frac(v.x), frac(v.y)}; }
-
-static inline sgVec2 add2(sgVec2 a, sgVec2 b) { return (sgVec2){a.x + b.x, a.y + b.y}; }
-static inline sgVec2 scale2(sgVec2 v, R s) { return (sgVec2){v.x * s, v.y * s}; }
-
+static inline sgVec3 add3(sgVec3 a, sgVec3 b) { return (sgVec3){a.x + b.x, a.y + b.y, a.z + b.z}; }
 static inline sgVec3 neg3(sgVec3 v) { return (sgVec3){-v.x, -v.y, -v.z}; }
 static inline sgVec3 sub3(sgVec3 a, sgVec3 b) { return (sgVec3){a.x - b.x, a.y - b.y, a.z - b.z}; }
 static inline sgVec3 div3(sgVec3 a, sgVec3 b) { return (sgVec3){a.x / b.x, a.y / b.y, a.z / b.z}; }
+static inline sgVec3 scale3(sgVec3 v, R s) { return (sgVec3){v.x * s, v.y * s, v.z * s}; }
 static inline R dot3(sgVec3 a, sgVec3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
 static inline R normsq3(sgVec3 v) { return v.x * v.x + v.y * v.y + v.z * v.z; }
 static inline R norm3  (sgVec3 v) { return (R)sqrt(normsq3(v)); }
 
+static inline sgVec4 add4(sgVec4 a, sgVec4 b) { return (sgVec4){a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w}; }
+static inline sgVec4 sub4(sgVec4 a, sgVec4 b) { return (sgVec4){a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w}; }
+static inline sgVec4 scale4(sgVec4 v, R s) { return (sgVec4){v.x * s, v.y * s, v.z * s, v.w * s}; }
 static inline sgVec4 lerp4(sgVec4 a, sgVec4 b, R t) {
   return (sgVec4){
     .x = a.x + t * (b.x - a.x),
@@ -85,6 +90,14 @@ static inline sgVec4 lerp4(sgVec4 a, sgVec4 b, R t) {
     .w = a.w + t * (b.w - a.w)};
 }
 
+/// Return the curl of 'a' towards 'b', which is defined as the z-coordinate of
+/// the cross product a x b, or as the determinant det(a,b).
+///
+/// The curl of 'a' towards 'b' is positive if 'a' curls towards 'b' like the
+/// positive x-axis curls towards the positive y-axis.
+static inline R curl2(sgVec2 a, sgVec2 b) {
+  return (a.x * b.y) - (a.y * b.x);
+}
 static inline sgVec3 cross3(sgVec3 a, sgVec3 b) {
   return (sgVec3) {
     a.y * b.z - a.z * b.y,
@@ -97,6 +110,7 @@ static inline sgVec3 normalize3(sgVec3 v) {
   return (n > 0) ? (sgVec3){v.x / n, v.y / n, v.z / n} : (sgVec3){0, 0, 0};
 }
 
+static inline sgVec2 Vec2FromVec4(sgVec4 v) { return (sgVec2){v.x, v.y}; }
 static inline sgVec4 Vec4FromVec3(sgVec3 v, R w) { return (sgVec4){v.x, v.y, v.z, w}; }
 
 static inline sgMat4 Mat4(
@@ -240,15 +254,15 @@ static inline sgMat4 Mat4Look(sgVec3 position, sgVec3 forward, sgVec3 up) {
 }
 
 static inline sgMat4 Mat4Perspective(R fovy, R aspect, R near, R far) {
-  R f = (R)tan(fovy / 2.0);
-  assert(f > 0.0);
-  f = 1.f / f;
+  assert(fovy > 0.f);
+  assert(near < far);
+  const R f = 1.f / tanf(fovy / 2.f);
   const R a = near - far;
   return Mat4(
-    f / aspect, 0,                 0,                    0,
-    0,          f,                 0,                    0,
+    f / aspect, 0,                 0,                   0,
+    0,          f,                 0,                   0,
     0,          0,  (far + near) / a, (2 * far * near / a),
-    0,          0,                 -1,                  0);
+    0,          0,                -1,                   0);
 }
 
 static inline sgVec4 PixelToVec4(sgPixel p) {
@@ -341,8 +355,10 @@ static inline sgPixel SampleBilinear(const sgImage* texture, sgVec2 uv) {
 // TODO: Clamping and other addressing strategies.
 static inline sgPixel Sample(const sgImage* texture, sgVec2 uv) {
   // TODO: Add a member to sgImage that determines how it should be filtered.
-  //return SampleNearest(texture, uv);
-  return SampleBilinear(texture, uv);
+  return SampleNearest(texture, uv);
+  // TODO: Debug bilinear. It's producing random colours on the ground, likely
+  //  out-of-bounds memory accesses.
+  //return SampleBilinear(texture, uv);
 }
 
 static inline sgAABB2 TriangleAabb2(sgVec2 p0, sgVec2 p1, sgVec2 p2) {
@@ -437,29 +453,157 @@ static inline sgVec4 ViewportToWindow(sgViewport_t vp, sgVec4 p) {
   return (sgVec4){p.x, (R)vp.height - p.y, p.z, p.w};
 }
 
-static inline sgVec4 TransformPosition(const swgfx* gfx, sgVec3 p) {
+/// Line segment-plane intersection special-case for the near camera plane.
+/// All quantities assumed to be in camera space.
+/// outP = a + outT*(b-a)
+static inline R IntersectSegmentPlane(R near, const sgVec3* const a, const sgVec3* const b) {
+  // D = near plane distance = perpendicular distance from the origin to the plane.
+  // o = line origin = a
+  // d = line direction = b-a
+  // Plane normal = (0, 0, +1) --- Could be -1, need to be consistent with D.
+  // Point in plane: p=(0, 0, -near)
+  //   <=> p dot n     + D = 0
+  //   === -near * n.z + D = 0
+  //   === -near * 1   + D = 0
+  //   === D = near
+  // Denominator = n dot d = (0,0,1) dot d = d.z = (b.z - a.z)
+  const R t = (-near - a->z) / (b->z - a->z);
+  assert(t >= 0.f);
+  assert(t <= 1.f);
+  return t;
+}
+
+/// Interpolate depth and vertex attributes at the in/out vertex 'out'.
+static void InterpolateAttributes(const sgVert4* const a, const sgVert4* const b, R t, sgVert4* out) {
+  assert(a);
+  assert(b);
+  assert(out);
+  assert(t >= 0.f);
+  assert(t <= 1.f);
+  const sgVec4 d = sub4(b->pos, a->pos); // Line direction.
+  out->pos = add4(a->pos, scale4(d, t));
+  out->uv  = lerp2(a->uv, b->uv, t);
+}
+
+/// Clip a triangle, vertices in clip space. Return the number of output
+/// triangles.
+///
+/// 4 possible cases:
+/// 1. All vertices in front of the camera near plane => draw.
+/// 2. All vertices behind   => discard.
+/// 3. One vertex in front   => draw 1 clipped triangle.
+/// 4. Two vertices in front => draw 2 clipped triangles.
+static inline int ClipTriangle(R near, const sgTri4* const tri, sgTri4 out[2]) {
+#define VERTEX(IDX) (&tri->p0)[IDX]
+#define VALID(X) ((0 <= (X)) && ((X) < 3))
+#define IN_FRONT(P) (P.z >= -near) // +Z points into the screen in clip space.
+  const bool f[3] = {IN_FRONT(tri->p0.pos), IN_FRONT(tri->p1.pos), IN_FRONT(tri->p2.pos)};
+  const int numFront = f[0] + f[1] + f[2];
+  int numTris;
+  if (numFront == 3) {
+    numTris = 1;
+    out[0] = *tri;
+  } else if (numFront == 2) {
+    numTris = 2;
+    int back = 0;
+    for (; f[back] && (back < 3); ++back) {}
+    assert(VALID(back));
+    assert(!f[back]);
+    int front[2] = {(back+1)%3, (back+2)%3};
+    assert(VALID(front[0]));
+    assert(VALID(front[1]));
+    const sgVert4* const backVert = &VERTEX(back);
+    sgVert4 p[2];
+    for (int i = 0; i < 2; ++i) {
+      const R t = IntersectSegmentPlane(near, (const sgVec3*)&backVert->pos, (const sgVec3*)&VERTEX(front[i]).pos);
+      InterpolateAttributes(backVert, &VERTEX(front[i]), t, &p[i]);
+    }
+    // We must preserve the winding order here for culling.
+    // Note that p[i] corresponds to front[i] = back+(i+1).
+    out[0] = (sgTri4){p[1], p[0], VERTEX(front[1])};
+    out[1] = (sgTri4){p[0], VERTEX(front[0]), VERTEX(front[1])};
+  } else if (numFront == 1) {
+    numTris = 1;
+    int front = 0;
+    for (; !f[front] && (front < 3); ++front){}
+    assert(VALID(front));
+    assert(f[front]);
+    int back[2] = {(front+1)%3, (front+2)%3};
+    assert(VALID(back[0]));
+    assert(VALID(back[1]));
+    const sgVert4* const frontVert = &VERTEX(front);
+    sgVert4 p[2];
+    for (int i = 0; i < 2; ++i) {
+      const R t = IntersectSegmentPlane(near, (const sgVec3*)&frontVert->pos, (const sgVec3*)&VERTEX(back[i]).pos);
+      InterpolateAttributes(frontVert, &VERTEX(back[i]), t, &p[i]);
+    }
+    // We must preserve the winding order here for culling.
+    // Note that p[i] corresponds to back[i] = front+(i+1).
+    out[0] = (sgTri4){*frontVert, p[0], p[1]};
+  } else {
+    numTris = 0;
+  }
+  return numTris;
+#undef IN_FRONT
+#undef VALID
+#undef VERTEX
+}
+
+static inline int TransformTri(const swgfx* gfx, const sgTri3* const tri, sgTri4 out[2]) {
   assert(gfx);
+  assert(tri);
   // Model to clip space.
-  const sgVec4 p_clip = Mat4MulVec4(gfx->mvp, Vec4FromVec3(p, 1));
-  // TODO: Backface culling.
-  // Perspective divide.
-  const sgVec4 p_ndc = PerspDivide(p_clip);
-  // TODO: Clip.
-  const sgVec4 p_vp = ViewportTransform(gfx->viewport, p_ndc);
-  return ViewportToWindow(gfx->viewport, p_vp);
+  const sgVec4 p0_clip = Mat4MulVec4(gfx->mvp, Vec4FromVec3(tri->p0.pos, 1));
+  const sgVec4 p1_clip = Mat4MulVec4(gfx->mvp, Vec4FromVec3(tri->p1.pos, 1));
+  const sgVec4 p2_clip = Mat4MulVec4(gfx->mvp, Vec4FromVec3(tri->p2.pos, 1));
+  const sgTri4 tri_clip = {
+    (sgVert4){ p0_clip, tri->p0.uv },
+    (sgVert4){ p1_clip, tri->p1.uv },
+    (sgVert4){ p2_clip, tri->p2.uv }};
+  // Clip.
+  // Our perspective matrix maps the near plane to z=-1 in clip space.
+  constexpr R near_clip = -1.f;
+  const int numTris = ClipTriangle(near_clip, &tri_clip, out);
+  assert((0 <= numTris) && (numTris <= 2));
+  for (int i = 0; i < numTris; ++i) {
+    sgTri4* const tri4 = &out[i];
+    // Perspective divide.
+    const sgVec4 p0_ndc = PerspDivide(tri4->p0.pos);
+    const sgVec4 p1_ndc = PerspDivide(tri4->p1.pos);
+    const sgVec4 p2_ndc = PerspDivide(tri4->p2.pos);
+    // To viewport.
+    const sgVec4 p0_vp = ViewportTransform(gfx->viewport, p0_ndc);
+    const sgVec4 p1_vp = ViewportTransform(gfx->viewport, p1_ndc);
+    const sgVec4 p2_vp = ViewportTransform(gfx->viewport, p2_ndc);
+    // To window.
+    const sgVec4 p0_wn = ViewportToWindow(gfx->viewport, p0_vp);
+    const sgVec4 p1_wn = ViewportToWindow(gfx->viewport, p1_vp);
+    const sgVec4 p2_wn = ViewportToWindow(gfx->viewport, p2_vp);
+    // Output.
+    tri4->p0.pos = p0_wn;
+    tri4->p1.pos = p1_wn;
+    tri4->p2.pos = p2_wn;
+  }
+  return numTris;
 }
 
-static void DrawTriangle3(swgfx* gfx, const sgTri3* const tri) {
+static void DrawTriangle3PostClip(swgfx* gfx, const sgTri4* const tri) {
   assert(gfx);
   assert(tri);
-  // TODO: Inline the transform here and interleave its operations to perform
-  //  backface culling and clipping as early as possible.
-  const sgVec4 p0    = TransformPosition(gfx, tri->p0.pos);
-  const sgVec4 p1    = TransformPosition(gfx, tri->p1.pos);
-  const sgVec4 p2    = TransformPosition(gfx, tri->p2.pos);
+  const sgVec4 p0 = tri->p0.pos;
+  const sgVec4 p1 = tri->p1.pos;
+  const sgVec4 p2 = tri->p2.pos;
   const sgVec2 p0_2d = (sgVec2){p0.x, p0.y};
   const sgVec2 p1_2d = (sgVec2){p1.x, p1.y};
   const sgVec2 p2_2d = (sgVec2){p2.x, p2.y};
+  // Backface culling, assume front face = ccw.
+  // In screen space, +Y goes down.
+  // p0p1p2 is ccw <=> p0p1 curls negatively towards p0p2. If the curl is
+  // positive (cw winding), cull.
+  if (curl2(sub2(p1_2d, p0_2d),
+            sub2(p2_2d, p0_2d)) > 0.f) {
+    return;
+  }
   const sgAABB2 bbox = TriangleAabb2(p0_2d, p1_2d, p2_2d);
   // We consider (x,y) to be the pixel center.
   // Draw all pixels touched by the bounding box. TODO: Multi-sampling.
@@ -497,13 +641,12 @@ static void DrawTriangle3(swgfx* gfx, const sgTri3* const tri) {
         if ((0.f <= p_depth) && (p_depth <= 1.f) && (p_depth <= *depth)) {
           *depth = p_depth;
           const sgPixel colour = Sample(gfx->texture, uv);
-          // TODO: When doing lighting, need to tone-map here.
-          /*const int d = (int)(z*255.f);
-          const sgPixel colour = (sgPixel){d,d,d,255};*/
-          //const sgPixel colour = (sgPixel){255, 0, 255, 255};
-          /*const int r = (int)(uv.x * 255.f);
-          const int g = (int)(uv.y * 255.f);
-          const sgPixel colour = (sgPixel){r, g, 255, 255};*/
+          // TODO: When doing lighting, need to tone-map here and apply inverse
+          //  gamma here.
+          //const sgPixel colour = {(uint8_t)(bar.x*255.f), (uint8_t)(bar.y*255.f), (uint8_t)(bar.z*255.f), 255};
+          //const sgPixel colour = {(int)(z*255.f), (int)(z*255.f), (int)(z*255.f), 255};
+          //const sgPixel colour = {255, 0, 255, 255};
+          //const sgPixel colour = {(int)(uv.x * 255.f), (int)(uv.y * 255.f), 255, 255};
           SetPixel(gfx, (sgVec2i){x,y}, colour);
         }
       }
@@ -511,6 +654,17 @@ static void DrawTriangle3(swgfx* gfx, const sgTri3* const tri) {
   }
 }
 
+static void DrawTriangle3(swgfx* gfx, const sgTri3* const tri) {
+  assert(gfx);
+  assert(tri);
+  sgTri4 tris[2];
+  const int numTris = TransformTri(gfx, tri, tris);
+  assert((0 <= numTris) && (numTris <= 2));
+  for (int i = 0; i < numTris; ++i) {
+    DrawTriangle3PostClip(gfx, &tris[i]);
+  }
+}
+
 #define is_pow2_or_0(X) ((X & (X - 1)) == 0)
 #define SG_ALIGN 64
 #define SG_ALLOC(PP_MEM, COUNT, TYPE) (TYPE*)Alloc(PP_MEM, COUNT, sizeof(TYPE))
@@ -557,6 +711,7 @@ swgfx* sgNew(int width, int height, void* mem) {
     .height = 1,
     .pixels = &gfx->defaultPixel,
   };
+  gfx->texture = &gfx->defaultTexture;
   return gfx;
 }