Wenqi He · Research Software Engineer, NCSA, University of Illinois
Extracted from my presentation slides by Claude (Anthropic).
A vector graphics renderer takes a list of shapes and produces a pixel buffer — a table whose rows are pixels and whose columns are shapes, with each cell holding a coverage value. The problem looks deceptively simple. In practice it spans a surprisingly wide slice of computer science: computational geometry, numerical algorithms, GPU architecture, and signal processing. This is an algorithmic tour from first principles to the state of the art.
The fundamental question in vector rasterization is: given a closed curve and a point, is the point inside or outside? The standard answer is the winding number — the number of times the curve wraps counterclockwise around the point. For a closed curve \(\gamma\) and a point \(p = (x_p, y_p)\):
$$w(\gamma, p) = \frac{1}{2\pi} \oint_\gamma \frac{(x - x_p)\,dy - (y - y_p)\,dx}{(x - x_p)^2 + (y - y_p)^2}$$Two fill rules are in common use. The non-zero rule considers a point inside if \(w \ne 0\); the even-odd rule if \(w \bmod 2 \ne 0\). They agree whenever the curve has no self-intersections.
In practice, curves are composed of cubic polynomial segments \(\gamma_i\):
$$x_i(t) = a_1 t^3 + b_1 t^2 + c_1 t + d_1, \quad y_i(t) = a_2 t^3 + b_2 t^2 + c_2 t + d_2$$To compute the winding number, cast a ray leftward from \(p\) and count signed crossings. Solve \(y_i(t) = y_p\) for roots \(t_{ij}\), then sum the sign of each crossing that lies to the left of \(p\):
$$w(\gamma, p) = \sum_{\substack{(i,j) \\ 0 \le t_{ij} < 1 \\ x_i(t_{ij}) < x_p}} -\operatorname{sgn}(y_i'(t_{ij}))$$In practice the curve segments are Bézier curves. In the standard (power) basis, a cubic is just a cubic polynomial:
$$x(t) = a_1 t^3 + b_1 t^2 + c_1 t + d_1, \quad y(t) = a_2 t^3 + b_2 t^2 + c_2 t + d_2$$The more useful representation is the Bernstein basis. The degree-\(n\) Bernstein polynomials are:
$$b_{i,n}(t) = \binom{n}{i} t^i (1-t)^{n-i}$$A cubic Bézier is then a convex combination of four control points \(p_0, \ldots, p_3\):
$$B(t) = p_0\, b_{0,3}(t) + p_1\, b_{1,3}(t) + p_2\, b_{2,3}(t) + p_3\, b_{3,3}(t)$$The Bernstein form makes geometric properties — convex hull containment, endpoint interpolation, affine invariance — immediately apparent.
The most elegant way to evaluate a Bézier curve is de Casteljau’s algorithm, which replaces the polynomial with repeated linear interpolation:
def bezier(t, p0, p1, p2, p3):
p4 = lerp(p0, p1, t)
p5 = lerp(p1, p2, t)
p6 = lerp(p2, p3, t)
p7 = lerp(p4, p5, t)
p8 = lerp(p5, p6, t)
p9 = lerp(p7, p8, t)
return p9
def lerp(a, b, t):
return (1 - t) * a + t * b
Beyond evaluation, the intermediate points at parameter \(t\) are the control points of the two sub-curves on \([0, t]\) and \([t, 1]\). This underlies adaptive subdivision.
Adaptive subdivision recursively splits a curve until each piece is flat enough to be treated as a line segment:
def flatten(c):
if flatness(c) < eps:
vertices.push(c.p3)
else:
l, r = subdivide(c)
flatten(l)
flatten(r)
Once curves are flattened to polygons, scanline rasterization fills them. For each scanline \(y\), find all edge crossings, track the winding number left to right, and fill spans where the fill rule is satisfied:
def rasterize(poly, buf):
for y in poly.y_range:
xs = poly.intersect(y)
wn = 0
for i in range(len(xs) - 1):
wn += xs[i].sign
if not inside(wn): continue
for x in range(xs[i].x, xs[i+1].x):
buf[y, x] = over(poly.color(x, y), buf[y, x])
The naive approach recomputes all intersections for each scanline. The active edge list maintains only the edges currently crossing the scanline and updates them incrementally as \(y\) advances, reducing the work per scanline from \(O(n)\) to \(O(k)\) where \(k\) is the number of active edges:
def rasterize(polygon, buf):
xs = []
for y in polygon.y_range:
xs = update(xs, polygon.pending, y)
wn = 0
for i in range(len(xs) - 1):
wn += xs[i].sign
if not inside(wn): continue
for x in range(xs[i].x, xs[i+1].x):
buf[y, x] = over(polygon.color(x, y), buf[y, x])
Modern rendering engines — Chrome’s RenderingNG, Mozilla’s WebRender — off-load rasterization to the GPU. The pipeline proceeds as follows: the CPU uploads data and shader programs over PCIe to GPU memory; a vertex shader transforms each vertex in parallel; hardware rasterization interpolates attributes across triangles; a fragment shader computes the final color per pixel in parallel, taking interpolated vertex outputs as input; testing and blending composite the results into the framebuffer.
Two broad strategies have emerged for vector graphics on this pipeline, depending on whether one starts from shapes or from pixels.
Starting from shapes. Approximate each shape with triangles. The fragment shader is trivial — every fragment inside a triangle receives the fill color. The difficulty lies in tessellation. For a circle \(x(t) = \cos t,\ y(t) = \sin t\), sample vertices \((x_i, y_i) = (\cos t_i, \sin t_i)\) and fan-triangulate. Curves are approximated to a fixed tolerance, which breaks down at high zoom levels.
Instead of approximating, encode the curve as an implicit function evaluated per fragment. For a circle of radius \(r\), \(f(x, y) = x^2 + y^2 - r^2 \le 0\):
# in parallel
def fragment(x, y):
if x * x + y * y - r * r <= 0:
return color
Loop and Blinn (2005) generalized this to arbitrary cubic curves using a GPU-friendly implicit formulation, enabling resolution-independent rendering without tessellation. Kokojima et al. (2006) extended this to deformable vector objects.
NVIDIA’s NV_path_rendering (NVPR) exposes path rendering as a first-class GPU API, using stencil-buffer operations to accumulate winding numbers and cover passes to fill resolved pixels.
Starting from pixels. A third approach works entirely in compute kernels. Each pixel independently computes its winding number by iterating over curve segments:
# in parallel
def pixel_kernel(x, y):
for c in curves:
compute winding number w
compute signed distance d
blend based on w, d
Two acceleration structures prune the candidate curves per pixel: a regular spatial grid (Nehab and Hoppe, 2008) and a shortcut tree (Ganacim et al., 2014).
An additional optimization precomputes the winding parity of each tile from its boundary alone, allowing most pixels within the tile to skip the per-curve iteration entirely.
Rasterization is sampling: the rasterizer converts a continuous 2D signal (the filled shape) to discrete pixel values at a finite sampling frequency. The Nyquist–Shannon theorem applies: if the spatial frequency of the shape boundary reaches half the pixel density, aliasing results — for instance, tiles whose width equals half a pixel appear not as stripes but as a solid color flickering between the two tile colors. The remedy is bandlimiting: convolving the signal with a filter \(g\) before sampling:
$$I(x, y;\, \Theta) = \iint_U g(u, v)\, f(x - u,\, y - v;\, \Theta)\, du\, dv$$In the rasterization context this means computing anti-aliased coverage rather than a binary inside/outside value.
The same integral formulation \(I(x, y; \Theta)\) makes differentiable rendering possible. Li et al. (2020) show how to propagate gradients through a rasterizer with respect to shape parameters \(\Theta\), enabling vector graphics to be optimized by gradient descent. The key challenge is that the rasterization function is discontinuous at shape boundaries; automatic differentiation handles the smooth interior but requires special treatment at edges.