2.2 常见 BRDF 的定点化

#!/usr/bin/env python3
# file: lantern/sim/cg_shade_02_2_brdf_fixed.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_shade_02_2_brdf_fixed.py
"""2.2 常见 BRDF 的定点化（S 层）。

在 2.1 的同一个 512 方向半球采样集上，把三种 BRDF 写成**全整数**形式，
并逐个实跑它们的**能量守恒残差**：

  - **Lambert**：f_r = ρ/π，1/π 的整数近似 INV_PI = round(ONE/π) = 20861。
  - **Phong**：归一化系数 (n+2)/(2π)，n=16 时 = 9/π = 9·INV_PI = 187749；
    指数 (R·V)^n 用**整数平方链**实现——n=16 = 2⁴，故 4 次平方
    c→c²→c⁴→c⁸→c¹⁶，每次都是 fmul（右移 16），无浮点、无 pow。
  - **Blinn-Phong**：归一化系数 (n+8)/(8π)，n=16 时 = 3/π = 3·INV_PI = 62583；
    半角向量 N·H = (1+cosθ)/√(2+2cosθ) = cos(θ/2)，整数里用 fsqrt 求分母。

入射方向取**法线**，故 Phong 的反射向量 = 法线、R·V = cosθ；
Blinn-Phong 的半角向量 N·H = cos(θ/2)。两种情形都被积函数各向同性，
8 个方位精确求和。

**参照值**用**闭式积分的有理数**给出（无求积），故「半球 512 方向估计」与
「闭式参照」之差就是**求积误差**（采样密度造成的），而**归一化误差**（BRDF
本身是否守恒）则看闭式参照与 ONE 的关系——两者能分开，结论才不含糊。

另给出**未归一化 Phong**（f = k_s·(R·V)^n）作为负例：它的能量远小于 1，
说明「不归一化就丢能量」。

产物：out/cg_shade_02_2_brdf.json。
"""

import json
import os
import sys

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

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))

ONE = cglib.ONE
INV_PI = base.INV_PI          # 20861
N_EXP = 16                    # Phong / Blinn-Phong 指数
PHONG_NORM = 9 * INV_PI       # (n+2)/(2π)·ONE，n=16 -> 9/π
BLINN_NORM = 3 * INV_PI       # (n+8)/(8π)·ONE，n=16 -> 3/π

def pow_chain(c, sq):
    """整数平方链：c^(2^sq)，每次一步 fmul。sq=4 -> c^16。"""
    for _ in range(sq):
        c = cglib.fmul(c, c)
    return c


def half_cos(cz):
    """Blinn-Phong 的 N·H = cos(θ/2) = (1+cosθ)/√(2+2cosθ)，整数。"""
    num = ONE + cz
    den = cglib.fsqrt(2 * num)
    return 0 if den == 0 else cglib.fdiv(num, den)


def closed_form_reference():
    """闭式参照（全整数有理数，无求积）：

    - Lambert：∫ (ρ/π)·cosθ dω = ρ = ONE。
    - Phong 归一化：(n+2)/(2π)·∫cosⁿθ cosθ dω = (n+2)/(2π)·2π/(n+2) = ONE。
    - Blinn 归一化：(n+8)/(8π)·2π·I_b，其中
      I_b = ∫₀^{π/2} cosⁿ(θ/2)·cosθ·sinθ dθ，换元 u = cosθ 得
      I_b = (1/256)∫₀¹ u(1+u)^8 du = (1/256)(1023/10 − 511/9) = 4097/23040。
      故 E = (n+8)/4·I_b = 6·4097/23040 = 4097/3840。
    - Phong 未归一化：∫cosⁿθ cosθ dω = 2π/(n+2) = π/9。

    返回四个 Q16.16 参照值。
    """
    e_blinn = (4097 * 6 << 16) // 23040          # 24582/23040 的 Q16.16
    e_raw = base.TWO_PI // (N_EXP + 2)           # π/9 的 Q16.16
    return {"lambert": ONE, "phong_normalized": ONE,
            "blinn_normalized": e_blinn, "phong_raw": e_raw}


def main():
    samples = base.build_hemisphere()

    lambert = lambda dx, dy, cz: base.lambert(ONE)                       # ρ/π
    phong = lambda dx, dy, cz: cglib.fmul(PHONG_NORM, pow_chain(cz, 4))
    blinn = lambda dx, dy, cz: cglib.fmul(BLINN_NORM, pow_chain(half_cos(cz), 4))
    phong_raw = lambda dx, dy, cz: pow_chain(cz, 4)                      # 负例：未归一化

    energy = {
        "lambert": base.hemisphere_integral(samples, lambert),
        "phong_normalized": base.hemisphere_integral(samples, phong),
        "blinn_normalized": base.hemisphere_integral(samples, blinn),
        "phong_raw": base.hemisphere_integral(samples, phong_raw),
    }

    analytic = closed_form_reference()

    c = samples[32 * base.N_PHI][2]  # 取中间一带（k=32，cosθ≈0.51）展示平方链
    chain = [c]
    for _ in range(4):
        chain.append(cglib.fmul(chain[-1], chain[-1]))

    record = {
        "constants": {
            "inv_pi": INV_PI,
            "phong_norm_q16": PHONG_NORM,
            "blinn_norm_q16": BLINN_NORM,
            "phong_norm_expr": "(n+2)/(2π) = 9/π，n=16",
            "blinn_norm_expr": "(n+8)/(8π) = 3/π，n=16",
        },
        "pow_chain": {
            "n": N_EXP, "squarings": 4,
            "c0_q16": chain[0], "c2_q16": chain[1], "c4_q16": chain[2],
            "c8_q16": chain[3], "c16_q16": chain[4],
        },
        "closed_form": {
            "blinn_fraction": "4097/3840",
            "blinn_q16": analytic["blinn_normalized"],
            "phong_raw_q16": analytic["phong_raw"],
        },
        "energy_q16": energy,
        "analytic_q16": analytic,
        "residual_q16": {k: ONE - v for k, v in energy.items()},
        "quadrature_gap_q16": {k: analytic[k] - energy[k] for k in energy},
        "permille_of_one": {k: round(v * 1000.0 / ONE, 4) for k, v in energy.items()},
        "analytic_permille": {k: round(v * 1000.0 / ONE, 4) for k, v in analytic.items()},
        "assert": {
            "lambert_le_one": energy["lambert"] <= ONE,
            "phong_normalized_le_one": energy["phong_normalized"] <= ONE,
            "blinn_normalized_over_one": energy["blinn_normalized"] > ONE,
            "phong_raw_loses_energy": energy["phong_raw"] < energy["phong_normalized"] // 2,
        },
    }

    out = os.path.join(ROOT, "out")
    os.makedirs(out, exist_ok=True)
    with open(os.path.join(out, "cg_shade_02_2_brdf.json"), "w", encoding="utf-8") as fh:
        json.dump(record, fh, ensure_ascii=False, indent=2)

    print(f"1/π / Phong 归一化 / Blinn 归一化 : {INV_PI} / {PHONG_NORM} / {BLINN_NORM}")
    print(f"平方链 c^16 (4 次平方)           : {chain[0]} -> {chain[4]}")
    for name in energy:
        print(f"{name:18s} 估计 {energy[name]:6d} ({record['permille_of_one'][name]:8.4f}‰)  "
              f"参照 {analytic[name]:6d} ({record['analytic_permille'][name]:8.4f}‰)  "
              f"求积差 {record['quadrature_gap_q16'][name]:5d}")
    print(f"断言全过                         : {all(record['assert'].values())}")
    return 0


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

1/π / Phong 归一化 / Blinn 归一化 : 20861 / 187749 / 62583
平方链 c^16 (4 次平方)           : 33280 -> 1
lambert            估计  65533 (999.9542‰)  参照  65536 (1000.0000‰)  求积差     3
phong_normalized   估计  65321 (996.7194‰)  参照  65536 (1000.0000‰)  求积差   215
blinn_normalized   估计  69886 (1066.3757‰)  参照  69922 (1066.9250‰)  求积差    36
phong_raw          估计  22800 (347.9004‰)  参照  22876 (349.0601‰)  求积差    76
断言全过                         : True

// file: lantern/core/test/cg_shade_02_2_brdf_fixed.test.ts
// tier: P
// verify: node --test 'core/test/cg_shade_02_2_brdf_fixed.test.ts'
// 2.2 的 P 层：TypeScript 侧重算三种 BRDF 的定点能量，断言与 S 侧
// out/cg_shade_02_2_brdf.json 逐位相同；Blinn 归一化超 1 是**真实结论**，不是误差。

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 { ONE, fmul } from "../src/lantern/fixed.ts";
import {
  BLINN_NORM, INV_PI, PHONG_NORM, buildHemisphere, halfCos, hemisphereIntegral, powChain,
} from "../src/lantern/cg_shade_02_material.ts";

const root = join(dirname(fileURLToPath(import.meta.url)), "..", "..");
const g = JSON.parse(
  readFileSync(join(root, "out", "cg_shade_02_2_brdf.json"), "utf8"),
) as {
  constants: { inv_pi: number; phong_norm_q16: number; blinn_norm_q16: number };
  pow_chain: { c0_q16: number; c16_q16: number };
  energy_q16: Record<string, number>;
  analytic_q16: Record<string, number>;
  assert: Record<string, boolean>;
};

const samples = buildHemisphere();

test("归一化常量与 S 侧逐位相同（1/π、9/π、3/π）", () => {
  assert.equal(INV_PI, BigInt(g.constants.inv_pi));
  assert.equal(PHONG_NORM, BigInt(g.constants.phong_norm_q16));
  assert.equal(BLINN_NORM, BigInt(g.constants.blinn_norm_q16));
  assert.equal(INV_PI, 20861n);
  assert.equal(PHONG_NORM, 187749n);
  assert.equal(BLINN_NORM, 62583n);
});

test("平方链 c^16：4 次平方与 S 侧逐位相同", () => {
  const c0 = samples[32 * 8][2];
  assert.equal(c0, BigInt(g.pow_chain.c0_q16));
  assert.equal(powChain(c0, 4), BigInt(g.pow_chain.c16_q16));
});

test("【核心】三种 BRDF 的能量估计与 S 侧逐位相同", () => {
  const energy = {
    lambert: hemisphereIntegral(samples, () => fmul(ONE, INV_PI)),
    phong_normalized: hemisphereIntegral(samples, (_x, _y, cz) => fmul(PHONG_NORM, powChain(cz, 4))),
    blinn_normalized: hemisphereIntegral(samples, (_x, _y, cz) => fmul(BLINN_NORM, powChain(halfCos(cz), 4))),
    phong_raw: hemisphereIntegral(samples, (_x, _y, cz) => powChain(cz, 4)),
  };
  for (const k of Object.keys(energy)) {
    assert.equal(energy[k as keyof typeof energy], BigInt(g.energy_q16[k]), `${k} 与 S 侧不符`);
  }
});

test("【核心】Lambert 与归一化 Phong ≤ ONE；归一化 Blinn **超过** ONE", () => {
  assert.ok(BigInt(g.energy_q16.lambert) <= ONE, "Lambert 必须 ≤ ONE");
  assert.ok(BigInt(g.energy_q16.phong_normalized) <= ONE, "归一化 Phong 必须 ≤ ONE");
  assert.ok(BigInt(g.energy_q16.blinn_normalized) > ONE, "归一化 Blinn 的能量必须 > ONE（真实结论）");
  assert.ok(g.assert.blinn_normalized_over_one);
});

test("负例：未归一化 Phong 严重丢能量（不到归一化的一半）", () => {
  assert.ok(
    BigInt(g.energy_q16.phong_raw) < BigInt(g.energy_q16.phong_normalized) / 2n,
    "未归一化 Phong 的能量必须远小于归一化版本",
  );
});

✔ 归一化常量与 S 侧逐位相同（1/π、9/π、3/π）
✔ 平方链 c^16：4 次平方与 S 侧逐位相同
✔ 【核心】三种 BRDF 的能量估计与 S 侧逐位相同
✔ 【核心】Lambert 与归一化 Phong ≤ ONE；归一化 Blinn **超过** ONE
✔ 负例：未归一化 Phong 严重丢能量（不到归一化的一半）
ℹ tests 5  ℹ pass 5  ℹ fail 0