Part IV · Path tracing · Chapter 12
Dice and averages
The light reaching one pixel comes from every direction, bounced off everything in the scene. Nobody can add all of that up exactly. A path tracer guesses instead: it follows one random path, then another, and averages what they bring back. This chapter is about why an average of random guesses gets the right answer, how fast it gets there, how the renderer knows how far it still is, and what happens when the dice are loaded.
1An answer from random guesses
Throw darts at a square one unit on a side, anywhere at random. A quarter circle of radius one covers π/4 of the square, so that is the chance a dart lands inside it. Score 4 for a dart inside and 0 for one outside: the average score is then π. Twenty darts give 3.4. That is a Monte Carlo estimate: an answer worked out as the average of random tries.
Throw more (the button) and the estimate wanders less and settles towards π. The plot shows it after every few darts, on a log scale (each label along the bottom four or sixteen times the one before, so the first darts get as much room as the last thousands), with an error bar either side: the faint lines. The bar says how far from the true answer the estimate typically is, and it is worked out from the darts themselves, without knowing π. At the very start it cannot be trusted: the first twelve darts here all landed inside, they agree perfectly, and the bar says ±0 while the estimate, 4, is far from π. A bar worked out from the tries is only as good as the tries, and it settles after a few dozen.
After 81,920 darts the estimate is 3.1440 ± 0.0057. The true value is 0.4 error bars away. An error bar is not a wall: about one estimate in three lands more than one bar away, one in twenty more than two. (Those two shares come from the bell-shaped curve that averages of many tries follow, a result to take on trust here; the figure lets you check it, below.) The readout counts the bars each time; Start again throws new darts. And the bar only measures how much the tries disagree with each other. Dice loaded the same way every time (section 4) make tries that agree on a wrong answer, and the bar shrinks round it just the same: it can see noise, never bias.
A pixel is the same game. A path is a ray from the eye that bounces off surfaces until it finds light (chapter 13 follows one). Each sample sends one random path through the pixel and brings back the light it found, the pixel's value is the average of what its samples brought back, and the noise in a young picture is each pixel's estimate still wandering.
core/measure.js).2Four times the work, half the error
How fast does it settle? Here are 240 separate estimates of π for each number of darts N, each dot one estimate. With 4 darts they scatter over the whole range, and sit only on 0, 1, 2, 3 or 4 (0 is rare; each dart adds 0 or 1 to the estimate); with 16, half as far; with 64, half again. The solid bar beside each column is how far they spread, measured from the dots (their standard deviation); the dashed one is the prediction, the spread of one dart divided by √N. Four times the darts, half the spread: that rule holds for every Monte Carlo estimate, and it is why a noisy picture needs four times the samples to look twice as clean. The measured spread is itself worked out from 240 dots, so it wanders too, typically by about 5% and now and then by twice that (New darts shows it): at N = 256 it may read 0.093 against the predicted 0.103. With 960 dots a column it wanders half as much: the same √N rule, applied to the spread itself. Point at a column to read its numbers.
The mathsWhy √N
Each dart scores X = 4 with probability p = π/4, else 0, so its average is 4p = π. How far it lands from that average, squared and averaged, is its variance: a score of 4 is 4 − π away, a score of 0 is π away, and
σ² = p(4 − π)² + (1 − p) π² = 16 p(1 − p)
(put π = 4p in: 16p(1 − p)² + 16p²(1 − p) = 16p(1 − p)). Its square root, the standard deviation, is σ = 4√(p(1 − p)) = 1.642. Now the √N. For tries that do not depend on each other, variances add. For two tries, the squared distance of their sum from its average is the first's squared distance, plus the second's, plus twice the first's distance times the second's; and when the tries do not depend on each other, that last product is as often positive as negative, and averages to 0. So the sum of N darts has variance Nσ², and spreads by √N σ. The average is that sum divided by N, and so is its spread:
σmean = √N σN = σ√N
- σ
- how far one try typically lands from the answer
- N
- how many tries are averaged
The error bar in section 1 is this, with σ worked out from the darts so far: σmean ≈ s/√N, where s is the spread of the scores about their own average. It is called the standard error. Nothing in it depends on how many random numbers each try is made of. A dart is two; a path through a scene is dozens (a direction at every bounce, a point of a lamp, a point of the lens). Spacing tries evenly instead, ten values for each of a dozen numbers, would take 10¹² tries; random tries settle as 1/√N however many numbers each one draws. That is why Monte Carlo is the way to add up light arriving from every direction.
3An error bar for every pixel
After each sample, a kernel (a program the GPU runs once for every pixel) called
resolve adds the sample to each pixel's average. It keeps no sum. It moves the average towards the new
sample by a share w = 1/cnt of the way, where cnt is how many samples the pixel
now holds: the first sample is the average (w = 1), the second moves it halfway, the third a third of
the way, and the result is always the plain average of all of them. With values 2, 4 and 9: 2, then 3, then
3 + (9 − 3)/3 = 5, which is (2 + 4 + 9)/3.
Beside the average it keeps two more numbers per pixel, so that the renderer knows each pixel's error bar too: the cursor readout in ShadyTracer shows it, and the test bench reads pictures by it. V is how much the pixel's samples scatter about their average (the average squared distance from it; for 2, 4 and 9, 8.67). S is the sum of the squares of the weights each sample has in the average. In a true average of n samples each weighs 1/n, so S = n · (1/n)² = 1/n, and the error bar the code works out, √(S · V / (1 − S)), becomes √(V / (n − 1)): the darts' s/√N exactly, with s worked out the way statisticians do, dividing by n − 1 rather than n (it corrects for measuring the scatter about the samples' own average, which sits closer to them than the true answer does). For 2, 4 and 9 that is √(8.67/2) = 2.08.
Both are updated a sample at a time, with d the new sample's distance from the old average. For 2, 4 and 9: after 2, V = 0 and S = 1. After 4 (w = ½, d = 4 − 2 = 2): V = (1 − ½)(0 + ½ × 2²) = 1, S = (1 − ½)² × 1 + ½² = ½. After 9 (w = ⅓, d = 9 − 3 = 6): V = ⅔ × (1 + ⅓ × 6²) = 8.67, S = (⅔)² × ½ + (⅓)² = ⅓. The same 8.67 and ⅓ as from all three at once. S is kept, rather than worked out from n, because the average can also forget old samples (section 6), and then the weights are not equal.
The host's side, and the darts figure above, use the same formula from a JavaScript copy of the update: .
4Loaded dice
All of this assumes the random numbers are fair: every value between 0 and 1 equally likely, every small stretch
of [0, 1) hit as often as its width says. A GPU has no dice; it has a pseudo-random number
generator, a formula that turns one number into the next in a way that looks random. ShadyTracer's is PCG4D
(section 5). Its old one was fract(f32(state) / 671040), and the comment above the new one says why it
went.
Here is the share of all possible states that land in each of 24 small windows of [0, 1), worked out exactly, for the old formula (left bars) and for PCG4D's output (right). The line at 1 is the fair share. In windows 0.001 wide, the old one is within about 4% of fair.
In windows of 0.0001 it is not: some get 0.95 of their share and some 1.83. The reason
starts with the state, the number the generator keeps and steps on: a whole number from 0 to about 4.3
billion (2³², 32 binary digits). The old formula divided it by 671,040 and kept what is after the point
(fract), a number in [0, 1): from 0 up to but not including 1. But first it turned the state into the
GPU's usual kind of number, an f32, which keeps 24 significant binary digits, the way a
calculator that shows seven digits turns 12,345,678 into 12,345,680. Every state below 2²⁴ (about 16.8 million)
fits whole. Larger ones are rounded to the nearest multiple of 2, 4, 8 …, and every state from 2²⁹ up, which is
seven-eighths of all of them, becomes a multiple of 64 or more: from 2²⁹ a state has 30 binary digits, six more
than f32 keeps, and dropping six binary digits leaves multiples of 2⁶ = 64. (From 2³⁰, multiples of 128; from
2³¹, of 256.) Divided by 671,040, which is 64 × 10,485, a multiple of 64 lands on a multiple of 1/10,485. A
multiple of 256 would seem to land only on a coarser grid, 256/671,040, about 1/2600, but fract keeps
only what is after the point, and as the numbers wrap round past each whole number, multiples of 256 fall on every
point of the 1/10,485 grid in turn. So seven-eighths of all states land on that grid, and the rest, below 2²⁹, on
finer ones. A window 0.0001 wide holds one of those points or two: on average 1.05. One point carries about
1/1.05 of seven-eighths of the window's fair share, plus a little of the finer eighth, which is why the bars read
0.95 and 1.83 rather than exactly 1 and 2.
In windows of 0.00001, narrower than the grid, some get almost nothing (0.01 of their share), some a little (the low bars: windows between the coarse grid's points, reached only by finer grids of the states below 2²⁹), and some nearly nine times (a window that holds a coarse point). PCG4D's output, the top 24 bits of its state divided by 2²⁴, is a grid 1/16,777,216 apart with every point equally likely. A window of 0.00001 holds 167 or 168 of those points, so its share is within half a percent of fair (the readout's 0.995×), and the wider the window, the closer.
Why would a window matter? A small, bright sun is a small window: of all the random numbers a bounce can draw, the few that send it into the sun. Which numbers those are depends on the surface the bounce leaves (the numbers are turned into a direction about that surface's own normal), so, in the comment's words, "each pixel sees it at a different place on the grid": one pixel's window holds an extra grid point, its neighbour's misses one. No number of samples averages that away: each pixel settled, confidently, on its own wrong answer, and its error bar shrank round it.
5How ShadyTracer does it
The generator, and the story of the one before it:
The comment's numbers, corrected. Until version 0.67.3 the comment said the old output's share of a 0.0001 window "ranged from 0.08x to 3.2x", at 0.001 "from 0.96x to 1.29x", and called the grid "about 1/2600", the multiples of 256 without the wrap of section 4. Counting every one of the 2³² states through the formula it gives (this tutorial's test does, in node, and once checked it against 63 million random states) gives 0.95× to 1.83×, 0.96× to 1.04×, and a grid of 1/10,485, as in the figure; this chapter found the difference, and the comment now gives these numbers and says what it used to. The old code is not in the repository's history, which starts at version 0.46.0, so whether it did more than the formula says cannot be checked. Either way the comment's point stands: a bias in small windows that no number of samples removes.
Each pass of a sample hands the kernels a fresh random seed, and each ray mixes it with its pixel and its bounce (), so no two pixels, bounces or samples start their dice from the same place.
6Tricks and trade-offs
Each is stamped with how it is known; point at a stamp to see what it means.
An error bar from a count. With aiming at the light (next-event
estimation, chapter 14) switched off, a sun 5° across is found only when a bounce happens to go into it:
about 50,000 times in 37 million paths (1188 × 969 pixels, 32 samples each, on the user's graphics card, an
NVIDIA Pascal; docs/TESTING.md, the section "BSDF parity"). The document puts the standard error near
0.45%; 1/√50,000 is 0.45%. Light that arrives through rare events is as noisy as the count of those events: the
sun's light comes only through the hits, a count of rare hits wanders by about its own square root (50,000 ±
224), and so its share of the picture wanders by 1/√(hits), whatever the number of paths around them.
A true mean, not a window. The average can forget old samples, which keeps a moving picture fresh: cnt stops growing at the window's size plus one, so each new sample keeps a fixed weight and the old ones fade. It used to, by default: a window of 10, each new sample weighing 1/11, which holds as much as a true average of about 21 samples (the code's comment: S settles near 1/(2N + 1)). A still picture then never converged. Through the usual tone curve (chapter 7) nobody saw it; in the false-colour view, which paints each band of brightness its own colour, every pixel kept flickering between bands. Moving the camera restarts the picture anyway, so the default is now the true mean ().
No sum of squares. The textbook variance, the mean of the squares less the square of the mean, subtracts two nearly equal numbers wherever the noise is small beside the light, and in 32-bit numbers that leaves nothing. With a calculator that keeps seven digits: samples 1000.1 and 1000.3 have squares 1,000,200 and 1,000,600 (their last digits gone), whose mean is 1,000,400; their mean, 1000.2, squared is 1,000,400 too. The difference is 0, and the true variance, 0.01, is lost. The update in section 3 never forms either sum.
7What went wrong before
- Loaded dice: the old generator, section 4.
- A light that leaked into the next sample. Each pixel has a slot where a sample's light is
added up before
resolvetakes it.resolveused to stop at once for a pixel "whose camera path then did not count" in a sample (the comment's words), before clearing its slot, so light a shadow ray had added there was counted in the next sample instead. It now takes and clears every slot first (the comments inprog_resolve). - Large pictures that would not run. The GPU runs a kernel in groups of threads, here 64 pixels
a group, and allows at most 65,535 groups along one direction.
resolveonce numbered pixels in one long row of groups, so it stopped at 64 × 65,535 = 4,194,240 pixels: a 2048 × 2048 picture was over by one group, and a 4K one far over. It is now one row of the picture per row of groups. - A window nobody needed: the averaging window of 10, section 6.
8Try it
- In the darts figure, press Throw 4 times as many once (80 darts, past the few dozen a bar needs) and then Start again a dozen times, throwing once more each time, and count how often π is more than one error bar away. About a third?
- In the dice figure, choose windows of 0.0001 and move where. Do the old bars' heights ever change pattern?
- In ShadyTracer, point at the picture with the cursor readout on (Render stats) and watch a pixel's ± shrink: how many samples does it take to halve?