Rasterization: a Practical Implementation

News (August, 31): We are working on Scratchapixel 3.0 at the moment (current version of 2). The idea is to make the project open source by storing the content of the website on GitHub as Markdown files. In practice, that means you and the rest of the community will be able to edit the content of the pages if you want to contribute (typos and bug fixes, rewording sentences). You will also be able to contribute by translating pages to different languages if you want to. Then when we publish the site we will translate the Markdown files to HTML. That means new design as well.

That's what we are busy with right now and why there won't be a lot of updates in the weeks to come. More news about SaP 3.0 soon.

We are looking for native Engxish (yes we know there's a typo here) speakers that will be willing to readproof a few lessons. If you are interested please get in touch on Discord, in the #scratchapixel3-0 channel. Also looking for at least one experienced full dev stack dev that would be willing to give us a hand with the next design.

Feel free to send us your requests, suggestions, etc. (on Discord) to help us improve the website.

And you can also donate). Donations go directly back into the development of the project. The more donation we get the more content you will get and the quicker we will be able to deliver it to you.

15 mns read.

Interpolating Vertex Attributes

All we need to get a basis rasterizer working is to know how to project triangles onto the screen, convert the projected coordinates to raster space, then rasterize triangles, and potentially use a depth-buffer to solve the visibility problem. That would already be enough to create images of 3D scenes, which would be both perspective correct, and in which the visibility problem would be solved (objects that are hidden by others, are indeed not appearing in front of objects that are supposed to occlude them). This is already a good result. The code we have presented to do this is functional, but could be optimized greatly; however, optimizing the rasterization algorithm is not something we will look into in this lesson.

Perspective Correct Vertex Attribute Interpolation

So, what are we going to talk about in this chapter? In the chapter devoted to rasterization, we talked about using barycentric coordinates to interpolate vertex data or vertex attribute (which is the most common name). We can use barycentric coordinates to interpolate depth values of course, and this is one of their main function, but they can also be used to interpolate vertex attributes; vertex attributes play a very important role in rendering, especially when it comes to lighting and shading. We will provide more information on how vertex attributes are used in shading when we get to shading in this section. But you don't need to know anything about shading to actually understand the concept of perspective correct interpolation and vertex attributes.

You might wonder why we don't talk about this topic in shading if it then relates to shading more specifically. In fact vertex attributes do relate to shading more specifically but the the topic of perspective correct interpolation relates more specifically to the topic of rasterization.

As you already know, we need to store the z-coordinates of the original vertices (camera space vertices) in the z-coordinate of our projected vertices (screen space vertices). This is required so that we can compute the depth of points lying across the surface of projected triangle. Depth, as explained in the last chapter, is required to solve the visibility problem and is computed by linearly interpolating the reciprocal of the triangle vertices z-coordinates using barycentric coordinates. Though the exact same technique can be used to interpolate any other variable we want as long as the value of this variable is defined at the triangle's vertices, similarly to the way we've stored in the projected point z-coordinates, the vertices original z-coordinates. For example it is very common to store at the triangle vertices a color. The other two most common variables or attributes (which is the proper term in CG) stored at the triangle vertices are texture coordinates and normals. Texture coordinates are 2D coordinates used for texturing (a technique we will study in this section). Normals are used in shading and define the orientation of the surface (check the lessons on shading to learn more about normals and smooth shading in particular). In this lesson we will more specifically use color and texture coordinates to illustrate the problem of perspective correct interpolation.

Figure 1: Finding the value of Cp requires to interpolate the colors defined at the triangle vertices using P's barycentric coordinates.

As mentioned in the chapter on the rasterization stage, we can specify colors or anything else we want for the triangle vertices. These attributes can be interpolated using barycentric coordinates to find what the value of these attribute should be for any point inside the triangle. In other words, vertex attributes must be interpolated across the surface of a triangle when it is rasterized. The process is as follows:

Figure 2: in 3D space the point P is in the middle of the quad, though in the perspective view, this point clearly doesn't seem to be in the middle of the geometry. In 3D space, point P on the green triangle's edge is exactly in the middle of V1 and V2. But in the bottom image, we can see that P is closer to V1 than it is to V2.

This technique though just doesn't work. To understand why let's see what happens to a point located in the middle of a 3D quad. As you can see in the top view of figure 2, we clearly have a quad and point P, is clearly in the middle of that quad (P is located at the intersection of the quad's diagonals). Though, when we look at this good from a random viewpoint it is easy to see that depending on the quad's orientation with respect to the camera, P doesn't appear to be in the centre of the quad anymore. This is due to perspective projection which as mentioned before, preserves lines but doesn't preserve distances. Though remember that barycentric coordinates are computed in screen space. Imagine that the quad is made of two triangles. In 3D space, P is located at equal distance between V1-V2 thus somehow it's barycentric coordinates in 3D space are (0,0.5,0.5). Though, in screen space, since P is closer to V1 than it is to V2, then \(\lambda_1\) is greater than \(\lambda_2\) (and \lambda_0 is equal to 0). The problem though is that these are the coordinates that are used to interpolate the triangle's vertex attributes. If V1 is white and V2 is black then the color at P should be 0.5. But if \(\lambda_1\) is greater than \(\lambda_2\) then we will get a value greater than 0.5. Thus clearly the technique we are using to interpolate vertex attributes doesn't work. Let's assume like in figure 1 that \(\lambda_1\) and \(\lambda_2\) are equal to 0.666 and 0.334 respectively. If we interpolate the triangle's vertex colors, we get:

$$C_P = \lambda_0 * C_0 + \lambda_1 * C_1 + \lambda_2 * C_2 = 0 * C_0 + 0.666 * 1 + 0.334 * 0 = 0.666.$$

We get 0.666 for the color of P and not 0.5 as we should. There is a problem and this problem relates to some extent to what we learned in the previous chapter regarding the interpolation of the vertices z-coordinates.

Figure 3: the ratio \((Z-Z_0)/(Z_1-Z_0)\) matches the ratio \((C-C_0)/(C_1-C_0)\). We already know from the previous chapter how to compute \(Z\). We can use these two equations to solve for \(C\).

Hopefully finding the right solution is not hard. Let's imagine that we have a triangle with two z-coordinates \(Z_0\) and \(Z_1\) on each side of the triangle as shown in figure 3. If we connect these two points we can interpolate the z-coordinate of a point on this line using linear interpolation. We can do the same thing with two values of a vertex attributes \(C_0\) and \(C_1\) defined at the same positions on the triangle than \(Z_0\) and \(Z_1\) respectively. Technically, because both \(Z\) and \(C\) are computed using linear interpolation, we can write the following equality (equation 1):

$$\dfrac{Z - Z_0}{Z_1 - Z_0} = \dfrac {C-C_0}{C_1 - C_0}.$$

We also know from the last chapter that (equation 2):

$$Z = \dfrac{1}{\dfrac{1}{Z_0}(1-q) + \dfrac{1}{Z_1}q}.$$

The first thing we are going to do is substitute the equation for Z (equation 2) in the left-hand side of equation 1. The trick to simplify the resulting equation is to multiply the numerator and denominator of equation 2 by \(Z_0Z_1\) to get rid of the \(1/Z_0\) and \(1/Z_1\) terms (this is of course not explained anywhere but here on Scratchapixel):

$$ \begin{align} \dfrac{\dfrac{1}{\dfrac{1}{Z_0}(1-q)+\dfrac{1}{Z_1}q} - Z_0}{(Z_1 - Z_0)} & = \dfrac{\dfrac{Z_0Z_1}{Z_1(1-q)+Z_0q} - Z_0}{(Z_1 - Z_0)}\\ & = \dfrac{\dfrac{Z_0Z_1 - Z_0(Z_1(1-q) + Z_0q)}{Z_1(1-q) + Z_0q}}{Z_1 - Z_0}\\ & = \dfrac{\dfrac{Z_0Z_1q - Z_0^2q}{Z_1(1-q) + Z_0q}}{Z_1 - Z_0}\\ & = \dfrac{\dfrac{Z_0q(Z_1 - Z_0)}{Z_1(1-q) + Z_0q}}{Z_1 - Z_0}\\ & = \dfrac{Z_0q}{Z_1(1-q) + Z_0q}\\ & = \dfrac{Z_0q}{q(Z_0 - Z_1) + Z_1}.\\ \end{align}$$

We can now solve for \(C\):

$$ \begin{array}{l} \dfrac{Z_0q}{q(Z_0 - Z_1) + Z_1} = \dfrac {C-C_0}{C_1 - C_0},\\ C = \dfrac{Z_0q}{q(Z_0 - Z_1) + Z_1}(C_1 - C_0) + C_0,\\ C = \dfrac{Z_0q(C_1 - C_0) + C_0(q(Z_0 - Z_1) + Z_1)}{q(Z_0 - Z_1) + Z_1},\\ C = \dfrac{Z_0qC_1 - Z_0qC_0 + C_0qZ_0 - C_0qZ_1 + C_0Z_1}{q(Z_0 - Z_1) + Z_1},\\ C = \dfrac{Z_0qC_1 - C_0qZ_1 + C_0Z_1}{q(Z_0 - Z_1) + Z_1},\\ C = \dfrac{C_0Z_1(1-q) + C_1Z_0q}{Z_1(1 - q) + Z_0q}.\\ \end{array}$$

If we now multiply the numerator and the denominator by \(1/Z_0Z_1\), we can then extract a factor of \(Z\) from the right-hand side of the equation:

$$ \begin{array}{l} C = \dfrac{\dfrac{1}{Z_0Z_1}(C_0Z_1(1-q) + C_1Z_0q)}{\dfrac{1}{Z_0Z_1}(Z_1(1 - q) + Z_0q)},\\ C = \dfrac{\dfrac{C_0}{Z_0}(1-q) + \dfrac{C_1}{Z_1}q}{\dfrac{1}{Z_0}(1 - q) + \dfrac{1}{Z_1}q},\\ C = Z \left [\dfrac{C_0}{Z_0}(1-q) + \dfrac{C_1}{Z_1}q \right ]. \end{array}$$

It took a while to get to that result, but this a very fundamental equation in rasterization (which by the way is almost never explained anywhere. It is sometimes explained but the steps to get to that result are almost never provided) because it is used to interpolate vertex attributes which is a very important and common feature in rendering. What the equation says is that to interpolate a vertex attribute correctly, we first need to divide by the vertex attribute value by the z-coordinate of the vertex it is defined to, then linearly interpolate them using q (which in our case, will be the barycentric coordinates of the pixel on the 2D triangle), and then finally multiply the result by \(Z\), which is the depth of the point on the triangle, that the pixel overlaps (the depth of the point in camera space where the vertex attribute is being interpolated). Here is another version of the code we used in chapter three, that shows an example of perspective correct vertex attribute interpolation:

// compile with: // c++ -o raster3d raster3d.cpp // for naive vertex attribute interpolation and: // c++ -o raster3d raster3d.cpp -D PERSP_CORRECT // for perspective correct interpolation // (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 Vec3 &a, const Vec3 &b, const Vec3 &c) { return (c[0] - a[0]) * (b[1] - a[1]) - (c[1] - a[1]) * (b[0] - a[0]); } int main(int argc, char **argv) { Vec3 v2 = { -48, -10, 82}; Vec3 v1 = { 29, -15, 44}; Vec3 v0 = { 13, 34, 114}; Vec3 c2 = {1, 0, 0}; Vec3 c1 = {0, 1, 0}; Vec3 c0 = {0, 0, 1}; const uint32_t w = 512; const uint32_t h = 512; // project triangle onto the screen v0[0] /= v0[2], v0[1] /= v0[2]; v1[0] /= v1[2], v1[1] /= v1[2]; v2[0] /= v2[2], v2[1] /= v2[2]; // convert from screen space to NDC then raster (in one go) v0[0] = (1 + v0[0]) * 0.5 * w, v0[1] = (1 + v0[1]) * 0.5 * h; v1[0] = (1 + v1[0]) * 0.5 * w, v1[1] = (1 + v1[1]) * 0.5 * h; v2[0] = (1 + v2[0]) * 0.5 * w, v2[1] = (1 + v2[1]) * 0.5 * h; #ifdef PERSP_CORRECT // divide vertex-attribute by the vertex z-coordinate c0[0] /= v0[2], c0[1] /= v0[2], c0[2] /= v0[2]; c1[0] /= v1[2], c1[1] /= v1[2], c1[2] /= v1[2]; c2[0] /= v2[2], c2[1] /= v2[2], c2[2] /= v2[2]; // pre-compute 1 over z v0[2] = 1 / v0[2], v1[2] = 1 / v1[2], v2[2] = 1 / v2[2]; #endif 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) { Vec3 p = {i + 0.5, h - j + 0.5, 0}; 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]; #ifdef PERSP_CORRECT float z = 1 / (w0 * v0[2] + w1 * v1[2] + w2 * v2[2]); // if we use perspective correct interpolation we need to // multiply the result of this interpolation by z, the depth // of the point on the 3D triangle that the pixel overlaps. r *= z, g *= z, b *= z; #endif 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; }

Computing the sample death requires to use the reciprocal of the vertices z-coordinates. For this reason, we can pre-compute these values before we loop over all pixels (line 52). If we decide to use perspective correct interpolation, then the vertex attribute values are divided by the z-coordinate of the vertex they are associated to (lines 48-50). The following image shows on the left, an image computed without perspective correct interpolation, an image with (middle) and the content of the z-buffer (as a greyscale image. The closer the object is to the screen, the brighter). The difference is subtle though you can see in the left image how each color seems to roughly fill up the same area. This is due to the fact that colors in this case are interpolated within the "space" of the 2D triangle (as if the triangle was a flat surface parallel to the plane of the screen). However if you inspect the triangle vertices (and the depth buffer), you will notice that the triangle is not at all parallel to the screen (but oriented with a certain angle). In fact, because the vertex "painted" in green is closer to the camera than the other two, this part of the triangle fills up a larger part of the screen which is visible in the middle image (the green area is larger than the blue or red area). The image in the middle shows the correct interpolation and what you will get if you render this triangle with a graphics API such as OpenGL or Direct3D.

The difference between correct and incorrect perspective interpolation is even more visible when applied to texturing. In the next example, we assigned texture coordinates to the triangle vertices as vertex attributes, and use these coordinates to create a checkerboard patter to the triangle. Rendering the triangle with or without perspective correct interpolation is left as an exercise. In the image below you can see the result with. As you can see, it also matches an image of the same triangle with the same pattern rendered in Maya. Hopefully, our code so far seems to do the right thing. As with color, all you need to do (and this is true of all vertex attributes) is to divide the texture coordinates (which are generally denoted ST coordinates) by the z-coordinate of the vertex they are associated to, and then later in the code, multiply the texture coordinate interpolated value by Z. Here are the changes we made to the code:

... int main(int argc, char **argv) { Vec3 v2 = { -48, -10, 82}; Vec3 v1 = { 29, -15, 44}; Vec3 v0 = { 13, 34, 114}; ... Vec2 st2 = { 0, 0 }; Vec2 st1 = { 1, 0 }; Vec2 st0 = { 0, 1 }; ... #ifdef PERSP_CORRECT // divide vertex-attribute by the vertex z-coordinate c0[0] /= v0[2], c0[1] /= v0[2], c0[2] /= v0[2]; c1[0] /= v1[2], c1[1] /= v1[2], c1[2] /= v1[2]; c2[0] /= v2[2], c2[1] /= v2[2], c2[2] /= v2[2]; st0[0] /= v0[2], st0[1] /= v0[2]; st1[0] /= v1[2], st1[1] /= v1[2]; st2[0] /= v2[2], st2[1] /= v2[2]; // pre-compute 1 over z v0[2] = 1 / v0[2], v1[2] = 1 / v1[2], v2[2] = 1 / v2[2]; #endif ... for (uint32_t j = 0; j < h; ++j) { for (uint32_t i = 0; i < w; ++i) { ... if (w0 >= 0 && w1 >= 0 && w2 >= 0) { ... float s = w0 * st0[0] + w1 * st1[0] + w2 * st2[0]; float t = w0 * st0[1] + w1 * st1[1] + w2 * st2[1]; #ifdef PERSP_CORRECT float z = 1 / (w0 * v0[2] + w1 * v1[2] + w2 * v2[2]); // if we use perspective correct interpolation we need to // multiply the result of this interpolation by z, the depth // of the point on the 3D triangle that the pixel overlaps. s *= z, t *= z; #endif const int M = 10; // checkerboard pattern float p = (fmod(s * M, 1.0) > 0.5) ^ (fmod(t * M, 1.0) < 0.5); framebuffer[j * w + i][0] = (unsigned char)(p * 255); framebuffer[j * w + i][1] = (unsigned char)(p * 255); framebuffer[j * w + i][2] = (unsigned char)(p * 255); ... } } } ... return 0; }

Here is the result you should get:

What's Next?

In the final chapter of this lesson, we will briefly speak about ways of improving the rasterization algorithm (though we won't implement these techniques specifically) and explain how the final code of our lesson works.