Triangle rasterization

A GPU draws a frame by turning triangles into pixels billions of times per second, and the algorithm it uses fits in thirty lines: for each pixel near the triangle, ask three edge functions whether the pixel's centre is inside; if so, use the same three numbers to blend the vertices' attributes and to test depth. This lesson builds that loop from scratch — including the parts that are easy to get almost right: shared edges, perspective, depth and the viewport.

It matters for neural rendering too. 3D Gaussian Splatting (Kerbl et al. 2023) ships its own tile-based software rasterizer in CUDA, nvdiffrast (Laine et al. 2020) makes the hardware one differentiable, and every "render a mesh to supervise a network" pipeline relies on exactly these rules.

1. Where the rasterizer sits

Each vertex arrives with a screen position (x, y) in pixels (y down), a depth z after the perspective divide, 1/w (kept from the divide, needed for perspective-correct interpolation), and attributes: colour, UVs, normals.

2. Edge functions

For an edge from a to b, the edge function

E_ab(p) = (b.x − a.x)·(p.y − a.y) − (b.y − a.y)·(p.x − a.x)

is the 2D cross product of b − a and p − a: twice the signed area of the triangle (a, b, p). It is zero on the line through a and b, positive on one side, negative on the other, and affine in p — moving one pixel right adds the constant −(b.y − a.y). Pineda (1988) built the first edge-function rasterizer on exactly that property, and every GPU since does the same.

For a triangle (v0, v1, v2) the total E_01(v2) is twice its signed area. When it's positive, the inside is where all three edge functions are positive. And the three values are the barycentric coordinates, up to the area:

λ0 = E_12(p) / area     λ1 = E_20(p) / area     λ2 = E_01(p) / area

— each vertex's weight is the area of the sub-triangle opposite it.

#include <print>

struct P { float x, y; };

float edge(P a, P b, P p) { return (b.x - a.x) * (p.y - a.y) - (b.y - a.y) * (p.x - a.x); }

int main() {
    P v0{0, 0}, v1{8, 0}, v2{0, 8};
    float area = edge(v0, v1, v2);
    std::println("2 x area = {}", area);
    for (P p : {P{2, 4}, P{0.5f, 0.5f}, P{4, 4}, P{6, 6}}) {
        float w0 = edge(v1, v2, p), w1 = edge(v2, v0, p), w2 = edge(v0, v1, p);
        std::println("p = ({}, {}): w = {:>4} {:>4} {:>4}   barycentrics {:.3f} {:.3f} {:.3f}   {}", p.x, p.y,
                     w0, w1, w2, w0 / area, w1 / area, w2 / area,
                     w0 > 0 && w1 > 0 && w2 > 0 ? "inside" : (w0 == 0 || w1 == 0 || w2 == 0) ? "on an edge" : "outside");
    }
}

With y pointing down, positive area means the vertices go clockwise on screen — which is counter-clockwise in OpenGL's y-up NDC, OpenGL's default front face. Every axis flip between spaces flips that sign; that's why winding bugs are so common.

3. Bounding box and pixel centres

Pixel (i, j) is covered when its centre (i + 0.5, j + 0.5) is inside. Only centres in the triangle's bounding box can be, so loop from ceil(min_x − 0.5) to floor(max_x − 0.5) (and the same for y), clamped to the viewport. Clamp in float before converting to int: after guard-band clipping a triangle may span ±10⁵ pixels, and an out-of-range float-to-int conversion is undefined behaviour.

4. The top-left fill rule

Meshes share edges, and edges often pass exactly through pixel centres (any axis-aligned edge at a half-integer, any diagonal through the lattice). A centre with E = 0 must belong to exactly one of the two triangles:

• E ≥ 0 for both → drawn twice: seams with blending, double counting in coverage, splats or particle passes.
• E > 0 for both → drawn never: one-pixel cracks.

Direct3D's and OpenGL's answer is the top-left rule: a centre on an edge is inside only if that edge is a top edge (exactly horizontal, the triangle below it) or a left edge (the triangle to its right). Two neighbours see their shared edge in opposite directions, so it's a left edge for exactly one of them. With y down and positive area, for the edge a → b:

left edge:  b.y − a.y < 0            top edge:  b.y − a.y == 0  and  b.x − a.x > 0

The == 0 test only means something with exact arithmetic. Hardware snaps vertices to a fixed-point grid (typically 8 sub-pixel bits) and evaluates the edge functions in integers, so the rule is applied consistently. With floats it works when vertices are dyadic (as in the tests here); a production software rasterizer snaps, or adds a bias of −1 in fixed point to non-top-left edges.

#include <print>
#include <vector>

struct V { float x, y; };
float edge(V a, V b, float px, float py) { return (b.x - a.x) * (py - a.y) - (b.y - a.y) * (px - a.x); }
bool top_left(V a, V b) { return b.y - a.y < 0 || (b.y == a.y && b.x > a.x); }

// Adds 1 to every pixel the triangle covers (positive area assumed).
void raster(V v0, V v1, V v2, int W, int H, std::vector<int>& count, bool fill_rule) {
    bool tl[3] = {top_left(v1, v2), top_left(v2, v0), top_left(v0, v1)};
    for (int y = 0; y < H; ++y)
        for (int x = 0; x < W; ++x) {
            float px = x + 0.5f, py = y + 0.5f;
            float w[3] = {edge(v1, v2, px, py), edge(v2, v0, px, py), edge(v0, v1, px, py)};
            bool in = true;
            for (int e = 0; e < 3; ++e)
                in = in && (fill_rule ? (w[e] > 0 || (w[e] == 0 && tl[e])) : w[e] >= 0);
            count[y * W + x] += in;
        }
}

int main() {
    for (bool rule : {false, true}) {
        std::vector<int> c(36, 0);
        raster({0, 0}, {6, 0}, {6, 6}, 6, 6, c, rule);   // A
        raster({0, 0}, {6, 6}, {0, 6}, 6, 6, c, rule);   // B
        std::println("{}:", rule ? "top-left rule" : "inclusive edges (w >= 0)");
        for (int y = 0; y < 6; ++y) {
            for (int x = 0; x < 6; ++x) std::print("{} ", c[y * 6 + x]);
            std::println();
        }
    }
}

5. Barycentrics, depth and perspective-correct attributes

Depth. After the divide, NDC z is an affine function of screen position for a planar triangle, so z = λ0·z0 + λ1·z1 + λ2·z2 is exact. That's what the z-buffer stores and compares.

Attributes are not. A colour or UV is affine on the triangle in 3D; the perspective divide squeezes the far part of the triangle, so equal steps on screen are unequal steps on the surface. What is affine on screen is attr / w and 1 / w (divide the 3D-affine function by the same w that the projection divides by). So interpolate both and divide per pixel:

q_i = λ_i / w_i          attr(p) = Σ q_i·attr_i / Σ q_i

#include <print>

int main() {
    // A checkerboard texture along an edge: u from 0 (near, w = 1) to 1 (far, w = 4).
    float uA = 0, wA = 1, uB = 1, wB = 4;
    std::println("  t    affine u   correct u   checker(affine) checker(correct)");
    for (int k = 0; k <= 8; ++k) {
        float t = k / 8.0f;
        float affine = (1 - t) * uA + t * uB;
        float q0 = (1 - t) / wA, q1 = t / wB;
        float correct = (q0 * uA + q1 * uB) / (q0 + q1);
        auto checker = [](float u) { return static_cast<int>(u * 4) % 2 ? "##" : ".."; };
        std::println("{:5.3f}  {:8.3f}   {:9.3f}        {}              {}", t, affine, correct, checker(affine),
                     checker(correct));
    }
}

The correct column bunches the far squares together, as perspective should. Affine interpolation is the wobbling, bending textures of the original PlayStation, which had no per-pixel divide. (Screen-aligned primitives — UI, Gaussian splats rendered as camera-facing quads — have constant w and get away with affine interpolation.)

6. The z-buffer

Clear depth to the far value (+∞ here, 1.0 in a normalized buffer), then per fragment: if z < depth[p], write the colour and the depth. Forgetting the depth write lets every later fragment pass; using <= instead of < lets a coplanar second pass overwrite (sometimes wanted: decals, multi-pass with LESS_EQUAL). With reverse-Z (see the camera lesson) clear to 0 and use GREATER.

The z-buffer resolves opaque visibility in any order. Transparency needs order: blending is not commutative — which is why 3D Gaussian Splatting sorts its splats per tile by depth and alpha-blends front to back instead of using a depth test.

Back-face culling is free once you have the signed area: skip triangles whose area has the back-facing sign. It halves the work for closed meshes, and is wrong for open/two-sided geometry, and for any transform with negative determinant (a mirror), which reverses every winding.

7. Clipping: why the rasterizer may assume w > 0

A vertex behind the camera has w < 0; dividing by it mirrors the vertex through the centre of the screen, and a triangle with one vertex behind the camera would be set up as a huge wrong sheet. So triangles are clipped against the near plane in clip space, before the divide (Sutherland–Hodgman against z ≥ −w, yielding up to two triangles). Clipping against the other planes is optional: GPUs use a guard band — they rasterize triangles that extend far off-screen and rely on the clamped bounding box, clipping only when coordinates would overflow the fixed-point range. That's why your bounding box clamp is not an optimization but a correctness requirement.

8. How GPUs actually do it

• Fixed point. Vertices snap to 1/256 pixel; edge functions are exact integer arithmetic, so the fill rule is watertight.
• Incremental and hierarchical. Edge functions are affine, so stepping one pixel is one add. Traversal is hierarchical: evaluate the three edges at the corners of an 8×8 tile and reject (all outside one edge) or accept (all inside all edges) whole tiles at once.
• Quads. Pixels are shaded in 2×2 quads so the shader can take screen-space derivatives; partially covered quads run helper lanes. A one-pixel triangle costs a whole quad — the "micropolygon" problem that led Unreal's Nanite and 3DGS to software rasterizers built around tiles and compute shaders.
• Early-Z, MSAA. The depth test runs before shading when the shader doesn't write depth; MSAA evaluates coverage at several samples per pixel but shades once.

Fabian Giesen's "A trip through the Graphics Pipeline" (resources) walks through all of this on real hardware.

In an interview. "Write a function that rasterizes a triangle." Start with the bounding box and three edge functions, then volunteer the three things they're listening for: pixel centres at +0.5, the top-left rule for shared edges, and perspective-correct interpolation through 1/w. Bonus: incremental evaluation and tiles.

In Python (NumPy / PyTorch). Vectorize over the bounding box: xs, ys = meshgrid(arange(x0, x1+1) + 0.5, arange(y0, y1+1) + 0.5), compute the three edge arrays at once, build the inside mask, and torch.where it into the image. Differentiable rasterizers (SoftRas, nvdiffrast) replace the hard inside test or add analytic antialiasing so gradients flow to vertex positions.

Summary

• E_ab(p) = 2D cross product; three edge values = coverage test and barycentrics × area.
• Sample at (i + 0.5, j + 0.5); clamp the bounding box to the viewport, in float.
• Top-left rule: zero counts only for top (horizontal, interior below) or left (y-down, positive area: dy < 0) edges — every pixel exactly once.
• Depth interpolates linearly; attributes via attr/w and 1/w.
• Depth test writes depth too. Cull by the sign of the area. Clip at the near plane before the divide.