5.3 视锥与近远平面

近平面取 0 时（depth 恒为 ONE，深度缓冲失效）：
  D= 1  depth=65536
  D= 4  depth=65536
  D=16  depth=65536
  D=64  depth=65536

#!/usr/bin/env python3
# file: lantern/sim/cg_raster_05_3_depth_precision.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_raster_05_3_depth_precision.py
"""5.3 视锥与近远平面：视锥六平面的解析构造、近平面不能取 0 的原因，
以及深度缓冲的精度分布 depth(D) = far/(far-near)*(1 - near/D)。

近平面越靠近 0，越多深度量程被挤到相机附近 —— 第 10 章「深度冲突」的伏笔。

产物（SVG 全部由标准库 xml.etree 手写并回读校验，本机无 matplotlib）：
  out/cg_raster_05_3_depth_precision.svg  左：depth-D 曲线；右：各倍程占用量程
  out/cg_raster_05_3_depth_curve.csv / _octaves.csv / _frustum.json
"""

import json
import math
import os
import sys
import xml.etree.ElementTree as ET

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib  # noqa: E402

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
ONE = cglib.ONE
SVG_NS = "http://www.w3.org/2000/svg"


def frustum_planes(mvp):
    """Gribb–Hartmann：从 MVP 直接抽出视锥六平面（Q16.16）。

    世界点 p 在视锥内 <=> 对六平面都有 plane·(p,1) >= 0。
    """
    r0, r1, r2, r3 = mvp[0:4], mvp[4:8], mvp[8:12], mvp[12:16]
    return {
        "left": [r3[i] + r0[i] for i in range(4)],
        "right": [r3[i] - r0[i] for i in range(4)],
        "bottom": [r3[i] + r1[i] for i in range(4)],
        "top": [r3[i] - r1[i] for i in range(4)],
        "near": [r3[i] + r2[i] for i in range(4)],
        "far": [r3[i] - r2[i] for i in range(4)],
    }


def plane_eval(pl, p):
    """Q16.16 平面与 Q16.16 世界点求值（齐次 w = ONE），保留符号。"""
    return pl[0] * p[0] + pl[1] * p[1] + pl[2] * p[2] + pl[3] * ONE


def inside(planes, p):
    return all(plane_eval(pl, p) >= 0 for pl in planes.values())


def project_inside(mvp, p):
    """用 project() 独立判定：非 None 且 |ndc|<=1 且 depth in [0,ONE]。"""
    r = cglib.project(mvp, p)
    if r is None:
        return False
    nx, ny, d, _w = r
    return -ONE <= nx <= ONE and -ONE <= ny <= ONE and 0 <= d <= ONE


def make_mvp(tan_half, aspect, near, far):
    """规范针孔：眼在原点、看 -z。世界点 (0,0,-D) 距相机 D。"""
    view0 = cglib.mat_look_at((0, 0, 0), (0, 0, -ONE), (0, ONE, 0))
    return cglib.mat_mul(cglib.mat_perspective(tan_half, aspect, near, far), view0)


def depth_at(mvp, distance_q16):
    """世界距离（Q16）-> 归一化深度（Q16，[0, ONE]）；相机后方返回 None。"""
    r = cglib.project(mvp, (0, 0, -distance_q16))
    return None if r is None else r[2]


def build_svg(path, curve, octaves, near_f, far_f):
    """纯 xml.etree 手写 SVG：左 depth-D 曲线（对数 x），右各倍程占用量程。"""
    W, H = 760, 380
    svg = ET.Element("svg", {"xmlns": SVG_NS, "width": str(W), "height": str(H),
                             "viewBox": f"0 0 {W} {H}"})

    def text(x, y, s, size=12, anchor="start", weight="normal"):
        e = ET.SubElement(svg, "text", {"x": f"{x:.1f}", "y": f"{y:.1f}",
                                        "font-family": "sans-serif", "font-size": str(size),
                                        "text-anchor": anchor, "font-weight": weight, "fill": "#333"})
        e.text = s

    text(W / 2, 24, "深度精度分布：depth(D) = far/(far-near)·(1 - near/D)", 16, "middle", "bold")

    # 左：depth-D 曲线，x 轴对数
    ax0, ax1, ay0, ay1 = 70, 370, 50, 320
    lx0, lx1 = math.log10(near_f), math.log10(far_f)
    ax = lambda d: ax0 + (math.log10(d) - lx0) / (lx1 - lx0) * (ax1 - ax0)
    ay = lambda v: ay1 - v * (ay1 - ay0)
    ET.SubElement(svg, "rect", {"x": str(ax0), "y": str(ay0), "width": str(ax1 - ax0),
                                "height": str(ay1 - ay0), "fill": "none", "stroke": "#999"})
    for k in range(5):
        v = k / 4
        ET.SubElement(svg, "line", {"x1": str(ax0), "y1": f"{ay(v):.1f}", "x2": str(ax1),
                                    "y2": f"{ay(v):.1f}", "stroke": "#eee"})
        text(ax0 - 8, ay(v) + 4, f"{v:.2f}", 10, "end")
    pts = " ".join(f"{ax(d):.1f},{ay(dep / ONE):.1f}" for d, dep in curve)
    ET.SubElement(svg, "polyline", {"points": pts, "fill": "none", "stroke": "#c0392b", "stroke-width": "2"})
    for d in (0.5, 1, 2, 4, 8, 16, 32, 64):
        if near_f <= d <= far_f:
            text(ax(d), ay1 + 16, f"{d:g}", 10, "middle")
            ET.SubElement(svg, "line", {"x1": f"{ax(d):.1f}", "y1": str(ay1),
                                        "x2": f"{ax(d):.1f}", "y2": f"{ay1 + 4:.1f}", "stroke": "#999"})
    text((ax0 + ax1) / 2, ay1 + 36, "世界距离 D（对数轴）", 12, "middle")
    text(20, (ay0 + ay1) / 2, "归一化深度", 12, "middle")
    text((ax0 + ax1) / 2, ay0 - 12, "depth(D) 曲线", 13, "middle", "bold")

    # 右：各倍程占用的深度量程
    bx0, bx1, by0, by1 = 450, 730, 50, 320
    ET.SubElement(svg, "rect", {"x": str(bx0), "y": str(by0), "width": str(bx1 - bx0),
                                "height": str(by1 - by0), "fill": "none", "stroke": "#999"})
    bw = (bx1 - bx0) / len(octaves)
    maxv = max(o["share_permille"] for o in octaves)
    for i, o in enumerate(octaves):
        h = (o["share_permille"] / maxv) * (by1 - by0 - 20)
        x = bx0 + i * bw + 3
        ET.SubElement(svg, "rect", {"x": f"{x:.1f}", "y": f"{by1 - h:.1f}", "width": f"{bw - 6:.1f}",
                                    "height": f"{h:.1f}", "fill": "#2980b9"})
        text(x + (bw - 6) / 2, by1 - h - 4, str(o["share_permille"]), 9, "middle", "bold")
        text(x + (bw - 6) / 2, by1 + 14, o["label"], 8, "middle")
    text((bx0 + bx1) / 2, ay0 - 12, "各倍程占用深度量程（‰）", 13, "middle", "bold")
    text((bx0 + bx1) / 2, by1 + 34, "距离倍程 [near·2^k, near·2^(k+1))", 11, "middle")

    ET.indent(svg, space="  ")
    ET.ElementTree(svg).write(path, encoding="utf-8", xml_declaration=True)
    root = ET.parse(path).getroot()  # 回读校验
    assert root.tag.endswith("svg"), root.tag
    return os.path.getsize(path), len(list(root.iter()))


def main():
    spec = cglib.load_spec()
    cam = spec["camera"]
    tan_half, near, far = cam["tanHalfFovYQ16"], cam["nearQ16"], cam["farQ16"]
    near_f, far_f = near / ONE, far / ONE
    aspect = cglib.fdiv(cglib.f_from_int(spec["width"]), cglib.f_from_int(spec["height"]))
    mvp = make_mvp(tan_half, aspect, near, far)
    out = os.path.join(ROOT, "out")

    # 1. 六平面：解析抽取 vs project() 独立判定
    planes = frustum_planes(mvp)
    agree = total = 0
    samples = []
    for dq in (ONE, 2 * ONE, 8 * ONE, 32 * ONE, 63 * ONE):
        for xq in (0, 8 * ONE, 32 * ONE):
            for yq in (0, 6 * ONE, 24 * ONE):
                p = (xq, yq, -dq)
                total += 1
                if inside(planes, p) == project_inside(mvp, p):
                    agree += 1
                else:
                    samples.append({"p": p, "planes": inside(planes, p), "project": project_inside(mvp, p)})

    # 2. 深度映射：曲线采样 + 各倍程占用
    curve = [(near_f * (far_f / near_f) ** (i / 240),
              depth_at(mvp, round(near_f * (far_f / near_f) ** (i / 240) * ONE))) for i in range(241)]
    octaves = []
    for k in range(7):
        d0, d1 = near_f * (2 ** k), near_f * (2 ** (k + 1))
        if d1 > far_f:
            break
        dep0, dep1 = depth_at(mvp, round(d0 * ONE)), depth_at(mvp, round(d1 * ONE))
        octaves.append({"label": f"{d0:g}–{d1:g}", "d0": d0, "d1": d1,
                        "depth0": dep0, "depth1": dep1, "share_permille": (dep1 - dep0) * 1000 // ONE})

    # 3. 世界距离分辨率：一个深度 LSB 对应多少世界距离（数值差分，步长随距离放大）
    def dd_at(dq):
        step = max(ONE // 4, dq // 64)
        return (depth_at(mvp, dq + step) - depth_at(mvp, dq - step)) / (2.0 * step)

    reso = []
    for d in (1, 2, 4, 8, 16, 32, 64):
        dd = dd_at(d * ONE)
        reso.append({"distance": d, "world_per_depth_lsb": round((1.0 / dd) / ONE, 8) if dd > 0 else None})

    # 4. 近平面取 0 的演示
    mvp0 = make_mvp(tan_half, aspect, 0, far)
    near_zero = [{"distance": d, "depth": depth_at(mvp0, d * ONE)} for d in (1, 4, 16, 64)]

    cglib.write_csv(os.path.join(out, "cg_raster_05_3_depth_curve.csv"),
                    ["distance", "depth_q16"], [[f"{d:.6f}", dep] for d, dep in curve])
    cglib.write_csv(os.path.join(out, "cg_raster_05_3_octaves.csv"),
                    ["octave", "d0", "d1", "depth0", "depth1", "share_permille"],
                    [[o["label"], f"{o['d0']:g}", f"{o['d1']:g}", o["depth0"], o["depth1"],
                      o["share_permille"]] for o in octaves])
    svg_bytes, svg_nodes = build_svg(os.path.join(out, "cg_raster_05_3_depth_precision.svg"),
                                     curve, octaves, near_f, far_f)
    with open(os.path.join(out, "cg_raster_05_3_frustum.json"), "w", encoding="utf-8") as fh:
        json.dump({
            "near": near_f, "far": far_f, "depthBits": spec["depthBits"],
            "frustum_planes": planes,
            "plane_vs_project_agreement": {"agree": agree, "total": total, "mismatch": samples},
            "octaves": octaves, "resolution": reso, "near_zero_demo": near_zero,
            "svg": {"bytes": svg_bytes, "nodes": svg_nodes},
        }, fh, ensure_ascii=False, indent=2)

    print(f"near/far      : {near_f} / {far_f}   depthBits={spec['depthBits']}")
    print("视锥六平面（Q16.16，[a b c d]，内部满足 a·x+b·y+c·z+d >= 0）：")
    for k, pl in planes.items():
        print(f"  {k:<7}: {pl}")
    print(f"平面判定 vs project() 一致: {agree}/{total}  （不一致 {len(samples)}）")
    print("各倍程占用的深度量程（‰）：")
    for o in octaves:
        print(f"  {o['label']:>10}  depth {o['depth0']} -> {o['depth1']}  占 {o['share_permille']}‰")
    print("世界距离分辨率（1 个深度 LSB 对应多少世界单位）：")
    for r in reso:
        print(f"  D={r['distance']:>2}  {r['world_per_depth_lsb']} 世界单位")
    print("近平面取 0 时（depth 恒为 ONE，深度缓冲失效）：")
    for r in near_zero:
        print(f"  D={r['distance']:>2}  depth={r['depth']}")
    print(f"SVG           : {svg_bytes} bytes / {svg_nodes} 节点（已回读校验）")
    return 0


if __name__ == "__main__":
    sys.exit(main())

// snippet: 5.3 P 层核心断言
const GOLDEN_OCTAVE_PERMILLE = [503, 251, 125, 62, 31, 15, 7];

test("各倍程占用的深度量程与 Python 侧逐位相同", () => {
  const shares: number[] = [];
  for (let k = 0; k < 7; k += 1) {
    const d0 = NEAR * (1n << BigInt(k));
    const d1 = NEAR * (1n << BigInt(k + 1));
    shares.push(Number(((depthAt(d1) - depthAt(d0)) * 1000n) / ONE));
  }
  assert.deepEqual(shares, GOLDEN_OCTAVE_PERMILLE);
});

test("近平面取 0：所有距离的深度都塌到 ONE，深度缓冲失效", () => {
  const m0 = makeMvp(0n);
  for (const d of [1, 4, 16, 64]) assert.equal(depthAt(fFromInt(d), m0), ONE);
});

node --test 'core/test/cg_raster_05_3_depth_precision.test.ts'
ℹ tests 6   ℹ pass 6   ℹ fail 0