Simple watertight triangle rasterizer

rast.c

// ---- triangle rasterizer

#define SUBPIXEL_SHIFT  8
#define SUBPIXEL_SCALE  (1 << SUBPIXEL_SHIFT)

static RADINLINE S64 det2x2(S32 a, S32 b, S32 c, S32 d)
{
   S64 r = (S64) a*d - (S64) b*c;
   return r >> SUBPIXEL_SHIFT;
}

static S64 det2x2_fill_convention(S32 a, S32 b, S32 c, S32 d)
{
   S64 r = (S64) a*d - (S64) b*c;         // the determinant
   if (c > 0 || c == 0 && a <= 0) r--;    // this implements the top-left fill convention
   return r >> SUBPIXEL_SHIFT;
}

static RADINLINE S32 fixed(F32 x)
{
   // -0.5f to place pixel centers at integer coords + 0.5
   // +0.5f afterwards is rounding factor
   return (S32) ((x - 0.5f) * SUBPIXEL_SCALE + 0.5f);
}

static RADINLINE S32 fixed_ceil(S32 x)
{
   return (x + SUBPIXEL_SCALE - 1) >> SUBPIXEL_SHIFT;
}

typedef S32 EdgeDist; // switch to S64 if there are overflows

static void point_to_fixed_xform(S32 *ox, S32 *oy, gswf_point *in, gswf_matrix *m)
{
   F32 x = in->x * m->m00 + (in->y *  m->m01 + m->trans[0]);
   F32 y = in->x * m->m10 + (in->y *  m->m11 + m->trans[1]);
   *ox = fixed(x);
   *oy = fixed(y);
}

static void draw_triangle(image *tgt, gswf_point *a, gswf_point *b, gswf_point *c, gswf_matrix *m, int allow_neg_winding)
{
   S32 x0,x1,x2;
   S32 y0,y1,y2;
   S32 minx,miny,maxx,maxy;
   S32 dx10,dx21,dx02;
   S32 dy10,dy21,dy02;
   EdgeDist e1,e2,e3;
   S32 t;
   S32 x,y;
   S64 det;
   int incr;

   // convert coordinates to fixed point
   point_to_fixed_xform(&x0, &y0, a, m);
   point_to_fixed_xform(&x1, &y1, b, m);
   point_to_fixed_xform(&x2, &y2, c, m);

   // check triangle winding order
   det = det2x2(x1-x0, x2-x0, y1-y0, y2-y0);
   if (det > 0)
      incr = 1;
   else if (det < 0) {
      incr = allow_neg_winding ? -1 : 1;
      t = x0; x0 = x1; x1 = t;
      t = y0; y0 = y1; y1 = t; 
   } else // zero-area triangle
      return;

   // bounding box / clipping
   minx = max(fixed_ceil(min3(x0,x1,x2)), 0);
   miny = max(fixed_ceil(min3(y0,y1,y2)), 0);
   maxx = min(fixed_ceil(max3(x0,x1,x2)), tgt->w);
   maxy = min(fixed_ceil(max3(y0,y1,y2)), tgt->h);
   if (minx >= maxx || miny >= maxy)
      return;

   // edge vectors
   dx10 = x1 - x0; dy10 = y1 - y0;
   dx21 = x2 - x1; dy21 = y2 - y1;
   dx02 = x0 - x2; dy02 = y0 - y2;

   // edge functions
   e1 = (EdgeDist) det2x2_fill_convention(dx10, (minx << SUBPIXEL_SHIFT) - x0, dy10, (miny << SUBPIXEL_SHIFT) - y0);
   e2 = (EdgeDist) det2x2_fill_convention(dx21, (minx << SUBPIXEL_SHIFT) - x1, dy21, (miny << SUBPIXEL_SHIFT) - y1);
   e3 = (EdgeDist) det2x2_fill_convention(dx02, (minx << SUBPIXEL_SHIFT) - x2, dy02, (miny << SUBPIXEL_SHIFT) - y2);

   // rasterize
   for (y=miny; y<maxy; y++) {
      U8 *line = tgt->data + y*tgt->w;
      EdgeDist ei1 = e1, ei2 = e2, ei3 = e3;

      for (x=minx; x<maxx; x++) {
         if ((ei1 | ei2 | ei3) >= 0) // pixel in triangle
            line[x] = (U8) (line[x] + incr);

         ei1 -= dy10;
         ei2 -= dy21;
         ei3 -= dy02;
      }

      e1 += dx10;
      e2 += dx21;
      e3 += dx02;
   }
}