Rasterization

The Rasterization Stage

Rasterization: What Are We Trying to Solve?

Rasterization is the process by which a primitive is converted to a two dimensional image. Each point of this image contains such information as color and depth Thus rasterizing a primitive consists of two parts The first is to determine which squares of an integer grid in window coordinates are occupied by the primitive The second is assigning a color and a depth value to each such square.

Quote from Chapter 3 of the OpenGL Specification (Version 1.0) document, written by Mark Segal, Kurt Akeley, and Chris Frazier in 1994.

In the previous chapter, we learned how to perform the first step of the rasterization algorithm, which is to project the triangle from 3D space onto the canvas. This description is not entirely accurate, in fact, since what we did was transform the triangle from camera space to screen space, which, as mentioned in the previous chapter, is also a three-dimensional space. However, the x- and y-coordinates of the vertices in screen-space correspond to the positions of the triangle vertices on the canvas, and by converting them from screen-space to NDC (Normalized Device Coordinates) space and then finally from NDC space to raster space, what we obtain in the end are the vertices' 2D coordinates in raster space. Finally, we also know that the z-coordinates of the vertices in screen-space retain the original z-coordinate of the vertices in camera space, inverted so that we deal with positive numbers rather than negative ones.

What we need to do next is to loop over the pixels in the image and find out if any of these pixels overlap the "projected image of the triangle" (Figure 1). In graphics APIs specifications, this test is sometimes called the inside-outside test or the coverage test. If they do, we then set the pixel in the image to the triangle's color. The idea is simple, but of course, we now need to come up with a method to find if a given pixel overlaps a triangle. This is essentially what we will study in this chapter. We will learn about the method that is typically used in rasterization to solve this problem, which uses a technique known as the edge function. We are now going to describe and study this edge function, which will also provide valuable information about the position of the pixel within the projected image of the triangle, known as barycentric coordinates. Barycentric coordinates play an essential role in computing the actual depth (or the z-coordinate) of the point on the surface of the triangle that the pixel overlaps. We will also explain what barycentric coordinates are in this chapter and how they are computed.

At the end of this chapter, you will be able to produce a very basic rasterizer. In the next chapter, we will look into the possible issues with this very naive implementation of the rasterization algorithm. We will list what these issues are as well as study how they are typically addressed.

A lot of research has been done to optimize the algorithm. The goal of this lesson is not to teach you how to write or develop an optimized and efficient renderer based on the rasterization algorithm. The goal of this lesson is to teach the basic principles of the rendering technique. It's important to note that the techniques we present in these chapters are indeed used, but how they are implemented, either on the GPU or in a CPU version of a production renderer, is likely to be a highly optimized version of the same idea. What is truly important is to understand the principle and how it works in general. From there, you can study on your own the different techniques used to speed up the algorithm. But the techniques presented in this lesson are generic and make up the foundations of any rasterizer.

Keep in mind that drawing a triangle (since the triangle is a primitive we will use in this case) is a two-step problem:

• We first need to find which pixels overlap the triangle.

• We then need to define the colors that the pixels overlapping the triangle should be set to, a process called shading.

The rasterization stage deals essentially with the first step. The reason we say "essentially" rather than "exclusively" is that at the rasterization stage, we will also compute something called barycentric coordinates, which to some extent, are used in the second step.

The Edge Function

As mentioned previously, there are several possible methods to determine if a pixel overlaps a triangle. While it would be beneficial to document older techniques, this lesson will only present the method that is generally used today. This method was introduced by Juan Pineda in a 1988 paper titled "A Parallel Algorithm for Polygon Rasterization."

Before delving into Pineda's technique itself, let's first understand the principle behind his method. Imagine the edge of a triangle as a line that splits the 2D plane (the image plane) into two (as shown in Figure 2). The principle of Pineda's method is to find a function, which he called the edge function, so that when we test on which side of this line a given point is (the point P in Figure 2), the function returns a positive number when it is to the right of the line, a negative number when it is to the left of this line, and zero when the point is exactly on the line.

In Figure 2, we applied this method to the first edge of the triangle (defined by the vertices v0 and v1; note that the order is important). If we now apply the same method to the other two edges (v1-v2 and v2-v0), we can see that there is an area (the white triangle) within which all points yield positive results (Figure 3). If we take a point within this area, then we will find that this point is to the right of all three edges of the triangle. If P is a point in the center of a pixel, we can then use this method to determine if the pixel overlaps the triangle. If for this point, the edge function returns a positive number for all three edges, then the pixel is contained within the triangle (or may lie on one of its edges). The function Pineda uses also happens to be linear, which means that it can be computed incrementally, but we will revisit this point later.

Now that we understand the principle, let's examine the function itself. The edge function is defined as follows (for the edge defined by vertices V0 and V1):

As the paper mentions, this function has the useful property that its value is related to the position of the point \(P\) relative to the edge defined by the points \(V0\) and \(V1\):

• \(E(P) > 0\) if \(P\) is to the "right" side

• \(E(P) = 0\) if \(P\) is exactly on the line

• \(E(P) < 0\) if \(P\) is to the "left" side

This function is mathematically equivalent to the magnitude of the cross products between the vector \((V1-V0)\) and \((P-V0)\). We can also represent these vectors in matrix form. Presenting them as a matrix serves no other purpose than to neatly display the two vectors:

If we denote \(A = (P-V0)\) and \(B = (V1 - V0)\), then we can represent the vectors \(A\) and \(B\) as a 2x2 matrix:

The determinant of this matrix can be calculated as:

If we substitute the vectors \(A\) and \(B\) back with \((P-V0)\) and \((V1-V0)\), respectively, we obtain:

As you can see, this is similar to the edge function we have defined above. In other words, the edge function can be interpreted either as the determinant of the 2x2 matrix defined by the components of the 2D vectors \((P-V0)\) and \((V1-V0)\), or as the magnitude of the cross product of the vectors \((P-V0)\) and \((V1-V0)\). Both the determinant and the magnitude of the cross product of two vectors have the same geometric interpretation. Let's explain.

Understanding what happens is easier when we examine the result of a cross-product between two 3D vectors (Figure 4). In 3D, the cross-product returns another 3D vector that is perpendicular (or orthonormal) to the two original vectors. However, as shown in Figure 4, the magnitude of that orthonormal vector also changes with the orientation of the two vectors with respect to each other. Assuming a right-hand coordinate system, when vectors A (red) and B (blue) are either pointing exactly in the same direction or in opposite directions, the magnitude of the third vector C (green) is zero. Vector A has coordinates (1,0,0) and is fixed. When vector B has coordinates (0,0,-1), then the green vector, vector C, has coordinates (0,-1,0). If we were to find its "signed" magnitude, we would discover that it is equal to -1. Conversely, when vector B has coordinates (0,0,1), vector C has coordinates (0,1,0), and its signed magnitude is equal to 1. In one case, the "signed" magnitude is negative, and in the second case, it is positive. In fact, in 3D, the magnitude of a vector can be interpreted as the area of the parallelogram having A and B as sides, as shown in Figure 5 (refer to the Wikipedia article on the cross product for more details on this interpretation):

An area should always be positive, though the sign of the above equation provides an indication of the orientation of vectors A and B with respect to each other. When B is within the half-plane defined by vector A and a vector orthogonal to A (let's call this vector D; note that A and D form a 2D Cartesian coordinate system), then the result of the equation is positive. When B is within the opposite half-plane, the result is negative (Figure 6). Another way to explain this result is that the result is positive when the angle \(\theta\) is in the range \([0,\pi]\) and negative when \(\theta\) is in the range \([\pi, 2\pi]\). Note that when \(\theta\) is exactly equal to 0 or \(\pi\), then the cross-product or the edge function returns 0.

To determine if a point is inside a triangle, what really matters is the sign of the function used to compute the area of the parallelogram. However, the area itself also plays a crucial role in the rasterization algorithm; it is used to compute the barycentric coordinates of the point within the triangle, a technique we will explore next. The cross-product in 3D and 2D shares the same geometric interpretation, thus the cross-product between two 2D vectors also returns the "signed" area of the parallelogram defined by the two vectors. The only difference is that in 3D, to compute the area of the parallelogram, you need to use this equation:

while in 2D, this area is given by the cross-product itself (which, as mentioned before, can also be interpreted as the determinant of a 2x2 matrix):

From a practical standpoint, all we need to do now is test the sign of the edge function computed for each edge of the triangle and another vector defined by a point and the first vertex of the edge (Figure 7).

If all three tests are positive or equal to 0, then the point is inside the triangle (or lying on one of its edges). If any one of the tests is negative, then the point is outside the triangle. In code, we have:

cpp
bool edgeFunction(const Vec2f &a, const Vec3f &b, const Vec2f &c) {
    return ((c.x - a.x) * (b.y - a.y) - (c.y - a.y) * (b.x - a.x) >= 0);
}

bool inside = true;
inside &= edgeFunction(V0, V1, p);
inside &= edgeFunction(V1, V2, p);
inside &= edgeFunction(V2, V0, p);

if (inside) {
    // point p is inside the triangle defined by vertices V0, V1, V2
    ...
}

Alternative to the Edge Function

There are alternative methods to the edge function for determining if pixels overlap triangles. However, as mentioned in the introduction of this chapter, we won't explore them in this lesson. For reference, another common technique is scanline rasterization, based on the Bresenham algorithm, generally used for drawing lines. GPUs mostly use the edge method because it is more generic and easier to parallelize than the scanline approach, but further details on this topic will not be provided in this lesson.

Be Careful! Winding Order Matters

One aspect we haven't discussed yet, but which holds great importance in CG, is the order in which the vertices that make up the triangles are declared. There are two possible conventions, illustrated in Figure 8: clockwise or counter-clockwise ordering, also known as winding. Winding is crucial because it essentially defines one important property of the triangle: the orientation of its normal. Remember, the normal of the triangle can be computed from the cross product of two vectors, \(A=(V2-V0)\) and \(B=(V1-V0)\). Suppose \(V0=\{0,0,0\}\), \(V1=\{1,0,0\}\), and \(V2=\{0,-1,0\}\), then \((V1-V0)=\{1,0,0\}\) and \((V2-V0)=\{0,-1,0\}\). Now, let's compute the cross product of these two vectors:

However, if you declare the vertices in counter-clockwise order, with \(V0=\{0,0,0\}\), \(V1=\{0,-1,0\}\), and \(V2=\{1,0,0\}\), \((V1-V0)=\{0,-1,0\}\) and \((V2-V0)=\{1,0,0\}\). Let's compute the cross-product of these two vectors again:

As expected, the two normals point in opposite directions. The orientation of the normal is crucial for several reasons, one of the most important being face culling. Most rasterizers, and even ray tracers, may not render triangles whose normals face away from the camera, a process known as back-face culling. Most rendering APIs, such as OpenGL or DirectX, allow back-face culling to be turned off. However, it's still important to be aware that vertex ordering affects rendering outcomes, among many other factors. Surprisingly, the edge function is also affected by vertex ordering. It's worth noting that there's no universal rule for choosing the order. Many details in a renderer's implementation can change the normal's orientation, so you can't assume that declaring vertices in a certain order guarantees the normal will face a specific direction. For instance, using vectors \(A=(V0-V1)\) and \(B=(V2-V1)\) instead could produce the same normal but flipped. Even if using \(A=(V1-V0)\) and \(B=(V2-V0)\), the order of vectors in the cross-product) (\(A \times B = -B \times A\)) also affects the normal's direction. For these reasons, assuming a specific outcome based on vertex order is unreliable. What's important is adhering to the chosen convention. Generally, graphics APIs like OpenGL and DirectX expect triangles to be declared in counter-clockwise order, which we will also adopt. Now, let's explore how ordering impacts the edge function.

Why does winding matter when it comes to the edge function? You may have noticed that, since the beginning of this chapter, we have drawn the triangle vertices in clockwise order. We have also defined the edge function as:

If we adhere to this convention, then points to the right of the line defined by vertices A and B will be positive. For example, a point to the right of V0V1, V1V2, or V2V0 would be positive. However, if we were to declare the vertices in counter-clockwise order, points to the right of an edge defined by vertices A and B would still be positive, but then they would be outside the triangle. In other words, points overlapping the triangle would not be positive but negative (Figure 10). You can potentially still make the code work with positive numbers with a small change to the edge function:

In conclusion, depending on the winding convention you use, you may need to employ one version of the edge function or the other.

Barycentric Coordinates

Computing barycentric coordinates is not necessary to make the rasterization algorithm work. For a naive implementation of the rendering technique, all you need is to project the vertices and use a technique like the edge function that we described above to find if pixels are inside triangles. These are the only two necessary steps to produce an image. However, the result of the edge function, which, as we explained above, can be interpreted as the area of the parallelogram defined by vectors A and B, can directly be used to compute these barycentric coordinates. Thus, it makes sense to study the edge function and barycentric coordinates simultaneously.

Before we go any further, let's explain what barycentric coordinates are. First, they come in a set of three floating-point numbers which, in this lesson, we will denote \(\lambda_0\), \(\lambda_1\), and \(\lambda_2\). Many different conventions exist, but Wikipedia uses the Greek letter lambda (\(\lambda\)), which is also used by other authors (the Greek letter omega (\(\omega\)) is also sometimes used). This doesn't matter; you can call them whatever you want. In short, the coordinates can be used to define any point on the triangle in the following manner:

Where, as usual, V0, V1, and V2 are the vertices of a triangle. These coordinates can take on any value, but for points that are inside the triangle (or lying on one of its edges), they can only be in the range [0,1], and the sum of the three coordinates is equal to 1. In other words:

This is a form of interpolation, if you will. They are also sometimes defined as weights for the triangle's vertices (which is why in the code we will denote them with the letter w). A point overlapping the triangle can be defined as "a little bit of V0 plus a little bit of V1 plus a little bit of V2". Note that when any of the coordinates is 1 (which means that the others in this case are necessarily 0), then the point P is equal to one of the triangle's vertices. For instance, if \(\lambda_2 = 1\), then P is equal to V2. Interpolating the triangle's vertices to find the position of a point inside the triangle is not that useful. But the method can also be used to interpolate across the surface of the triangle any quantity or variable that has been defined at the triangle's vertices. Imagine, for instance, that you have defined a color at each vertex of the triangle. Say V0 is red, V1 is green, and V2 is blue (Figure 12). What you want to do is find how these three colors are interpolated across the surface of the triangle. If you know the barycentric coordinates of a point P on the triangle, then its color \(C_P\) (which is a combination of the triangle vertices' colors) is defined as:

This technique proves very useful for shading triangles. Data associated with the vertices of triangles, called vertex attributes, is a common and important technique in computer graphics (CG). The most common vertex attributes include colors, normals, and texture coordinates. Practically, this means that when you define a triangle, you don't only pass the triangle's vertices to the renderer but also its associated vertex attributes. For example, to shade the triangle, you may need color and normal vertex attributes, which means each triangle will be defined by 3 points (the triangle vertex positions), 3 colors (the colors of the triangle vertices), and 3 normals (the normals of the triangle vertices). Normals too can be interpolated across the surface of the triangle. Interpolated normals are used in a technique called smooth shading, first introduced by Henri Gouraud. We will explore this technique later when we discuss shading.

How do we find these barycentric coordinates? It turns out to be straightforward. As mentioned previously when presenting the edge function, the result of the edge function can be interpreted as the area of the parallelogram defined by vectors A and B. If you examine Figure 8, it's easy to see that the area of the triangle defined by vertices V0, V1, and V2 is just half the area of the parallelogram defined by vectors A and B. The area of the triangle is thus half the area of the parallelogram, which we know can be computed by the cross-product of the two 2D vectors A and B:

If point P is inside the triangle, then, as illustrated in Figure 3, we can draw three sub-triangles: V0-V1-P (green), V1-V2-P (magenta), and V2-V0-P (cyan). It is evident that the sum of these three sub-triangle areas is equal to the area of the triangle V0-V1-V2:

Let's first intuitively understand how they work, which might be easier if you look at Figure 14. Each image in the series shows what happens to the sub-triangle as point P, originally on the edge defined by vertices V1-V2, moves towards V0. Initially, P lies exactly on the edge V1-V2. This scenario is akin to basic linear interpolation between two points. In other words, we could write:

With \(\lambda_1 + \lambda_2 = 1\), thus \(\lambda_2 = 1 - \lambda_1\). More interestingly, in this specific case, if the generic equation for computing the position of P using barycentric coordinates is:

it clearly shows that, in this particular case, \(\lambda_0\) equals 0.

This is relatively straightforward. Also, note that in the first image, the red triangle is not visible, indicating that P is closer to V1 than to V2. Therefore, \(\lambda_1\) is necessarily greater than \(\lambda_2\). Furthermore, the green triangle is larger than the blue triangle in the first image. To summarize: when the red triangle is not visible, \(\lambda_0\) equals 0, \(\lambda_1\) is greater than \(\lambda_2\), and the green triangle is larger than the blue triangle. This suggests a relationship between the area of the triangles and the barycentric coordinates. Moreover, the red triangle appears to be associated with \(\lambda_0\), the green triangle with \(\lambda_1\), and the blue triangle with \(\lambda_2\):

• \(\lambda_0\) is proportional to the area of the red triangle,

• \(\lambda_1\) is proportional to the area of the green triangle,

• \(\lambda_2\) is proportional to the area of the blue triangle.

Jumping directly to the last image, where P equals V0, this situation can only occur if \(\lambda_0\) equals 1 and the other two coordinates equal 0:

In this specific scenario, the blue and green triangles have disappeared, and the area of the triangle V0-V1-V2 is the same as the area of the red triangle. This confirms our intuition about a relationship between the areas of the sub-triangles and the barycentric coordinates. From the observations above, we can also deduce that each barycentric coordinate is related to the area of the sub-triangle defined by the edge directly opposite the associated vertex and the point P. In other words (Figure 15):

• \(\color{red}{\lambda_0}\) is associated with V0. The edge opposite V0 is V1-V2, defining the red triangle with V1-V2-P.

• \(\color{green}{\lambda_1}\) is associated with V1. The edge opposite V1 is V2-V0, defining the green triangle with V2-V0-P.

• \(\color{blue}{\lambda_2}\) is associated with V2. The edge opposite V2 is V0-V1, defining the blue triangle with V0-V1-P.

The areas of the red, green, and blue triangles are given by the respective edge functions we have used before to determine if P is inside the triangle, divided by 2 (recall that the edge function provides the "signed" area of the parallelogram defined by vectors A and B, where A and B can be any of the triangle's three edges):

The barycentric coordinates can be computed as the ratio between the area of the sub-triangles and the area of the triangle V0V1V2:

Dividing by the triangle's area essentially normalizes the coordinates. For example, when P is positioned at V0, the area of the triangle V2V1P (the red triangle) is the same as the area of the triangle V0V1V2. Thus, dividing one by the other yields 1, which is the value of the coordinate \(\lambda_0\). Since the green and blue triangles' areas are 0 in this case, \(\lambda_1\) and \(\lambda_2\) equal 0, and we get:

which is the expected result.

To compute the area of a triangle, we can use the edge function as mentioned earlier. This method works for the sub-triangles as well as the main triangle V0V1V2. However, the edge function returns the area of the parallelogram instead of the triangle's area (Figure 8). But since the barycentric coordinates are computed as the ratio between the sub-triangle area and the main triangle area, we can ignore the division by 2 (this division, present in both the numerator and denominator, cancels out):

Note that \(E_{01}(V2) = E_{12}(V0) = E_{20}(V1) = 2 \cdot Area_{tri}(V0,V1,V2)\).

Let's examine the code implementation. Previously, we were computing the edge functions to test if points were inside triangles, simply returning true or false depending on whether the result of the function was positive or negative. To compute the barycentric coordinates, we need the actual result of the edge function. We can also use the edge function to compute the area (multiplied by 2) of the triangle. Here is an implementation that tests if a point P is inside a triangle and, if so, computes its barycentric coordinates:

float edgeFunction(const Vec2f &a, const Vec3f &b, const Vec2f &c)
{
    return (c.x - a.x) * (b.y - a.y) - (c.y - a.y) * (b.x - a.x);
}

float area = edgeFunction(v0, v1, v2); // Area of the triangle multiplied by 2
float w0 = edgeFunction(v1, v2, p); // Signed area of the triangle v1v2p multiplied by 2
float w1 = edgeFunction(v2, v0, p); // Signed area of the triangle v2v0p multiplied by 2
float w2 = edgeFunction(v0, v1, p); // Signed area of the triangle v0v1p multiplied by 2

// If point p is inside triangles defined by vertices v0, v1, v2
if (w0 >= 0 && w1 >= 0 && w2 >= 0) {
    // Barycentric coordinates are the areas of the sub-triangles divided by the area of the main triangle
    w0 /= area;
    w1 /= area;
    w2 /= area;
}

Let's apply this code to generate an actual image.

Interpolate vs. Extrapolate

It's important to note that the computation of barycentric coordinates works regardless of the point's position relative to the triangle. In other words, the coordinates are valid whether the point is inside or outside the triangle. When the point is inside, using barycentric coordinates to evaluate the value of a vertex attribute is referred to as interpolation, and when the point is outside, it's called extrapolation. This distinction is crucial because, in some instances, we may need to evaluate the value of a given vertex attribute for points that do not overlap with triangles. Specifically, this is necessary to compute the derivatives of the triangle texture coordinates, which are used for proper texture filtering. If you're interested in delving deeper into this topic, we recommend reading the lesson on Texture Mapping. For now, remember that barycentric coordinates are valid even when the point doesn't overlap with the triangle, and it's important to understand the difference between vertex attribute extrapolation and interpolation.

Rasterization Rules

In certain cases, a pixel may overlap more than one triangle, especially when it lies precisely on an edge shared by two triangles, as illustrated in Figure 17. Such a pixel would pass the coverage test for both triangles. If the triangles are semi-transparent, a dark edge may appear where the pixels overlap the two triangles due to the way semi-transparent objects are combined (imagine two superimposed semi-transparent sheets of plastic. The surface becomes more opaque and appears darker than the individual sheets). This phenomenon results in something similar to what is depicted in Figure 18, where a darker line appears where the two triangles share an edge.

To address this issue, a rule is needed to ensure that a pixel can never overlap two triangles sharing an edge twice. Most graphics APIs, such as OpenGL and DirectX, define what they call the top-left rule. This rule stipulates that a pixel or point is considered to overlap a triangle if it is either inside the triangle or lies on a triangle's top edge or any edge considered to be a left edge. What constitutes a top and left edge? Figure 19 clearly demonstrates the definition of top and left edges.

• A top edge is an edge that is perfectly horizontal and whose defining vertices are positioned above the third one. Technically, this means that the y-coordinates of the vector V[(X+1)%3]-V[X] are equal to 0, and its x-coordinates are positive (greater than 0).

• A left edge is essentially an edge that is ascending. Given that vertices are defined in clockwise order in our context, an edge is considered ascending if its respective vector V[(X+1)%3]-V[X] (where X can be 0, 1, 2) has a positive y-coordinate.

In pseudo-code, we have:

// Does it pass the top-left rule?
Vec2f v0 = { ... };
Vec2f v1 = { ... };
Vec2f v2 = { ... };

float w0 = edgeFunction(v1, v2, p); 
float w1 = edgeFunction(v2, v0, p); 
float w2 = edgeFunction(v0, v1, p); 

Vec2f edge0 = v2 - v1;
Vec2f edge1 = v0 - v2;
Vec2f edge2 = v1 - v0;

bool overlaps = true;

// If the point is on the edge, test if it is a top or left edge, 
// otherwise test if the edge function is positive
overlaps &= (w0 == 0 ? ((edge0.y == 0 && edge0.x > 0) || edge0.y > 0) : (w0 > 0));
overlaps &= (w1 == 0 ? ((edge1.y == 0 && edge1.x > 0) || edge1.y > 0) : (w1 > 0));
overlaps &= (w2 == 0 ? ((edge2.y == 0 && edge2.x > 0) || edge2.y > 0) : (w2 > 0));

if (overlaps) {
    // The pixel overlaps the triangle.
    ...
}

This version serves as a proof of concept but is highly unoptimized. The fundamental idea is to first check whether any of the values returned by the edge function equals 0, indicating that the point lies on the edge. In such cases, we test if the edge in question is a top-left edge. If it is, it returns true. If the value returned by the edge function is not 0, we then return true if the value is greater than 0. The top-left rule won't be implemented in the program provided with this lesson.

Putting Things Together: Finding if a Pixel Overlaps a Triangle

Let's apply the techniques we've learned in this chapter to a program that produces an actual image. We'll assume that the triangle has already been projected (refer to the last chapter of this lesson for a complete implementation of the rasterization algorithm). We will also assign a color to each vertex of the triangle. The process is as follows: we will loop over all the pixels in the image and test if they overlap the triangle using the edge function method. If the edge function returns a positive number for all the edges, then the pixel overlaps the triangle. We can then compute the pixel's barycentric coordinates and use these coordinates to shade the pixel by interpolating the color defined at each vertex of the triangle. The result of the framebuffer is saved to a PPM file (which you can open with Photoshop). The output of the program is shown in Figure 20.

An optimization for this program could be to loop over the pixels contained within the bounding box of the triangle. This optimization is not implemented in this version of the program, but you are encouraged to try it yourself (using the code from the previous chapters). The source code of this lesson is available in the last chapter.

Note that in this version of the program, we move point P to the center of each pixel. Alternatively, you could use the integer coordinates of the pixel. More details on this topic will be discussed in the next chapter.

// c++ -o raster2d raster2d.cpp
// (c) www.scratchapixel.com

#include <cstdio>
#include <cstdlib>
#include <fstream>

typedef float Vec2[2];
typedef float Vec3[3];
typedef unsigned char Rgb[3];

inline float edgeFunction(const Vec2 &a, const Vec2 &b, const Vec2 &c) {
    return (c[0] - a[0]) * (b[1] - a[1]) - (c[1] - a[1]) * (b[0] - a[0]);
}

int main(int argc, char **argv) {
    Vec2 v0 = {491.407, 411.407};
    Vec2 v1 = {148.593, 68.5928};
    Vec2 v2 = {148.593, 411.407};
    Vec3 c0 = {1, 0, 0};
    Vec3 c1 = {0, 1, 0};
    Vec3 c2 = {0, 0, 1};
    
    const uint32_t w = 512;
    const uint32_t h = 512;
    
    Rgb *framebuffer = new Rgb[w * h];
    memset(framebuffer, 0x0, w * h * 3);
    
    float area = edgeFunction(v0, v1, v2);
    
    for (uint32_t j = 0; j < h; ++j) {
        for (uint32_t i = 0; i < w; ++i) {
            Vec2 p = {i + 0.5f, j + 0.5f};
            float w0 = edgeFunction(v1, v2, p);
            float w1 = edgeFunction(v2, v0, p);
            float w2 = edgeFunction(v0, v1, p);
            if (w0 >= 0 && w1 >= 0 && w2 >= 0) {
                w0 /= area;
                w1 /= area;
                w2 /= area;
                float r = w0 * c0[0] + w1 * c1[0] + w2 * c2[0];
                float g = w0 * c0[1] + w1 * c1[1] + w2 * c2[1];
                float b = w0 * c0[2] + w1 * c1[2] + w2 * c2[2];
                framebuffer[j * w + i][0] = (unsigned char)(r * 255);
                framebuffer[j * w + i][1] = (unsigned char)(g * 255);
                framebuffer[j * w + i][2] = (unsigned char)(b * 255);
            }
        }
    }
    
    std::ofstream ofs;
    ofs.open("./raster2d.ppm");
    ofs << "P6\n" << w << " " << h << "\n255\n";
    ofs.write((char*)framebuffer, w * h * 3);
    ofs.close();
    
    delete [] framebuffer;
    
    return 0; 
}

As you can see, and in conclusion, we can state that the rasterization algorithm is quite straightforward, and its basic implementation is relatively easy as well.

Conclusion and What's Next?

Many fascinating techniques and pieces of trivia are related to the topic of barycentric coordinates, but this lesson serves merely as an introduction to the rasterization algorithm, so we won't delve deeper. However, one interesting piece of trivia to note is that barycentric coordinates remain constant along lines parallel to an edge, as illustrated in Figure 21.

In this lesson, we've explored two important methods and various concepts:

• Firstly, we learned about the edge function and its application in determining whether a point \(P\) overlaps a triangle. The edge function is calculated for each edge of the triangle, with a second vector defined by the edge's first vertex and the point \(P\). If the function yields a positive result for all three edges, then point \(P\) overlaps the triangle.

• Additionally, we discovered that the outcome of the edge function can be used to compute the barycentric coordinates of point \(P\). These coordinates allow for the interpolation of vertex data or attributes across the triangle's surface, acting as weights for the triangle's vertices. The most commonly interpolated vertex attributes include color, normals, and texture coordinates.