# The Noise Falls Like 1/√n

*Task 5 of 5 · [Progressive Path Tracing: Noise Melting Into an Image](https://gpu.rocks/learn/path-tracing-c99efc67.md) · GPU.js Learn*

You have watched the noise go. Now measure it, because the rate it goes at is the
single most important number in rendering — and it is the same `1/√n` that
*Monte Carlo Methods* derives for a dart-throwing estimate of π. A path tracer is
that estimator with a much more interesting integrand.

Measuring it needs a reference, and the honest trick is to build one out of thin air:
run the accumulation **twice**, with two different slices of the random stream.
Neither run is the truth, but their errors are independent, so the RMS gap between them is
√2 times either one's own error. Halve the square and you have measured your own noise with
no ground truth anywhere in sight.

Plot that against `errs[0] / √n` on a log axis. If the two curves lie on top
of each other, the estimator is behaving exactly as theory says — and the shape of them is
the bad news in the theory: halving the noise costs *four times* the samples. That
single fact is why production renderers spend their effort on importance sampling and
denoisers instead of simply waiting longer.

## Goal

**Goal:** finish `sqDiff` — half the squared difference between
the two runs at this pixel — and build the `1/√n` reference curve to plot
against it.

## Requirements

- In `sqDiff`, return `0.5 * d * d` for `d = a - b` at this pixel
- Push `errs[0] / Math.sqrt(f + 1)` into `ideal` each pass
- Both curves go to `plot` on a log axis (already wired)

## Hint 1 — where the half comes from

Run A is off by `eₐ` and run B by `e_b`, independently.
The gap between them has variance `var(eₐ) + var(e_b)` = twice one run's.
So the square of the gap, halved, estimates the square of one run's error — and the
square root of the mean of that is the RMS error you want.

## Hint 2 — the kernel

```js
const d = a[this.thread.y][this.thread.x] - b[this.thread.y][this.thread.x];
return 0.5 * d * d;
```

Squaring in the kernel and square-rooting once in JavaScript is the usual split: the
per-pixel work is parallel, the final scalar is not.

## Hint 3 — the reference curve

Anchor it to the measurement you already have, so the two lines start together
and any drift between them is real:

```js
ideal.push(errs[0] / Math.sqrt(f + 1));
```

## Same idea elsewhere

√n is why every serious renderer stopped brute-forcing: importance sampling,
next-event estimation and multiple importance sampling all attack the *constant* in
front of the `1/√n` because the exponent itself cannot be moved; low-discrepancy
(quasi-Monte Carlo) sequences buy a slightly better exponent on smooth integrands; and
OptiX's and OIDN's neural denoisers give up unbiasedness altogether to get the picture
twenty times sooner. The measurement you just made — two independent runs, no reference
image — is also how those denoisers are evaluated in practice.

## Starter code

```js
const gpu = new GPU({ mode });

// The whole random-number generator. It is a pure function: same state in,
// same number out, on every backend and every run.
gpu.addFunction(function nextRandom(state) {
  let p = state * 78.233 + 0.7213;
  p = p - Math.floor(p);
  p = p * (p + 61.7);
  return p - Math.floor(p);
});

// Nearest hit of a unit-direction ray against a sphere, or -1 for a miss.
gpu.addFunction(function sphereT(ox, oy, oz, dx, dy, dz, cx, cy, cz, r) {
  const ex = ox - cx;
  const ey = oy - cy;
  const ez = oz - cz;
  const b = ex * dx + ey * dy + ez * dz;
  const c = ex * ex + ey * ey + ez * ez - r * r;
  const disc = b * b - c;
  if (disc <= 0) {
    return -1;
  }
  const t = -b - Math.sqrt(disc);
  if (t <= 0.001) {
    return -1;
  }
  return t;
});

// The floor: one division, no quadratic.
gpu.addFunction(function planeT(oy, dy, planeY) {
  if (dy > -0.0001) {
    return -1;
  }
  const t = (planeY - oy) / dy;
  if (t <= 0.001) {
    return -1;
  }
  return t;
});

// The only light in the scene is the sky itself.
gpu.addFunction(function skyLight(dy) {
  if (dy < 0) {
    return 0.3;
  }
  return 0.3 + 0.7 * dy;
});

const trace = gpu.createKernel(function (frame, spp) {
  // This thread's own stream: start from the frame, stir in the column,
  // then the row. Every stir goes through nextRandom, so two threads one
  // pixel apart come out completely unrelated.
  let seed = 0.1237 + frame * 0.0173;
  seed = nextRandom(seed + this.thread.x * 0.0713);
  seed = nextRandom(seed + this.thread.y * 0.0917);

  let total = 0;
  for (let s = 0; s < spp; s++) {
    // One camera ray, jittered inside this pixel — so the edges antialias
    // themselves for free as the frames pile up.
    seed = nextRandom(seed);
    const u = ((this.thread.x + seed) / 64) * 2 - 1;
    seed = nextRandom(seed);
    const v = ((this.thread.y + seed) / 64) * 2 - 1;
    let dx = u;
    let dy = v;
    let dz = -1.7;
    let len = Math.sqrt(dx * dx + dy * dy + dz * dz);
    dx = dx / len;
    dy = dy / len;
    dz = dz / len;
    let ox = 0;
    let oy = 0.3;
    let oz = 1.5;
    let through = 1;

    // Up to three segments: camera → surface → surface → sky.
    for (let d = 0; d < 3; d++) {
      const tA = sphereT(ox, oy, oz, dx, dy, dz, -0.45, 0, -0.4, 0.5);
      const tB = sphereT(ox, oy, oz, dx, dy, dz, 0.5, -0.2, -0.1, 0.3);
      const tP = planeT(oy, dy, -0.5);
      let t = 999;
      let which = -1;
      if (tA > 0 && tA < t) {
        t = tA;
        which = 0;
      }
      if (tB > 0 && tB < t) {
        t = tB;
        which = 1;
      }
      if (tP > 0 && tP < t) {
        t = tP;
        which = 2;
      }
      if (which < 0) {
        // Nothing in the way — the path ends in the sky, which is the light.
        total = total + through * skyLight(dy);
        break;
      }

      const hx = ox + t * dx;
      const hy = oy + t * dy;
      const hz = oz + t * dz;
      let nx = 0;
      let ny = 1;
      let nz = 0;
      let albedo = 0.55;
      if (which === 0) {
        nx = (hx + 0.45) / 0.5;
        ny = hy / 0.5;
        nz = (hz + 0.4) / 0.5;
        albedo = 0.85;
      }
      if (which === 1) {
        nx = (hx - 0.5) / 0.3;
        ny = (hy + 0.2) / 0.3;
        nz = (hz + 0.1) / 0.3;
        albedo = 0.3;
      }

      // A uniformly random point on the unit sphere, out of this thread's
      // stream: z is uniform on [-1, 1] and phi uniform around the axis.
      seed = nextRandom(seed);
      const z = 2 * seed - 1;
      seed = nextRandom(seed);
      const phi = 6.2831853 * seed;
      const rr = Math.sqrt(Math.max(1 - z * z, 0));
      const sx = rr * Math.cos(phi);
      const sy = rr * Math.sin(phi);
      const sz = z;

      through = through * albedo;
      dx = nx + sx;
      dy = ny + sy;
      dz = nz + sz;
      len = Math.sqrt(dx * dx + dy * dy + dz * dz);
      dx = dx / len;
      dy = dy / len;
      dz = dz / len;

      // Step off the surface so the next ray does not hit it at t = 0.
      ox = hx + 0.001 * nx;
      oy = hy + 0.001 * ny;
      oz = hz + 0.001 * nz;
    }
  }
  return total / spp;
}, { output: [64, 64], loopMaxIterations: 8, pipeline: true });

// A buffer of zeros to start the average from — the same trick Pipelines &
// Textures uses to get its seed onto the card as a texture.
const black = gpu.createKernel(function () {
  return 0;
}, { output: [64, 64], pipeline: true });

const accumulate = gpu.createKernel(function (previous, sample, n) {
  const old = previous[this.thread.y][this.thread.x];
  const now = sample[this.thread.y][this.thread.x];
  return (old * n + now) / (n + 1);
}, { output: [64, 64], pipeline: true, immutable: true });

const sqDiff = gpu.createKernel(function (a, b) {
  // TODO: half the squared difference between the two runs at this pixel.
  // Half, because two independent runs disagree by √2 times either one's
  // own error.
  return 0;
}, { output: [64, 64] });

// Radiance in, pixels out. The square root is a gamma curve: it is what
// makes a 0.2 and a 0.4 look as different as they are.
const paint = gpu.createKernel(function (image) {
  const g = Math.sqrt(Math.max(image[this.thread.y][this.thread.x], 0));
  this.color(g, g, g, 1);
}, { output: [64, 64], graphical: true });

const FRAMES = 24;
let accA = await black();
let accB = await black();
const errs = [];
const ideal = [];

for (let f = 0; f < FRAMES; f++) {
  // The same accumulation twice, on two different slices of the stream.
  accA = await accumulate(accA, await trace(f, 1), f);
  accB = await accumulate(accB, await trace(f + 500, 1), f);

  const d = await sqDiff(accA, accB);
  let sum = 0;
  for (let y = 0; y < 64; y++) {
    for (let x = 0; x < 64; x++) sum += d[y][x];
  }
  errs.push(Math.sqrt(sum / 4096));
  // TODO: the theory line — errs[0] / Math.sqrt(f + 1)
  ideal.push(errs[0]);
}

await paint(accA);
render(paint.canvas);
console.log('RMS noise per frame:', errs);
plot({ measured: errs, 'errs[0] / √n': ideal }, {
  title: 'RMS noise vs frames accumulated',
  log: true,
});
```

---

Interactive version: https://gpu.rocks/learn/path-tracing-c99efc67/5

[Previous task](https://gpu.rocks/learn/path-tracing-c99efc67/4.md)