11.2 ray marching 求交

# snippet: sim/cg_trace_11_sdf.py（节选：sphere tracing 主循环）
# verify: /opt/homebrew/bin/python3 sim/cg_trace_11_2_ray_marching.py
def march(scene, orig, dirv, cfg, scale=ONE, eps=None, max_steps=None, step_scale=ONE):
    """sphere tracing。返回 `(t, steps, status, min_f)`。

    `status` ∈ {"hit", "miss", "maxsteps"}；`min_f` 是过程中遇到过的最小场值
    （`< 0` 说明某一步落进了表面内部——即**穿透**）。
    """
    mc = cfg["marching"]
    eps = mc["epsQ16"] if eps is None else eps
    max_steps = mc["maxSteps"] if max_steps is None else max_steps
    tmax = mc["tMaxQ16"]
    t = 0
    min_f = 1 << 60
    for i in range(max_steps):
        p = v_add(orig, tuple(fmul(dirv[k], t) for k in range(3)))
        f = scene(p, scale)
        if f < min_f:
            min_f = f
        if f < eps:
            return t, i + 1, "hit", min_f
        t += fmul(f, step_scale)
        if t > tmax:
            return None, i + 1, "miss", min_f
    return None, max_steps, "maxsteps", min_f

// file: lantern/core/test/cg_trace_11_2_ray_marching.test.ts
// tier: P
// verify: node --test 'core/test/cg_trace_11_2_ray_marching.test.ts'
// 卷三第 11.2 节《ray marching 求交》的 P 层（上）：基线一遍与容差扫描。
//
// 断言对象：
//   1. 基线一遍 —— 主光线数、解析命中、marching 命中、真漏交、步数耗尽、最大命中误差；
//   2. 容差—漏交率表 —— eps 从 16 扫到 512 的漏交率 / 步数 / 穿透条数。
// 场放大导致的真漏交与不变量见 cg_trace_11_2_miss_vs_overestimate.test.ts。
// 金标准由 S 侧产出：
//   /opt/homebrew/bin/python3 sim/cg_trace_11_2_ray_marching.py

import assert from "node:assert/strict";
import { readFileSync } from "node:fs";
import { dirname, join } from "node:path";
import { test } from "node:test";
import { fileURLToPath } from "node:url";

import {
  type SdfConfig,
  analyticHit,
  cameraBasis,
  makeScene,
  march,
  primaryDir,
} from "../src/lantern/cg_trace_11_sdf.ts";
import { ONE, fdiv, fmul } from "../src/lantern/fixed.ts";
import { type Vec3 } from "../src/lantern/vec.ts";

const root = join(dirname(fileURLToPath(import.meta.url)), "..", "..");
const cfg = JSON.parse(readFileSync(join(root, "config", "cg_trace_11_sdf.json"), "utf8")) as SdfConfig;
const golden = JSON.parse(readFileSync(join(root, "out", "cg_trace_11_2_march.json"), "utf8"));

const scene = makeScene(cfg);

type Ray = { px: number; py: number; d: Vec3 };

function pixelRays(): { eye: Vec3; rays: Ray[] } {
  const cam = cfg.camera;
  const w = cfg.render.width;
  const h = cfg.render.height;
  const tanHalf = BigInt(cam.tanHalfFovYQ16);
  const tanX = fmul(tanHalf, fdiv(BigInt(w), BigInt(h)));
  const [right, up, back] = cameraBasis(cfg);
  const eye: Vec3 = [BigInt(cam.eyeQ16[0]), BigInt(cam.eyeQ16[1]), BigInt(cam.eyeQ16[2])];
  const rays: Ray[] = [];
  for (let py = 0; py < h; py += 1) {
    for (let px = 0; px < w; px += 1) {
      rays.push({ px, py, d: primaryDir(cfg, px, py, right, up, back, tanHalf, tanX) });
    }
  }
  return { eye, rays };
}

const { eye, rays } = pixelRays();
const half = rays.filter((r) => r.px % 2 === 0 && r.py % 2 === 0);

test("11.2 基线：命中数、漏交数、最大命中误差与 S 层一致", () => {
  let hit = 0;
  let miss = 0;
  let exhausted = 0;
  let analyticHits = 0;
  let tErrMax = 0n;
  let minFAll = 1n << 60n;
  for (const r of rays) {
    const ta = analyticHit(cfg, eye, r.d);
    const res = march(scene, eye, r.d, cfg);
    if (res.minF < minFAll) minFAll = res.minF;
    if (res.status === "hit") hit += 1;
    else if (res.status === "maxsteps") exhausted += 1;
    if (ta !== null) {
      analyticHits += 1;
      if (res.status === "miss") miss += 1;
      else if (res.t !== null) {
        const e = res.t > ta ? res.t - ta : ta - res.t;
        if (e > tErrMax) tErrMax = e;
      }
    }
  }
  assert.equal(rays.length, golden.primaryRays);
  assert.equal(analyticHits, golden.baseline.analyticHits);
  assert.equal(hit, golden.baseline.hit);
  assert.equal(miss, golden.baseline.trueMiss);
  assert.equal(exhausted, golden.baseline.exhausted);
  assert.equal(tErrMax, BigInt(golden.baseline.tErrMaxQ16));
  assert.equal(minFAll, BigInt(golden.baseline.minFQ16));
  // 原场是精确距离场：真漏交必须是 0
  assert.equal(miss, 0);
});

test("11.2 容差—漏交率：eps 扫描与 S 层逐档一致", () => {
  const near = (a: number, b: number, what: string): void => {
    assert.ok(Math.abs(a - b) < 1e-3, `${what}: ${a} vs ${b}`);
  };
  const sweep = (scale: bigint, table: typeof golden.toleranceExact): void => {
    for (const g of table) {
      const eps = BigInt(g.epsQ16);
      let analytic = 0;
      let misses = 0;
      let exhausted = 0;
      let pen = 0;
      let terr = 0n;
      let stepsSum = 0;
      for (const r of half) {
        const ta = analyticHit(cfg, eye, r.d);
        const res = march(scene, eye, r.d, cfg, scale, eps);
        stepsSum += res.steps;
        if (res.minF < 0n) pen += 1;
        if (res.status === "maxsteps") exhausted += 1;
        if (ta !== null) {
          analytic += 1;
          if (res.status === "miss") misses += 1;
          else if (res.t !== null) {
            const e = res.t > ta ? res.t - ta : ta - res.t;
            if (e > terr) terr = e;
          }
        }
      }
      const permille = (1000 * misses) / analytic;
      near(permille, g.missPermille, `eps=${g.epsQ16} 漏交率`);
      near(stepsSum / half.length, g.meanSteps, `eps=${g.epsQ16} 平均步数`);
      assert.equal(misses, g.misses, `eps=${g.epsQ16} 漏交数`);
      assert.equal(exhausted, g.exhausted, `eps=${g.epsQ16} 耗尽数`);
      assert.equal(pen, g.penetrations, `eps=${g.epsQ16} 穿透数`);
      assert.equal(terr, BigInt(g.maxTErrQ16), `eps=${g.epsQ16} 最大命中误差`);
    }
  };
  sweep(ONE, golden.toleranceExact);
  sweep((5n * ONE) / 4n, golden.toleranceOver);
});

# snippet: 本机实跑（工作区 lantern/）
/opt/homebrew/bin/python3 sim/cg_trace_11_2_ray_marching.py
node --test 'core/test/cg_trace_11_2_*.test.ts'

== 11.2 ray marching 求交 ==
主光线 3072（64x48），解析命中 2048，marching 命中 1995，真漏交 0，步数耗尽 53
基线 max|t - t_analytic| = 2046 Q16.16，过程最小场值 2（< 0 即穿透）
迭代次数分箱（每 8 步一箱），maxSteps = 256：
    1..  8 : 151
    9.. 16 : 1245
   17.. 24 : 667
   25.. 32 : 288
   33.. 40 : 160
   41.. 48 : 90
   49.. 56 : 105
   57.. 64 : 31
   65.. 72 : 36
   73.. 80 : 10
   81.. 88 : 32
   89.. 96 : 13
  105..112 : 48
  113..120 : 36
  121..128 : 29
  129..136 : 1
  137..144 : 1
  145..152 : 21
  153..160 : 21
  161..168 : 13
  169..176 : 8
  177..184 : 2
  249..256 : 64
容差—漏交率（原场 / 放大 1.25；漏交 = 真越过表面，耗尽 = 步数用完）:
  eps= 16  原场 漏交     0.0‰ 耗尽   32 步数 37.457   |  放大 漏交     0.0‰ 穿透   49 步数  31.57
  eps= 32  原场 漏交     0.0‰ 耗尽   32 步数 36.072   |  放大 漏交     0.0‰ 穿透   49 步数 30.053
  eps= 64  原场 漏交     0.0‰ 耗尽   26 步数 34.721   |  放大 漏交     0.0‰ 穿透   49 步数 28.367
  eps=128  原场 漏交     0.0‰ 耗尽    9 步数 32.919   |  放大 漏交     0.0‰ 穿透   49 步数 26.632
  eps=256  原场 漏交     0.0‰ 耗尽    0 步数 30.767   |  放大 漏交     0.0‰ 穿透   49 步数 24.882
  eps=512  原场 漏交     0.0‰ 耗尽    0 步数 28.484   |  放大 漏交     0.0‰ 穿透   49 步数 23.121
步长放大—穿透（原场 / 放大 1.25）:
  x0.5  原场 穿透    0 漏交     0.0‰   |  放大 穿透    0 漏交     0.0‰ max|Δt| 810
  x0.75 原场 穿透    0 漏交     0.0‰   |  放大 穿透    0 漏交     0.0‰ max|Δt| 808
  x1.0  原场 穿透    0 漏交     0.0‰   |  放大 穿透   49 漏交     0.0‰ max|Δt| 14577
  x1.25 原场 穿透   49 漏交     0.0‰   |  放大 穿透   64 漏交     0.0‰ max|Δt| 62151
  x1.5  原场 穿透   57 漏交     0.0‰   |  放大 穿透  143 漏交     0.0‰ max|Δt| 1703079
  x2.0  原场 穿透  178 漏交     0.0‰   |  放大 穿透  268 漏交     0.0‰ max|Δt| 180504
无地面（球+盒）：场放大倍数 -> 真漏交（Lipschitz 被违反时）:
  x1.0   解析命中  231 真漏交    0（     0.0‰）穿透    0 max|Δt| 1251
  x1.1   解析命中  231 真漏交    0（     0.0‰）穿透  212 max|Δt| 40263
  x1.25  解析命中  231 真漏交   70（303.0303‰）穿透  161 max|Δt| 100811
  x1.5   解析命中  231 真漏交  231（  1000.0‰）穿透    0 max|Δt| 0
  x2.0   解析命中  231 真漏交  231（  1000.0‰）穿透    0 max|Δt| 0
csv/json/png      : out/cg_trace_11_2_{march.json,iterations,tolerance,stepscale}.csv