16.1 光栅化管线的装配

#!/usr/bin/env python3
# file: lantern/sim/cg_raster_16_1_assembly.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_raster_16_1_assembly.py
"""16.1 光栅化管线的装配：把前 15 章的环节串成一条完整管线。

逐级核对契约与不变量，并断言串联后的成片 sha256 与卷一基线一致：
  ① 顶点变换的 w 恒等于 (-z_view_full)>>16（透视矩阵第三行为 (0,0,-ONE,0)）
  ② 重心权重 w0+w1+w2 恒等于两倍三角形面积 area（与求值点无关）
  ③ 覆盖像素的深度落在 [0, ONE]，未覆盖处落回哨兵
产物：out/cg_raster_16_1_assembly.json。
"""

import hashlib
import json
import os
import sys
import time

import numpy as np

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

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
BASE_SHA = "bac659770178037ad3f6a6ad184945a96ff37a7fe4d754ef1ebc5b0aaa538ced"
STAGES = ["规格->MVP", "顶点变换", "近平面剔除", "屏幕映射", "三角形装配",
          "背面剔除", "边函数覆盖", "重心插值", "深度测试", "像素输出"]
SENTINEL = int(np.iinfo(np.int64).max)


def clip_stats(mvp, view, verts):
    """第 2/3 级：w == (-z_view_full)>>16 不变量 + 近平面剔除计数。

    透视矩阵第三行是 (0, 0, -ONE, 0)，故 w 完全由**未截断**的视图空间 z 点积决定。
    """
    kept = w_ok = 0
    for v in verts:
        w = cglib.mat_vec4(mvp, (v[0], v[1], v[2], cglib.ONE))[3]
        z_full = view[8] * v[0] + view[9] * v[1] + view[10] * v[2] + view[11] * cglib.ONE
        w_ok += w == ((-z_full) >> 16)
        kept += w > 0
    return kept, len(verts) - kept, w_ok


def coverage_invariant(rtris, width, height):
    """第 7/8 级：对每个覆盖判定为真的像素断言 w0+w1+w2 == area（向量化）。"""
    half = cglib.SUB1 // 2
    checked = bad = 0
    for t in rtris:
        x0, y0, x1, y1, x2, y2 = t[0], t[1], t[6], t[7], t[12], t[13]
        area = (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0)
        if area <= 0:
            continue
        ax = np.arange(max(0, min(x0, x1, x2) // cglib.SUB1),
                       min(width - 1, max(x0, x1, x2) // cglib.SUB1) + 1, dtype=np.int64)
        ay = np.arange(max(0, min(y0, y1, y2) // cglib.SUB1),
                       min(height - 1, max(y0, y1, y2) // cglib.SUB1) + 1, dtype=np.int64)
        if len(ax) == 0 or len(ay) == 0:
            continue
        PX, PY = np.meshgrid(ax * cglib.SUB1 + half, ay * cglib.SUB1 + half)
        w0 = (x2 - x1) * (PY - y1) - (y2 - y1) * (PX - x1)
        w1 = (x0 - x2) * (PY - y2) - (y0 - y2) * (PX - x2)
        w2 = (x1 - x0) * (PY - y0) - (y1 - y0) * (PX - x0)
        mask = (w0 >= 0) & (w1 >= 0) & (w2 >= 0)
        checked += int(mask.sum())
        bad += int((mask & ((w0 + w1 + w2) != area)).sum())
    return checked, bad


def main():
    t0 = time.perf_counter()
    spec = cglib.load_spec()
    width, height = spec["width"], spec["height"]
    cam = spec["camera"]

    # ---- 1 规格 -> MVP；2/3 顶点变换与近平面剔除
    aspect = cglib.fdiv(cglib.f_from_int(width), cglib.f_from_int(height))
    view = cglib.mat_look_at(tuple(cam["eyeQ16"]), tuple(cam["targetQ16"]), tuple(cam["upQ16"]))
    proj = cglib.mat_perspective(cam["tanHalfFovYQ16"], aspect, cam["nearQ16"], cam["farQ16"])
    mvp = cglib.mat_mul(proj, view)
    verts, tris, vcolors = render_raster.build_scene(spec)
    kept, dropped_v, w_ok = clip_stats(mvp, view, verts)

    # ---- 4 屏幕映射
    screen, lo_x, hi_x, lo_y, hi_y = [], 1 << 40, -(1 << 40), 1 << 40, -(1 << 40)
    for v in verts:
        p = cglib.project(mvp, v)
        if p is None:
            screen.append(None)
            continue
        sx, sy = cglib.to_screen(p[0], p[1], width, height)
        lo_x, hi_x, lo_y, hi_y = min(lo_x, sx), max(hi_x, sx), min(lo_y, sy), max(hi_y, sy)
        screen.append((sx, sy, p[2]))

    # ---- 5 三角形装配
    rtris, dropped_t = [], 0
    for a, b, c in tris:
        sa, sb, sc = screen[a], screen[b], screen[c]
        if sa is None or sb is None or sc is None:
            dropped_t += 1
            continue
        rtris.append((sa[0], sa[1], sa[2], vcolors[a][0], vcolors[a][1], vcolors[a][2],
                      sb[0], sb[1], sb[2], vcolors[b][0], vcolors[b][1], vcolors[b][2],
                      sc[0], sc[1], sc[2], vcolors[c][0], vcolors[c][1], vcolors[c][2]))

    # ---- 6..9 背面剔除 -> 边函数覆盖 -> 重心插值 -> 深度测试
    color, depth, covered, culled = cglib.rasterize(rtris, width, height, spec["backgroundQ8"])
    checked, bad = coverage_invariant(rtris, width, height)

    # ---- 10 像素输出
    px = np.clip(color, 0, 255).astype(np.uint8)
    digest = hashlib.sha256(px.tobytes()).hexdigest()
    fin = depth[covered]
    d_lo, d_hi = int(fin.min()), int(fin.max())
    uncov_ok = bool((depth[~covered] == SENTINEL).all())
    assert digest == BASE_SHA, "串联后的成片 sha256 与基线不一致"
    assert w_ok == len(verts), "w == (-z_view_full)>>16 不变量被破坏"
    assert bad == 0, "重心权重和不等于面积"
    assert 0 <= d_lo <= d_hi <= cglib.ONE, "覆盖像素深度越界"
    assert uncov_ok, "未覆盖像素未落回哨兵"

    ms = (time.perf_counter() - t0) * 1000
    out = os.path.join(ROOT, "out")
    os.makedirs(out, exist_ok=True)
    stats = {"canvas": [width, height], "stages": STAGES,
             "vertices": len(verts), "triangles": len(tris),
             "vertices_kept": kept, "vertices_clipped": dropped_v, "w_invariant_ok": w_ok,
             "screen_x_range_q6": [lo_x, hi_x], "screen_y_range_q6": [lo_y, hi_y],
             "triangles_assembled": len(rtris), "triangles_dropped": dropped_t,
             "triangles_culled": culled, "coverage_pixels": checked,
             "coverage_invariant_violations": bad, "pixels_covered": int(covered.sum()),
             "depth_range_q16": [d_lo, d_hi], "depth_sentinel_uncovered": uncov_ok,
             "sha256": digest, "sha256_matches_baseline": digest == BASE_SHA, "ms": round(ms, 1)}
    with open(os.path.join(out, "cg_raster_16_1_assembly.json"), "w", encoding="utf-8") as fh:
        json.dump(stats, fh, ensure_ascii=False, indent=2)

    print(f"顶点变换        : {len(verts)} 顶点，w == (-z_view_full)>>16 成立 {w_ok}/{len(verts)}")
    print(f"近平面剔除      : 保留 {kept}，丢弃 {dropped_v}（w<=0 即相机后方）")
    print(f"屏幕映射        : sx in [{lo_x}, {hi_x}]  sy in [{lo_y}, {hi_y}] Q6")
    print(f"三角形装配      : {len(rtris)} 装配 / {len(tris)} 总数（相机后方丢弃 {dropped_t}）")
    print(f"背面剔除        : 剔除 {culled}，保留 {len(rtris) - culled}")
    print(f"覆盖判定        : {checked} 像素，w0+w1+w2 == area 违反 {bad} 次")
    print(f"深度测试        : 覆盖 {int(covered.sum())} 像素，深度 in [{d_lo}, {d_hi}] Q16.16")
    print(f"未覆盖深度      : 落回哨兵 {SENTINEL} -> {uncov_ok}")
    print(f"sha256(rgb)     : {digest}")
    print(f"与基线一致      : {digest == BASE_SHA}   装配耗时 {ms:.1f} ms")
    return 0


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

顶点变换        : 487 顶点，w == (-z_view_full)>>16 成立 487/487
近平面剔除      : 保留 487，丢弃 0（w<=0 即相机后方）
屏幕映射        : sx in [-96501, 117959]  sy in [5218, 30345] Q6
三角形装配      : 896 装配 / 896 总数（相机后方丢弃 0）
背面剔除        : 剔除 290，保留 606
覆盖判定        : 47673 像素，w0+w1+w2 == area 违反 0 次
深度测试        : 覆盖 43731 像素，深度 in [60516, 65335] Q16.16
未覆盖深度      : 落回哨兵 9223372036854775807 -> True
sha256(rgb)     : bac659770178037ad3f6a6ad184945a96ff37a7fe4d754ef1ebc5b0aaa538ced
与基线一致      : True   装配耗时 51.1 ms

// snippet: 16.1 的 P 层断言（摘自 lantern/core/test/cg_raster_16_1_assembly.test.ts）
// 16.1 的装配断言：TypeScript 侧把同一份规格串成完整管线，逐级核对十级契约与不变量，
// 并断言成片 sha256 与 S 层基线**逐字节一致**；负例（顶点顺序交换）必须让摘要分叉。
// 契约与计数由 S 侧 out/cg_raster_16_1_assembly.json 独立产出，本文件只做对账。

test("装配成片 + 十级契约对账：复现基线 sha256 与逐级计数", () => {
  const frame = renderScene(spec);
  assert.equal(sha(frame.result.rgb), BASE);
  assert.equal(frame.scene.verts.length, 487);
  assert.equal(frame.scene.tris.length, 896);
  assert.equal(frame.dropped, 0, "相机后方丢弃必须为 0");
  assert.equal(frame.result.culled, 290, "背面剔除必须为 290");
  assert.equal(frame.result.coveredCount, 43731);
  assert.equal(frame.screen.filter((s) => s !== null).length, 487);
  assert.equal(g.stages.length, 10, "契约表必须是十级");
  assert.equal(g.vertices, 487);
  assert.equal(g.triangles, 896);
  assert.equal(g.triangles_culled, 290);
  assert.equal(g.pixels_covered, 43731);
  assert.equal(g.w_invariant_ok, 487, "w 不变量必须对全部顶点成立");
  assert.equal(g.coverage_invariant_violations, 0);
  assert.equal(g.sha256, BASE);
});

test("不变量：重心权重和 == 面积（对全部覆盖判定像素）", () => {
  const { rtris } = prepare(spec);
  let checked = 0;
  for (const t of rtris) {
    const [x0, y0] = t.xy[0], [x1, y1] = t.xy[1], [x2, y2] = t.xy[2];
    ...
        if (w0 < 0 || w1 < 0 || w2 < 0) continue;
        assert.equal(w0 + w1 + w2, area, "重心权重和不等于面积");
        checked += 1;
      }
    }
  }
  assert.equal(checked, g.coverage_pixels, "覆盖判定像素数与 S 层不一致");
});

test("负例：三角形顶点顺序交换（b<->c）必须让成片摘要分叉", () => {
  const { rtris } = prepare(spec);
  const swapped = rtris.map((t) => ({
    xy: [t.xy[0], t.xy[2], t.xy[1]],
    d: [t.d[0], t.d[2], t.d[1]],
    c: [t.c[0], t.c[2], t.c[1]],
  }));
  const r = rasterize(swapped, W, H, spec.backgroundQ8);
  assert.notEqual(sha(r.rgb), BASE, "顶点顺序交换后摘要必须分叉");
});

✔ 装配成片 + 十级契约对账：复现基线 sha256 与逐级计数
✔ 不变量：重心权重和 == 面积（对全部覆盖判定像素）
✔ 负例：三角形顶点顺序交换（b<->c）必须让成片摘要分叉
ℹ tests 3  ℹ pass 3  ℹ fail 0