Implementing a tiny CPU rasterizer

Part 3: Interpolating colors


The code for this part is available in this commit.

Last time we finally drew a triangle on screen, but it was filled with one fixed color. Today our goal is to fill it with a gradient of colors instead. And for that, we need something called attribute interpolation.

Contents

Attributes

In graphics, a vertex attribute is just any data that is stored per-vertex. Right now our only vertex attribute is 3D position; other examples of attributes include colors, normals, texture coordinates, joint IDs and weights, etc.

To fill our triangle with a gradient, we need to set the colors of our vertices to some values, and then our engine will interpolate the colors in the pixels, thus forming a gradient. On the GPU this happens somewhere between the vertex and the fragment shaders.

Later we'll need even more attributes, so let's define a convenience wrapper for it:

#pragma once

#include <cstdint>

namespace rasterizer
{

    template <typename T>
    struct attribute
    {
        void const * pointer = nullptr;
        std::uint32_t stride = sizeof(T);

        // Get attribute by vertex index
        T const & operator[] (std::uint32_t i) const
        {
            return *(T const *)((char const *)(pointer) + stride * i);
        }
    };

}

Here we have the pointer to the attribute data, and something called a stride. Stride is the distance in bytes between two values of this attribute corresponding to two consecutive vertices. This makes vertex storage more flexible:

This is how it is typically done in GPUs, and this is the convention we will adopt. Note that in our attribute class stride is equal to the size of the attribute by default, i.e. it assumes all values of this particular attribute are packed in a separate array. Also note that in case of stride = 0 we effectively can pass a single value as the common attribute value for all vertices, which can be really useful. Some graphics APIs like OpenGL don't allow a stride value of 0, and instead interpret it as our default value of stride = sizeof(T).

Now, let's tweak our mesh class to add the color attribute:

#pragma once

#include <rasterizer/types.hpp>
#include <rasterizer/attribute.hpp>

namespace rasterizer
{

    struct mesh
    {
        attribute<vector3f> positions = {};
        attribute<vector4f> colors = {};
        std::uint32_t vertex_count = 0;
    };

}

The syntax = {} means the field will be initialized with default values. I've chosen to use vectors of 4 floating-points for colors for simplicity, but you can instead use color4ub for them and add the conversion in the rasterizer, it doesn't matter that much for us.

Barycentric coordinates: theory

Who the hell is Bary?

Now we need to figure our how to actually compute interpolated colors. Luckily, in a triangle there is a unique way of doing this! If we want a linear color gradient, we could assume that the function mapping pixel positions to colors is linear, or, technically, affine, i.e. of the form \(Ax+By+C\). Then, if we know the value of this function at 3 different points, we can make a 3x3 system of linear equations for the unknowns \(A,B,C\), which by general linear algebra almost always has a solution (unless the three points lie on a straight line). So, our color-for-pixel function is uniquely determined by the values at 3 vertices!

Though, the usual way of doing this doens't solve this system, but instead computes something called the barycentric coordinates. Given three vertices \(V_0,V_1,V_2\), the barycentric coordinates of some point \(P\) are special numbers \(\lambda_0,\lambda_1,\lambda_2\) such that

\[ \lambda_0+\lambda_1+\lambda_2=1 \] \[ \lambda_0V_0+\lambda_1V_1+\lambda_2V_2=P \]

That is, point \(P\) is an affine combination of points \(V_0,V_1,V_2\) with coefficients \(\lambda_0,\lambda_1,\lambda_2\). They are also unique, which follows from a similar argument of constructing a 3x3 linear system for them. The condition that their sum is equal to 1 is essential, because otherwise the sum \(\lambda_0V_0+\lambda_1V_1+\lambda_2V_2\) is not well-defined: it depends on a choice of coordinate system origin.

We can plug \(P=V_0\), and since \(V_0=1\cdot V_0 + 0\cdot V_1 + 0\cdot V_2\), and since the barycentric coordinates are unique, we see that for this point the coordinates must be simply \((1,0,0)\). Similarly, they are \((0,1,0)\) for \(V_1\) and \((0,0,1)\) for \(V_2\).

Now, if we have some value at the three vertices \(V_0,V_1,V_2\), for example the colors \(C_0,C_1,C_2\), we can interpolate them with the formula

\[ C = \lambda_0C_0+\lambda_1C_1+\lambda_2C_2 \]

This is precisely what's called linear interpolation inside a triangle: it is interpolation using barycentric coordinates. Notice that for the case \(P=V_0\) the coordinates are \((1,0,0)\), and the resulting color will be \(C=1\cdot C_0+0\cdot C_1+0\cdot C_2=C_0\), the original color of this vertex, and similarly for \(V_1\) and \(V_2\).

Now, how do we compute these coordinates for the current pixel we're rasterizing? There are many different ways of doing it. One is purely algebraic: take the definition of the barycentric coordinates, and subtract, for example, \(V_0\) from both sides:

\[ \lambda_0V_0+\lambda_1V_1+\lambda_2V_2\color{blue}{-V_0}=P\color{blue}{-V_0} \]

Now, since, \(\lambda_0+\lambda_1+\lambda_2=1\), we can rewrite it as

\[ \lambda_0V_0+\lambda_1V_1+\lambda_2V_2-\color{blue}{(\lambda_0+\lambda_1+\lambda_2)}\cdot V_0=P-V_0 \] \[ \lambda_0V_0+\lambda_1V_1+\lambda_2V_2\color{blue}{-\lambda_0V_0-\lambda_1V_0-\lambda_2V_0}=P-V_0 \] \[ \lambda_1(V_1\color{blue}{-V_0})+\lambda_2(V_2\color{blue}{-V_0})=P-V_0 \]

This is an equation on 2D points, so it is effectively 2 scalar equations, with 2 unknowns \(\lambda_1,\lambda_2\). So it's a 2x2 linear system, and we can directly use Cramer's rule to solve it. This rule gives explicit formulas for the unknowns in terms of some determinants. Let's get the formula for \(\lambda_2\):

\[ \lambda_2 = \frac{\det(V_1-V_0,P-V_0)}{\det(V_1-V_0,V_2-V_0)} \]

Do you recognize this determinants? The numerator is exactly what we used to figure out if the point \(P\) is to the left of an edge \(V_0V_1\), while the denominator is exactly what we used to figure out the orientation of our triangle!

Because none of the three triangle vertices are more special than the others, it must be that the formulas for \(\lambda_0\) and \(\lambda_1\) can be found by making a cyclic permutation of the vertices \(V_0V_1V_2\), and we'll get

\[ \lambda_1 = \frac{\det(V_0-V_2,P-V_2)}{\det(V_0-V_2,V_1-V_2)} \] \[ \lambda_0 = \frac{\det(V_2-V_1,P-V_1)}{\det(V_2-V_1,V_0-V_1)} \]

Again, the numerators look like what we've already get, but the denominators look different, but they're actually still the same. This follows from various properties of determinants, or more simply from the fact that all 3 denominators equal 2x of the area of our triangle.

Barycentric coordinates: code

This was quite a bit of formulas, but the change to our rasterizer code is really tiny. When rasterizing a triangle, grab its colors together with vertex positions:

auto v0 = command.transform * as_point(command.mesh.positions[vertex_index + 0]);
auto v1 = command.transform * as_point(command.mesh.positions[vertex_index + 1]);
auto v2 = command.transform * as_point(command.mesh.positions[vertex_index + 2]);

auto c0 = command.mesh.colors[vertex_index + 0];
auto c1 = command.mesh.colors[vertex_index + 1];
auto c2 = command.mesh.colors[vertex_index + 2];

Then, when drawing a pixel, compute its barycentric coordinates and interpolate the colors. Remember: the barycentric coordinates are just fractions composed of determinants that we have already at hand:

if (det01p >= 0.f && det12p >= 0.f && det20p >= 0.f)
{
    float l0 = det12p / det012;
    float l1 = det20p / det012;
    float l2 = det01p / det012;

    color_buffer.at(x, y) = to_color4ub(l0 * c0 + l1 * c1 + l2 * c2);
}

Now let's tweak our mesh definition in main.cpp:

vector3f positions[] =
{
    {  0.f,   0.f, 0.f},
    {100.f,   0.f, 0.f},
    {  0.f, 100.f, 0.f},
};

vector4f colors[] =
{
    {1.f, 0.f, 0.f, 1.f},
    {0.f, 1.f, 0.f, 1.f},
    {0.f, 0.f, 1.f, 1.f},
};

for (int i = 0; i < 100; ++i)
    draw(color_buffer,
        draw_command{
            .mesh = {
                .positions = {positions},
                .colors = {colors},
                .vertex_count = 3,
            },
            .transform = {
                1.f, 0.f, 0.f, mouse_x + 100.f * (i % 10),
                0.f, 1.f, 0.f, mouse_y + 100.f * (i / 10),
                0.f, 0.f, 1.f, 0.f,
                0.f, 0.f, 0.f, 1.f,
            },
        }
    );

Today's article was more formulas as less code than usual. This will happen again next time, when we start working on actual 3D graphics and tackle perspective projection!

\(\leftarrow\) Part 2: Drawing a triangle Source code for this part Part 4: Changing perspective \(\rightarrow\)