4.2 分辨率、偏差与走样

# snippet: lantern/sim/cg_shade_04_2_bias.py（深度量化步长）
span_world_q16 = far_q16 - near_q16        # 5242880
lsb_world_q16  = span_world_q16 // ONE     # 80

# snippet: lantern/sim/cg_shade_04_2_bias.py（每纹素覆盖的世界尺寸）
texel_world_q16 = (2 * half_world_q16) // RES     # 20480

# snippet: lantern/sim/cg_shade_04_2_bias.py（每纹素深度步长的解析上界）
ndl_flat      = abs(light[1])                                # 37837  (0.5774)
tan_theta     = fsqrt(max(0, ONE - fmul(ndl_flat, ndl_flat)))  # 92682  (1.4142)
tan_theta     = fdiv(tan_theta, max(1, ndl_flat))
step_analytic = fmul(texel_world_q16, tan_theta)             # 28963

need_bias_p50           : 1
need_bias_p95           : 184
need_bias_p99           : 1012
need_bias_max           : 5620
step_measured_median    : 188
step_measured_p95       : 576
step_measured_max       : 5845

#!/usr/bin/env python3
# file: lantern/sim/cg_shade_04_2_bias.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_shade_04_2_bias.py
"""4.2《分辨率、偏差与走样》：偏差的解析推导 + 实测的「自遮挡 ↔ 阴影脱离」权衡曲线。

三个量都直接量出来，不引用任何文献常数：

1. **深度量化步长**：深度图把世界深度区间 `[near, far]` 线性映到 Q16.16 的 `[0, ONE]`，
   故 1 个 LSB 对应多少世界单位是纯算术。
2. **每纹素深度步长**：在深度图上量相邻纹素的深度差，得到「表面斜率 × 纹素尺寸」的
   实测分布；它的高分位数就是偏差必须盖住的自遮挡量。
3. **真值**：用 `shadowlib.ground_truth_mask` 在光源屏幕空间做竖直射线求交，等价于
   **无限分辨率**的深度图，从而把「偏差误差」与「分辨率走样」两个误差源分开。

然后扫偏差，逐档数出 `acne`（真值可见却被判遮挡）与 `peter_panning`（真值被遮挡却判可见），
得到权衡曲线，并落一张 SVG。
"""

import hashlib
import json
import os
import sys

import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
ROOT = os.path.dirname(HERE)
sys.path.insert(0, HERE)
import cglib  # noqa: E402
import cgshade  # noqa: E402
import cg_shade_04_shadowlib as lib  # noqa: E402

OUT = os.path.join(ROOT, "out")
os.makedirs(OUT, exist_ok=True)

ONE = cglib.ONE
ORACLE_EPS = 8  # 真值判定的自遮挡护栏（Q16.16 的 8 个 LSB）

spec = cglib.load_spec(os.path.join(ROOT, "config", "scene_spec.json"))
sh = cglib.load_spec(os.path.join(ROOT, "config", "cg_shade_04_shadow.json"))

W, H = spec["width"], spec["height"]
RES = sh["lightRes"]


def run(res, bias):
    """一次两趟渲染，返回结果字典 + 真值掩码 + 每像素「所需偏差」分布。"""
    r = lib.render_shadow_frame(spec, sh, res, bias, 1)
    truth = lib.ground_truth_mask(r["shadow_tris"], r["light_uvd"], r["covered_mask"], ORACLE_EPS)
    cov = r["covered_mask"]
    return r, truth, cov


# ---------------------------------------------------------------- 解析量

near_q16 = sh["nearQ16"]
far_q16 = sh["farQ16"]
span_world_q16 = far_q16 - near_q16          # 光源正交投影覆盖的世界深度区间
lsb_world_q16 = span_world_q16 // ONE        # 深度 1 LSB 对应的世界单位（Q16.16）
half_world_q16 = sh["orthoHalfQ16"]
texel_world_q16 = (2 * half_world_q16) // RES

light = cglib.v_norm(tuple(sh["lightDirQ16"]))
# 平坦地形的法线 ~ +Y，故 n·L 给出表面与光源方向的夹角余弦
ndl_flat = abs(light[1])
tan_theta = cglib.fsqrt(max(0, ONE - cglib.fmul(ndl_flat, ndl_flat)))
tan_theta = cglib.fdiv(tan_theta, max(1, ndl_flat))
# 每纹素深度步长（解析上界）：纹素尺寸 × tan(theta)
step_analytic = cglib.fmul(texel_world_q16, tan_theta)

# ---------------------------------------------------------------- 实测：深度图相邻纹素差
r0, truth0, cov0 = run(RES, 0)
dmap = r0["shadow_depth"]
valid = dmap != lib.DEPTH_MAX
dx = np.zeros_like(dmap)
dx[:, :-1] = dmap[:, :-1]
dh = np.where(valid[:, :-1] & valid[:, 1:], np.abs(dmap[:, 1:] - dmap[:, :-1]), 0)
dv = np.where(valid[:-1, :] & valid[1:, :], np.abs(dmap[1:, :] - dmap[:-1, :]), 0)
steps = np.concatenate([dh[dh > 0], dv[dv > 0]])
step_median = int(np.median(steps))
step_p95 = int(np.percentile(steps, 95))
step_max = int(steps.max())

# ---------------------------------------------------------------- 实测：每像素「所需偏差」
shift = lib.res_shift(RES)
lu = r0["light_uvd"][:, :, 0]
lv = r0["light_uvd"][:, :, 1]
ld = r0["light_uvd"][:, :, 2]
tx = (lu >> shift) >> 6
ty = (lv >> shift) >> 6
inb = cov0 & (tx >= 0) & (tx < RES) & (ty >= 0) & (ty < RES)
md = np.where(inb, dmap[np.clip(ty, 0, RES - 1), np.clip(tx, 0, RES - 1)], lib.DEPTH_MAX)
delta = np.where(inb & (md != lib.DEPTH_MAX), md - ld, 0)
lit_by_truth = inb & (~truth0) & (md != lib.DEPTH_MAX)
need = delta[lit_by_truth]
need_p50 = int(np.percentile(need, 50))
need_p95 = int(np.percentile(need, 95))
need_p99 = int(np.percentile(need, 99))
need_max = int(need.max())

# ---------------------------------------------------------------- 扫偏差：权衡曲线
curve = []
for b in sh["biasSweep"]:
    r, truth, cov = run(RES, b)
    m = r["shadow_mask"] & cov
    t = truth & cov
    acne = int((m & ~t).sum())
    pan = int((t & ~m).sum())
    hit = int((m & t).sum())
    curve.append({"bias": int(b), "acne": acne, "peter_panning": pan, "hit": hit})

# ---------------------------------------------------------------- 扫分辨率
res_curve = []
for res in sh["resolutionSweep"]:
    r, truth, cov = run(res, sh["pcfBias"])
    m = r["shadow_mask"] & cov
    t = truth & cov
    res_curve.append({
        "res": int(res),
        "acne": int((m & ~t).sum()),
        "peter_panning": int((t & ~m).sum()),
        "map_covered": int(r["map_covered"]),
    })

truth_px = int((truth0 & cov0).sum())

report = {
    "canvas": [W, H],
    "lightRes": RES,
    "oracle_eps": ORACLE_EPS,
    "oracle_shadow_pixels": truth_px,
    "oracle_shadow_permille": round(1000.0 * truth_px / max(1, int(cov0.sum())), 1),
    "depth_span_world_q16": int(span_world_q16),
    "depth_lsb_world_q16": int(lsb_world_q16),
    "texel_world_q16": int(texel_world_q16),
    "ndl_flat_q16": int(ndl_flat),
    "tan_theta_q16": int(tan_theta),
    "step_analytic_q16": int(step_analytic),
    "step_measured_median": step_median,
    "step_measured_p95": step_p95,
    "step_measured_max": step_max,
    "need_bias_p50": need_p50,
    "need_bias_p95": need_p95,
    "need_bias_p99": need_p99,
    "need_bias_max": need_max,
    "bias_curve": curve,
    "resolution_curve": res_curve,
    "per_pixel_sha256_bias0": hashlib.sha256(
        np.clip(r0["color"], 0, 255).astype(np.uint8).tobytes()
    ).hexdigest(),
}

with open(os.path.join(OUT, "cg_shade_04_2_bias.json"), "w") as f:
    json.dump(report, f, indent=2)

# ---------------------------------------------------------------- SVG 曲线
Wp, Hp = 640, 400
pad = 56
xs = [c["bias"] for c in curve]
acne = [c["acne"] for c in curve]
pan = [c["peter_panning"] for c in curve]
xmax = max(xs)
ymax = max(max(acne), max(pan)) or 1


def sx(v):
    return pad + (Wp - 2 * pad) * (v / xmax)


def sy(v):
    return Hp - pad - (Hp - 2 * pad) * (v / ymax)


def poly(vals):
    return " ".join(f"{sx(x):.1f},{sy(y):.1f}" for x, y in zip(xs, vals))


svg = [
    f'<svg xmlns="http://www.w3.org/2000/svg" width="{Wp}" height="{Hp}" viewBox="0 0 {Wp} {Hp}">',
    f'<rect width="{Wp}" height="{Hp}" fill="#0f1116"/>',
    f'<line x1="{pad}" y1="{Hp-pad}" x2="{Wp-pad}" y2="{Hp-pad}" stroke="#4a5568"/>',
    f'<line x1="{pad}" y1="{pad}" x2="{pad}" y2="{Hp-pad}" stroke="#4a5568"/>',
    f'<polyline points="{poly(acne)}" fill="none" stroke="#ff6b6b" stroke-width="2"/>',
    f'<polyline points="{poly(pan)}" fill="none" stroke="#4dabf7" stroke-width="2"/>',
    f'<text x="{Wp/2}" y="{Hp-16}" fill="#cbd5e0" font-size="13" text-anchor="middle">bias (Q16.16 LSB)</text>',
    f'<text x="{pad-8}" y="{pad-16}" fill="#ff6b6b" font-size="12">acne</text>',
    f'<text x="{pad+52}" y="{pad-16}" fill="#4dabf7" font-size="12">peter-panning</text>',
    "</svg>",
]
with open(os.path.join(OUT, "cg_shade_04_2_bias.svg"), "w") as f:
    f.write("\n".join(svg))

for k, v in report.items():
    if k.endswith("curve"):
        print(f"{k}:")
        for row in v:
            print("   ", row)
    else:
        print(f"{k:24s}: {v}")