Path tracer scene, light and sky buffers

/** Every mesh item in world space (matrix applied, drawRange and index honoured), cut by the planes (keep n·p ≤ d). */
export function bakeTriangles(items, planes = [], ids = null) {
  planes = planes || [];
  const P = new Grow(), N = new Grow(), U = new Grow(), look = [], obj = [];
  // a clipped polygon has at most 3 + planes vertices: fixed scratch, reused
  const cap = 3 + planes.length * 2, pp = new Float64Array(cap * 3), pn = new Float64Array(cap * 3), pu = new Float64Array(cap * 2);
  const qp = new Float64Array(cap * 3), qn = new Float64Array(cap * 3), qu = new Float64Array(cap * 2);
  for (let ii = 0; ii < items.length; ii++) {
    const it = items[ii];
    if (!it || it.lines || !it.geometry) continue;
    const oid = ids ? ids[ii] : ii;
    const g = it.geometry, pa = g.getAttribute('position'), na = g.getAttribute('normal'), ua = g.getAttribute('uv');
    if (!pa) continue;
    const m = it.matrix || null;
    // normal matrix = inverse transpose of the upper 3 × 3 (column-major m)
    let nm = null;
    if (m) {
      const a = m[0], b = m[4], c = m[8], d = m[1], e = m[5], f = m[9], gg = m[2], h = m[6], i = m[10];
      const A = e * i - f * h, B = -(d * i - f * gg), C = d * h - e * gg, det = a * A + b * B + c * C || 1;
      // inverse transpose, row-major rows
      nm = [A / det, B / det, C / det, -(b * i - c * h) / det, (a * i - c * gg) / det, -(a * h - b * gg) / det, (b * f - c * e) / det, -(a * f - c * d) / det, (a * e - b * d) / det];
    }
    const idx = g.index ? g.index.array : null, start = (g.drawRange && g.drawRange.start) || 0;
    const total = idx ? idx.length : pa.count;
    const rc = g.drawRange && Number.isFinite(g.drawRange.count) ? g.drawRange.count : Infinity;
    const end = Math.min(total, start + rc);
    const key = it.look || 'default';
    for (let t = start; t + 2 < end; t += 3) {
      let nv = 3;
      for (let k = 0; k < 3; k++) {
        const vi = idx ? idx[t + k] : t + k;
        const x = pa.getX(vi), y = pa.getY(vi), z = pa.getZ(vi);
        if (m) { pp[k * 3] = m[0] * x + m[4] * y + m[8] * z + m[12]; pp[k * 3 + 1] = m[1] * x + m[5] * y + m[9] * z + m[13]; pp[k * 3 + 2] = m[2] * x + m[6] * y + m[10] * z + m[14]; }
        else { pp[k * 3] = x; pp[k * 3 + 1] = y; pp[k * 3 + 2] = z; }
        let nx = na ? na.getX(vi) : 0, ny = na ? na.getY(vi) : 0, nz = na ? na.getZ(vi) : 0;
        if (nm) { const X = nm[0] * nx + nm[1] * ny + nm[2] * nz, Y = nm[3] * nx + nm[4] * ny + nm[5] * nz, Z = nm[6] * nx + nm[7] * ny + nm[8] * nz; nx = X; ny = Y; nz = Z; }
        const l = Math.hypot(nx, ny, nz);
        pn[k * 3] = l ? nx / l : 0; pn[k * 3 + 1] = l ? ny / l : 0; pn[k * 3 + 2] = l ? nz / l : 0;
        pu[k * 2] = ua ? ua.getX(vi) : 0; pu[k * 2 + 1] = ua ? ua.getY(vi) : 0;
      }
      // section cuts (Sutherland–Hodgman against each plane)
      let src = [pp, pn, pu], dst = [qp, qn, qu];
      for (const pl of planes) {
        if (nv < 3) break;
        const [sp, sn, su] = src, [dp, dn, du] = dst;
        let out = 0;
        for (let i = 0; i < nv; i++) {
          const j = (i + 1) % nv;
          const di = pl.n[0] * sp[i * 3] + pl.n[1] * sp[i * 3 + 1] + pl.n[2] * sp[i * 3 + 2] - pl.d;
          const dj = pl.n[0] * sp[j * 3] + pl.n[1] * sp[j * 3 + 1] + pl.n[2] * sp[j * 3 + 2] - pl.d;
          if (di <= 0) { for (let k = 0; k < 3; k++) { dp[out * 3 + k] = sp[i * 3 + k]; dn[out * 3 + k] = sn[i * 3 + k]; } du[out * 2] = su[i * 2]; du[out * 2 + 1] = su[i * 2 + 1]; out++; }
          if ((di <= 0) !== (dj <= 0)) {
            const f = di / (di - dj);
            for (let k = 0; k < 3; k++) { dp[out * 3 + k] = sp[i * 3 + k] + (sp[j * 3 + k] - sp[i * 3 + k]) * f; dn[out * 3 + k] = sn[i * 3 + k] + (sn[j * 3 + k] - sn[i * 3 + k]) * f; }
            du[out * 2] = su[i * 2] + (su[j * 2] - su[i * 2]) * f; du[out * 2 + 1] = su[i * 2 + 1] + (su[j * 2 + 1] - su[i * 2 + 1]) * f;
            out++;
          }
        }
        nv = out;
        const tmp = src; src = dst; dst = tmp;
      }
      if (nv < 3) continue;
      const [sp, sn, su] = src;
      for (let i = 1; i + 1 < nv; i++) {
        for (const k of [0, i, i + 1]) { P.push3(sp[k * 3], sp[k * 3 + 1], sp[k * 3 + 2]); N.push3(sn[k * 3], sn[k * 3 + 1], sn[k * 3 + 2]); U.push2(su[k * 2], su[k * 2 + 1]); }
        look.push(key); obj.push(oid);
      }
    }
  }
  return { n: look.length, pos: P.done(), nor: N.done(), uv: U.done(), look, obj };
}

/**
 * The derived words of a record (§ 6.4 vec4 14.zw, 15): the base interface's relative IOR after OpenPBR's coat rule (η_s =
 * lerp(n_b, n_b/n_c or n_c/n_b when n_c > n_b, C)) and its specular_weight modulation; the coat's hemispherical Fresnel
 * average E_F (for the darkening); the medium's extinction μt = −ln(transmission_color) / transmission_depth (0 = none: the
 * colour then tints at the surface); Cauchy's B for dispersion (n(λ) = A + B/λ², λ in µm, n(λ_d) = η, V_d = abbe / scale).
 */
export function finishRecord(r) {
  const nb = Math.max(1e-4, r[9]), nc = Math.max(1e-4, r[29]), C = r[27];
  const ebc = nc > nb ? nc / nb : nb / nc;
  const etaS = nb + (ebc - nb) * C;
  r[58] = modulatedEta(etaS, r[7]);
  let ef = 0;
  for (let i = 0; i < 32; i++) { const mu = (i + 0.5) / 32; ef += 2 * frDielectric(mu, nc) * mu / 32; }
  r[59] = ef;
  const depth = r[51], T = r[11];
  for (let k = 0; k < 3; k++) r[60 + k] = depth > 0 && T > 0 ? -Math.log(Math.max(1e-6, Math.min(1, r[12 + k]))) / depth : 0;
  const scale = r[53], abbe = r[54];
  r[63] = scale > 0 && T > 0 ? (r[58] - 1) / ((abbe / scale) * (1 / 0.4861 ** 2 - 1 / 0.6563 ** 2)) : 0;
  return r;
}

/** Preetham's sky for a sun direction (toward the sun) and turbidity (2 = very clear … 10 = hazy). */
export function preetham(sunDir, up, turbidity = 3) {
  const T = Math.min(10, Math.max(1.7, turbidity)), s = vnorm(sunDir), u = vnorm(up);
  const thetaS = Math.acos(Math.min(1, Math.max(0, dot(s, u))));   // the sky model wants the sun at or above the horizon
  const co = (tab) => tab.map(([a, b]) => a * T + b);
  const cY = co(PERZ.Y), cx = co(PERZ.x), cy = co(PERZ.y);
  const chi = (4 / 9 - T / 120) * (Math.PI - 2 * thetaS);
  const Yz = Math.max(0, (4.0453 * T - 4.9710) * Math.tan(chi) - 0.2155 * T + 2.4192);
  const ts = [thetaS ** 3, thetaS ** 2, thetaS, 1], tv = [T * T, T, 1];
  const mat = (M) => tv.reduce((acc, t, r) => acc + t * M[r].reduce((a, m, k) => a + m * ts[k], 0), 0);
  const xz = mat(XZ), yz = mat(YZ);
  const nY = perez(cY, 0, thetaS), nx = perez(cx, 0, thetaS), ny = perez(cy, 0, thetaS);
  return {
    zenithY: Yz, turbidity: T, thetaS,
    radiance(dir) {
      const d = vnorm(dir), cz = dot(d, u);
      const theta = Math.acos(Math.min(1, Math.max(0, cz)));   // below the horizon: the horizon's colour (the ground covers it)
      const gamma = Math.acos(Math.min(1, Math.max(-1, dot(d, s))));
      const Y = (Yz * perez(cY, theta, gamma)) / nY, x = (xz * perez(cx, theta, gamma)) / nx, y = (yz * perez(cy, theta, gamma)) / ny;
      const X = (x / Math.max(1e-6, y)) * Y, Z = ((1 - x - y) / Math.max(1e-6, y)) * Y;
      return mul3(XYZ_TO_REC709, [X, Y, Z]).map((v) => Math.max(0, v));
    },
  };
}

export function profileGrid(ies, nTheta = 64, nPhi = 32) {
  const g = new Float32Array(nTheta * nPhi), k = (4 * Math.PI) / ies.flux;
  for (let a = 0; a < nTheta; a++) for (let b = 0; b < nPhi; b++) g[a * nPhi + b] = candelaAt(ies, ((a + 0.5) * 180) / nTheta, ((b + 0.5) * 360) / nPhi) * k;
  return g;
}