4.3 软阴影与 PCF

#!/usr/bin/env python3
# file: lantern/sim/cg_shade_04_3_pcf.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_shade_04_3_pcf.py
"""4.3《软阴影与 PCF》：N×N 邻域比对、半影半径与核大小的权衡、以及 PCF 的开销曲线。

场景换成一块**地面 + 一面悬空的墙**（8 顶点 / 4 三角形）：地形那种起伏会把半影
淹没在自遮挡里，看不出 PCF 在做什么。方向光刻意调成让墙的影子朝相机方向落。

PCF 判定仍是纯整数：N² 个纹素各做一次 `ld - map > bias` 比较，`occ` 累加，可见度
= `(N² - occ) / N²`。除法只在**非负**分子上做，两侧同为 floor 语义。

质量用「与 N=15 参考的可见度 MAE」度量（把 N=15 当作半影的数值参考），
开销用整帧墙钟毫秒。
"""

import hashlib
import json
import os
import sys
import time

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 cg_shade_04_shadowlib as lib  # noqa: E402

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

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

RES = sh["lightRes"]
BIAS = 2048  # 4.2 的曲线说明：这个斜率下 ~2048 才能让 N≥5 的邻域不误伤
REF_N = 15


def lit_fraction(r, pcf_n):
    """可见度 = (N² - occ) / N²，整数量化到 Q16.16。"""
    occ = r["occ_count"].astype(np.int64)
    den = pcf_n * pcf_n
    return ((den - occ) << 16) // den


def run(pcf_n, bias=BIAS):
    t0 = time.perf_counter()
    r = lib.render_explicit_scene(scene, sh, RES, bias, pcf_n)
    ms = (time.perf_counter() - t0) * 1000.0
    return r, ms


ref, _ = run(REF_N)
ref_lit = lit_fraction(ref, REF_N)
cov = ref["covered_mask"]

rows = []
for n in sh["pcfSizes"]:
    r, ms = run(n)
    lit = lit_fraction(r, n)
    occ = r["occ_count"].astype(np.int64)
    den = n * n
    penumbra = int(((occ > 0) & (occ < den) & cov).sum())
    full = int(((occ == den) & cov).sum())
    diff = np.abs(lit[cov] - ref_lit[cov])
    mae = float(diff.mean()) / 65536.0
    rows.append({
        "n": int(n),
        "samples": int(den),
        "penumbra_pixels": penumbra,
        "full_shadow_pixels": full,
        "shadow_pixels": int(r["shadow_pixels"]),
        "mae_vs_ref_q16": int(diff.mean()),
        "mae_vs_ref": round(mae, 5),
        "ms": round(ms, 1),
        "per_pixel_us": round(ms * 1000.0 / max(1, int(cov.sum())), 3),
        "sha256": hashlib.sha256(np.clip(r["color"], 0, 255).astype(np.uint8).tobytes()).hexdigest(),
    })

# 半影半径：在影子中心那一行上量「从全亮到全暗」跨越多少像素
ys, xs = np.nonzero(cov)
mid_y = int(np.median(ys[ys > 0])) if len(ys) else 0
widths = {}
for n in sh["pcfSizes"]:
    r, _ = run(n)
    occ = r["occ_count"].astype(np.int64)
    den = n * n
    row = cov[mid_y]
    soft = (occ[mid_y] > 0) & (occ[mid_y] < den) & row
    widths[str(n)] = int(soft.sum())

cglib.save_png(os.path.join(OUT, "cg_shade_04_3_pcf1.png"), lib.render_explicit_scene(scene, sh, RES, BIAS, 1)["color"])
cglib.save_png(os.path.join(OUT, "cg_shade_04_3_pcf7.png"), lib.render_explicit_scene(scene, sh, RES, BIAS, 7)["color"])

report = {
    "canvas": [scene["width"], scene["height"]],
    "lightRes": RES,
    "bias": BIAS,
    "reference_n": REF_N,
    "covered_pixels": int(cov.sum()),
    "scanline_y": mid_y,
    "penumbra_width_on_scanline": widths,
    "pcf_curve": rows,
}
with open(os.path.join(OUT, "cg_shade_04_3_pcf.json"), "w") as f:
    json.dump(report, f, indent=2)

for k, v in report.items():
    if k == "pcf_curve":
        print("pcf_curve:")
        for row in v:
            print("   ", row)
    else:
        print(f"{k:28s}: {v}")

canvas                      : [320, 240]
lightRes                    : 256
bias                        : 2048
reference_n                 : 15
covered_pixels              : 22704
scanline_y                  : 173
penumbra_width_on_scanline  : {'1': 0, '3': 76, '5': 88, '7': 94}

// snippet: core/src/lantern/cg_shade_04_shadow.ts（PCF 邻域累加）
// verify: node --test 'core/test/cg_shade_04_shadow.test.ts'
export function shadowLit(
  ld: number, lu: number, lv: number,
  shadowDepth: Float64Array, res: number, bias: number, pcfN: number,
): [number, number] {
  const shift = resShift(res);
  if (pcfN <= 1) {                                   // 硬阴影：单个纹素
    const tx = (lu >> shift) >> 6;
    const ty = (lv >> shift) >> 6;
    if (tx >= 0 && tx < res && ty >= 0 && ty < res) {
      const md = shadowDepth[ty * res + tx];
      if (md !== DEPTH_MAX && ld - md > bias) return [0, 1];
    }
    return [1, 1];
  }
  const half = (pcfN - 1) >> 1;                      // N=3 -> 1, N=5 -> 2, N=7 -> 3
  let occ = 0;
  let total = 0;
  for (let dy = -half; dy <= half; dy += 1) {
    const ty = ((lv >> shift) >> 6) + dy;
    for (let dx = -half; dx <= half; dx += 1) {
      const tx = ((lu >> shift) >> 6) + dx;
      total += 1;
      if (tx >= 0 && tx < res && ty >= 0 && ty < res) {
        const md = shadowDepth[ty * res + tx];
        if (md !== DEPTH_MAX && ld - md > bias) occ += 1;
      }
    }
  }
  return [total - occ, total];                       // 分子分母都非负
}