6.1 微表面理论

# snippet: lantern/sim/cg_shade_06_pbr.py 的 ggx_d 主体
def ggx_d(a2, nh):
    t = fmul(fmul(nh, nh), a2 - ONE) + ONE
    return fdiv(fdiv(a2, t), fmul(t, PI))

#!/usr/bin/env python3
# file: lantern/sim/cg_shade_06_pbr.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_shade_06_pbr.py
"""6.1 微表面理论：GGX 分布 D、Smith 几何项 G、Schlick 菲涅尔 F 的整数实现。

全部 Q16.16 定点（Python int 任意精度），只用 fmul / fdiv / fsqrt——没有三角函数、
没有浮点常量（π 取整数 205887）。两个实验：

  1. D 的归一化 ∫ D(h)(n·h) dω_h ≈ 1：在**投影圆盘**上均匀取样，度量天然正确
     （余弦因子与 1/cosθ 相消：∫ f cosθ dω = ∫∫_{圆盘} f dr ds）。
  2. 镜面波瓣 f_s = D G F / (4 (n·l)(n·v)) 随入射角的扫描与半高全宽。

产物：out/cg_shade_06_1_lobe.csv、out/cg_shade_06_1_lobe.svg、out/cg_shade_06_1_microfacet.json。
"""

import json
import os
import sys

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib  # noqa: E402
from cg_shade_06_util import asin_deg, svg  # noqa: E402
from cglib import ONE, fdiv, fmul, fsqrt  # noqa: E402

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
PI = 205887            # round(pi * 2^16)
GRID = 64              # 投影圆盘每边的格数
A2, F0 = 4096, 2621    # 粗糙度 0.5 -> alpha = 0.25、alpha^2 = 0.0625；绝缘体 F0 = 0.04


def ggx_d(a2, nh):
    """GGX / Trowbridge-Reitz：D = a2 / (pi t^2)，t = (n·h)^2 (a2-1) + 1。

    拆成两次 fdiv（a2/t 再 /(pi t)）而不是一次除以 pi*t^2：t >= a2，两次商都不过界；
    若先算 t^2，小 alpha 下 (t*t)>>16 会下溢成 0，fdiv 直接除零。
    """
    t = fmul(fmul(nh, nh), a2 - ONE) + ONE
    return fdiv(fdiv(a2, t), fmul(t, PI))


def smith_g1(a2, nu):
    """Smith-GGX 遮蔽：G1(u) = 2 nu / (nu + sqrt(a2 + (1-a2) nu^2))。"""
    return fdiv(2 * nu, nu + fsqrt(a2 + fmul(ONE - a2, fmul(nu, nu))))


def schlick_f(f0, vh):
    """Schlick 菲涅尔：F = f0 + (1-f0)(1 - v·h)^5。"""
    c = ONE - vh
    c = 0 if c < 0 else c
    c2 = fmul(c, c)
    return f0 + fmul(ONE - f0, fmul(fmul(c2, c2), c))


def spec_brdf(a2, nl, nv, nh, vh, f0):
    """镜面项 f_s = D G F / (4 (n·l)(n·v))。"""
    g = fmul(smith_g1(a2, nv), smith_g1(a2, nl))
    return fdiv(fmul(fmul(ggx_d(a2, nh), g), schlick_f(f0, vh)), 4 * fmul(nl, nv))


def disk_dirs(n=GRID):
    """投影圆盘上的确定性方向集：l = (r, s, sqrt(1 - r^2 - s^2))，r、s 取均匀格。"""
    step, half = 2 * ONE // n, ONE // n
    out = []
    for i in range(n):
        r = -ONE + half + i * step
        for j in range(n):
            s = -ONE + half + j * step
            q = ONE - fmul(r, r) - fmul(s, s)
            if q >= 0:
                out.append((r, s, fsqrt(q)))
    return out


def cell_of(n=GRID):
    """圆盘上一个格子的面积（Q16.16）；n 个格子铺满单位圆盘时 Σ = pi。"""
    return fmul(2 * ONE // n, 2 * ONE // n)


def d_norm(a2, dirs, n=GRID):
    """∫ D(h)(n·h) dω_h ≈ cell * Σ D(h)：h 在圆盘上均匀取样，n·h 就是 h 的 z 分量。"""
    return (sum(ggx_d(a2, h[2]) for h in dirs) * cell_of(n)) >> 16


def lobe(a2, f0, rv, n=128):
    """平面内波瓣：v 固定在 (rv, 0, sqrt(1-rv^2))，l 扫过整个平面，r 从 -1 到 1。

    镜面方向 l = reflect(v, n) = (-rv, 0, zv) 在法线的**另一侧**，那里 h = n、n·h = 1、
    D 取峰；但 f_s 还带 1/(4 (n·l)(n·v))，故 BRDF 的峰比镜面方向略偏掠射。
    """
    zv = fsqrt(ONE - fmul(rv, rv))
    step = 2 * ONE // n
    rows = []
    for k in range(n + 1):
        r = min(max(-ONE + k * step, 1 - ONE), ONE - 1)   # 避开 n·l = 0 的掠射端点
        z = fsqrt(ONE - fmul(r, r))
        hx, hz = r + rv, z + zv
        hl = fsqrt(fmul(hx, hx) + fmul(hz, hz))
        rows.append((k, r, z, spec_brdf(a2, z, zv, fdiv(hz, hl),
                                       fdiv(fmul(rv, hx) + fmul(zv, hz), hl), f0)))
    return rows


def main():
    out = os.path.join(ROOT, "out")
    os.makedirs(out, exist_ok=True)
    dirs = disk_dirs()
    table = [(lab, d_norm(a2, dirs)) for lab, a2 in
             [("0.05", 164), ("0.25", 4096), ("0.50", 16384), ("1.00", 65536)]]
    assert all(abs(n - ONE) < 4000 for _, n in table), "D 未归一化"

    rows = lobe(A2, F0, ONE // 2)
    ip = max(range(len(rows)), key=lambda i: rows[i][3])
    half, lo, hi = rows[ip][3] // 2, ip, ip
    while lo > 0 and rows[lo - 1][3] >= half:
        lo -= 1
    while hi < len(rows) - 1 and rows[hi + 1][3] >= half:
        hi += 1
    deg = asin_deg(rows[hi][1]) - asin_deg(rows[lo][1])
    assert ggx_d(A2, ONE) > ggx_d(A2, ONE - 1024) > ggx_d(A2, ONE - 4096), "D 未在 n·h = 1 取峰"

    cglib.write_csv(os.path.join(out, "cg_shade_06_1_lobe.csv"),
                    ["k", "r_q16", "z_q16", "fs_q16"], rows)
    size = svg(os.path.join(out, "cg_shade_06_1_lobe.svg"),
               [("f_s(r)  峰值 %.4f" % (rows[ip][3] / ONE), "#ff8a5c", max(t[3] for t in rows),
                 [((t[1] + ONE) // 2, t[3]) for t in rows])],
               "6.1  镜面波瓣 f_s（v 在 30 度，alpha = 0.25）")
    data = {"alpha2": A2, "f0": F0, "grid": GRID, "cell_q16": cell_of(GRID), "directions": len(dirs),
            "d_g1_f_at_nh1": {"d_q16": ggx_d(A2, ONE), "g1_q16": smith_g1(A2, ONE),
                              "f_q16": schlick_f(F0, ONE)},
            "d_norm": [{"alpha": lab, "value_q16": n} for lab, n in table],
            "peak": {"r_q16": rows[ip][1], "fs_q16": rows[ip][3]},
            "fwhm_samples": hi - lo + 1, "fwhm_deg_q16": deg, "svg_bytes": size}
    with open(os.path.join(out, "cg_shade_06_1_microfacet.json"), "w", encoding="utf-8") as fh:
        json.dump(data, fh, ensure_ascii=False, indent=2)

    print(f"投影圆盘取样   : {len(dirs)} 个方向（{GRID}x{GRID} 格，格面积 {cell_of(GRID)}/65536）")
    print("D 归一化       : " + "  ".join(f"a={lab}: {n / ONE:.4f}" for lab, n in table))
    print(f"波瓣峰值       : r={rows[ip][1] / ONE:.4f}（sin）  f_s={rows[ip][3] / ONE:.4f}")
    print(f"半高全宽       : {hi - lo + 1} / {len(rows)} 个采样 = {deg / ONE:.2f} 度")
    print(f"D/G1/F(n·h=1)  : {ggx_d(A2, ONE) / ONE:.4f} / {smith_g1(A2, ONE) / ONE:.4f} / "
          f"{schlick_f(F0, ONE) / ONE:.4f}")
    print(f"svg            : out/cg_shade_06_1_lobe.svg ({size} bytes)")
    return 0


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

投影圆盘取样   : 3228 个方向（64x64 格，格面积 64/65536）
D 归一化       : a=0.05: 0.9927  a=0.25: 1.0000  a=0.50: 1.0008  a=1.00: 1.0034
波瓣峰值       : r=-0.5312（sin）  f_s=0.0680
半高全宽       : 38 / 129 个采样 = 40.74 度
D/G1/F(n·h=1)  : 5.0933 / 1.0000 / 0.0400
svg            : out/cg_shade_06_1_lobe.svg (1788 bytes)

// file: lantern/core/src/lantern/cg_shade_06_pbr.ts
// tier: P
// verify: node --test 'core/test/cg_shade_06_*.test.ts'
// 6.1 微表面理论的 TypeScript 侧：GGX 分布 D、Smith 几何项 G1、Schlick 菲涅尔 F、
// 镜面项 f_s，以及投影圆盘上的确定性方向集。与 sim/cg_shade_06_pbr.py 逐字对齐。
//
// 为什么必须 BigInt：小 alpha 下 D 在 Q16.16 里能到 5.3e7，乘上格面积再累加 3228 项
// 会到 1e12 量级；Number 的 `>>` 只对 32 位有效，超了会静默截断成错误结果。
// 舍入口径与 Python 侧一致：降精度用 `>>`（向下取整），除法一律走 fdiv（向零截断）。

import { ONE, fdiv, fmul, fsqrt } from "./fixed.ts";

export const PI = 205887n; // round(pi * 2^16)
export const GRID = 64; // 投影圆盘每边的格数

export type Dir = [bigint, bigint, bigint];

/** GGX / Trowbridge-Reitz：D = a2 / (pi t^2)，t = (n·h)^2 (a2-1) + 1。 */
export function ggxD(a2: bigint, nh: bigint): bigint {
  const t = fmul(fmul(nh, nh), a2 - ONE) + ONE;
  return fdiv(fdiv(a2, t), fmul(t, PI));
}

/** Smith-GGX 遮蔽：G1(u) = 2 nu / (nu + sqrt(a2 + (1-a2) nu^2))。 */
export function smithG1(a2: bigint, nu: bigint): bigint {
  return fdiv(2n * nu, nu + fsqrt(a2 + fmul(ONE - a2, fmul(nu, nu))));
}

/** Schlick 菲涅尔：F = f0 + (1-f0)(1 - v·h)^5。 */
export function schlickF(f0: bigint, vh: bigint): bigint {
  const raw = ONE - vh;
  const c = raw < 0n ? 0n : raw;
  const c2 = fmul(c, c);
  return f0 + fmul(ONE - f0, fmul(fmul(c2, c2), c));
}

/** 镜面项 f_s = D G F / (4 (n·l)(n·v))。 */
export function specBrdf(
  a2: bigint, nl: bigint, nv: bigint, nh: bigint, vh: bigint, f0: bigint,
): bigint {
  const g = fmul(smithG1(a2, nv), smithG1(a2, nl));
  return fdiv(fmul(fmul(ggxD(a2, nh), g), schlickF(f0, vh)), 4n * fmul(nl, nv));
}

/** 投影圆盘上的确定性方向集：l = (r, s, sqrt(1 - r^2 - s^2))，r、s 取均匀格。 */
export function diskDirs(n: number = GRID): Dir[] {
  const step = (2n * ONE) / BigInt(n);
  const half = ONE / BigInt(n);
  const out: Dir[] = [];
  for (let i = 0; i < n; i += 1) {
    const r = -ONE + half + BigInt(i) * step;
    for (let j = 0; j < n; j += 1) {
      const s = -ONE + half + BigInt(j) * step;
      const q = ONE - fmul(r, r) - fmul(s, s);
      if (q >= 0n) out.push([r, s, fsqrt(q)]);
    }
  }
  return out;
}

/** 圆盘上一个格子的面积（Q16.16）；n 个格子铺满单位圆盘时总和 = pi。 */
export function cellOf(n: number = GRID): bigint {
  return fmul((2n * ONE) / BigInt(n), (2n * ONE) / BigInt(n));
}

/** ∫ D(h)(n·h) dω_h ≈ cell * Σ D(h)：圆盘均匀取样，n·h 就是 h 的 z 分量。 */
export function dNorm(a2: bigint, dirs: Dir[], n: number = GRID): bigint {
  let sum = 0n;
  for (const h of dirs) sum += ggxD(a2, h[2]);
  return (sum * cellOf(n)) >> 16n;
}