Path tracer kernels

function render({ spp = 1, s0 = 0, accum = null, aux = null, x0 = 0, y0 = 0, x1 = w, y1 = h, active = null } = {}) {   // active (RN.7, pt-adapt.js): 0 = this pixel gets no sample
  accum = accum || new Float64Array(w * h * 4);
  if (STRIDE) aux = aux || new Float64Array(w * h * STRIDE * 4);
  for (let s = s0; s < s0 + spp; s++) {
    for (let py = y0; py < y1; py++) for (let px = x0; px < x1; px++) {
      if (active && !active[py * w + px]) continue;
      const c = radiance(px, py, s), o = (py * w + px) * 4;
      const fine = Number.isFinite(c[0]) && Number.isFinite(c[1]) && Number.isFinite(c[2]);
      if (fine) { accum[o] += Math.max(0, c[0]); accum[o + 1] += Math.max(0, c[1]); accum[o + 2] += Math.max(0, c[2]); }
      accum[o + 3] += 1;
      if (STRIDE && fine) {
        const b = (py * w + px) * STRIDE * 4;
        for (let i = 0; i < STRIDE * 4; i++) {
          if (lit && i >= EL.IDS * 4 && i < EL.IDS * 4 + 4) { if (s === 0) aux[b + i] = AX[i]; continue; }
          aux[b + i] += Number.isFinite(AX[i]) ? AX[i] : 0;
        }
      }
    }
  }
  return STRIDE ? Object.assign(accum, { aux }) : accum;
}

export function frDielectric(cosI, eta) {
  let c = Math.min(1, Math.max(-1, cosI));
  if (c < 0) { eta = 1 / eta; c = -c; }
  const s2 = 1 - c * c, t2 = s2 / (eta * eta);
  if (t2 >= 1) return 1;
  const ct = Math.sqrt(1 - t2);
  const rpa = (eta * c - ct) / (eta * c + ct), rpe = (c - eta * ct) / (c + eta * ct);
  return (rpa * rpa + rpe * rpe) / 2;
}

/**
 * The sample for pixel (px, py), sample index s, dimension d (RENDER-SPEC § 4). d < BN_LIMIT: an additive recurrence (R2 for
 * the 2-D dimensions, the golden ratio for the 1-D ones: alphaFor) rotated by blue noise, in 32-bit fixed point; else a PCG hash.
 */
export function sample(mask, px, py, s, d) {
  if (d < DIMS.BN_LIMIT) {
    const q = d >> 2, ox = (13 * q) & 63, oy = (29 * q) & 63;
    const bn = mask[((((py + oy) & 63) * 64 + ((px + ox) & 63)) << 2) + (d & 3)];
    const off = ((bn << 24) | (hash4(px, py, d, 0x5bd1e995) & 0xffffff)) >>> 0;
    return u32ToUnit((Math.imul(s >>> 0, alphaFor(d)) + off) >>> 0);
  }
  return u32ToUnit(hash4(px, py, s, d));
}

export function makeAdaptive(w, h, { noise = 0.02, min = 16 } = {}) {
  const n = w * h, active = new Uint8Array(n).fill(1), S = new Float64Array(n * 4);   // S1, S2, Y before, N before
  const j0 = minPasses(min);
  let count = n;
  return {
    active, noise, min, w, h,
    get count() { return count; },
    update(accum) {
      for (let i = 0; i < n; i++) {
        if (!active[i]) continue;
        const o = i * 4, N = accum[o + 3], Y = 0.2722 * accum[o] + 0.6741 * accum[o + 1] + 0.0537 * accum[o + 2], dn = N - S[o + 3];
        if (!(dn > 0)) continue;
        const m = (Y - S[o + 2]) / dn;
        S[o] += m; S[o + 1] += m * m; S[o + 2] = Y; S[o + 3] = N;
        const J = Math.round(N / ADAPT.PASS);
        if (J >= j0) {
          const mu = S[o] / J, v = Math.max(0, (S[o + 1] - (S[o] * S[o]) / J) / (J - 1)), se = Math.sqrt(v / J);
          if (se <= noise * (mu + ADAPT.FLOOR)) { active[i] = 0; count--; }
        }
      }
      return count;
    },
  };
}