Rasterisation: Edge Functions and the Z-Buffer

Rasterisation takes the first of the two possible loops. For each triangle, find the pixels it covers, and write to them. The outer loop is over geometry and the inner loop is over the pixels of one triangle.

That choice fixes the cost model immediately. Total work is proportional to the number of triangles plus the total area they cover, and it never depends on how many triangles are in the scene behind the one being drawn, except through the covered area.

It also creates the problem the rest of this lesson solves. Because triangles are processed independently and in arbitrary order, a triangle can be drawn into a pixel that a nearer triangle has already claimed, or will claim later. Something must arbitrate, and doing it by sorting triangles turns out to be both slow and, for interpenetrating geometry, impossible.

The inner loop needs to answer one question per pixel: is this point inside the triangle? Juan Pineda gave the formulation that hardware still uses, in "A Parallel Algorithm for Polygon Rasterization", ACM SIGGRAPH Computer Graphics, volume 22, issue 4, 1988, pages 17 to 20.

For a directed edge from A to B, define E(P) = (B.x - A.x) * (P.y - A.y) - (B.y - A.y) * (P.x - A.x). This is the cross product's z component. Its sign says which side of the edge P lies on, and it is zero exactly on the line.

A point is inside the triangle when all three edge functions share a sign. Three multiplications and some subtractions, no division, no branching on geometry.

The property that made it a hardware algorithm is that E is linear in the pixel coordinates. Stepping one pixel right adds a constant. Coverage for a whole block of pixels can be evaluated by addition alone, in parallel, from one setup per triangle.

The whole idea fits in one function, including interpolation.

def edge(a, b, p):
    return (b[0] - a[0]) * (p[1] - a[1]) - (b[1] - a[1]) * (p[0] - a[0])

def raster(v0, v1, v2, depth, colour_of, w, h):
    area = edge(v0, v1, v2)
    if area == 0:
        return                            # degenerate triangle
    xs = [v0[0], v1[0], v2[0]]
    ys = [v0[1], v1[1], v2[1]]
    for y in range(max(0, int(min(ys))), min(h, int(max(ys)) + 1)):
        for x in range(max(0, int(min(xs))), min(w, int(max(xs)) + 1)):
            p = (x + 0.5, y + 0.5)        # sample at the pixel centre
            w0, w1, w2 = edge(v1, v2, p), edge(v2, v0, p), edge(v0, v1, p)
            if (w0 >= 0) == (w1 >= 0) == (w2 >= 0):
                b0, b1, b2 = w0 / area, w1 / area, w2 / area   # barycentric
                z = b0 * v0[2] + b1 * v1[2] + b2 * v2[2]
                if z < depth[y][x]:       # nearer than what is there
                    depth[y][x] = z
                    colour_of(x, y, b0, b1, b2)

The bounding box keeps the loop off the rest of the screen. The three edge values, divided by the total area, are the barycentric weights, so the same numbers that decided coverage also interpolate depth and every vertex attribute.

The three-line depth test in that code is the second key idea, and it is doing more work than it looks.

Keep a buffer of one depth value per pixel, initialised to infinity. When a fragment is produced, compare its interpolated depth to the stored one. If it is nearer, overwrite both depth and colour. If not, discard it.

Edwin Catmull described the technique in his 1974 University of Utah doctoral thesis, A Subdivision Algorithm for Computer Display of Curved Surfaces. Wolfgang Strasser had described the same idea slightly earlier in his own thesis, and the two arrived at it independently.

What makes it structural rather than merely convenient: the result is order independent. Triangles may arrive in any sequence and the image is the same. There is no sort, no global data structure, and no communication between triangles. That is precisely what allows thousands of them to be processed at once.

Two stages are programmable and the rest are fixed function, and the split is not arbitrary. The vertex stage runs once per vertex; the fragment stage runs once per covered sample. Everything between them is coverage and interpolation, which is the same arithmetic for every scene and is therefore worth casting into silicon.

The ratio between those two rates is the thing to hold on to. A full-screen triangle at 4K produces around eight million fragments from three vertices. Fragment work dominates almost everything, which is why shader cost is usually discussed per pixel and why anything that reduces covered area pays immediately.

Order independence has a cost. If a distant wall is drawn first and a near one second, the wall's fragments were shaded and then thrown away. Shading a fragment that never reaches the screen is overdraw, and in a scene with depth complexity of four it can waste three quarters of the fragment budget.

The fix is to move the depth test before shading rather than after. Coverage and interpolated depth are known at triangle setup, so a fragment that already loses the depth comparison can be discarded before its shader runs. This is early depth testing, and it is automatic on modern hardware, with one important condition: it is disabled when the fragment shader itself writes depth or discards fragments, because then the hardware cannot know the outcome in advance.

Drawing roughly front to back therefore matters, and so does keeping shaders free of discard when it is not needed.

The z-buffer resolves visibility for opaque surfaces exactly. For transparency it does not, and the reason is worth stating precisely.

Blending a translucent surface over what is behind it is not a commutative operation: red over blue is a different colour from blue over red. So the correct result depends on the order the surfaces are composited, and the z-buffer's whole virtue was that order does not matter.

The usual remedy is to draw opaque geometry first with depth writing on, then sort the translucent geometry back to front and draw it with depth writing off. That is a per-object sort, so it fails exactly where sorting always fails: two translucent surfaces that intersect have no correct ordering, and neither does a cycle of three overlapping quads.

Transparency is the standing exception in a rasteriser. It is handled by convention and approximation rather than by the depth buffer.

The last row is the one that shapes the next lesson. A rasteriser answers one question extremely well, and it is a narrow question: what is nearest along the rays leaving this viewpoint.

A shadow needs to know what is visible from the light. A mirror needs to know what is visible from a point on the surface, in a direction that depends on the surface. Neither is expressible in the loop above without setting up an entirely new render from a new viewpoint. That limitation, not raw speed, is the reason the other architecture exists.

Pull the properties together and the hardware follows almost forcibly.

Coverage is three linear functions, so it is addition in a loop. Interpolation reuses the same three numbers. Visibility is a local comparison against one memory location. Nothing in the inner loop needs to know about any other triangle, any other pixel, or the scene as a whole.

An algorithm with no global data structure and no inter-element communication is exactly what massively parallel hardware executes well, and that is the real reason rasterisation dominated real-time graphics for thirty years. It was not that it approximated better. It was that its work decomposes into millions of independent pieces.

The next lesson takes the opposite loop, and the first thing that appears is the global data structure this one so carefully avoided.