5618 字

28 分钟

计算机图形学笔记(四)：光线追踪与全局光照

2025-07-23

计算机图形学

C++

/

编程

/

计算机图形学

/

OpenGL

NOTE
这是全系列最硬核也最美的一部分。光栅化追求”快得能实时交互”，光线追踪追求”正确到能当参考图”。本章作为渲染方程 / 辐射度量学 / BRDF / 蒙特卡洛积分 / 路径追踪的 canonical 位置——Part 2 光栅化里写到的光照模型、Part 6 里的 PBR/RTX/NeRF，都回链到这里。

目录#

辐射度量学基础 — 立体角 / 辐射通量 / 辐射度 / 辐照度（整个体系的物理单位）
渲染方程 — 积分形式 / 递归形式 / Neumann 级数 / 光传输算子
BRDF 与材质模型 — 定义 / 性质 / Lambert / Phong / Cook-Torrance / GGX
光线-几何体相交 — Ray 参数化 / 球 / 三角形（Möller-Trumbore）/ AABB Slab
空间加速结构 — BVH（canonical）/ SAH / KD-Tree / 八叉树对比
蒙特卡洛积分 — 估计量方差 / 重要性采样 / 分层 / MIS
路径追踪 — Whitted / 基础路径追踪 / 俄罗斯轮盘赌 / 直接+间接分离
高级技术概述 — 双向路径追踪 / 光子映射 / MLT（Part 6 前导）

十六、辐射度量学基础#

WARNING
为什么先讲辐射度量学：渲染方程里每一个符号都有严格的物理单位。Part 2 的 Blinn-Phong 用 “光强 × cos” 的朴素模型可以写出看起来对的图，但要做 PBR、对比真实世界、做反演（逆渲染），就必须回到辐射度量学的严格定义。

16.1 立体角#

16.1.1 定义#

立体角是平面角在三维空间的推广：一个锥体在单位球上截下的面积。

\Omega = \frac{A}{r^2} \quad \text{[sr, 球面度]}

整个球面的立体角是 $4\pi$ sr
半球是 $2\pi$ sr

16.1.2 球面坐标下的微分立体角#

用 $(\theta, \phi)$ 参数化（ $\theta$ 为极角/与法线夹角， $\phi$ 为方位角）：

d\omega = \sin\theta \, d\theta \, d\phi

NOTE
这个 $\sin\theta$ 因子很关键：接近两极时相同的 $d\theta \, d\phi$ 对应的球面面积缩小。后面半球积分 $\int_\Omega (\cdot) \, d\omega = \int_0^{2\pi} \int_0^{\pi/2} (\cdot) \sin\theta \, d\theta \, d\phi$ 都要记得这一点。

16.1.3 验证半球面积#

\int_0^{2\pi}\int_0^{\pi/2} \sin\theta \, d\theta \, d\phi = 2\pi \cdot [-\cos\theta]_0^{\pi/2} = 2\pi \quad \checkmark

16.2 辐射度量四件套#

量	符号/单位	定义	直觉
辐射通量	$\Phi$ [W]	单位时间通过表面/区域的能量	”每秒多少焦耳的光”
辐射强度	$I = d\Phi/d\omega$ [W/sr]	单位立体角上的辐射通量	点光源的各向强度分布
辐照度	$E = d\Phi/dA$ [W/m²]	表面单位面积接收到的辐射通量	墙上”每平米多亮”，与观察方向无关
辐射度	$L$ [W/(m²·sr)]	每单位投影面积、单位立体角的辐射通量	”某方向上看某点的亮度”，相机直接感知这个量

16.2.1 辐射度的严格定义#

L(\mathbf{p}, \omega) = \frac{d^2 \Phi}{dA \cos\theta \, d\omega}

其中 $\theta$ 是 $\omega$ 与表面法线的夹角。投影面积 $dA \cos\theta$ 是辐射度比辐照度更基本的原因——辐射度沿真空传播不变（能量守恒的简洁表达）。

16.2.2 辐照度与辐射度的关系#

E(\mathbf{p}) = \int_\Omega L_i(\mathbf{p}, \omega_i) \cos\theta_i \, d\omega_i

这就是半球辐照度积分，是渲染方程的右侧的雏形。

16.3 Lambert 余弦定律#

对一个面元接收平行光，单位法向面积上的功率 = 入射功率 × $\cos\theta_i$ 。这就是为什么渲染方程里到处有 $\cos\theta_i$ 。

十七、渲染方程（canonical）#

WARNING
Kajiya 1986 的渲染方程是整个物理渲染的根基。这一章是全系列渲染方程的规范位置，后面 Part 6 的 PBR、RTX、NeRF 都会回链这里。

17.1 渲染方程的积分形式#

在表面点 $\mathbf{p}$ 、出射方向 $\omega_o$ 上，出射辐射度满足：

L_o(\mathbf{p}, \omega_o) = L_e(\mathbf{p}, \omega_o) + \int_\Omega f_r(\mathbf{p}, \omega_i, \omega_o) \, L_i(\mathbf{p}, \omega_i) \cos\theta_i \, d\omega_i

各项含义：

$L_o$ :出射辐射度（我们要算的目标）
$L_e$ :自发光项（光源的内禀辐射）
$f_r$ ，下一章详述
$L_i$ :入射辐射度
$\cos\theta_i = \vec{n} \cdot \omega_i$ 余弦
$\Omega$ :上半球，积分微元 $d\omega_i = \sin\theta_i \, d\theta_i \, d\phi_i$

17.2 递归形式#

$L_i(\mathbf{p}, \omega_i)$ 本身就是从 $-\omega_i$ 方向射向 $\mathbf{p}$ 的光——如果场景中沿 $-\omega_i$ 走遇到点 $\mathbf{p}'$ ，那么：

L_i(\mathbf{p}, \omega_i) = L_o(\mathbf{p}', -\omega_i)

代入原方程得到递归形式:

L_o(\mathbf{p}, \omega_o) = L_e(\mathbf{p}, \omega_o) + \int_\Omega f_r(\mathbf{p}, \omega_i, \omega_o) \, L_o(\mathbf{p}', -\omega_i) \cos\theta_i \, d\omega_i

这是一个第二类 Fredholm 积分方程——未知函数 $L_o$ 同时出现在等号两侧。

17.3 Neumann 级数展开#

定义光传输算子 $T$ :

(TL)(\mathbf{p}, \omega_o) = \int_\Omega f_r(\mathbf{p}, \omega_i, \omega_o) L(\mathbf{p}', -\omega_i) \cos\theta_i \, d\omega_i

渲染方程写成算子形式： $L = L_e + TL$ 。形式上求解：

L = (I - T)^{-1} L_e = \sum_{n=0}^{\infty} T^n L_e = L_e + T L_e + T^2 L_e + T^3 L_e + \cdots

📌 各项的物理意义#

$L_e$ :0 次弹射——直接看到光源本身
$T L_e$ :1 次弹射——光源发光 → 经 1 次反射到眼睛（Whitted 直接光照模型止步于此）
$T^2 L_e$ :2 次弹射——间接光照的第一层
$T^n L_e$ :n 次弹射——完整全局光照

路径追踪本质上就是用蒙特卡洛方法对这个无穷级数做无偏估计。

17.4 收敛性#

谱半径条件: $\rho(T) < 1$ 时级数收敛。

物理上对应能量守恒:每次反射能量损失一部分（BRDF 半球积分 $\leq 1$ ），无穷次弹射总能量有限。

TIP
工程上用俄罗斯轮盘赌把无穷级数无偏地截断成有限路径（见 §22.3）。

十八、BRDF 与材质模型#

18.1 BRDF 的定义#

双向反射分布函数（Bidirectional Reflectance Distribution Function）：

f_r(\mathbf{p}, \omega_i, \omega_o) = \frac{dL_o(\mathbf{p}, \omega_o)}{dE_i(\mathbf{p}, \omega_i)} = \frac{dL_o(\mathbf{p}, \omega_o)}{L_i(\mathbf{p}, \omega_i) \cos\theta_i \, d\omega_i} \quad \text{[sr}^{-1}\text{]}

直觉：单位入射辐照度带来的单位出射辐射度。注意它的单位是 $\text{sr}^{-1}$ ，不是无量纲。

18.2 BRDF 必须满足的三条性质#

性质 1 — 非负性: $f_r \geq 0$

性质 2 — Helmholtz 互易（时间反演对称）：

f_r(\omega_i, \omega_o) = f_r(\omega_o, \omega_i)

性质 3 — 能量守恒:对任意 $\omega_o$

\int_\Omega f_r(\omega_i, \omega_o) \cos\theta_i \, d\omega_i \leq 1

WARNING
Part 2 的 Blinn-Phong 模型不满足能量守恒（没归一化），PBR 里用 Blinn-Phong 前要乘归一化因子 $(n+8)/(8\pi)$ 才合法。

18.3 Lambert 漫反射（完美各向同性散射）#

定义:

f_r^{\text{Lambert}} = \frac{\rho_d}{\pi}, \qquad \rho_d \in [0, 1]^3

$\rho_d$ 称为漫反射率（albedo）。

📌 为什么分母是 $\pi$ （完整推导）#

要求能量守恒（取等号）：

\int_\Omega \frac{\rho_d}{\pi} \cos\theta_i \, d\omega_i = \frac{\rho_d}{\pi} \int_0^{2\pi} \int_0^{\pi/2} \cos\theta \sin\theta \, d\theta \, d\phi = \frac{\rho_d}{\pi} \cdot 2\pi \cdot \tfrac{1}{2} = \rho_d \leq 1

即 $\rho_d = 1$ 时是完美白漫反射——入多少出多少，无方向偏好。

18.4 Phong / Blinn-Phong（经验模型）#

f_r^{\text{Phong}} = \frac{\rho_d}{\pi} + \rho_s \frac{(n+2)}{2\pi} (\vec{r} \cdot \omega_o)^n

$\vec{r}$ 是入射方向关于法线的完美镜面反射方向。归一化因子 $(n+2)/(2\pi)$ 是满足能量守恒的必要条件（Part 2 §10 里省略了，这里补上）。

Blinn-Phong 用半向量 $\vec{h} = (\omega_i + \omega_o)/\|\omega_i + \omega_o\|$ 取代 $\vec{r}$ ,归一化因子为 $(n+8)/(8\pi)$ 。

18.5 Cook-Torrance（物理微表面模型）#

现代 PBR 的主流选择。假设表面由大量微小镜面（microfacet）组成：

f_r = \frac{D(\vec{h}) \, F(\omega_i, \vec{h}) \, G(\omega_i, \omega_o, \vec{h})}{4 \, (\omega_i \cdot \vec{n})(\omega_o \cdot \vec{n})}

三个分量各司其职：

法线分布 $D$ :微表面法线与宏观法线夹角的分布（GGX/Beckmann）
菲涅尔 $F$ :镜面反射率随入射角的变化（Schlick 近似）
几何项 $G$ :自阴影/自遮蔽（Smith 模型）

18.5.1 GGX（Trowbridge-Reitz）分布#

D_{\text{GGX}}(\vec{h}) = \frac{\alpha^2}{\pi \left[(\vec{n}\cdot\vec{h})^2 (\alpha^2 - 1) + 1\right]^2}

$\alpha = \text{roughness}^2$ 。GGX 相比 Beckmann 有更长的高光尾巴,匹配真实材质的观测。

18.5.2 Schlick 菲涅尔近似#

F(\theta) = F_0 + (1 - F_0)(1 - \cos\theta)^5

$F_0$ 是正入射时的反射率；金属用 RGB 三通道值，非金属典型 $F_0 \approx 0.04$ 。

18.5.3 Smith 几何项#

G(\omega_i, \omega_o) = G_1(\omega_i) \cdot G_1(\omega_o), \quad G_1(\omega) = \frac{2(\vec{n}\cdot\omega)}{(\vec{n}\cdot\omega) + \sqrt{\alpha^2 + (1-\alpha^2)(\vec{n}\cdot\omega)^2}}

Part 6 会给完整实现，这里先有个概念即可。

十九、光线-几何体相交#

19.1 光线参数化#

\mathbf{r}(t) = \mathbf{o} + t \vec{d}, \quad t \in [t_{\min}, t_{\max}]

$\mathbf{o}$ :起点
$\vec{d}$ :方向（通常取单位向量）
$t_{\min} > 0$ :避免自相交（典型 $10^{-4}$ ）

1
struct Ray {
2
    Eigen::Vector3f origin;
3
    Eigen::Vector3f direction;   // 单位向量
4
    float t_min = 1e-4f;
5
    float t_max = std::numeric_limits<float>::infinity();
6

7
    Eigen::Vector3f at(float t) const { return origin + t * direction; }
8
};

19.2 光线-球体求交#

球面方程 $\|\mathbf{p} - \mathbf{c}\|^2 = r^2$ ,代入 $\mathbf{p} = \mathbf{o} + t\vec{d}$ 得：

t^2 (\vec{d}\cdot\vec{d}) + 2t \, (\mathbf{o} - \mathbf{c})\cdot\vec{d} + (\mathbf{o} - \mathbf{c})\cdot(\mathbf{o} - \mathbf{c}) - r^2 = 0

取 $\vec{d}$ 为单位向量时 $a = 1$ ,判别式 $\Delta = b^2 - 4c$ （其中 $b = 2\vec{d}\cdot(\mathbf{o}-\mathbf{c})$ , $c = \|\mathbf{o}-\mathbf{c}\|^2 - r^2$ ）。

1
bool intersect_sphere(const Ray& r, const Eigen::Vector3f& c, float radius, float& t_hit) {
2
    Eigen::Vector3f oc = r.origin - c;
3
    float b = 2.0f * oc.dot(r.direction);
4
    float c_term = oc.squaredNorm() - radius * radius;
5
    float disc = b * b - 4.0f * c_term;
6
    if (disc < 0) return false;
7
    float sqrt_d = std::sqrt(disc);
8
    float t0 = (-b - sqrt_d) * 0.5f;
9
    float t1 = (-b + sqrt_d) * 0.5f;
10
    if (t0 > t1) std::swap(t0, t1);
11
    t_hit = (t0 >= r.t_min) ? t0 : t1;
12
    return t_hit >= r.t_min && t_hit <= r.t_max;
13
}

WARNING
当 $\mathbf{o}$ 远离球心、 $r$ 很小时 $b^2 \approx 4c$ ,裸减法会灾难性消去。工业实现（PBRT）用 $\mathbf{o} - (\mathbf{o}\cdot\vec{d})\vec{d}$ 先投影再算垂直距离，数值更稳。

19.3 光线-三角形：Möller-Trumbore#

核心思想：把交点用重心坐标表示（Part 2 §8.3 的定义回顾）， $(u, v)$ 既是求解量也是插值坐标，一次解完。

设三角形 $\mathbf{V}_0, \mathbf{V}_1, \mathbf{V}_2$ ,交点满足：

\mathbf{o} + t\vec{d} = (1-u-v)\mathbf{V}_0 + u \mathbf{V}_1 + v \mathbf{V}_2

移项：

\underbrace{[-\vec{d}, \mathbf{E}_1, \mathbf{E}_2]}_{3\times 3} \begin{pmatrix} t \\ u \\ v \end{pmatrix} = \mathbf{o} - \mathbf{V}_0

其中 $\mathbf{E}_1 = \mathbf{V}_1 - \mathbf{V}_0, \mathbf{E}_2 = \mathbf{V}_2 - \mathbf{V}_0$ 。

用 Cramer 法则 + 标量三重积恒等式 $(\vec{a}\times\vec{b})\cdot\vec{c} = (\vec{b}\times\vec{c})\cdot\vec{a}$ ,可不显式求逆：

t = \frac{(\mathbf{S}\times\mathbf{E}_1)\cdot\mathbf{E}_2}{\mathbf{P}\cdot\mathbf{E}_1}, \quad u = \frac{\mathbf{P}\cdot\mathbf{S}}{\mathbf{P}\cdot\mathbf{E}_1}, \quad v = \frac{(\mathbf{S}\times\mathbf{E}_1)\cdot\vec{d}}{\mathbf{P}\cdot\mathbf{E}_1}

其中 $\mathbf{S} = \mathbf{o} - \mathbf{V}_0, \mathbf{P} = \vec{d}\times\mathbf{E}_2$ 。

1
bool moller_trumbore(const Ray& r,
2
                     const Eigen::Vector3f& v0,
3
                     const Eigen::Vector3f& v1,
4
                     const Eigen::Vector3f& v2,
5
                     float& t, float& u, float& v) {
6
    const float EPS = 1e-8f;
7
    Eigen::Vector3f e1 = v1 - v0;
8
    Eigen::Vector3f e2 = v2 - v0;
9
    Eigen::Vector3f p  = r.direction.cross(e2);
10
    float det = p.dot(e1);
11
    if (std::abs(det) < EPS) return false;        // 光线平行于三角形
12
    float inv_det = 1.0f / det;
13
    Eigen::Vector3f s = r.origin - v0;
14
    u = p.dot(s) * inv_det;
15
    if (u < 0 || u > 1) return false;
16
    Eigen::Vector3f q = s.cross(e1);
17
    v = q.dot(r.direction) * inv_det;
18
    if (v < 0 || u + v > 1) return false;
19
    t = q.dot(e2) * inv_det;
20
    return t >= r.t_min && t <= r.t_max;
21
}

TIP
得到 $(u, v)$ 后，可以直接做 Part 2 §8.3 的重心坐标插值——法线、纹理坐标、顶点色全都在这一次求交里顺便算出。

19.4 光线-AABB：Slab 法#

AABB = 三对平行轴对齐平面。光线与每对平面求交得 $[t_{\text{near}, i}, t_{\text{far}, i}]$ ,取三组区间的交,非空即相交。

t_{\text{enter}} = \max_i t_{\text{near}, i}, \quad t_{\text{exit}} = \min_i t_{\text{far}, i}

相交条件 $t_{\text{enter}} \leq t_{\text{exit}}$ 且 $t_{\text{exit}} \geq 0$ 。

1
struct AABB {
2
    Eigen::Vector3f pmin, pmax;
3

4
    bool intersect(const Ray& r, float& t_enter, float& t_exit) const {
5
        t_enter = -std::numeric_limits<float>::infinity();
6
        t_exit  =  std::numeric_limits<float>::infinity();
7
        for (int i = 0; i < 3; ++i) {
8
            float d = r.direction[i];
9
            if (std::abs(d) < 1e-8f) {
10
                // 光线平行于第 i 轴
11
                if (r.origin[i] < pmin[i] || r.origin[i] > pmax[i]) return false;
12
                continue;
13
            }
14
            float inv_d = 1.0f / d;
15
            float t0 = (pmin[i] - r.origin[i]) * inv_d;
16
            float t1 = (pmax[i] - r.origin[i]) * inv_d;
17
            if (t0 > t1) std::swap(t0, t1);
18
            t_enter = std::max(t_enter, t0);
19
            t_exit  = std::min(t_exit,  t1);
20
            if (t_enter > t_exit) return false;
21
        }
22
        return t_exit >= std::max(0.0f, r.t_min);
23
    }
24
};

二十、空间加速结构（BVH canonical）#

WARNING
本章是全系列 BVH 的规范位置。Part 2 光栅化不用 BVH（光栅化按屏幕顺序扫描），但光线追踪、碰撞检测、剔除都会回链这里。

20.1 为什么需要加速结构#

朴素光线追踪对每条光线都与全场景 $n$ 个图元求交，复杂度 $O(n)$ 。一张 1080p 图像就是 $\approx 2\times 10^6$ 条主光线，Cornell Box 几十个三角形还行，换到百万三角形场景直接崩。加速结构把每条光线的期望代价降到 $O(\log n)$ 。

20.2 BVH 的基本思想#

Bounding Volume Hierarchy = 包围盒层次：用二叉树递归把场景切分，内部节点的 AABB 包含所有后代,叶子节点放少量图元。

光线遍历时：

与根 AABB 求交失败 → 整棵树跳过
与内部节点相交 → 递归进两个子节点
到达叶子 → 才与真正的三角形求交

20.2.1 数据结构#

1
struct BVHNode {
2
    AABB bounds;
3
    BVHNode* left = nullptr;
4
    BVHNode* right = nullptr;
5
    std::vector<const Object*> prims;   // 仅叶子非空
6
    bool is_leaf() const { return left == nullptr && right == nullptr; }
7
};

20.2.2 递归构建#

1
BVHNode* build(std::vector<const Object*>& objs) {
2
    auto* node = new BVHNode;
3
    for (auto* o : objs) node->bounds.expand(o->bounds());
4

5
    if (objs.size() <= 4) {                 // 叶子阈值
6
        node->prims = objs;
7
        return node;
8
    }
9

10
    // 沿最长轴划分：按质心中位数分两半（Naive SAH 的替代）
11
    int axis = node->bounds.longest_axis();
12
    std::nth_element(objs.begin(), objs.begin() + objs.size()/2, objs.end(),
13
        [axis](const Object* a, const Object* b) {
14
            return a->bounds().centroid()[axis] < b->bounds().centroid()[axis];
15
        });
16
    std::vector<const Object*> left(objs.begin(), objs.begin() + objs.size()/2);
17
    std::vector<const Object*> right(objs.begin() + objs.size()/2, objs.end());
18

19
    node->left  = build(left);
20
    node->right = build(right);
21
    return node;
22
}

20.2.3 光线遍历#

1
bool traverse(const BVHNode* n, const Ray& r, Intersection& best) {
2
    float t_enter, t_exit;
3
    if (!n->bounds.intersect(r, t_enter, t_exit)) return false;
4
    if (t_enter > best.t) return false;      // 已有更近的交点，整棵子树剪枝
5

6
    if (n->is_leaf()) {
7
        bool hit = false;
8
        for (auto* p : n->prims) {
9
            Intersection tmp;
10
            if (p->intersect(r, tmp) && tmp.t < best.t) { best = tmp; hit = true; }
11
        }
12
        return hit;
13
    }
14
    // 内部节点：先访问更近的一侧
15
    bool hl = traverse(n->left,  r, best);
16
    bool hr = traverse(n->right, r, best);
17
    return hl || hr;
18
}

20.3 SAH：表面积启发式#

目标:让划分方式最小化期望光线-节点相交代价。对一个内部节点，代价为：

C = C_{\text{trav}} + \frac{SA(L)}{SA(N)} \cdot n_L \cdot C_{\text{isect}} + \frac{SA(R)}{SA(N)} \cdot n_R \cdot C_{\text{isect}}

其中 $SA(\cdot)$ 是 AABB 表面积, $n_L, n_R$ 是左右子图元数。

核心假设：光线方向均匀 → 光线击中子盒的概率正比于表面积。

20.3.1 桶分（Bucketing）加速#

朴素 SAH 要试 $O(n)$ 个划分位置,每次 $O(n)$ 计算成本 → $O(n^2)$ 。实用做法分桶:把图元质心按最长轴分到 $K \approx 12$ 个桶里,只在 $K-1$ 个桶边界处试划分:

1
constexpr int NUM_BUCKETS = 12;
2
struct Bucket { int count = 0; AABB bounds; };
3

4
int best_split_sah(const std::vector<const Object*>& objs, int axis, const AABB& all) {
5
    Bucket buckets[NUM_BUCKETS];
6
    for (auto* o : objs) {
7
        Eigen::Vector3f c = o->bounds().centroid();
8
        float t = (c[axis] - all.pmin[axis]) / (all.pmax[axis] - all.pmin[axis]);
9
        int b = std::min(NUM_BUCKETS - 1, int(t * NUM_BUCKETS));
10
        buckets[b].count++;
11
        buckets[b].bounds.expand(o->bounds());
12
    }
13

14
    float best_cost = std::numeric_limits<float>::infinity();
15
    int best_i = -1;
16
    for (int i = 0; i < NUM_BUCKETS - 1; ++i) {
17
        AABB b0, b1;
18
        int n0 = 0, n1 = 0;
19
        for (int j = 0; j <= i; ++j) { b0.expand(buckets[j].bounds); n0 += buckets[j].count; }
20
        for (int j = i + 1; j < NUM_BUCKETS; ++j) { b1.expand(buckets[j].bounds); n1 += buckets[j].count; }
21
        if (n0 == 0 || n1 == 0) continue;
22
        float cost = 0.125f + (n0 * b0.surface_area() + n1 * b1.surface_area()) / all.surface_area();
23
        if (cost < best_cost) { best_cost = cost; best_i = i; }
24
    }
25
    return best_i;
26
}

20.4 与其他加速结构对比#

结构	划分方式	优势	劣势
BVH	对象空间,每节点一个 AABB	一个图元只进一个叶子;动态更新友好;实现简洁	兄弟 AABB 可能重叠
KD-Tree	空间二分,轴对齐分割面	遍历高效,经典理论最优	图元可能跨边界被复制;重建贵
八叉树	空间八分(2³ 均匀)	实现极简;体素/粒子友好	图元分布不均时浪费严重
均匀网格	空间均匀切成固定大小格子	$O(1)$ 遍历每格;构建极快	Teapot-in-stadium(尺寸悬殊)最坏情况极差

TIP
生产级光线追踪器(Embree、OptiX、PBRT)基本全用 SAH BVH。动态场景用 refit 或 LBVH(线性 BVH)。

二十一、蒙特卡洛积分#

NOTE
渲染方程是一个高维积分,维度随弹射次数指数增长。高维数值积分唯一不怕维度诅咒的方法就是蒙特卡洛:收敛速率 $O(N^{-1/2})$ 与维度无关。

21.1 基本蒙特卡洛估计量#

I = \int_\Omega f(x) \, dx \approx \hat{I} = \frac{1}{N} \sum_{i=1}^{N} \frac{f(X_i)}{p(X_i)}, \quad X_i \sim p

无偏性: $\mathbb{E}[\hat{I}] = \int f(x)/p(x) \cdot p(x) \, dx = \int f(x) \, dx = I$ 。

21.2 方差分析#

\text{Var}[\hat{I}] = \frac{1}{N} \left( \int \frac{f^2(x)}{p(x)} dx - I^2 \right)

📌 方差最小的采样分布(最优重要性采样)#

用拉格朗日乘子法最小化 $\int f^2/p$ 在约束 $\int p = 1, p \geq 0$ 下：

p^*(x) = \frac{|f(x)|}{\int |f(y)| \, dy}

此时方差为零(若 $f \geq 0$ )。但这需要知道 $\int f$ ——正是我们要算的量。实用策略：让 $p$ 尽量正比于 $|f|$ ,即重要性采样。

21.3 重要性采样(IS)#

渲染方程被积函数 $f_r L_i \cos\theta_i$ 里, $L_i$ 通常只在光源方向上大, $\cos\theta_i$ 在法线附近大, $f_r$ 在镜面方向附近大。沿这些”大值方向”采样能显著降方差。

21.3.1 余弦加权半球采样#

采样 $p(\omega) = \cos\theta/\pi$ （恰好和 Lambert 的 $\cos$ 项匹配）。用 Malley 方法：

在单位圆盘上均匀采样 $(x, y)$
取 $z = \sqrt{1 - x^2 - y^2}$

1
Eigen::Vector3f sample_cosine_hemisphere(float u1, float u2) {
2
    float r = std::sqrt(u1);
3
    float phi = 2.0f * MY_PI * u2;
4
    float x = r * std::cos(phi);
5
    float y = r * std::sin(phi);
6
    float z = std::sqrt(std::max(0.0f, 1.0f - u1));
7
    return {x, y, z};
8
}
9
float cosine_hemisphere_pdf(float cos_theta) { return cos_theta / MY_PI; }

21.3.2 均匀半球采样(对比基线)#

$p(\omega) = 1/(2\pi)$ :

\theta = \arccos(1 - u_1), \quad \phi = 2\pi u_2

Lambert 材质下余弦加权方差远小于均匀方差。

21.4 分层采样#

将积分域分成 $M$ 个等分块,每块独立采 $N/M$ 个样本:

\text{Var}_{\text{strat}} = \frac{1}{M^2} \sum_{k=1}^{M} \text{Var}_k \leq \text{Var}_{\text{uniform}}

工程上最简单的 $\sqrt{N}\times\sqrt{N}$ 网格 + 格内抖动,就能明显减噪。

21.5 多重重要性采样(MIS)#

当有多个采样策略(比如 BRDF 采样 + 光源采样)时,用加权组合:

\hat{I}_{\text{MIS}} = \sum_k \frac{1}{n_k} \sum_{j=1}^{n_k} w_k(X_{k,j}) \frac{f(X_{k,j})}{p_k(X_{k,j})}

权函数要求 $\sum_k w_k(x) = 1$ 在 $f(x) \neq 0$ 处。Veach 1997 的幂启发式( $\beta = 2$ ):

w_k(x) = \frac{(n_k p_k(x))^2}{\sum_j (n_j p_j(x))^2}

TIP
MIS 是”既采 BRDF 又采光源”能无偏结合的数学基础——没它就要在两种策略里二选一,光泽高光(BRDF 好采)和点光源(光源好采)只能兼顾其一。

二十二、路径追踪#

22.1 Whitted 风格(经典基线)#

Whitted 1980 的非物理但漂亮的模型:

直接光照:显式遍历光源 + Phong 模型
镜面反射 / 折射:递归追踪理想反射方向 / Snell 折射方向
终止:固定递归深度

1
Eigen::Vector3f whitted(const Ray& r, const Scene& s, int depth) {
2
    if (depth <= 0) return Eigen::Vector3f::Zero();
3
    Intersection hit = s.intersect(r);
4
    if (!hit.happened) return s.background;
5

6
    Eigen::Vector3f color = hit.material->emission;
7

8
    // 直接光照(所有光源累加)
9
    for (const auto& light : s.lights) {
10
        Eigen::Vector3f L_dir = (light.pos - hit.p).normalized();
11
        float dist = (light.pos - hit.p).norm();
12
        Ray shadow{hit.p + EPSILON * hit.n, L_dir, EPSILON, dist - EPSILON};
13
        if (!s.occluded(shadow)) {
14
            float cos_t = std::max(0.0f, hit.n.dot(L_dir));
15
            color += hit.material->brdf(-r.direction, L_dir, hit.n)
16
                   * light.intensity * cos_t / (dist * dist);
17
        }
18
    }
19
    // 镜面递归
20
    if (hit.material->specular) {
21
        Eigen::Vector3f R = reflect(r.direction, hit.n);
22
        Ray refl{hit.p + EPSILON * hit.n, R};
23
        color += hit.material->Ks * whitted(refl, s, depth - 1);
24
    }
25
    return color;
26
}

局限:只有”光源 → 镜面链 → 眼”这一种路径被显式算。diffuse 之间的相互反射(color bleeding)完全缺失。

22.2 基础路径追踪(Kajiya 1986)#

每次弹射用蒙特卡洛随机采一个方向,递归直到俄罗斯轮盘赌终止。本质是对 Neumann 级数 $\sum_n T^n L_e$ 的无偏估计。

22.2.1 朴素实现(只按 BRDF 采样)#

1
Eigen::Vector3f path_trace_naive(const Ray& r, const Scene& s, int depth) {
2
    Intersection hit = s.intersect(r);
3
    if (!hit.happened) return s.background;
4

5
    Eigen::Vector3f L = hit.material->emission;
6
    if (depth >= MAX_DEPTH) return L;
7

8
    // 俄罗斯轮盘赌
9
    const float p_rr = 0.8f;
10
    if (random_float() > p_rr) return L;
11

12
    // 按 BRDF 采样新方向
13
    Eigen::Vector3f wi;
14
    float pdf;
15
    Eigen::Vector3f f = hit.material->sample(-r.direction, hit.n, wi, pdf);
16
    if (pdf < 1e-8f) return L;
17

18
    float cos_t = std::max(0.0f, hit.n.dot(wi));
19
    Ray next{hit.p + EPSILON * hit.n, wi};
20
    Eigen::Vector3f Li = path_trace_naive(next, s, depth + 1);
21

22
    L += f.cwiseProduct(Li) * cos_t / pdf / p_rr;
23
    return L;
24
}

问题:小光源击中率低,噪点可怕。

22.3 分直接/间接的路径追踪(GAMES101 PA7 风格)#

关键改进:每次交点显式采光源算直接光照,只让递归负责”下一次弹射后的间接光”——这样光源采样的低方差被充分利用。

1
Eigen::Vector3f shade(const Intersection& hit, const Eigen::Vector3f& wo, const Scene& s) {
2
    if (hit.material->has_emission()) return hit.material->emission;
3

4
    // —— 直接光照:从光源上采一点 ——
5
    Eigen::Vector3f L_dir = Eigen::Vector3f::Zero();
6
    Intersection light_hit;
7
    float light_pdf;
8
    s.sample_light(light_hit, light_pdf);
9
    Eigen::Vector3f ws = (light_hit.p - hit.p).normalized();
10
    float dist2 = (light_hit.p - hit.p).squaredNorm();
11
    Ray to_light{hit.p + EPSILON * hit.n, ws};
12
    Intersection chk = s.intersect(to_light);
13
    if (chk.happened && (chk.p - light_hit.p).squaredNorm() < 1e-3f) {
14
        Eigen::Vector3f f = hit.material->eval(ws, wo, hit.n);
15
        float cos_t  = std::max(0.0f, hit.n.dot(ws));
16
        float cos_tl = std::max(0.0f, light_hit.n.dot(-ws));
17
        L_dir = light_hit.material->emission.cwiseProduct(f) * cos_t * cos_tl / dist2 / light_pdf;
18
    }
19

20
    // —— 间接光照:俄罗斯轮盘赌 + BRDF 采样 ——
21
    Eigen::Vector3f L_indir = Eigen::Vector3f::Zero();
22
    const float p_rr = 0.8f;
23
    if (random_float() < p_rr) {
24
        Eigen::Vector3f wi = hit.material->sample(wo, hit.n);
25
        float pdf = hit.material->pdf(wi, wo, hit.n);
26
        if (pdf > 1e-8f) {
27
            Ray next{hit.p + EPSILON * hit.n, wi};
28
            Intersection next_hit = s.intersect(next);
29
            if (next_hit.happened && !next_hit.material->has_emission()) {
30
                Eigen::Vector3f f = hit.material->eval(wi, wo, hit.n);
31
                float cos_t = std::max(0.0f, hit.n.dot(wi));
32
                L_indir = shade(next_hit, -wi, s).cwiseProduct(f) * cos_t / pdf / p_rr;
33
            }
34
        }
35
    }
36
    return L_dir + L_indir;
37
}

22.3.1 为什么要跳过命中光源的间接项#

直接光照里我们已经显式采了光源。如果间接递归再撞上光源,那条路径被算了两次,结果偏亮。判据 !next_hit.material->has_emission() 就是排除这种重复。

22.3.2 俄罗斯轮盘赌为什么无偏#

设继续概率 $p$ ,继续时把贡献除以 $p$ :

\mathbb{E}[X] = p \cdot \frac{f}{p} + (1 - p) \cdot 0 = f

与直接算 $f$ 等价。但方差被放大 $1/p - 1$ ,所以 $p$ 不能太小(典型 0.8)。

22.4 BRDF/光源 MIS 组合#

当表面是光泽(非漫反射、非纯镜面)时,BRDF 采样可能比光源采样更有效。MIS 把两种策略的样本加权合并——工业级路径追踪器的标配。完整实现见 PBRT 第 13 章。

二十三、高级技术概述(Part 6 前导)#

23.1 双向路径追踪(BDPT)#

从摄像机走一条路径 + 从光源走一条路径,然后在所有 $(s, t)$ 组合连接点上评估贡献,MIS 合并。

强项:焦散、镜面-漫反射-镜面(SDS)路径
代价:每像素评估连接次数爆炸

23.2 光子映射#

两阶段:(1) 从光源发射光子、存到 KD-Tree; (2) 相机视角光线交到表面时,对附近光子做密度估计。

强项:焦散(如玻璃杯在桌面上的光斑)
弱项:低频偏差,需要仔细选半径

23.3 Metropolis Light Transport(MLT)#

基于马尔科夫链的采样,对贡献大的路径局部扰动,自适应聚焦困难光路。对 SDS、细缝光极其有效,但收敛分析困难。

23.4 体渲染与参与介质#

当光穿过雾、烟、皮下组织时,Beer-Lambert 衰减 + phase function + in-scattering/out-scattering——渲染方程要扩展成体渲染方程(Volume Rendering Equation)。NeRF 的数学底层就是它,Part 6 展开。

小结#

主题	工具	canonical 位置
辐射度量学	$L, E, \Phi, I$ 四件套	§16 — 是 Part 2 光照、Part 6 PBR 的物理基础
渲染方程	Kajiya 1986 积分/递归/Neumann	§17 — 后面所有章节回链
BRDF	定义 / 互易 / 能量守恒 / Lambert-Phong-Cook-Torrance	§18 — Part 2 光照模型的物理正名
求交	Ray 参数化 / 球 / Möller-Trumbore / Slab	§19
BVH	递归构建 / SAH / 桶分	§20 — 光线追踪与碰撞检测的通用加速
蒙特卡洛积分	$O(N^{-1/2})$ / 重要性采样 / 分层 / MIS	§21 — 渲染方程唯一的数值工具
路径追踪	Whitted / 朴素 PT / 直接+间接分离 / RR 终止	§22 — 物理渲染的黄金算法
高级技术	BDPT / 光子映射 / MLT / 体渲染	§23 → 细节在 Part 6

下一部分进入 Part 5:动画与物理模拟——关键帧、样条插值、粒子系统、弹簧质点、刚体动力学、流体模拟。

计算机图形学笔记(四)：光线追踪与全局光照

https://kyc001.github.io/posts/计算机图形学笔记四/

作者

kyc001

发布于

2025-07-23

许可协议

CC BY-NC-SA 4.0

计算机图形学笔记(五)：动画与物理仿真

计算机图形学笔记(三)：几何建模与处理

18.3 Lambert 漫反射（完美各向同性散射）

18.4 Phong / Blinn-Phong（经验模型）

18.5 Cook-Torrance（物理微表面模型）

19.3 光线-三角形：Möller-Trumbore

19.4 光线-AABB：Slab 法

5

二十、空间加速结构（BVH canonical）

22.1 Whitted 风格(经典基线)

22.2 基础路径追踪(Kajiya 1986)

22.3 分直接/间接的路径追踪(GAMES101 PA7 风格)

22.4 BRDF/光源 MIS 组合

8

二十三、高级技术概述(Part 6 前导)

23.1 双向路径追踪(BDPT)

23.2 光子映射

23.3 Metropolis Light Transport(MLT)

23.4 体渲染与参与介质

小结

目录#

十六、辐射度量学基础#

16.1 立体角#

16.1.1 定义#

16.1.2 球面坐标下的微分立体角#

16.1.3 验证半球面积#

16.2 辐射度量四件套#

16.2.1 辐射度的严格定义#

16.2.2 辐照度与辐射度的关系#

16.3 Lambert 余弦定律#

十七、渲染方程（canonical）#

17.1 渲染方程的积分形式#

17.2 递归形式#

17.3 Neumann 级数展开#

📌 各项的物理意义#

17.4 收敛性#

十八、BRDF 与材质模型#

18.1 BRDF 的定义#

18.2 BRDF 必须满足的三条性质#

18.3 Lambert 漫反射（完美各向同性散射）#

📌 为什么分母是 π\piπ（完整推导）#

18.4 Phong / Blinn-Phong（经验模型）#

18.5 Cook-Torrance（物理微表面模型）#

18.5.1 GGX（Trowbridge-Reitz）分布#

18.5.2 Schlick 菲涅尔近似#

18.5.3 Smith 几何项#

十九、光线-几何体相交#

19.1 光线参数化#

19.2 光线-球体求交#

19.3 光线-三角形：Möller-Trumbore#

19.4 光线-AABB：Slab 法#

二十、空间加速结构（BVH canonical）#

20.1 为什么需要加速结构#

20.2 BVH 的基本思想#

20.2.1 数据结构#

20.2.2 递归构建#

20.2.3 光线遍历#

20.3 SAH：表面积启发式#

20.3.1 桶分（Bucketing）加速#

20.4 与其他加速结构对比#

二十一、蒙特卡洛积分#

21.1 基本蒙特卡洛估计量#

21.2 方差分析#

📌 方差最小的采样分布(最优重要性采样)#

21.3 重要性采样(IS)#

21.3.1 余弦加权半球采样#

21.3.2 均匀半球采样(对比基线)#

21.4 分层采样#

21.5 多重重要性采样(MIS)#

二十二、路径追踪#

22.1 Whitted 风格(经典基线)#

22.2 基础路径追踪(Kajiya 1986)#

22.2.1 朴素实现(只按 BRDF 采样)#

22.3 分直接/间接的路径追踪(GAMES101 PA7 风格)#

22.3.1 为什么要跳过命中光源的间接项#

22.3.2 俄罗斯轮盘赌为什么无偏#

22.4 BRDF/光源 MIS 组合#

二十三、高级技术概述(Part 6 前导)#

23.1 双向路径追踪(BDPT)#

23.2 光子映射#

23.3 Metropolis Light Transport(MLT)#

23.4 体渲染与参与介质#

小结#

📌 为什么分母是 $\pi$ （完整推导）#