计算机图形学-11：复杂光传输与外观

估计性质

渲染方程是随弹射次数增加而升维的积分，实际渲染只能从路径空间抽取有限样本。一个渲染算法首先要明确它的估计量 Estimator，再讨论它的视觉质量。

设真值为 $I$，用 $N$ 个样本得到的估计量为 $\hat I_N$。偏差是

$$\operatorname{Bias}(\hat I_N)=\mathbb E[\hat I_N]-I,$$

无偏表示对每个有限样本数，估计量的期望都等于 $I$；它不表示单次图像没有噪声。方差衡量单次结果围绕期望的波动。蒙特卡洛估计在独立同分布、方差有限时满足

$$\operatorname{Var}(\hat I_N)=\frac{\operatorname{Var}(X)}{N}, \qquad \operatorname{Std}(\hat I_N)\propto\frac1{\sqrt N}.$$

因此噪声减半通常需要约四倍样本。有限预算下还要比较单样本成本、误差在图像中的空间分布以及是否适合逐帧播放。均方误差把两类误差放在同一评价里：

$$\operatorname{MSE}(\hat I_N) =\mathbb E[(\hat I_N-I)^2] =\operatorname{Var}(\hat I_N)+\operatorname{Bias}(\hat I_N)^2.$$

有偏但能显著降低方差的方法在有限预算下可能更好；有偏且不一致的方法则会在样本增加后仍停在错误极限。

若

$$\hat I_N\xrightarrow[N\to\infty]{\mathsf P}I,$$

就称为一致估计 Consistent Estimator。一致性约束的是无限样本极限，无偏性约束的是期望；两者不能互换。

难采样光路

基础路径追踪的困难不只来自样本数量，还来自路径结构。某些光路对最终图像贡献很大，但从相机一侧随机生成它们的概率极低。一个典型场景是光源只照亮天花板的一小块区域，房间其他表面依靠间接光照明。从相机打到漫反射表面后，下一次随机弹射只有极少方向能回到那块高能量区域。

游泳池底的焦散是另一个例子。起伏水面像许多小透镜，把光线折射并聚焦到池底的网状亮纹。从相机端看，这条路径是 specular–diffuse–specular：先穿过水面到池底，再从池底经过水面连接光源。中间的漫反射点很容易采到，但两端的镜面方向由折射定律唯一确定，普通随机方向几乎不会命中。

图中上排为 BDPT、下排为 MLT，左列是室内茶壶，右列是水池焦散：BDPT 一侧仍有明显颗粒，MLT 一侧的高光和水面网状亮纹更完整。这里的画面对照只说明困难路径上的方差差异，不足以证明 MLT 在所有场景都更好。

BDPT 与 MLT

双向路径追踪 Bidirectional Path Tracing，BDPT 同时从光源和相机生成子路径，再连接两个子路径的端点。光源侧已经显式找到的路径不必依靠相机随机撞上，因此适合光源一侧结构复杂的场景。同一条完整路径可以由不同长度的相机、光源子路径组成，需要用 MIS 合并这些生成方式。连接会增加可见性测试与 PDF 计算；复杂光照下也可能因更容易找到有效路径而减少总耗时，不能只按单条路径成本判断效率。

理想镜面是 delta 分布，连接两个镜面点几乎不可能成功。遇到 delta 事件时应沿确定的反射或折射方向继续追踪，把连接留给可连接的非镜面端点。

Metropolis Light Transport，MLT 把完整路径当作马尔可夫链 Markov Chain 的状态：下一状态的生成只依赖当前路径。先找到一条贡献较高的路径，再对它做小幅变异，新的路径通常仍在高贡献区域附近。若当前路径为 $x$、候选为 $y$，一般提议分布下的 Metropolis–Hastings 接受率为

$$a(x\to y)=\min\left(1,\frac{\pi(y)q(y\to x)}{\pi(x)q(x\to y)}\right).$$

其中 $\pi$ 是与路径贡献成比例的目标分布，$q$ 是候选路径的提议分布；若提议对称，两个 $q$ 相消，接受率简化为 $\min(1,\pi(y)/\pi(x))$。拒绝的候选仍要把当前状态计入链，否则会改变停留时间比例。MLT 擅长焦散、狭缝等困难路径，但样本相关，像素间收敛速度不均，难以事先估计渲染多久才干净，动画中相邻帧还可能出现抖动。

光子与 VPL

Photon Mapping 分两个阶段。第一阶段从光源发射光子，使其反射或折射，直到到达漫反射表面，并保存位置、方向和通量。第二阶段从相机追踪路径，在着色点附近查询光子并做局部密度估计。

取着色点附近最近的 $N$ 个光子，估计它们覆盖的面积 $A$，密度近似为 $N/A$。以半径 $r$ 的圆邻域为例，一种简化形式是

$$L_o(\mathbf x,\boldsymbol\omega_o) \approx \frac1{\pi r^2} \sum_{k=1}^{N} f_r(\boldsymbol\omega_k,\boldsymbol\omega_o)\Delta\Phi_k.$$

真正的密度是 $\mathrm dN/\mathrm dA$，有限邻域计算的是 $\Delta N/\Delta A$，所以有限样本下有偏。发射光子越来越多时，固定邻居数 $N$ 覆盖的区域会变小，偏差会减小；但若 $N$ 始终固定，邻域内的随机样本数不增长，方差不一定消失，不能仅凭半径收缩就称为严格一致。若固定搜索半径而只增加光子，平滑范围不收缩，会保留系统性偏差。

邻居数小，结果噪声大；邻居数大，区域变宽，结果模糊。Vertex Connection and Merging，VCM 把 BDPT 的端点连接和光子式邻域合并放到同一框架中：端点能直接连接时使用连接，端点在同一表面上很近但无法精确连接时使用密度合并。

在局部二维表面上，若发射总光子数为 $M$、搜索半径为 $r_M$，一致收敛通常要求

$$r_M\to0, \qquad Mr_M^2\to\infty.$$

前一条件让空间平均造成的偏差消失，后一条件保证收缩的邻域中仍有越来越多的光子。固定半径只会收敛到被平滑过的结果；固定邻居数虽然让半径收缩，却还要随 $M$ 增长才能让方差也消失。

等价地，若用最近邻数 $N_M$ 表示核的样本量，需要同时满足

$$N_M\to\infty, \qquad \frac{N_M}{M}\to0.$$

前者控制统计方差，后者让邻域面积相对整个采样域缩小；固定 $N$ 只满足后者的直觉，不能保证估计量一致。

Instant Radiosity 或 many-light 方法从光源追踪短路径，把反弹点当作虚拟点光源（VPL），再用大量点光源的直接光照近似间接光。它速度快，但 VPL 接近着色点时几何项中的 $1/r^2$ 会造成缝隙亮斑，也不能很好处理窄的 glossy 高光。

光子映射把发射与查询分成两个阶段。光子是功率的带权样本；若 M 个等权样本代表总功率，每个初始权重与总功率/M 成比例：

发射阶段：
    按光源分布产生光子光线和初始通量权重
    沿光路寻找交点
    若命中漫反射表面，保存位置、入射方向和通量
    按材质继续反射或折射，更新权重
    通过吸收或 Russian Roulette 终止
查询阶段：
    从相机找到可见表面点
    查询表面邻域的光子
    累加 BRDF 加权的通量，并除以核面积

保存光子不等于立即停止传播，继续的光子才能记录后续间接光。使用面积核时还要避免把不同表面、墙的另一侧或不相容法线的光子混进同一个邻域。

参与介质

雾、云、烟和浑浊液体不是定义在表面上的材质，而是填充空间的参与介质 Participating Media。光在其中传播时可能被吸收，也可能被散射。吸收让能量消失，散射把能量导向新的方向，同时当前位置也会收到其他方向散射来的光。

吸收系数 Absorption Coefficient $\sigma_a$ 与散射系数 Scattering Coefficient $\sigma_s$ 都具有长度倒数的量纲，表示每单位传播距离的事件率。两者使光离开原传播方向，共同构成消光 Extinction，其系数为

$$\sigma_t=\sigma_a+\sigma_s.$$

沿距离 $s$ 的透射率满足

$$\frac{\mathrm dT}{\mathrm ds}=-\sigma_t(s)T(s),$$

解这个方程并令起点透射率为 1，得到 Beer–Lambert 定律：

$$T(\mathbf x,\mathbf y)= \exp\left(-\int_{\mathbf x}^{\mathbf y}\sigma_t(\mathbf s)\,\mathrm ds\right).$$

积分 $\tau=\int\sigma_t\,\mathrm ds$ 是无量纲的光学厚度 Optical Depth，$T=e^{-\tau}$。均匀介质中 $T(d)=e^{-\sigma_td}$，首次事件的累积分布为 $1-T(d)$，因此自由飞行距离 Free-Flight Distance 的概率密度和反演采样为

$$p(t)=\sigma_t e^{-\sigma_t t}, \qquad t=-\frac{\ln(1-\xi)}{\sigma_t}.$$

这里要求 $\sigma_t>0$；$\sigma_t=0$ 时没有介质事件，应直接穿过。若下一表面或介质边界距离为 $d_b$，样本 $t<d_b$ 才发生体事件，$t\ge d_b$ 表示无事件到达边界。无事件分支的概率正好是 $T(d_b)$；已经通过随机分支计算这一概率时，不应再把同一透射率乘一次。若确定性地穿过介质段，则需显式乘 $T(d_b)$。

Beer–Lambert 透射率来自一个局部概率平衡。传播到 $s$ 时尚未发生事件的比例为 $T(s)$，再前进 $\mathrm ds$ 后，损失比例约为 $\sigma_t(s)\mathrm ds$，所以

$$T(s+\mathrm ds)\approx T(s)[1-\sigma_t(s)\mathrm ds].$$

移项并取极限得到 $\mathrm dT/\mathrm ds=-\sigma_tT$，从起点积分才得到指数形式。第一次事件是吸收还是散射的条件概率分别为

$$P(\mathrm{scatter})=\frac{\sigma_s}{\sigma_t}, \qquad P(\mathrm{absorb})=\frac{\sigma_a}{\sigma_t}.$$

这两个概率和相函数承担不同职责：消光系数决定走多远，事件类型决定是否继续，phase function 决定继续时往哪走。

散射事件发生后，散射方向由相函数 $p(\boldsymbol\omega_i,\boldsymbol\omega_o)$ 给出：

$$p\ge0, \qquad \int_{S^2}p(\boldsymbol\omega_i,\boldsymbol\omega_o)\,\mathrm d\omega_o=1.$$

相函数与 BRDF 的分工相似，但它作用于体积内部，不带表面余弦因子。相函数可以是各向同性的 $1/(4\pi)$，也可以用 Henyey–Greenstein 等模型描述前向或后向散射。

Henyey–Greenstein 相函数用 $g\in(-1,1)$ 控制方向性。为避免与表面方向约定混淆，先用实际传播方向 $\mathbf d_i$、$\mathbf d_o$ 表示进入、离开散射点的光，令 $\mu=\mathbf d_i\cdot\mathbf d_o$。若入射参数 $\boldsymbol\omega_i$ 从点指向上一顶点，则应代入 $\mu=-\boldsymbol\omega_i\cdot\boldsymbol\omega_o$。常见形式为

$$p_{\mathrm{HG}}(\mu)= \frac{1-g^2}{4\pi(1+g^2-2g\mu)^{3/2}}.$$

$g=0$ 时是各向同性；$g>0$ 偏向前向散射，$g<0$ 偏向后向散射。采样时可以先在 $[0,1)$ 上抽取均匀随机数，根据相函数的累计分布反解 $\mu$，再均匀抽取方位角，并把局部方向旋转到当前入射方向的坐标系。采样坐标系的轴取传播方向 $\mathbf d_i$，此时 $g>0$ 才直接对应继续沿原方向传播。

对 $g\ne0$，若 $\xi_1\sim U[0,1)$，一种常见的反演采样写法是

$$\mu=\frac{1+g^2-\left(\frac{1-g^2}{1-g+2g\xi_1}\right)^2}{2g}, \qquad \phi=2\pi\xi_2,$$

其中 $\xi_2\sim U[0,1)$；$g=0$ 时改用 $\mu=1-2\xi_1$。先在以入射方向为轴的局部坐标中构造方向，再旋转回世界坐标。这里的“前向”必须以传播方向为准；若存储的是从散射点指向上一顶点的方向，就要相应取反，不能混用两个约定。

在位置 $\mathbf x$、方向 $\boldsymbol\omega_o$ 上，来自其他方向的 in-scattering 为

$$L_s(\mathbf x,\boldsymbol\omega_o)= \sigma_s(\mathbf x)\int_{S^2} p(\boldsymbol\omega_i,\boldsymbol\omega_o)L_i(\mathbf x,\boldsymbol\omega_i)\,\mathrm d\omega_i.$$

若介质自身发光，还要加入体发射项 $q_e$；沿光线的输运方程可写成

$$\frac{\mathrm dL}{\mathrm ds}=-\sigma_tL+q_e+ \sigma_s\int_{S^2}p(\boldsymbol\omega_i,\boldsymbol\omega_o)L_i(\boldsymbol\omega_i)\,\mathrm d\omega_i.$$

局部输运方程的衰减项包含吸收和 out-scattering，$q_e$ 表示发射，积分项表示 in-scattering。令 $q=q_e+L_s$，沿正向传播坐标的一段 $[0,d]$，解可写成

$$L(d)=T(0,d)L(0)+\int_0^d T(s,d)q(s)\,\mathrm ds.$$

第一项是到达终点的边界光，第二项把各位置产生的光乘以它到终点的透射率。相机反向追踪时则按相反方向评价对应线段，不能混用 $s$ 的朝向。

非均匀介质中 $\sigma_t(\mathbf x)$ 随位置变化，可以用 Woodcock tracking：以一个上界 $\bar\sigma_t$ 采样候选事件，再按 $\sigma_t(\mathbf x)/\bar\sigma_t$ 接受真实碰撞，拒绝时记作不改变路径的 null collision。

毛发与皮毛

一根飘散的头发更像曲线周围的细长圆柱，不能简单当作普通表面。Kajiya–Kay 是便宜的经验模型，把纤维切线周围的散射近似成高光和漫反射的组合。Marschner 模型把头发看成带鳞片表皮和含色素皮质的玻璃状圆柱，区分三条主要路径：

• $R$：在表面直接反射，基本保持光源颜色；
• $TT$：折射进入纤维，再折射穿出；
• $TRT$：进入后在内部反射一次，再穿出。

穿过皮质的路径会被色素吸收，因而产生带发色的高光。真实头发还会在许多根纤维之间发生多次散射，数十万根纤维使渲染成本迅速上升。

动物毛发的髓质（medulla）比人发大得多，内部结构会把光像参与介质一样散射。只用人发模型会得到过于透明、过于锐利的皮毛。双圆柱模型在外层皮质中加入内层髓质，并在 $TT$、$TRT$ 路径中增加经髓质散射的 $TT_s$、$TRT_s$，总共五个主要分量。

一个仓鼠渲染例子由约 60 万根毛发组成，使用 1024 samples per pixel 时每帧约需 36.9 分钟。双圆柱模型后来用于 2017 年上映的《War for the Planet of the Apes》，该片获 2018 年奥斯卡最佳视觉效果提名；也用于 2019 年上映的《The Lion King》，该片获 2020 年奥斯卡最佳视觉效果提名。数量级说明了毛发模型的困难：单根曲线的散射公式并不是主要瓶颈，纤维之间的多次传播和几何规模才是。

耳缘和身体轮廓可见细毛，内部整体亮度则来自大量纤维之间的传播。60 万根、1024 SPP 和 36.9 分钟每帧是这一场景的设置与耗时，说明细长几何和多次散射的共同代价，不是所有毛发渲染的固定成本。

颗粒材质

沙、盐、糖、香料和雪由许多离散颗粒堆积而成。近处可以看到单颗粒的轮廓、空隙和高光，远处却只看见沙丘或雪堆的整体形状。可以逐粒建模，也可以用颗粒大小、密度、成分比例和散射参数建立统计表示；随着观察距离变化，表示方式也要改变。

颗粒材质的通用高效表示仍然困难。Pixar 的短片《Piper》展示了沙子这类颗粒外观，但它不等同于一个简单表面 BRDF。

近处沙粒的轮廓、间隙和不同反光仍可分辨，远处逐渐呈现整体沙滩。表示方法必须随投影尺度改变，不能在所有距离上都简单使用一个平滑表面或逐粒几何。

次表面散射

BRDF 假设入射和出射发生在同一个表面点。玉石、皮肤、蜡、牛奶和耳廓中的光却可能从 $\mathbf x_i$ 进入，在内部散射后从另一个位置 $\mathbf x_o$ 离开。这个空间迁移需要 BSSRDF Bidirectional Surface Scattering Reflectance Distribution Function：

$$S(\mathbf x_i,\boldsymbol\omega_i,\mathbf x_o,\boldsymbol\omega_o).$$

出射辐亮度写成

$$L_o(\mathbf x_o,\boldsymbol\omega_o) =\int_A\int_{\Omega^+} S(\mathbf x_i,\boldsymbol\omega_i,\mathbf x_o,\boldsymbol\omega_o) L_i(\mathbf x_i,\boldsymbol\omega_i) \cos\theta_i\,\mathrm d\omega_i\,\mathrm dA_i.$$

当 $\mathbf x_i=\mathbf x_o$ 时，BSSRDF 退化为带有表面 Dirac delta 的 BRDF。BSSRDF 多出的面积积分正是内部空间传播的数学表达，其单位为 $\mathrm{m}^{-2}\mathrm{sr}^{-1}$；面积上的 delta 具有面积倒数单位。

更具体地，局部 BRDF 是

$$S(\mathbf x_i,\boldsymbol\omega_i,\mathbf x_o,\boldsymbol\omega_o) =f_r(\mathbf x_i,\boldsymbol\omega_i,\boldsymbol\omega_o) \delta_A(\mathbf x_i-\mathbf x_o),$$

其中 $\delta_A$ 是表面上的面积 Dirac delta。代入面积积分后，$\mathbf x_i$ 被压到 $\mathbf x_o$，就恢复普通 BRDF 反射方程。这说明 BSSRDF 的新增自由度是出射位置，而不是简单给 BRDF 加一层模糊。

偶极近似用一对点源近似强散射介质中的扩散：表面下方的真实源代表入射能量，表面上方的虚拟负源负责近似边界条件。只用 BRDF 的雕塑像石膏，加入 BSSRDF 后背光处也会透出光，才有大理石、玉石或皮肤的柔和过渡。

在扩散近似中，吸收和约化散射系数决定有效扩散距离。用 $\sigma_a$ 和 $\sigma_s'$ 表示吸收与约化散射，常写 $\sigma_t'=\sigma_a+\sigma_s'$、$\sigma_s'=(1-g)\sigma_s$，扩散系数为

$$D=\frac1{3\sigma_t'}.$$

点源到表面下方真实位置和表面上方虚拟位置的距离分别记为 $d_r,d_v$。偶极模型用两项扩散剖面叠加，并通过虚拟负源满足边界条件；材料的颜色和散射半径由各波长的 $\sigma_a,\sigma_s'$ 决定。薄片、分层和弱散射材料不能直接套用半无限介质的偶极假设。

布料

布料是纤维加捻、成股、成线后通过机织或针织形成的集合。若按观察尺度和预算区分，常见有三种表示：

1. 把织法的总体行为拟合成表面 BRDF，适合远景；
2. 在空间网格中记录纤维密度和方向，把布料当作参与介质；
3. 显式建模每根纤维，适合特写但计算量最大。

表面 BRDF 速度快，却不能真实表现绒毛轮廓和纱线空隙；参与介质保留统计透光；逐根纤维最完整但路径和求交数量巨大。The BFG 中的亚麻衣物说明最简单的表面模型在合适场景中仍有实用价值。

微观细节

标准微表面模型把 NDF 当作平滑的统计分布，因此车漆、鼠标外壳和金属可能显得过度规整。真实表面还有划痕、灰尘、磨痕和金属亮片。高分辨率法线图中，每个像素只覆盖极少数真正满足镜面反射条件的微面，直接路径追踪会严重 miss。一个约 $200\mathrm{k}\times200\mathrm{k}$ 的法线图例子在数小时内不能收敛，按样本量外推接近一个月。

对像素覆盖区域 $P$ 统计真实法线分布，得到 patch NDF：

$$D_P(\mathbf h)=\frac1{|P|}\int_P D_{\mathbf x}(\mathbf h)\,\mathrm d\mathbf x.$$

用 $D_P$ 替换光滑统计 NDF，可以在不逐个命中微面的情况下保留稀疏高光和 glint。footprint 变大时，许多微面被平均；footprint 变小时，离散峰和方向性更明显。

当微结构尺寸接近可见光波长，几何光学不再足够，需要波动光学。多条微观路径的复振幅为

$$E_j=A_j e^{i2\pi\ell_j/\lambda},$$

总强度必须先叠加复振幅再取模平方：

$$I=\left|\sum_jE_j\right|^2.$$

两条路径的强度展开为

$$I=\sum_jA_j^2+2\sum_{i<j}A_iA_j \cos\left(\frac{2\pi(\ell_i-\ell_j)}{\lambda}\right).$$

因此真正起作用的是光程差 $\ell_i-\ell_j$ 与波长的比值：相位接近时相长，相差约 $\pi$ 时相消。界面反射还可能加入额外相位；独立光源的相位在时间平均后可能消掉交叉项，不能把任意两束光都相干相加。光盘、金属薄膜、屏幕像素结构和餐具划痕都可能出现这种现象。

程序外观

程序化外观不预先存储完整纹理，而是在查询坐标 $(x,y,z)$ 时计算函数 $f(x,y,z)$。二维贴图适合表面，三维噪声则能让物体被切开后仍露出连续的内部结构。对归一化噪声 $n(\mathbf x)$ 做阈值，可得到二值材料：

$$m(\mathbf x)= \begin{cases} 1,&n(\mathbf x)>\tau,\\ 0,&n(\mathbf x)\le\tau. \end{cases}$$

Perlin、Worley、多频噪声和域扭曲可以生成木纹、石材、地形、水面和细胞状结构。重要的是区分两种工具链：shader 中按需查询的程序化外观不需要存储体积；Houdini 等工业工具也可以先生成程序化资产，再交给后续流程。

多频噪声常写成 fractal Brownian motion（fBm）：

$$n_{\mathrm{fbm}}(\mathbf x)=\sum_{k=0}^{K-1}a_k\,n(2^k\mathbf x), \qquad a_k\text{ 随频率递减}.$$

它把多个尺度的结构叠加起来；Worley noise 通过到随机特征点的距离形成细胞状纹理，域扭曲先改变输入坐标再查询原噪声。高频程序纹理若不按像素 footprint 过滤，远处会闪烁和走样，因此仍需限制最高频率或做预过滤。