Ray casting for rendering and the render bench

/**
 * Quantize the child boxes of one node into its words. nodeBox: 6 numbers (f32-exact). kidBox(i, out6) fills child i's box
 * and returns 1 inner / 2 leaf / 0 empty.
 */
function quantize(words, f, base, nodeBox, nb0, kidKind, kidBoxes) {
  // origin = node min (already f32, rounded down); scale per axis: smallest power of two with 255 · s ≥ extent
  const ox = nodeBox[nb0], oy = nodeBox[nb0 + 1], oz = nodeBox[nb0 + 2];
  let ebits = 0, l0 = 0, l1 = 0, l2 = 0, h0 = 0, h1 = 0, h2 = 0, innerMask = 0;
  const S = [0, 0, 0], O = [ox, oy, oz];
  for (let ax = 0; ax < 3; ax++) {
    const ext = nodeBox[nb0 + 3 + ax] - nodeBox[nb0 + ax];
    let e = ext > 0 ? Math.ceil(Math.log2(ext / 255)) : -126;
    if (e < -126) e = -126; if (e > 127) e = 127;
    while (e < 127 && 255 * 2 ** e < ext) e++;   // log2 can land one short through rounding
    S[ax] = 2 ** e;
    ebits |= ((e + 127) & 255) << (8 * ax);
  }
  f[base] = ox; f[base + 1] = oy; f[base + 2] = oz;
  for (let i = 0; i < 4; i++) {
    const kind = kidKind[i];
    if (!kind) continue;
    if (kind === 1) innerMask |= 1 << i;
    const kb = i * 6;
    for (let ax = 0; ax < 3; ax++) {
      const o = O[ax], s = S[ax];
      let ql = Math.floor((kidBoxes[kb + ax] - o) / s), qh = Math.ceil((kidBoxes[kb + 3 + ax] - o) / s);
      ql = ql < 0 ? 0 : ql > 255 ? 255 : ql; qh = qh < 0 ? 0 : qh > 255 ? 255 : qh;
      // conservative in f32 (o + q·s with one rounding, as the kernels compute it)
      while (ql > 0 && f32(o + ql * s) > kidBoxes[kb + ax]) ql--;
      while (qh < 255 && f32(o + qh * s) < kidBoxes[kb + 3 + ax]) qh++;
      const sh = 8 * i;
      if (ax === 0) { l0 |= ql << sh; h0 |= qh << sh; } else if (ax === 1) { l1 |= ql << sh; h1 |= qh << sh; } else { l2 |= ql << sh; h2 |= qh << sh; }
    }
  }
  words[base + 3] = (ebits | (innerMask << 24)) >>> 0;
  words[base + 8] = l0 >>> 0; words[base + 9] = l1 >>> 0; words[base + 10] = l2 >>> 0;
  words[base + 11] = h0 >>> 0; words[base + 12] = h1 >>> 0; words[base + 13] = h2 >>> 0;
}

/** stk: a stack from newStack(); a traversal nested in another (instances: TLAS → BLAS) passes its own. */
export function traverseWide(words, f, root, ray, tmax, leafHit, any = false, st = st0) {
  const { ox, oy, oz, dx, dy, dz } = ray;
  const ix = SAFE(dx), iy = SAFE(dy), iz = SAFE(dz);
  let best = tmax;
  let sp = 0;
  const S = st.s, T = st.t;
  S[sp] = root; T[sp] = 0; sp++;
  const ct = [0, 0, 0, 0], ck = [0, 0, 0, 0];
  while (sp) {
    sp--;
    if (T[sp] > best) continue;
    const node = S[sp], base = node * WORDS;
    const e = words[base + 3], mask = e >>> 24, counts = words[base + 6];
    const sx = 2 ** ((e & 255) - 127), sy = 2 ** (((e >>> 8) & 255) - 127), sz = 2 ** (((e >>> 16) & 255) - 127);
    const nx = f[base], ny = f[base + 1], nz = f[base + 2];
    let ci = words[base + 4], slot = words[base + 5], nk = 0;
    for (let i = 0; i < 4; i++) {
      const sh = 8 * i;
      const inner = (mask >> i) & 1;
      const c = (counts >>> sh) & 255;
      if (!inner && !c) continue;
      const x0 = f32(nx + ((words[base + 8] >>> sh) & 255) * sx), x1 = f32(nx + ((words[base + 11] >>> sh) & 255) * sx);
      const y0 = f32(ny + ((words[base + 9] >>> sh) & 255) * sy), y1 = f32(ny + ((words[base + 12] >>> sh) & 255) * sy);
      const z0 = f32(nz + ((words[base + 10] >>> sh) & 255) * sz), z1 = f32(nz + ((words[base + 13] >>> sh) & 255) * sz);
      let a0 = (x0 - ox) * ix, a1 = (x1 - ox) * ix;
      let lo = a0 < a1 ? a0 : a1, hi = a0 < a1 ? a1 : a0;
      a0 = (y0 - oy) * iy; a1 = (y1 - oy) * iy;
      lo = Math.max(lo, a0 < a1 ? a0 : a1); hi = Math.min(hi, a0 < a1 ? a1 : a0);
      a0 = (z0 - oz) * iz; a1 = (z1 - oz) * iz;
      lo = Math.max(lo, a0 < a1 ? a0 : a1); hi = Math.min(hi, a0 < a1 ? a1 : a0);
      hi *= 1 + 4e-7;
      const hit = lo <= hi && hi >= 0 && lo <= best;
      if (inner) {
        if (hit) { ct[nk] = Math.max(lo, 0); ck[nk] = ci; nk++; }
        ci++;
      } else {
        if (hit) {
          const t = leafHit(slot, c, ray, best);
          if (t < best) { best = t; if (any) return best; }
        }
        slot += c;
      }
    }
    // push far to near so the nearest pops first (insertion sort of ≤ 4)
    for (let i = 1; i < nk; i++) { const t = ct[i], k = ck[i]; let j = i - 1; while (j >= 0 && ct[j] < t) { ct[j + 1] = ct[j]; ck[j + 1] = ck[j]; j--; } ct[j + 1] = t; ck[j + 1] = k; }
    if (sp + nk > S.length) { const s2 = new Uint32Array(S.length * 2), t2 = new Float64Array(S.length * 2); s2.set(S); t2.set(T); st.s = s2; st.t = t2; return traverseWide(words, f, root, ray, tmax, leafHit, any, st); }
    for (let i = 0; i < nk; i++) { S[sp] = ck[i]; T[sp] = ct[i]; sp++; }
  }
  return best;
}

export function renderTile({ x0, y0, tw, th, sample0, spp, W, H, world, proj, sky = [0.6, 0.7, 0.85] }, ctx) {
  const nodes = ctx.shared['bvh.nodes'];
  if (cachedNodes !== nodes) { cached = sceneView({ nodes, tris: ctx.shared['bvh.tris'], inst: ctx.shared['bvh.inst'], meta: ctx.shared['bvh.meta'] }); cachedNodes = nodes; }
  const out = new Float32Array(tw * th * 4);
  const o = [world[12], world[13], world[14]], d = [0, 0, 0];
  const T = cached.tris, I = cached.inst;
  for (let y = 0; y < th; y++) for (let x = 0; x < tw; x++) {
    const px = x0 + x, py = y0 + y, k = (y * tw + x) * 4;
    for (let s = sample0; s < sample0 + spp; s++) {
      const jx = hash(px, py, s * 2), jy = hash(py, px, s * 2 + 1);
      const nx = ((2 * (px + jx)) / W - 1) / proj[0], ny = (1 - (2 * (py + jy)) / H) / proj[5];
      d[0] = world[0] * nx + world[4] * ny - world[8]; d[1] = world[1] * nx + world[5] * ny - world[9]; d[2] = world[2] * nx + world[6] * ny - world[10];
      const l = Math.hypot(d[0], d[1], d[2]); d[0] /= l; d[1] /= l; d[2] /= l;
      const h = sceneHit(cached, o, d, Infinity, hit);
      if (!h) { out[k] += sky[0]; out[k + 1] += sky[1]; out[k + 2] += sky[2]; out[k + 3] += 1; continue; }
      // geometric normal (object space), into world space through the instance rows when instanced
      const q = h.slot * 12;
      let n0 = T[q + 5] * T[q + 10] - T[q + 6] * T[q + 9], n1 = T[q + 6] * T[q + 8] - T[q + 4] * T[q + 10], n2 = T[q + 4] * T[q + 9] - T[q + 5] * T[q + 8];
      if (I && hit.instSlot >= 0) {
        // world normal = rowsᵀ · n (rows = world → object)
        const b = hit.instSlot * 16;
        const a0 = I[b] * n0 + I[b + 4] * n1 + I[b + 8] * n2, a1 = I[b + 1] * n0 + I[b + 5] * n1 + I[b + 9] * n2, a2 = I[b + 2] * n0 + I[b + 6] * n1 + I[b + 10] * n2;
        n0 = a0; n1 = a1; n2 = a2;
      }
      const nl = Math.hypot(n0, n1, n2) || 1;
      const c = Math.abs(n0 * d[0] + n1 * d[1] + n2 * d[2]) / nl;
      const v = 0.08 + 0.72 * c;
      out[k] += v; out[k + 1] += v; out[k + 2] += v; out[k + 3] += 1;
    }
  }
  return out;
}

export function castTile({ x0, y0, tw, th, W, H, world, proj, spp = 1, seed = 0 }, ctx) {
  const nodes = ctx.shared['bvh.nodes'];
  if (!nodes) throw new Error('raycast: no scene shared');
  if (cachedNodes !== nodes) { cached = sceneView({ nodes, tris: ctx.shared['bvh.tris'], inst: ctx.shared['bvh.inst'], meta: ctx.shared['bvh.meta'] }); cachedNodes = nodes; }
  const o = [world[12], world[13], world[14]], d = [0, 0, 0];
  let rays = 0, hits = 0, sumT = 0;
  for (let y = y0; y < y0 + th; y++) {
    if ((y & 7) === 0 && ctx.cancelled()) break;
    for (let x = x0; x < x0 + tw; x++) for (let s = 0; s < spp; s++) {
      const jx = hash(x, y, s, seed), jy = hash(y, x, s, seed + 1);
      const nx = ((2 * (x + jx)) / W - 1) / proj[0], ny = (1 - (2 * (y + jy)) / H) / proj[5];
      d[0] = world[0] * nx + world[4] * ny - world[8]; d[1] = world[1] * nx + world[5] * ny - world[9]; d[2] = world[2] * nx + world[6] * ny - world[10];
      const l = Math.hypot(d[0], d[1], d[2]); d[0] /= l; d[1] /= l; d[2] /= l;
      rays++;
      const h = sceneHit(cached, o, d, Infinity, out);
      if (h) { hits++; sumT += h.t; }
    }
  }
  return { rays, hits, sumT };
}