Stratified Pixel Sampling and the Aliasing of Regular Grids

A pixel is not a point. It is an area, and the colour that belongs to it is the average of the radiance arriving across that area, weighted by the reconstruction filter. A renderer that fires one ray through the pixel centre computes the radiance at one point and reports it as the average, which is a sampling approximation. Everything that follows — the moire on a checkerboard floor receding into the distance, the stair-stepped edges of a silhouette, the crawling speckle on a moving surface — is the consequence of that approximation, and the shape of the sample pattern determines which errors survive.

Aliasing is a spectrum problem, not a resolution problem

Take a plane covered in a square-wave pattern of alternating black and white stripes of spatial period $p$. Project it to the image plane, where the stripe period becomes $q$ pixels. Sample the image on a regular grid of spacing 1 pixel and reconstruct by treating each sample as the pixel’s value. When $q > 2$ the stripes are reproduced faithfully. As $q$ falls below 2, the sampled pattern’s apparent period diverges from the true period and eventually reverses — the classic wagon-wheel effect — because the sampled signal contains a component at a frequency that is an alias of the true one. With a pure sine wave of frequency $f$ cycles per pixel, sampling at unit spacing produces a discrete spectrum with energy at every $f + k$ for integer $k$, and any of those that fall inside the interval $[-1/2, 1/2]$ cycles per pixel are indistinguishable from the true signal.

The reconstruction step is where the error becomes visible, because a box filter — treating each sample as its own pixel — has a frequency response that is not zero outside the passband. It passes the aliases through with substantial weight. Crow’s 1977 analysis of the problem framed it exactly this way: the renderer is sampling a function whose spectrum is not band-limited, so no sampling rate eliminates aliasing entirely, and the only options are to pre-filter the function (compute the pixel’s true integral) or to post-filter the samples with a better reconstruction kernel.

Pre-filtering is the physically motivated answer and the expensive one: computing the exact average of radiance over a pixel requires integrating over an area of the scene, which for a ray tracer means either analytic area coverage — available for simple primitives and impractical in general — or an approximation by taking many samples. The number of samples needed to make the residual error small is governed by the variance of the estimator, and reducing variance is the entire purpose of the sampling pattern.

What stratification buys and what it does not

Stratified sampling divides the pixel’s area into $N$ cells of equal area and places exactly one sample in each, with the sample’s position within its cell random or pseudo-random. The estimator remains unbiased — each sample’s expected contribution is still the pixel’s average — but the samples are now negatively correlated in a useful way: a sample that lands in a bright part of the cell pushes its neighbours’ cells to be sampled elsewhere, so the errors across cells partially cancel instead of adding.

The benefit is quantifiable for a smooth integrand. With $N$ independent uniform samples the error falls as $N^{-1/2}$, the standard Monte Carlo rate, because the estimator’s variance is the integrand’s variance divided by $N$. With $N = n^2$ stratified samples in two dimensions the error falls as approximately $N^{-3/4}$ for a function with bounded variation, because the variance of the stratified estimator is bounded by a term proportional to $N^{-3/2}$. The exponent improves because stratification removes the low-frequency component of the sampling noise: the within-cell mean is exact by construction, so only the variation inside each cell remains.

$$\text{variance} ;\propto; \frac{1}{N^{2}} \quad \text{for a partition of } N \text{ cells, versus} \quad \frac{1}{N} \text{ for independent samples.}$$

That reasoning is the standard one and it has a sharp precondition: the integrand must be smooth relative to the cell size. Across a silhouette edge the integrand jumps from full radiance to zero, and no stratification helps inside a cell that straddles the edge, because the cell’s content is not well approximated by a smooth function. The variance reduction from stratification on a hard edge is therefore whatever the edge’s area fraction produces, and the visible result is that stratified sampling cleans up smooth gradients — the shading across a sphere, the falloff of a light — much better than it cleans up a stair-stepped silhouette. The two artifacts have different causes and different remedies, which is why a renderer with heavy stratification still shows jaggies and needs either more samples or a reconstruction filter that widens across the edge.

The pattern is the filter

A regular grid of samples, one per pixel centre, is the worst pattern for aliasing because all its energy is concentrated at one spatial frequency. Jittering each sample uniformly within its pixel moves that energy into a broad spectrum — the jitter’s own spectrum — which converts coherent moire into incoherent noise. This is the standard trade: aliasing is a structured error that the eye interprets as pattern, noise is a random error that the eye interprets as grain. At equal variance, noise is far less objectionable, and it also converges away with more samples while a moire pattern does not.

import math, random

def stratified_pixel(n, rng):
    """n x n stratified samples in [0,1)^2, one per cell of an n x n grid."""
    jitter = []
    inv = 1.0 / n
    for j in range(n):
        for i in range(n):
            jitter.append(((i + rng.random()) * inv,
                           (j + rng.random()) * inv))
    return jitter

def energy_at(freq, samples):
    """Normalised magnitude of the sample pattern's Fourier component at `freq`,
       a crude measure of how much a grid frequency will alias."""
    re = sum(math.cos(2 * math.pi * freq * (x + y)) for x, y in samples)
    im = sum(math.sin(2 * math.pi * freq * (x + y)) for x, y in samples)
    return math.hypot(re, im) / len(samples)

The function above computes a single Fourier coefficient of the sample set, which is what the reconstruction filter sees: a pattern with a large coefficient at some frequency will imprint that frequency on the output whenever the scene contains content near it. A regular grid with $n = 1$ has coefficient 1 at every integer frequency, which is the algebraic statement of “a single centred sample reproduces every alias at full strength”. Stratification spreads the coefficient across frequencies, and the coefficient’s magnitude falls as $1/N$ for the stratified set, matching the variance argument.

Because the eye is sensitive to oriented structure and insensitive to uncorrelated noise, the choice between patterns is not symmetric even when their variances are equal. A pattern with residual structure at low frequencies — visible as blotches a few pixels across — is worse than a pattern with the same energy at high frequencies, which reads as fine grain. This is why a stratified grid with one sample per cell is not the end of the story: the grid’s regularity leaves energy at the cell-boundary frequencies, and at high sample counts those become visible as a faint weave.

Scrambling, and why rotation is not free

A stratified grid can be decorrelated between pixels by applying a per-pixel offset. Offsetting the whole pattern by a random amount in $[0,1)^2$ and wrapping it modulo one preserves the stratification exactly — every cell still holds exactly one sample — and randomises the pattern’s phase across pixels, which breaks the long-range regularity that produces structured artifacts. This is the standard jitter, and it costs two random numbers per pixel rather than two per sample.

A per-pixel rotation does not preserve the stratification. Rotating a square grid by an angle that is not a multiple of 90 degrees aligns the cells with a different lattice, and a sample that was inside its cell before rotation may fall outside after. Rotating and wrapping does preserve the count per cell in the sense that the total number of samples in the pixel is unchanged, but not the one-per-cell property, and the variance benefit degrades. A rotation of the grid is a valid sampling strategy — it is a form of Latin hypercube in the angular sense — but it is not a stratification of the pixel’s square, and the two should not be conflated.

The scrumbling technique that preserves stratification at every level is the one worth knowing: for a $2^k \times 2^k$ grid, apply a random permutation of the cells at each level of a binary subdivision, with a nesting property so that each quadrant of four cells contains one cell from each of the four sub-quadrants. The result is a pattern that is stratified at 1, 4, 16, and every intermediate count up to $4^k$, so a renderer using a progressive sample count can stop at any power of four and still have a stratified set. This is the mechanism behind progressive rendering that refines smoothly.

Where stratification meets the rest of the kernel

Pixel sampling decides where rays start; it does not decide how many, nor what they do afterwards. A path tracer draws further samples at each bounce — for the direction of the next ray, for the choice of a light — and the same stratification argument applies to each of those dimensions. Combining dimensions is where naive stratification fails: stratifying two dimensions separately and pairing the samples by index correlates them, so the pattern is stratified in each projection but the joint distribution has structure. The remedy is to stratify in a higher dimension, which costs $N^d$ samples for $d$ dimensions, or to use a sequence that is well distributed in every projection — the subject of low-discrepancy sequences and the correlation they leave behind.

The interaction with the intersection kernel is through coherence. Jittered samples within a pixel produce rays that are nearly parallel, so they hit the same leaves of the hierarchy, which is why ordered traversal works well with adaptive sampling and poorly with samples scattered across the whole image. Stratification at the pixel level is therefore the cheapest form of ray coherence available: it keeps the sample set dense in space, which helps both the sampling variance and the traversal’s cache behaviour.

Failure modes worth knowing

The first is a fixed pattern with a fixed sample count and a filter that is too narrow. Samples at exactly the cell centres reproduce a regular grid, whose aliases are strong and whose appearance depends on the scene’s spatial frequencies in a way that changes discontinuously as geometry moves. Any renderer that jitters for stills must jitter for animation, or the aliasing will crawl.

The second is a random number generator whose low-order bits are poor. Stratification consumes random numbers in a fixed pattern — two per cell, traversed in cell order — and a generator with a short period or with correlation between consecutive outputs will place the samples in a pattern that is not random within cells. The failure is subtle: the samples are still one per cell, so the stratification’s variance bound holds, but the within-cell positions correlate with the cell index, which reintroduces structure at the cell frequency. A generator with a long period and good equidistribution in two dimensions is required; a linear congruential generator with a modulus of $2^{31}$ and two consecutive outputs paired as a point is not adequate, because the points lie on a small number of planes in two dimensions — the spectral test’s concern.

The third is treating the pixel’s area as the integration domain when the reconstruction filter is wider than one pixel. A filter with a radius of 2 pixels overlaps its neighbours’ areas, so each sample contributes to several pixels and the per-pixel sample count is not the number of samples in that pixel’s square. Stratifying over the pixel square then does not stratify the domain the filter actually integrates over, and the variance reduction is weaker than the count suggests. The correct domain for stratification is the filter’s support, and a renderer that uses a wide filter should stratify over that wider region, which costs more samples per pixel for the same quality per pixel but avoids the mismatch.

The fourth is adaptive sampling that terminates early on the basis of a variance estimate computed from the samples themselves. A stratified set has artificially low variance between neighbouring samples because they are arranged to cancel, so a variance estimate from a small stratified set underestimates the true error and the renderer stops too soon. Adaptive schemes need to estimate error from independent replicates, not from the stratified pattern, which usually means computing the estimate over a coarser grouping than the sample set itself.