真实感渲染：BSDF

前段时间写了一个 GPU BVH 构建器，顺手跑了个漫反射光追测了下效果。

虽然看着挺像那么一回事的，但还是有很多令人不太舒服的伪影，也不知道是代码的问题还是没有用重要性采样的原因，反正是懒得 Debug 了。

不过既然光追的框架已经完成了，也是时候升级成 BSDF 看看效果了。

这篇首先会简单的介绍一下 BSDF，随后回顾一下渲染方程和微表面理论，这篇文章不会涉及到具体的代码实现。

光的行为与BSDF

在图形学中，光的行为通常可以概括为三种基本形式：吸收、散射和发光。

• 吸收（Absorption）：指光在传播过程中，其能量逐渐被介质吸收而减弱的现象。在一些均匀介质中（例如常见的空气），光通常沿直线传播，但其强度会随着传播距离的增加而逐渐衰减。

• 散射（Scattering）：通常发生在参与介质中（例如牛奶，云朵等）。当光进入此类介质后，会不断与其中的微小颗粒发生相互作用，使光被重新分配到不同方向传播。该过程满足能量守恒，即散射后的总能量不超过散射前的能量。

• 发光（Emission）：发生在自发光物体中。这类物体能够向各个方向发射新的光，从而成为光的来源。

然而，BSDF（Bidirectional Scattering Distribution Function，双向散射分布函数） 中的“散射”与上述在参与介质中的散射有着本质区别。BSDF 实际上是一种用于描述表面光照行为的模型，它统一了 BRDF（Bi-directional Transmit Distribution Function，双向反射分布函数） 和 BTDF（Bi-directional Scattering Distribution Function， 双向透射分布函数）。因此，这里的“散射”并非指光在介质中被随机改变方向的过程，而是特指光在物体表面发生的反射与透射等相互作用。

我自己曾经对 BSDF 就存在过很长一段时间的误解，一度以为它描述的是参与介质中的散射过程，因此始终对其抱有一种“高深莫测”的敬畏之心。实际上，BSDF 仅用于刻画光在物体表面的反射与透射行为。

而真正涉及参与介质中散射效应时，则需要引入另一个更加复杂的主题——体积渲染（Volume Rendering）。在这一框架下，我们所使用的也不再是经典的渲染方程，而是更为一般的辐射传输方程（Radiative Transfer Equation, RTE）。当然，这篇文章并不会涉及这部分内容。

渲染方程

回顾一下渲染方程：

$$L_o(p, \omega_o) = L_e(p,\omega_o) + \int_{\Omega}L_i(p,\omega_i)f_s(p,\omega_i,\omega_o)(n\cdot \omega_i) \mathrm d \omega_i \tag{1}$$

其中，$L_o$​ 和 $L_i$​ 分别表示出射辐射亮度与入射辐射亮度，$L_e$​ 表示表面自身的发射项；$p$ 表示空间中的位置，$\omega_i$ 与 $\omega_o$​ 分别表示入射方向和出射方向。$f_s$​ 即 BSDF，用于刻画光从入射方向 $\omega_i$到出射方向 $\omega_o$​ 的能量分布关系，通俗来讲就是描述入射方向 $\omega_i$ 的光线有多少能量散射到出射方向 $\omega_o$ 上。

从物理意义上看，渲染方程描述的是：在点 $p$ 处，沿出射方向 $\omega_o$ 的辐射亮度，由自发光项与所有入射方向光经过表面相互作用后贡献的总和构成。换言之，它刻画了光在场景中传播与能量交换的基本规律。

需要注意的是，这里的 $f_s$​ 描述 宏观物体表面 的散射行为，由两部分组成，分别是 $f_r$​ 和 $f_t$​，并满足如下关系：

$$f_s = f_r + f_t$$

也就是说，BSDF 可以看作是 BRDF 与 BTDF 的统一表示形式：前者描述光在表面的反射行为，后者则描述光的透射行为。

为了推导 $f_s$ 的具体形式，我们需要借助 微表面模型（Microfacet Model） 来进行分析。

微表面模型

正如前文所述，在实际求解渲染方程时，我们通常是在 宏观表面 上进行计算。无论是现实世界中的物体表面，还是 3D 模型中的几何表面，在我们的建模与计算中，往往都被近似为 理想的平滑表面 。然而在微观尺度下，不同材质的表面实际上具有复杂的结构特征，呈现出各种细微的凹凸起伏。

显然，直接在微观层面求解渲染方程几乎是不可行的，这太复杂了。如果能够将这些复杂的微观结构行为，转化为一种简洁而有效的宏观表达形式，将极大地方便渲染方程的求解。

这正是微表面模型（Microfacet Model）的核心思想：将光在表面上的微观相互作用，映射为宏观层面的反射与透射行为 。

如图所示，微表面模型的目标，是将由黑色曲线表示的微观表面（Microsurface），等效为由青色曲线表示的宏观表面（Macrosurface），从而在保持物理合理性的同时，大幅简化计算过程。

显然，这种从微观表面到宏观表面的转化需要借助积分来完成。通常我们会假设微观表面由无数个极小的理想平面组成，而每一个这样的平面元都会发生理想的镜面反射与折射。

为了便于后续表述，我们先约定如下符号：

• $\mathbf{v}$：归一化后的向量；
• $\mathbf{n}$：宏观表面的法向；
• $\mathbf{m}$：微观表面的法向；
• $f_s(\mathbf{i},\mathbf{o},\mathbf{n})$：宏观表面的 BSDF，用于描述在法向为 $\mathbf{n}$ 的表面上，来自入射方向 $\mathbf{i}$ 的光，有多少能量被散射到出射方向 $\mathbf{o}$。需要注意的是，这里的表示方式与前文渲染方程略有不同：渲染方程中通常使用方向参数（如球坐标或方向向量）并显式依赖位置 $p$，而此处我们采用三维方向向量表示，并用法向 $\mathbf{n}$ 代替 $p$，因为在计算 $f_s$​ 时，主要关心的是该点处的表面方向；
• $f_s^m(\mathbf{i},\mathbf{o},\mathbf{n})$：微观表面的BSDF；
• $f_r,f_r^m,f_t,f_t^m$：分别表示宏观表面的 BRDF、微观表面的 BRDF，以及宏观与微观层面的 BTDF。

描述微观表面通常需要以下三个关键要素：

1. 微表面法线分布函数 $D(\mathbf{m})$：用于描述不同法向 $\mathbf{m}$ 的微表面在宏观表面上的分布情况；
2. 几何项（阴影-遮蔽函数） $G(\mathbf{i},\mathbf{o},\mathbf{m})$：用于刻画由于微表面之间的相互遮挡所产生的阴影与遮蔽效应；
3. 微观表面 BSDF $f_s^m$：描述单个微表面上光的反射与透射行为。

微表面法线分布函数 D

微表面法线分布函数 $D(\mathbf{m})$ 本质上是一个定义在单位球面上的概率密度函数，其单位为 $sr^{-1}$（即单位立体角的倒数）。它用于描述在宏观表面上，各个取向为 $\mathbf{m}$ 的微表面法线的分布情况。

具体而言，对于一块总面积为 $A$ 的微观表面，法向落在立体角 $\mathrm{d}\omega_m$​ 范围内、且方向为 $\mathbf{m}$ 的微表面所占的比例为 $D(\mathbf{m})\,\mathrm{d}\omega_m$。因此，对应的实际面积为：

$$A \cdot D(\mathbf{m})\, \mathrm{d}\omega_m$$

也就是说，$D(\mathbf{m})$ 刻画了不同法线方向在微观表面中的分布密度。

需要注意的是，微观表面的总面积通常大于宏观表面的面积，因此有：

$$1 \leq \int D(\mathbf{m})\, \mathrm{d}\omega_m$$

这一点其实是直观可理解的：在示意图中，用黑色曲线表示的微观表面由于存在起伏，其“实际长度”（或面积）显然大于用青色曲线表示的宏观表面。

但这似乎与概率密度函数的定义相矛盾——按照常规定义，概率密度函数在其定义域上的积分应当等于 1。这里的关键在于：$D(\mathbf{m})$ 并不是直接以“宏观面积”为基准定义的概率密度，而是描述微观表面法线在立体角上的分布，其度量基准与宏观表面并不一致。

为了得到与宏观表面相一致的概率分布，我们需要在积分中引入一个投影因子，将微观表面的面积投影到宏观表面上。通常这个投影项为 $(\mathbf{n} \cdot \mathbf{m})$，因此满足归一化条件的是：

$$\int_\Omega D(\mathbf{m})(\mathbf{n} \cdot \mathbf{m})\, \mathrm{d}\omega_m = 1$$

这样处理后，积分结果才对应于宏观表面上的面积比例分布，也符合概率密度函数的归一化要求。

换句话说，前面积分结果大于等于 1，本质上是因为我们是在宏观表面的参考系下，对微观表面（实际起伏的表面）上进行度量（即青色曲线为单位长度）。如果我们是在微观表面的参考系下对微观表面进行度量，那$\int D(\mathbf{m})\, \mathrm{d}\omega_m$ 显然就是1（即黑色曲线为单位长度），而用青色曲线表示的宏观表面此时长度应该小于1。

几何项（阴影-遮蔽函数） G

几何项（阴影-遮蔽函数）$G(\mathbf{i},\mathbf{o},\mathbf{m})$ 用于刻画微表面之间由于相互遮挡而产生的阴影与遮蔽效应。通俗来说，它描述的是：从入射方向$\mathbf{i}$ 到达微表面的光，有多少能够顺利地从出射方向 $\mathbf{o}$ 离开，而不是被表面的凹凸结构所阻挡。

如图所示，微表面之间可能发生多种遮挡情况，例如入射光被其他微表面遮挡（阴影），或者反射后的光在出射路径上被阻挡（遮蔽）。因此，并不是所有入射光都能有效参与最终的出射贡献。

本质上来看，几何项的核心作用在于对这类可见性问题进行建模：即在给定入射与出射方向的情况下，只有一部分光线能够真正完成从入射到出射的传播过程，其余部分则会被微观结构所遮挡。

计算宏观 BSDF

有了微表面法线分布函数 $D(\mathbf{m})$和几何项（阴影-遮蔽函数）$G(\mathbf{i},\mathbf{o},\mathbf{m})$，我们只需要对微观表面的无数微小理想镜面积分，就能得到渲染方程中的宏观BSDF $f_s$，即：

$$f_s(\mathbf{i},\mathbf{o},\mathbf{n})=\int\left|\frac{\mathbf{i}\cdot\mathbf{m}}{\mathbf{i}\cdot\mathbf{n}}\right|f_s^m(\mathbf{i},\mathbf{o},\mathbf{m})\left|\frac{\mathbf{o}\cdot\mathbf{m}}{\mathbf{o}\cdot\mathbf{n}}\right|G(\mathbf{i},\mathbf{o},\mathbf{m})D(\mathbf{m}) \mathrm d\mathbf{\omega}_m \tag{2}$$

其中$f_s^m$ 表示光在单个理想微表面上的相互作用，即理想的镜面反射与透射（遵循斯涅尔定律，Snell’s Law）。我们通常通过菲涅尔项（Fresnel Term） 来描述光在该微表面上发生反射与透射的比例：当入射角接近掠射角时，反射分量会显著增加，而透射分量相应减小。

在上述积分表达式中，还需要引入两个修正因子：$\left|\frac{\mathbf{i}\cdot\mathbf{m}}{\mathbf{i}\cdot\mathbf{n}}\right|$和$\left|\frac{\mathbf{o}\cdot\mathbf{m}}{\mathbf{o}\cdot\mathbf{n}}\right|$。其中，第一个因子用于将宏观表面接收到的光能量转换为微观表面上的入射能量；第二个因子则用于将微观层面的出射能量重新映射回宏观表面。

试想一束辐射亮度为$L$的光线打到宏观物体表面，根据基本光照关系，宏观表面单位面积所接收到的能量与$(\mathbf{i}\cdot\mathbf{n})$成正比，即为$L\times (\mathbf{i}\cdot\mathbf{n})$。

然而，在微表面模型中，我们关注的是具有法向 $\mathbf{n}$ 的微小平面元，因此对应的入射能量应当与 $(\mathbf{i}\cdot\mathbf{m})$ 成正比，即为$L\times (\mathbf{i}\cdot\mathbf{m})$

因此，从宏观到微观的转换过程中，就需要乘以修正因子 $\left|\frac{\mathbf{i}\cdot\mathbf{m}}{\mathbf{i}\cdot\mathbf{n}}\right|$​​，以完成入射能量的正确映射。

从微观层面将出射量重新转换回宏观层面时，也需要一个投影修正因子。其原因在于，辐射亮度（radiance）的定义依赖于单位投影面积：

$$L=\frac{\mathrm d^2 \Phi}{\mathrm d \omega \mathrm d A \cos\theta}$$

因此，对于同一束沿 $\mathbf o$ 方向出射的光，相对于宏观法向 $\mathbf n$ 与微观法向 $\mathbf m$ 而言，其对应的辐射亮度分别可写为

$$L_{marco}=\frac{\mathrm d^2 \Phi}{\mathrm d \omega \,\mathrm d A \left|\mathbf{o}\cdot\mathbf{n}\right|}$$

$$L_{mirco}=\frac{\mathrm d^2 \Phi}{\mathrm d \omega \,\mathrm d A \left|\mathbf{o}\cdot\mathbf{m}\right|}$$

于是两者满足关系

$$\frac{L_{marco}}{L_{mirco}}= \left|\frac{\mathbf{o}\cdot\mathbf{m}}{\mathbf{o}\cdot\mathbf{n}}\right|$$

因此，当我们把以微表面为参考定义的出射辐射亮度换算为以宏观表面为参考的出射辐射亮度时，需要乘上修正因子$\left|\frac{\mathbf{o}\cdot\mathbf{m}}{\mathbf{o}\cdot\mathbf{n}}\right|$。

微表面镜面 BSDF

正如前文所述，在微表面模型中，每一个微小表面都被视为一块理想镜面，因此只会发生理想的镜面反射与折射。这意味着：对于给定的入射方向 $\mathbf{i}$，只有唯一确定的反射方向或折射方向具有非零的辐射亮度，而其他所有方向的贡献均为 0。

为了刻画这种“只在特定方向上有贡献”的分布，我们通常引入 狄拉克函数（Dirac Delta Function） 进行描述：

$$\left. \delta(t) = \left\{ \begin{array} {cc}\infty & t = 0\\ 0 & \mathrm{otherwise} \end{array}\right.\right.$$

该函数满足如下归一化性质：

$$\int\delta(t)dt = 1$$

因此，它非常适合用于表示理想镜面中“能量集中在唯一方向”的特性。我们可以定义 义$\delta_{\omega_o}(\mathbf{s},\mathbf{o})$，表示只有当镜面反射（或折射）方向 $\mathbf{s}$ 与实际出射方向 $\mathbf{o}$ 完全一致时，该项才有贡献。

进一步地，引入反射（或透射）系数 $\rho$（即出射辐射亮度与入射辐射亮度之比），即可得到微表面镜面 BSDF 的一般形式：

$$f_s^m(\mathbf{i},\mathbf{o},\mathbf{n}) = \rho\frac{\delta_{\omega_o}(\mathbf{s},\mathbf{o})}{\left|\mathbf{o} \cdot\mathbf{m}\right|}$$

那么这里为什么又要除以 $\left|\mathbf{o} \cdot\mathbf{m}\right|$ 呢？其根本原因仍然来源于辐射亮度的定义：辐射亮度是相对于单位投影面积定义的，而 BSDF 的定义需要满足能量守恒关系。

因此我们这里除以 $\left|\mathbf{o} \cdot\mathbf{m}\right|$ 就是为了将表面接受的辐射亮度恢复到入射光线实际的辐射亮度。

上述表达还可以进一步变形。在原始形式中，我们需要通过判断$\mathbf{s}=\mathbf{o}$ 来确定某一出射方向是否有贡献；然而，直接用微表面法向 $\mathbf{m}$ 与入射方向 $\mathbf{i}$ 来描述这一关系相对比较复杂。因此，通常我们会引入半程向量（half-vector） $\mathbf{h}$ 来简化表达。

在镜面反射情况下，$\mathbf{h}$ 恰好等于微表面法向 $\mathbf{m}$，且可以表示为入射方向 $\mathbf{i}$ 和 $\mathbf{o}$ 的函数（例如镜面反射半程向量$\mathbf{h} = \frac{\mathbf{i} + \mathbf{o}}{\|\mathbf{i} + \mathbf{o}\|}$​）。于是，我们可以通过判断 $\mathbf{h}=\mathbf{m}$ 来替代原先的判断条件。

基于这一变换，微表面镜面 BSDF 可改写为：

$$f_s^m(\mathbf{i},\mathbf{o},\mathbf{n}) = \rho\frac{\delta_{\omega_o}(\mathbf{h(\mathbf{i},\mathbf{o})},\mathbf{m})}{\left|\mathbf{o} \cdot\mathbf{m}\right|}\left\|\frac{\partial\omega_{\mathbf{h}}}{\partial\omega_{\mathbf{o}}}\right\| \tag{3}$$

公式（3）引入了一个雅可比行列式（Jacobian determinant） $\left\|\frac{\partial \omega_h}{\partial \omega_o}\right\|$，因为我们对积分变量进行了变换：

$$\begin{align*} &\int \rho\frac{\delta_{\omega_o}(\mathbf{h(\mathbf{i},\mathbf{o})},\mathbf{m})}{\left|\mathbf{o} \cdot\mathbf{m}\right|}\left\|\frac{\partial\omega_{\mathbf{h}}}{\partial\omega_{\mathbf{o}}}\right\| \,\mathrm d\omega_o\\ = & \int \rho\frac{\delta_{\omega_o}(\mathbf{h(\mathbf{i},\mathbf{o})},\mathbf{m})}{\left|\mathbf{o} \cdot\mathbf{m}\right|}\ \mathrm d\omega_h \end{align*}$$

现在只要求解雅可比行列式 $\left\|\frac{\partial\omega_{\mathbf{h}}}{\partial\omega_{\mathbf{o}}}\right\|$，我们就能得到完整的$f_s^m$了。在这里反射项 $f_r^m$ 与透射项 $f_t^m$ 需要分别计算，因为它们的半程向量计算公式并不相同。

雅可比行列式的求解有点麻烦（还是太菜了^ _ ^），所以就直接给出结果了。

反射项的雅可比行列式为：

$$\left\|\frac{\partial\omega_{\mathbf{h_r}}}{\partial\omega_{\mathbf{o}}}\right\| =\frac{1}{4\left\|\mathbf{o} \cdot\mathbf{h_r}\right\|}$$

透射项的雅可比行列式为：

$$\left\|\frac{\partial\omega_{\mathbf{h_t}}}{\partial\omega_{\mathbf{o}}}\right\| =\frac{\eta_{o}^{2}\left|\mathbf{o}\cdot\mathbf{h_{t}}\right|}{\left(\eta_{i}(\mathbf{i}\cdot\mathbf{h_{t}})+\eta_{o}(\mathbf{o}\cdot\mathbf{h_{t}})\right)^{2}}$$

分别代入公式（3），并将 $\rho$ 替换为菲涅尔项 $F$ ，可以得到反射项 $f_r^m$ 与透射项 $f_t^m$ ：

$$\begin{align*} f_r^m(\mathbf{i},\mathbf{o},\mathbf{m}) & =F(\mathbf{i},\mathbf{m})\frac{\delta_{\mathbf{\omega}_m}(\mathbf{h}_\mathbf{r},\mathbf{m})}{4(\mathbf{i}\cdot\mathbf{h}_\mathbf{r})^2} \\ f_t^m(\mathbf{i},\mathbf{o},\mathbf{m})&=(1-F(\mathbf{i},\mathbf{m}))\frac{\mathbf{\delta}_{\mathbf{\omega}_m}(\mathbf{h}_\mathrm{t},\mathbf{m})\mathbf{\eta}_o^2}{(\mathbf{\eta}_i(\mathbf{i}\cdot\mathbf{h}_\mathrm{t})+\mathbf{\eta}_o(\mathbf{o}\cdot\mathbf{h}_\mathrm{t}))^2} \end{align*}$$

将微表面的反射项 $f_r^m$ 与透射项 $f_t^m$ 回代到公式（2），积分后就能得到宏观表面的 BSDF 了:

$$\begin{align*} f_s(\mathbf{i},\mathbf{o},\mathbf{m}) &=f_r(\mathbf{i},\mathbf{o},\mathbf{m})+f_t(\mathbf{i},\mathbf{o},\mathbf{m}) \\ f_r(\mathbf{i},\mathbf{o},\mathbf{n})&=\frac{F(\mathbf{i},\mathbf{h}_\mathbf{r})G(\mathbf{i},\mathbf{o},\mathbf{h}_\mathbf{r})D(\mathbf{h}_\mathbf{r})}{4\left|\mathbf{i}\cdot\mathbf{n}\right|\left|\mathbf{o}\cdot\mathbf{n}\right|} \\ f_t\left(\mathbf{i},\mathbf{o},\mathbf{n}\right)&=\frac{\left|\mathbf{i}\cdot\mathbf{h}_\mathrm{t}\right|\left|\mathbf{o}\cdot\mathbf{h}_\mathrm{t}\right|}{\left|\mathbf{i}\cdot\mathbf{n}\right|\left|\mathbf{o}\cdot\mathbf{n}\right|}\frac{\mathbf{\eta}_o^2\left(1-F(\mathbf{i},\mathbf{h}_\mathrm{t})\right)G(\mathbf{i},\mathbf{o},\mathbf{h}_\mathrm{t})D(\mathbf{h}_\mathrm{t})}{\left(\mathbf{\eta}_i(\mathbf{i}\cdot\mathbf{h}_\mathrm{t})+\mathbf{\eta}_o(\mathbf{o}\cdot\mathbf{h}_\mathrm{t})\right)^2} \end{align*}$$

其中 $\eta_i$ 和 $\eta_o$ 分别是入射光线和透射光线的折射率。

选择 D、G、F 项

D、G 和 F 有很多可选形式，比如LearnOpenGL中介绍的经验公式。

---

• 微表面法线分布函数 D：

  $$NDF_{GGXTR}(n,h,\alpha)=\frac{\alpha^2}{\pi((n\cdot h)^2(\alpha^2-1)+1)^2}$$

  其中$\alpha$是粗糙度。

  float D_GGX_TR(vec3 N, vec3 H, float a) 
  { 
   float a2 = a*a; 
   float NdotH = max(dot(N, H), 0.0); 
   float NdotH2 = NdotH*NdotH; 
   float nom = a2; 
   float denom = (NdotH2 * (a2 - 1.0) + 1.0); 
   denom = PI * denom * denom; 
   return nom / denom; 
  }

---

• 几何项（阴影-遮蔽函数） G：

  $$\begin{align*} G_{\mathrm{SchlickGG}X}(n,v,k) & =\frac{n\cdot v}{(n\cdot v)(1-k)+k} \\ k & =\frac{(\alpha+1)^2}{8} \\ G(n,v,l,k)&=G_{sub}(n,v,k)G_{sub}(n,l,k) \end{align*}$$

  其中$\alpha$是粗糙度。

  float GeometrySchlickGGX(float NdotV, float k)
  {
  	float nom   = NdotV;
  	float denom = NdotV * (1.0 - k) + k;
  	return nom / denom;
  }
  
  float GeometrySmith(vec3 N, vec3 V, vec3 L, float k)
  {
  	float NdotV = max(dot(N, V), 0.0);
  	float NdotL = max(dot(N, L), 0.0);
  	float ggx1 = GeometrySchlickGGX(NdotV, k);
  	float ggx2 = GeometrySchlickGGX(NdotL, k);
  	return ggx1 * ggx2;
  }

---

• 菲涅尔方程 F：

  $$F_{Schlick}(h,v,F_0)=F_0+(1-F_0)(1-(h\cdot v))^5$$

  其中$F_0$是基础反射率。

  vec3 F0 = vec3(0.04); 
  F0 = mix(F0, surfaceColor.rgb, metalness);
  
  vec3 fresnelSchlick(float cosTheta, vec3 F0) 
  { 
   return F0 + (1.0 - F0) * pow(1.0 - cosTheta, 5.0); 
  }

采样

这一部分涉及重要性采样（Importance Sampling） 的内容。对于给定的一条入射方向 $\mathbf{i}$，如果我们能够生成出射方向 $\mathbf{o}$，使其概率密度分布 $p_o(\mathbf{o})$ 与

$$f_s(\mathbf {i},\mathbf {o},\mathbf {n})(\mathbf {o} \cdot \mathbf {n})$$

成正比，就可以显著降低采样方差。

在路径追踪中，对应样本的权重为：

$$\mathrm{weight}(\mathbf{o}) = \frac{f_s(\mathbf{i}, \mathbf{o}, \mathbf{n})\, |\mathbf{o}\cdot\mathbf{n}|}{p_o(\mathbf{o})} \tag{4}$$

如果满足

$$p_o(\mathbf{o}) \propto f_s(\mathbf{i}, \mathbf{o}, \mathbf{n})\, |\mathbf{o} \cdot \mathbf{n}|$$

即

$$p_o(\mathbf{o}) = c \cdot f_s(\mathbf{i}, \mathbf{o}, \mathbf{n})\, |\mathbf{o} \cdot \mathbf{n}|$$

那么所有样本的权重将为常数，此时方差为 0（理想情况）。

然而，直接在 $\mathbf{o}$ 上进行采样通常是困难的，因为 $p_o(\mathbf{o})$ 的解析形式往往较为复杂。因此，更常见的做法是：先对微表面法向 $\mathbf{m}$ 进行采样，再通过反射或折射关系得到 $\mathbf{o}$。

微表面法向的概率密度分布为：

$$p_m(\mathbf{m}) = D(\mathbf{m})\, |\mathbf{m} \cdot \mathbf{n}|$$

随后，通过变量变换（由 $\mathbf{m} \rightarrow \mathbf{o}$ ）引入雅可比行列式，将 $p_m$​ 转换为 $p_o$ ​：

$$p_o(\mathbf{o}) = p_m(\mathbf{m})\, \left\|\frac{\partial \omega_h}{\partial \omega_o}\right\|$$

将反射项或透射项分别代入公式（4），最终可以得到各自的权重表达：

$$\mathrm{weight_r}(\mathbf{o})=\frac{f_r(\mathbf{i},\mathbf{o},\mathbf{n})\left|\mathbf{o}\cdot\mathbf{n}\right|}{p_o(\mathbf{o})}=F(\mathbf{i},\mathbf{h}_\mathbf{r})\frac{\left|\mathbf{i}\cdot\mathbf{m}\right|G(\mathbf{i},\mathbf{o},\mathbf{m})}{\left|\mathbf{i}\cdot\mathbf{n}\right|\left|\mathbf{m}\cdot\mathbf{n}\right|}$$

$$\mathrm{weight_t}(\mathbf{o})=\frac{f_t(\mathbf{i},\mathbf{o},\mathbf{n})\left|\mathbf{o}\cdot\mathbf{n}\right|}{p_o(\mathbf{o})}=(1 - F(\mathbf{i},\mathbf{h}_\mathbf{t}))\frac{\left|\mathbf{i}\cdot\mathbf{m}\right|G(\mathbf{i},\mathbf{o},\mathbf{m})}{\left|\mathbf{i}\cdot\mathbf{n}\right|\left|\mathbf{m}\cdot\mathbf{n}\right|}$$

虽然这种转换看着挺简单的，但是实际操作起来有着不少讲究。为了生成符合 $p_m$​ 分布的样本，通常需要先对 $p_m$​ 进行积分，得到其累积分布函数（CDF），再通过求逆变换 $cdf^{-1}(X)$ 来实现采样。此外，还需要一种能够生成均匀随机序列的方法，将其作为输入传入 $cdf^{-1}(X)$ ，从而得到最终的采样方向。

这部分想要一下子搞清楚还挺不容易，所以一开始我并没有采用准确的 $p_m$ 作为采样函数，而是用半球均匀采样替代。

参考文献

Microfacet Models for Refraction through Rough Surfaces
Physically Based Shading at Disney
Background: Physics and Math of Shading
PBR - LearnOpenGL CN
PBRT: Transforming between Distributions