1

I have been trying to implement the optimized rasterizer outlined in this blog: https://fgiesen.wordpress.com/2013/02/10/optimizing-the-basic-rasterizer/.

The naive approach outlined in his prior blog post https://fgiesen.wordpress.com/2013/02/08/triangle-rasterization-in-practice/, calculates the determinants (for barycentric weights) at each pixel. But his optimized version takes advantage of the fact that for three points a, b, c, the determinant function

(b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x) 

can be rewritten as

A * c.x + B * c.y + C 

for

A = (a.y - b.y), B = (b.x - a.x), C = a.x * b.y - a.y * b.x.  

Since during traversal, the points a, b are two points on the triangle, c is the only point whose values change.

For any edge function E which outputs the weight for the corresponding triangle partition made by it and the point c, we can express an x difference and y difference in discrete steps.

E(c.x, c.y) = A * c.x + B * c.y + C

so

E(c.x + 1, c.y) - E(c.x, c.y) = A, and E(c.x, c.y + 1) - E(c.x, c.y) = B

So at each iteration instead of recalculating the determinant, we can just find the three determinants for the first c, and then increment them the A or B which corresponds to their edge.

I'm currently trying to implement this in my own rasterizer, but I quickly noticed an issue with the formula. My triangles are fed to the draw function in screen space, so a low y value means high up whereas in vector space it means low down. I thought I would account for this by multiplying every y value in the formula by -1.

This gave me the following formula for the determinant of a, b, c:

(b.x - a.x) * (a.y - c.y) - (a.y - b.y) * (c.x - a.x)

From which I derived the following A, B, C:

A = b.y - a.y, B = a.x - b.x, C = b.x * a.y - a.x * b.y

In my tests, using this new determinant formula (calculated at every pixel) works fine. And for the first point c in traversal, it is equivalent to

A * c.x + B * c.y + C

But as it continues traversing along the triangle's bounding box, the step incremented determinant values go out of sync with the raw calculated determinant values. Somehow this means that the step sizes A and B are faulty-- which makes no sense to me.

The only two causes of this problem I can think of are either I calculated A B and C incorrectly, or I am not mapping from vector space to screen space in a way that preserves area or orientation.

But just in case, here is all of my code for the rasterizer:

helpers

typedef float* point_t;
typedef float* triangle_t;

typedef struct edge {
  point_t tail, tip;
  float step_x, step_y;
  int is_top_left;
} edge_t;

/* ... */

/* tail is the begining of the edge (i.e a), tip is the end of the edge (b) and c is variable */

static float init_edge(edge_t* edge, point_t tail, point_t tip, point_t origin) {

    edge->tail = tail;
    edge->tip = tip;
    edge->is_top_left = is_top_left(tail, tip);

    float A = tip[1] - tail[1];
    float B = tail[0] - tip[0];
    float C = tip[0] * tail[1] - tail[0] * tip[1];

    /* step sizes */
    edge->step_x = A;
    edge->step_y = B;

    /* edge function output at origin */
    return A * origin[0] + B * origin[1] + C;
 }

static float det(point_t a, point_t b, point_t c) {
   return (b[0] - a[0]) * (a[1] - c[1]) - (a[1] - b[1]) * (c[0] - a[0]);
}

draw_triangle

void draw_triangle(sr_pipeline_t* pipeline, triangle_t triangle) {

    /* orient triangle ccw */

    point_t v0 = (point_t)malloc(sizeof(float) * pipeline->num_attr);
    point_t v1 = (point_t)malloc(sizeof(float) * pipeline->num_attr);
    point_t v2 = (point_t)malloc(sizeof(float) * pipeline->num_attr);

    memcpy(v0, triangle, sizeof(float) * pipeline->num_attr);
    memcpy(v1, triangle + pipeline->num_attr, sizeof(float) * pipeline->num_attr);
    memcpy(v2, triangle + (2 * pipeline->num_attr), sizeof(float) * pipeline->num_attr);

    orient_ccw(&v0, &v1, &v2);

    /* find bounding box */

    float min_x = /* ... */;
    float min_y = /* ... */;
    float max_x = /* ... */;
    float max_y = /* ... */;

    /* store current point */

    point_t p = (point_t)calloc(pipeline->num_attr, sizeof(float));
    p[0] = min_x;
    p[1] = min_y;

    /* grab edge information */

   edge_t e01, e12, e20;
   float w0 = init_edge(&e12, v1, v2, p);
   float w1 = init_edge(&e20, v2, v0, p);
   float w2 = init_edge(&e01, v0, v1, p);

   /* rasterize */

  for (p[1] = min_y; p[1] <= max_y; p[1]++) {
    for (p[0] = min_x; p[0] <= max_x; p[0]++) {

      /* determinant calculated at every step (I suspect these are correct) */
      float s0 = det(v1, v2, p);
      float s1 = det(v2, v0, p);
      float s2 = det(v0, v1, p);

      if ( (s0 >= 0) && (s1 >= 0) && (s2 >= 0) ) {
        draw_point(pipeline, p);
      }

      w0 += e12.step_x;
      w1 += e20.step_x;
      w2 += e01.step_x;
   } 
   w0 += e12.step_y;
   w1 += e20.step_y;
   w2 += e01.step_y;
 }

 free(v0);
 free(v1);
 free(v2);
 free(p);
}

Code and functions that I have omitted I have verified work correctly.

To reiterate, my question is why are the values w0, w1, w2 not the same as s0, s1, s2 as they should be?

Any help is appreciated, thank you!

aurreco
  • 21
  • 3
  • just a sanity check Hope each triangle is rendered by multiple threads ... if not then this will be slow as barycentric overhead is far bigger (even with your optimization) then scan line based [convex polygon rasterizer](https://stackoverflow.com/a/19078088/2521214) as the BBOX passes through a lot of unused pixels while scan line rasterizer iterates only through used pixels... Also I would use [matrix approach](https://stackoverflow.com/a/37431636/2521214) and compute it just once per triangle, then use your optimization on the result instead of multiplication as you want to do now... – Spektre Jun 13 '22 at 18:24
  • No, for now I plan to have it single threaded just to get it working (threading in c is not something I've done outside of school projects). But I prefer this algorithm more since its a lot clearer to find interpolated attribute values per pixel. Also, in the future I may try to further optimize the bounding box traversal like mentioned in the ryg blog http://people.csail.mit.edu/ericchan/bib/pdf/p15-mccormack.pdf. As for the matrix approach, I'm not exactly sure what you mean by the multiplication-- Is it meant to save time from the multiply operations in the edge function? – aurreco Jun 13 '22 at 18:57
  • matrix form is just numerically more stable way of computing the same thing (without any if statements) once you obtain the inverse its the same as you current approach (if you expand the equation you will get the same stuff), by the multiplication I meant you can apply the same optimization like you do now so instead of multiply of consequent iteration you just add constant in a loop ... – Spektre Jun 13 '22 at 19:20

1 Answers1

0

Rather than multiply every y value in the formula by -1, replace them with (top - y).

greg-tumolo
  • 698
  • 1
  • 7
  • 30
  • I see how this is supposed to work, but it doesn't. I think because for example (c.y - a.y) term in the original determinant function becomes (TOP - c.y - (TOP - a.y)) = (TOP - c.y - TOP + a.y) = (a.y - c.y) which is what my new determinant function looks like already. – aurreco Jun 13 '22 at 19:45
  • I answered hastily…I was referring to the formula to transform from vector space to screen space. – greg-tumolo Jun 13 '22 at 20:06
  • Please, post `draw_point`. – greg-tumolo Jul 28 '22 at 15:19
  • 1
    I ended up solving the problem a while back. The code for init_edge was what was wrong-- I calculated the determinant incorrectly. – aurreco Jul 30 '22 at 17:38