萘Naie_Char

立体角#

定义

• A：锥体在球面上所截下的面积
• r：球半径

球面坐标下的微分立体角 微分立体角描述的是三维空间中一个无限小方向，用球面坐标表示为（这里的字母跟微积分教材有出入，在图形学中，一般 $\theta$ 是极角，$\phi$ 是方位角）：

NOTE

为什么需要 $\sin\theta$ 因子？因为在越靠近球顶部的位置，一个小的 $d\phi$ 实际覆盖的面积远小于在赤道附近的面积，所以需要一个 $\sin\theta$ 修正

辐射度量#

• 辐射通量：单位时间通过表面/区域的能量，相当于光的功率
• 辐射强度：单位立体角上的辐射通量
• 辐照度：表面单位面积接收到的辐射通量（收进来的光） $E=\frac{d\phi}{dA}$
• 辐射度：每单位投影面积、单位立体角的辐射通量（沿某个方向发出去/传过来的光）

辐射度定义

辐照度与辐射度的关系

这也就是为什么说：

• 辐照度：一个单位表面收到的所有方向光的总量。
• 辐射度：一个单位表面朝某个方向发出的光强度。

渲染方程的积分形式#

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

翻译：一个点朝某个方向发出的光 = 自己发出的光 + 从周围所有方向来的光经过材质反射后的总和。

各项含义（下标 $i$ 指的是 incoming 进入）：

• $L_o$：出射辐射度（在点 p 朝观察方向 $\omega_o$发出的光。）
• $L_e$：自发光项
• $f_r$：BRDF，材质把入射光变成出射光的比例
• $L_i$：从某方向 $\omega_i$ 射到 p 的入射辐射度
• $\cos\theta_i = \vec{n} \cdot \omega_i$：即 Lambert 余弦
• $\Omega$：上半球，积分微元 $d\omega_i = \sin\theta_i \ d\theta_i \ d\phi_i$

渲染方程递归#

$L_i(\mathbf{p}, \omega_i)$ 本身就是从 $-\omega_i$ 方向射向 $\mathbf{p}$ 的光——如果场景中沿 $-\omega_i$ 走遇到点 $\mathbf{p}'$，（$\mathbf{p}'$ 是该光线的源头，可能是反射面，也可能直接是光源）那么：

对于点 $\mathbf{p}'$ 又可以依照此法进行回溯，找到点 $\mathbf{p}''$，直到最后找到光源，这便是 全局光照。

Neumann 级数展开#

这一节便是把光无限反弹的过程展开，知道为什么计算机能够算出来。
定义光传输算子 $T$：

$T$：一次反弹操作，也就是说输入 $L$，经过 $T$，便得到 $TL$（让这束光在场景中传播一次）。
因此，渲染方程写成算子形式：$L = L_e + TL$
进一步形式上求解：

$T^n L_e$ 项表示第 n 次弹射

然而无限级数无法计算，于是我们通过

• Monte Carlo 随机采样（因为积分 $T$ 本身包含 $\int_\omega$，其方向无限，因此随机取几个方向进行估计积分）
• Russian Roulette 截断（把无限级数截断成有限级数进行计算，人为创造概率终止，最后为了保证平均结果不变，还要增加权重补偿。）

双向反射分布函数#

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

可以简单理解为：$BRDF = \frac{出射光}{入射光}$，根据单位为 ${sr}^{-1}$ 可知，BRDF 描述的是每单位立体角方向上的反射能力

BRDF三条性质：

• 非负性
• Helmholtz 互易性（交换入射出射方向，BRDF 一样）
• 能量守恒： 对任意 $\omega_o$

Lambert 漫反射#

定义:

$\rho_d$ 称为漫反射率（albedo）。
要求能量守恒：

当 $\rho_d = 1$ 时为完美白漫反射。
表面把入射光均匀散射到整个半球，但为了保证所有方向加起来反射能量等于 $\rho_d$，必须除以以 $cos\theta$ 加权后的半球积分面积 π

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

从 BRDF 角度看，Phong 本质是在设计一个 $f_r$，描述光如何从 $\omega_i$ 分布到 $\omega_o$

由于 $shininess$ 取不同值，高光总能量不同，因此我们需要乘一个系数（归一化因子）$\frac{(n + 2)}{2\pi}$，让整个半球积分恰巧为1

若是 Blinn-Phong 模型，则用半程向量 $\vec{h}$ 替换 $\vec{r}$，归一化因子也变化：

Cook-Torrance#

现代 PBR 核心数学模型：Cook-Torrance 认为真实物体表面不是一个平面，而是由无数微小镜子组成，最终看到的反射效果取决于这些微镜子的朝向分布、镜面反射强度以及微镜子之间的遮挡。

• $D$：有多少微表面朝向正确
• $F$：这些微表面反射多少光
• $G$：有多少光能真正到达/离开
• 分母：几何修正

光线参数化#

将光线翻译成高中解析几何里的直线方程：

• $\mathbf{o}$：起点位置
• $\vec{d}$：方向（通常取单位向量）
• $t$：光线上的第几个位置
• $t_{\min} > 0$：$t_{\min}$ 取 0 会因为自相交而产生 Shadow acne（阴影痤疮）（$t_{\min}$ 典型值取 $10^{-4}$，$t_{\max}$ 取无限大）

光线是光线追踪里最核心的数据结构：

1
struct Ray
2
{
3
    Eigen::Vector3f origin;
4
    Eigen::Vector3f direction;
5
    float t_min = 1e-4f;
6
    float t_max = std::numeric_limits<float>::infinity();   // 即 t_max 取正无穷
7

8
    Eigen::Vector3f at(float t) const
9
    {
10
        return origin + t * direction;
11
    }
12
}

光线-球体求交#

球面方程：$\|\mathbf{p} - \mathbf{c}\|^2 = r^2$ （其中，p 为球面上一点，c 为球心，R 为半径）
因为交点一定在光线上，所以将光线方程 $\mathbf{p} = \mathbf{o} + t\vec{d}$ 代入上式可得（关于 t 的一元二次方程）:

计算判别式 $\Delta = b^2 - 4c$ 判断交点情况：

• 当 $\vec{d}$ 为单位向量时， $a = \vec{d}\cdot\vec{d}=1$

代码实现：

1
bool intersect_sphere(const Ray& ray, const Eigen::Vector3f& c, float radius, float& t_hit)
2
{
3
    Eigen::Vector3f oc = ray.origin - c;
4
    float b = 2.0f * oc.dot(ray.direction);
5
    float c = oc.squaredNorm() - radius * radius;   // squaredNorm() 为 Eigen 函数，平方长度
6

7
    float disc = b * b - 4.0f * c;
8

9
    if (disc < 0)   return false;
10

11
    float sqrt_disc = std::sqrt(dist);  // 对判别式开方
12
    float t0 = (-b - sqrt_d) * 0.5f;
13
    float t1 = (-b + sqrt_d) * 0.5f;
14
    if (t0 > t1) std::swap(t0, t1);
15
    t_hit = (t0 >= r.t_min) ? t0 : t1;  // 如果 t0 合法，就选 t0 （更近）
16

17
    return t_hit >= r.t_min && t_hit <= r.t_max;
18
}

得到 $t_{hit}$ 后，将其代入光线求得交点，然后计算法线 $n=\frac{p - c}{R}$，再进入光照计算。

光线-三角形求交#

核心思想：把交点用重心坐标表示。（这里的重心坐标不要理解为以前几何中的重心，而是三个顶点对交点坐标贡献的权重）
由于交点既在光线上，也在三角形上：

为何用重心坐标（u, v）：

• 只有三个未知量，如果用重心（x, y, z）将有四个未知量（t 也是未知量）
• 计算好 （u, v）后，后续法线，uv，顶点颜色都可以直接进行插值，一举两得。

继续移项可得：

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

将上式变为矩阵形式可得：

这相当于就是三元一次方程组，$t, u, v$ 是未知量。此时理论上求解 $(t, u, v)$ 只用把最左边的矩阵变成逆矩阵移到右边即可，但求逆矩阵很慢，这里便使用一个数学技巧，Cramer 法则 + 标量三重积恒等式 $(\vec{a}\times\vec{b})\cdot\vec{c} = (\vec{b}\times\vec{c})\cdot\vec{a}$，计算结果全由叉乘和点乘表示。

最后结果可得：

其中 $\mathbf{S} = \mathbf{o} - \mathbf{V}_0$ ，$\mathbf{P} = \vec{d}\times\mathbf{E}_2$，$det = (\vec{d}\times\mathbf{E}_2)\cdot\mathbf{E}_1$（行列式）

代码实现：

1
bool moller_trumbore(const Ray& ray, float& t, float& u, float& v,
2
                     const Eigen::Vector3f& v0,
3
                     const Eigen::Vector3f& v1,
4
                     const Eigen::Vector3f& v2)
5
{
6
    const float EPS = 1e-8f;    // 图形学中经常出现 EPS，表示一个极小量，如果一个数小于它，那么就可以认为是0
7
    Eigen::Vector3f e1 = v1 - v0;
8
    Eigen::Vector3f e2 = v2 - v0;
9

10
    Eigen::Vector3f p  = ray.direction.cross(e2);
11
    Eigen::Vector3f s = ray.origin - v0;
12
    float det = p.dot(e1);
13
    float inv_det = 1.0f / det;
14
    if (std::abs(det) < EPS) return false;  // 光线平行于三角形，因为 det = 0，矩阵不可逆，三向量共面。
15

16
    // 求u
17
    u = p.dot(s) * inv_det;
18
    if (u < 0 || u > 1) return false;
19

20
    // 求v
21
    Eigen::Vector3f q = s.cross(e1);
22
    v = q.dot(r.direction) * inv_det;
23
    if (v < 0 || u + v > 1) return false;
24

25
    // 求t
26
    t = q.dot(e2) * inv_det;
27
    return t >= r.t_min && t <= r.t_max;
28
}

AABB 轴对齐包围盒#

目的：快速排除大量不可能相交的物体。
我们用 slab 法 得到空间包围盒。 slab 表示平行同一个轴的两平面之间的区域，这里会用到三个 slab，分别平行 x，y，z 轴，他们在三维空间中一定会围成一个立方体，这个立方体就是轴对齐包围盒。
下面计算一束光线进入和穿出包围盒的坐标：

先看下面这张图，先理解二维平面光线的行为。

AABB 光线求交的 slab 方法二维示意图

蓝色的竖直条带就是 “x方向的slab” —— 两条竖直虚线之间的无限长条形区域。光线只要在这个x范围内，就算 “在x的slab里”。
橙色的水平条带同理，是 “y方向的slab” —— 两条水平虚线之间的无限长条带。

对每个slab单独看，光线在该slab内部的t区间是 [t_near, t_far]。图中标注了:

• x方向：t_near.x ≈ 0.23, t_far.x ≈ 0.65
• y方向：t_near.y = 0.25, t_far.y ≈ 0.88

而光线真正”同时在两个slab里” (也就是在AABB里) 的那一段,是这两个区间的交集——图中绿色粗线部分，所以：

• t_enter = max(t_near_x, t_near_y) = 0.25
• t_exit = min(t_far_x, t_far_y) = 0.65

如果 t_enter > t_exit,说明光线在还没同时进入所有slab之前,就已经从某个slab跑出去了，即与包围盒不相交。

推广到三维同理：

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

1
struct AABB
2
{
3
    Eigen::Vector3f pmin, pmax; // AABB 顶点的最小、最大坐标，只需两个坐标即可确认一个 AABB
4

5
    bool intersect(const Ray& ray, float& t_enter, float& t_exit) const
6
    {
7
        const float EPS = 1e-8f;
8
        t_enter = -std::numeric_limits<float>::infinity();
9
        t_exit  =  std::numeric_limits<float>::infinity();
10

11
        // 三个轴分别计算
12
        for (int i = 0; i < 3; i++)
13
        {
14
            float d = ray.direction[i];   // Eigen 重载了[]。direction[0] 表示 x，direction[1] 表示 y
15
            if (std::abs(d) < EPS)  // i 分量几乎为0，则光线平行于 i 面
16
            {
17
                if (ray.origin[i] < pmin[i] || ray.origin[i] > pmax[i]) return false;
18
                continue;
19
            }
20
            float inv_d = 1.0f / d;
21
            // 求与盒子的交点，已知交点在 i 面上，pmin 和 pmax 知晓，代入光线方程求解 t
22
            float t0 = (pmin[i] - r.origin[i]) * inv_d;
23
            float t1 = (pmax[i] - r.origin[i]) * inv_d;
24
            if (t0 > t1) std::swap(t0, t1);
25
            // 更新进入与穿出坐标
26
            t_enter = std::max(t_enter, t0);
27
            t_exit  = std::min(t_exit,  t1);
28

29
            if (t_enter > t_exit) return false;
30
        }
31
        return t_exit >= std::max(0.0f, ray.t_min); // 作 max 是一个防御性写法，防止 t_min 被意外设置成负数。
32
    }
33
}

BVH 的核心思想#

BVH：Bounding Volume Hierarchy 即包围体层次结构。它将 Box 逐级分类直至形成一棵树，最后的叶子里才放真正的三角形。
光线射入过程：

• 光线与根 AABB 求交失败：直接整棵树跳过
• 若与根相交，判断光线与哪个子节点 AABB 相交，判断完后直接舍弃另外一半的子节点进入下一层。
• 如此一直递归到子叶与真正三角形求交

加速结构最后将复杂度变为了 $O(\log n)$

1
struct BVHNode
2
{
3
    AABB bounds;    // 每个节点都有一个自己的包围盒
4
    std::unique_ptr<BVHNode> left;
5
    std::unique_ptr<BVHNode> right;
6
    std::vector<const Object*> prims;   // 里面的元素是物体的指针，Object 是所有几何体的统一接口
7
    // 判断是否为叶子
8
    bool is_leaf() const
9
    {
10
        return !left && !right;
11
    }
12
}

1
// 保存交点信息的核心结构体
2
struct Intersection
3
{
4
    bool happened = false;        // 是否发生相交
5
    float t = std::numeric_limits<float>::infinity();  // 光线参数 t
6
    Eigen::Vector3f coords;     // 交点坐标
7
    Eigen::Vector3f normal;
8
    const Object* obj = nullptr;  // 被击中的物体
9
    Eigen::Vector2f uv;
10
};

1
std::unique_ptr<BVHNode> build(std::vector<const Object*>& objs)
2
{
3
    auto node = std::make_unique<BVHNode>();
4
    // 计算当前节点包围盒
5
    for (auto* o : objs)
6
        node->bounds.expand(o->get_bounds());
7

8
    // 叶子阈值
9
    if (objs.size() <= 4)
10
    {
11
        node->prims = objs;
12
        return node;
13
    }
14

15
    // 选择分割方向，一般分割最长的轴，比如返回0代表分割 x 轴
16
    int axis = node->bounds.longest_axis();
17
    // 找中位数但不完全排序
18
    std::nth_element(objs.begin(), objs.begin() + objs.size()/2, objs.end(),
19
                    [axis](const Object* a, const Object* b)
20
                    {
21
                        return a->get_bounds().centroid()[axis] < b->get_bounds().centroid()[axis];
22
                    }); // 匿名函数告诉算法如何比较两个物体
23
    // 分成左右数组 left 与 right
24
    std::vector<const Object*> left(objs.begin(), objs.begin() + objs.size()/2);
25
    std::vector<const Object*> right(objs.begin() + objs.size()/2, objs.end());
26

27
    node->left = build(left);
28
    node->right = build(right);
29
    return node;    // move 语义
30
}

1
bool  traverse(const BVHNode* n, const Ray& ray, Intersection& best)
2
{
3
    if (node == nullptr)  return false;
4

5
    float t_enter, t_exit;
6

7
    // 当前节点的 AABB 与光线不相交
8
    // 第一次裁剪。同时传入的是引用，间接返回 t_enter，t_exit
9
    if (!n->bounds.intersect(ray, t_enter, t_exit))   return false;
10

11
    // 第二次裁剪子树，已经有更近的交点
12
    if (t_enter > best.t)   return false;
13

14
    /*处理叶子节点*/
15
    if (n->is_leaf)
16
    {
17
        bool hit = false;
18
        for (auto* p : n->prims)    // p 就是一个图元
19
        {
20
            Intersection tmp;
21
            // 利用多态，intersect 根据不同图元调用不同图元的求交函数
22
            if (p->intersect(r, tmp) && tmp.t < bext.t)
23
            {
24
                // 如果碰到图元，且交点更近
25
                best = tmp;
26
                hit = true;
27
            }
28
        }
29
        return hit;
30
    }
31

32
    // 内部节点，先计算左右子节点 AABB 的进入距离 (t_enter)，然后优先遍历更近的一侧
33
    float left_enter, left_exit;
34
    float right_enter, right_exit;
35

36
    bool hit_left_box =
37
        node->left &&
38
        node->left->bounds.intersect(ray, left_enter, left_exit);
39

40
    bool hit_right_box =
41
        node->right &&
42
        node->right->bounds.intersect(ray, right_enter, right_exit);
43

44
    bool hit = false;
45

46
    // 两边都命中
47
    if (hit_left_box && hit_right_box)
48
    {
49
        if (left_enter < right_enter)
50
        {
51
            hit = traverse(node->left, ray, best) || hit;
52
            hit = traverse(node->right, ray, best) || hit;
53
        }
54
        else
55
        {
56
            hit = traverse(node->right, ray, best) || hit;
57
            hit = traverse(node->left, ray, best) || hit;
58
        }
59
    }
60
    else if (hit_left_box)
61
    {
62
        hit = traverse(node->left, ray, best) || hit;
63
    }
64
    else if (hit_right_box)
65
    {
66
        hit = traverse(node->right, ray, best) || hit;
67
    }
68

69
    return hit;
70
}

SAH：表面启发式#

目标：让划分方式最小化 期望光线-节点相交代价
前面我们建 BVH 用的是最长轴的中位数划分，这样简单且快，但查找不一定最快。

1
constexpr int NUM_BUCKETS = 12;
2
struct Bucket
3
{
4
    int count = 0;  // 桶中图元数
5
    AABB bounds;
6
}
7

8
int best_split_sah(const std::vector<const Object*>& objs, int axis, const AABB& all)
9
{
10
    Bucket buckets[NUM_BUCKETS];
11
    // 遍历所有图元
12
    for (auto* o : objs)
13
    {
14
        Eigen::Vector3f center = o->get_bounds().centroid();
15
        float t = (center[axis] - all.pmin[axis]) / (all.pmax[axis] - all.pmin[axis]);  // 比例
16
        int b = std::min(NUM_BUCKETS - 1, int(t * NUM_BUCKETS));    // 算桶编号
17
        buckets[b].count++;
18
        buckets[b].bounds.expand(o->bounds());
19
    }
20

21
    float best_cost = std::numeric_limits<float>::infinity();
22
    int best_i = -1;    // 最佳切分位置，之后比如得到 best_i = 5，表示在5,6桶之间切。
23
    // 开始尝试所有切法
24
    for (int i = 0; i < NUM_BUCKETS - 1; i++)
25
    {
26
        AABB b0, b1;    // 每种切法所分得的左右盒子
27
        int n0 = 0, n1 = 0; // 左右两边各多少图元
28

29
        // 计算左边
30
        for (int j = 0; j <= i; j++)
31
        {
32
            b0.expand(buckets[j].bounds);
33
            n0 += buckets[j].count;
34
        }
35
        // 计算右边
36
        for (int j = i + 1; j < NUM_BUCKETS; j++)
37
        {
38
            b1.expand(buckets[j].bounds);
39
            n1 += buckets[j].count;
40
        }
41
        if (n0 == 0 || n1 == 0) continue;   // 极端情况，图元全在0桶，导致右边没有图元，这样的切分是无意义的
42

43
        // 计算 SAH
44
        float cost = 0.125f + (n0 * b0.surface_area() + n1 * b1.surface_area()) / all.surface_area();
45
        // 更新最优
46
        if(cost < best_cost)
47
        {
48
            best_cost = cost;
49
            best_i = i;
50
        }
51
    }
52
    return best_i;
53
}

基本蒙特卡洛估计量#

这里先明确，我们要计算的目标是 $\int f(x) dx$，基本蒙特卡洛方法：从某个概率密度 $p(x)$ 里抽样 $X1,…,X_N$ 构造估计量：

之所以除以 $p(X_i)$，是为了补偿概率。如果是均匀采样，则 $p(X_i)=1$，如果是非均匀采样，概率小的权重应更大（概率小说明很难抽到，抽到一次便代表很多没抽到的点），概率大的权重应更小。

无偏性（期望正好为 $I$）：即平均意义下，答案一定正确，数学证明：

方差分析#

单个样本的方差是 $Var[Y] = E[Y^2] - E[Y]^2$，这里 $Y = \frac{f(X)}{p(X)}$

所以单个样本方差是 $\int \frac{f^2}{p} dx - I^2$。$N$ 个独立样本取平均，方差除以 $N$，就得到方差公式:

这里 $I^2$ 是与 $p$ 无关的常数，由于标准差与噪点息息相关，我们必须让方差最小，于是本质是让 $\int f^2 / p$ 最小。

我们可以用拉格朗日乘算法求 $\int f^2 / p$ 最小值，此时的约束条件是 $\int p(x)dx = 1$，最终得到最优 $p^*(x) \propto |f(x)|$：

还记得我们要求的是 $I$ 吗，但注意到最优 $p^*(x)$ 里的分母 $\int |f(y)| dy$ 恰巧就是 $I$，这便产生的悖论，我们只能理论知道方差最小能到多少（零），但是没法直接套公式使用，于是便引出了重要性采样。

重要性采样（IS）#

我们不能求出精确的 $p^*(x)$，那就只能退一步让 $p(x)$ 的形状大致跟 $|f(x)|$ 一样：

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

余弦加权采样#

Melly 方法与余弦加权采样直觉图

Melly 操作：

• 在一个单位圆盘上,均匀地(按面积均匀)撒点 $(x, y)$
• 把每个点竖直向上抬至半球面上,即补一个 $z = \sqrt{1 - x^2 - y^2}$

这就是余弦加权采样”点集中在法线附近”这个现象的几何来源，跟半球本身的曲率相关。所以，偏向法线方向（$\theta$ 小）采样的概率密度更大：

概率密度的严格推导：
关键恒等式：圆盘上的面积微元 $dA_{\text{disk}}$，和它投影对应的半球立体角微元 $d\omega$，满足

圆盘上是均匀采样，密度 $p_{\text{disk}} = \frac{1}{\pi}$ (单位圆盘面积为 $\pi$)。 根据概率密度在变量替换下必须保持”概率质量守恒”：

代入 $dA_{\text{disk}} = \cos \theta d\omega$ 即可得：

1
Eigen::Vector3f sample_cosine_hemisphere(float u1, float u2) // u1 和 u2 就是两个服从 [0,1) 均匀分布的随机数
2
{
3
    float r = std::sqrt(u1);    // 求半径，即此点在圆盘上对应的极坐标中的 r
4
    float phi = 2.0f * PI * u2; // 求极角
5
    float x = r * std::cos(phi);
6
    float y = r * std::sin(phi);
7
    float z = std::sqrt(std::max(0.0f, 1.0f - u1));
8
    return {x, y, z};
9
}

• 为什么求 r 要开方？我们取点的方法是按面积均匀采样的，这要求半径的概率满足 $r \propto \sqrt{u}$，（圆环面积随半径线性增长,越往外单位角度的面积越大）

均匀半球采样（作对比）#

这里换了个思路，不走投影，而是直接对半球的极角做反函数采样，半球上均匀分布的 $\theta$ 累积分布是 $F(\theta) = 1 - \cos\theta$：

这样得到的 $p(\omega) = \frac{1}{2\pi}$ 不会偏向任何方向，处处相同。

Lambert 材质下余弦加权方差远小于均匀采样#

渲染方程里 Lambert 材质的被积函数是 $f_r \cdot L_i(\omega) \cdot \cos \theta$，其中漫反射 BRDF 是常数 $f_r = \rho / \pi$( $\rho$ 是反照率)。 代入蒙特卡洛估计量 $\frac{1}{N} \sum \frac{f(\omega_i)}{p(\omega_i)}$：

用余弦加权 $p = \cos \theta / \pi$：

$\cos \theta$ 和 $\pi$ 完全约掉了，这说明每次采样只剩下 $\rho \cdot L_i(\omega)$，估计量的波动只取决于 $L_i$ (入射光本身在各方向的强弱差异)，而不再受 $\cos \theta$ 这个额外的、在掠射角处会把贡献压到接近零的因子影响。

分层采样#

纯随机(蒙特卡洛)采样虽然无偏，但不可避免会出现点随机扎堆，随机留出大片空白的情况，于是我们将积分域划分成 $M$ 哥等分块，每块独立采样 $N / M$ 个样本，这样保证每个格子都有样本，且样本在格子内是随机偏移的

• $\text{Var}_k$：每个格子内部的方差

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

多重重要性采样 (MIS)#

当被积函数是好几个”性格迥异”的因子相乘时，不存在一个单一的 $p(x)$ 能同时匹配所有因子的形状。

比如直接光照的被积函数为 $f_r(\omega)\cdot L_i(\omega)$：

• $f_r(\text{BRDF})$ 在光泽材质下会有一个又窄又高的镜面反射瓣，只有对着反射方向附近采样才能命中，这种情况下 “BRDF采样” 效率很高。
• $L_i(\text{入射光})$ 如果是一个小面积光源或点光源,只在光源所在的一小片方向上非零,这种情况下”光源采样”效率很高。

我们对两种方式进行加权组合：

多重重要性采样

幂启发式

最朴素的想法是”按 (p_k) 的大小线性分配权重”（这叫 balance heuristic，把上式的平方去掉就是），Veach 发现平方一下效果更好：如果某个策略在 (x) 点的 (p_k(x)) 相对其他策略很小（比如光源采样在一个高光溢出的方向上概率几乎为零），那 (f/p_k) 这一项本来就很容易因为分母极小而炸出一个离群的巨大数值（这就是渲染里常见的”萤火虫”噪点的来源之一）。平方之后，这个小 (p_k) 对应的权重会被压得比线性分配时更低，相当于在容易出现极端噪声的地方，更果断地把信任度让给 pdf 更大的那个策略，从源头上抑制了萤火虫噪点，这也是为什么工业界渲染器几乎都默认用 power heuristic 而不是最朴素的 balance heuristic。

TIP

BRDF采样天生擅长”追高光”，光源采样天生擅长”追亮点”，两者的peak常常根本不在一个方向上。MIS的价值就是让你两个都采，并且用一套有数学保证(无偏)的权重规则，自动让每个样本在它”更擅长”的地方发挥更大作用。

基础路径追踪（Kajiya 1986）#

整个函数的流程：

朴素实现（只按 BRDF 采样）

1
Eigen::Vector3f path_trace_naive(const Ray& r, const Scene& s, int depth)
2
{
3
    Intersection hit = s.intersect(r);
4
    if (!hit.happened) return s.background; // 光线未击中任何目标，返回背景颜色。
5

6
    Eigen::Vector3f L = hit.material->emission; // 自发光 L_e
7

8
    if (depth >= MAX_DEPTH) return L;
9

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

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

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

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

分直接与间接的路径追踪#

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

具体流程如图：

路径追踪 shader 函数的递归流程。

• 每个像素，相机只发出一条主光线，打到第一个交点 P1。
• 调用 shade(P1, ...)。这个函数内部会再发出两类新光线：一条固定射向光源的shadow ray (用来算P1这一点收到的直接光),一条按 BRDF 采样出来的弹射方向 (用来算 P1 这一点收到的间接光)。
• 弹射方向那条光线打到 P2，shade(P1) 内部递归调用自己 shade(P2, ...)，得到一个返回值。
• shade(P2) 内部又重复第2步<自己发一条shadow> ray,自己再决定要不要继续弹射到P3……
• 每一层的返回值，都乘上对应的 BRDF、cosθ、pdf 等权重，一路加回给上一层，最后回到 P1，得到 P1 这一个像素样本的最终颜色。

1
Eigen::Vector3f shade(const Intersection& hit, const Eigen::Vector3f& wo, const Scene& s)
2
{
3
    // 终止条件之一：交点是光源本身，直接返回自发光即可
4
    if (hit.material->has_emission()) return hit.material->emission;
5

6
    // —— 直接光照:从光源上采一点，不参与递归 ——
7
    Eigen::Vector3f L_dir = Eigen::Vector3f::Zero();
8
    Intersection light_hit;
9
    float light_pdf;
10
    s.sample_light(light_hit, light_pdf);   // 在场景所有光源(面积)上按面积均匀采一个点,拿到这个点和它的采样概率密度 light_pdf
11
    Eigen::Vector3f ws = (light_hit.p - hit.p).normalized();
12
    float dist2 = (light_hit.p - hit.p).squaredNorm();
13
    // shadow ray的可见性测试，若中间有遮挡，直接光贡献为0
14
    Ray to_light{hit.p + EPSILON * hit.n, ws};
15
    Intersection chk = s.intersect(to_light);
16
    if (chk.happened && (chk.p - light_hit.p).squaredNorm() < 1e-3f)
17
    {
18
        Eigen::Vector3f f = hit.material->eval(ws, wo, hit.n);
19
        float cos_t  = std::max(0.0f, hit.n.dot(ws));
20
        float cos_tl = std::max(0.0f, light_hit.n.dot(-ws));
21
        L_dir = light_hit.material->emission.cwiseProduct(f) * cos_t * cos_tl / dist2 / light_pdf;  // 把"面积域"的光源采样,转换成渲染方程需要的"立体角域"贡献,
22
    }
23

24
    // —— 间接光照:俄罗斯轮盘赌 + BRDF 采样 ——
25
    Eigen::Vector3f L_indir = Eigen::Vector3f::Zero();
26
    const float p_rr = 0.8f;    // 以 p_rr 的概率选择继续弹射
27
    if (random_float() < p_rr)
28
    {
29
        Eigen::Vector3f wi = hit.material->sample(wo, hit.n);
30
        float pdf = hit.material->pdf(wi, wo, hit.n);
31
        if (pdf > 1e-8f)
32
        {
33
            Ray next{hit.p + EPSILON * hit.n, wi};
34
            Intersection next_hit = s.intersect(next);
35
            if (next_hit.happened && !next_hit.material->has_emission())
36
            {
37
                Eigen::Vector3f f = hit.material->eval(wi, wo, hit.n);
38
                float cos_t = std::max(0.0f, hit.n.dot(wi));
39
                L_indir = shade(next_hit, -wi, s).cwiseProduct(f) * cos_t / pdf / p_rr;
40
            }
41
        }
42
    }
43
    return L;
44
}