光线传输方程

基本形式

之前我们已经见过了反射方程

\[ L_o(\mathbf{p},\boldsymbol{\omega}_o)=\int_{H^2}f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)L_i(\mathbf{p},\boldsymbol{\omega}_i)\cos\theta_i\:d\Omega_i. \]

如果我们将描述表面反射的双向反射分布函数 \(f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)\) 替换为双向散射分布函数 \(f\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)\)，并加入自发光辐照率 \(L_e(\mathbf{p},\boldsymbol{\omega}_o)\)，那么就得到了如下含自发光项的散射方程

\[ L_o(\mathbf{p},\boldsymbol{\omega}_o)=L_e(\mathbf{p},\boldsymbol{\omega}_o)+\int_{S^2}f\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)L_i(\mathbf{p},\boldsymbol{\omega}_i)\left|\cos\theta_i\right|d\Omega_i. \]

注意到在上式中，\(L_e(\mathbf{p},\boldsymbol{\omega}_o)\) 和 \(f\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)\) 都是已知量．如果我们能将函数 \(L_i\) 用函数 \(L_o\) 表示出来，那么我们将得到一个只含 \(L_o\) 的积分方程．为了实现这一点，我们需要引入 光线投射函数（ray-casting function）\(t(\mathbf{p},\boldsymbol{\omega})\)．假设在 \(\mathbf{p}\) 点向 \(\boldsymbol{\omega}\) 方向发出一条射线，射线第一次与场景中某个表面相交的位置就记作 \(t(\mathbf{p},\boldsymbol{\omega})\)．于是，我们得到

\[ L_i(\mathbf{p},\boldsymbol{\omega}_i)=L_o(t(\mathbf{p},\boldsymbol{\omega}_i),-\boldsymbol{\omega}_i). \]

为了处理射线不与场景中任何表面相交的特殊情况，我们规定此时 \(t(\mathbf{p},\boldsymbol{\omega})\) 取特殊值 \(\Lambda\)，并且对任意 \(\boldsymbol{\omega}\)，有 \(L_o(\Lambda,\boldsymbol{\omega})=0\)．

利用上面的关系式，并将 \(L_o(\mathbf{p},\boldsymbol{\omega}_o)\) 简记为 \(L(\mathbf{p},\boldsymbol{\omega}_o)\)，由此导出的方程就称为 光线传输方程（light transport equation）

\[ L(\mathbf{p},\boldsymbol{\omega}_o)=L_e(\mathbf{p},\boldsymbol{\omega}_o)+\int_{S^2}f\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)L(t(\mathbf{p},\boldsymbol{\omega}_i),-\boldsymbol{\omega}_i)\left|\cos\theta_i\right|d\Omega_i. \]

光线传输方程也叫做称渲染方程（rendering equation），它是关于函数 \(L(\mathbf{p},\boldsymbol{\omega}_o)\) 的积分方程，描述了均衡状态下物体表面任意一点、沿任意方向的辐照率．实际上，我们可以也可以得到场景中任意一点、沿任意方向的辐照率，因为借助前面导出的关系式

\[ L(\mathbf{p},\boldsymbol{\omega})=L(t(\mathbf{p},\boldsymbol{\omega}),-\boldsymbol{\omega}), \]

我们可以将场景中的辐照率转化为物体表面的辐照率．

表面形式

在上面的光线传输方程中，场景中物体表面的几何信息全部隐藏在了光线投射函数 \(t(\mathbf{p},\boldsymbol{\omega})\) 里．为了让这些信息显式地表达出来，我们将导出表面形式的光线传输方程．

之前在对某点各个方向上的入射光做积分时，我们是在单位球上做积分．考虑到某点接收的入射光实际上来自于其他物体表面各点的出射光，我们也可以把单位球上的积分转化为其他物体表面上的积分．如果我们要用蒙特卡洛积分法估计该积分，那么我们的采样区域将从单位球转到其他物体的表面．

下面引入一些记号．如果从 \(\mathbf{p}^\prime\) 点沿 \(\boldsymbol{\omega}\) 方向出发的射线经过了 \(\mathbf{p}\) 点，我们记

\[ L(\mathbf{p}^\prime\to\mathbf{p}):=L(\mathbf{p}^\prime,\boldsymbol{\omega})=L\left(\mathbf{p}^\prime,\frac{\mathbf{p}-\mathbf{p}^\prime}{\Vert\mathbf{p}-\mathbf{p}^\prime\Vert}\right). \]

如下图所示，

图1 - \(f(\mathbf{p}^{\prime\prime}\to\mathbf{p}^\prime\to\mathbf{p})\) 示意图

如果 \(t(\mathbf{p}^\prime,\boldsymbol{\omega}_o)=\mathbf{p}\)，\(t(\mathbf{p}^\prime,\boldsymbol{\omega}_i)=\mathbf{p}^{\prime\prime}\)，我们记

\[ f(\mathbf{p}^{\prime\prime}\to\mathbf{p}^\prime\to\mathbf{p}):=f\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)=f\left(\mathbf{p}^\prime,\frac{\mathbf{p}-\mathbf{p}^\prime}{\Vert\mathbf{p}-\mathbf{p}^\prime\Vert},\frac{\mathbf{p}^{\prime\prime}-\mathbf{p}^\prime}{\Vert\mathbf{p}^{\prime\prime}-\mathbf{p}^\prime\Vert}\right). \]

设 \(V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)\) 是可视性函数，当 \(\mathbf{p}\) 与 \(\mathbf{p}^{\prime}\) 两点之间没有遮挡、相互可见时，即

\[ t\left(\mathbf{p},\frac{\mathbf{p}^\prime-\mathbf{p}}{\Vert\mathbf{p}^\prime-\mathbf{p}\Vert}\right)=\mathbf{p}^\prime,\quad t\left(\mathbf{p}^\prime,\frac{\mathbf{p}-\mathbf{p}^\prime}{\Vert\mathbf{p}-\mathbf{p}^\prime\Vert}\right)=\mathbf{p} \]

时，\(V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)\) 取值为 \(1\)；否则 \(V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)\) 取值为 \(0\)．

设几何耦合项（geometric coupling term）为可视性函数与积分区域转换映射的 Jacobi 行列式的乘积，即

\[ G\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)=V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right) \frac{|\cos \theta|\left|\cos \theta^{\prime}\right|}{\left\|\mathbf{p}-\mathbf{p}^{\prime}\right\|^{2}}, \]

其中 \(\theta\) 为 \(\mathbf{p}\) 点处法向量与 \(\mathbf{p}^\prime-\mathbf{p}\) 的夹角， \(\theta^\prime\) 为 \(\mathbf{p}^\prime\) 点处法向量与 \(\mathbf{p}-\mathbf{p}^\prime\) 的夹角．

当 \(V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)=1\) 时，光线传输方程可以重写成场景中所有物体表面 \(A\) 上的积分

\[ L\left(\mathbf{p}^{\prime} \rightarrow \mathbf{p}\right)=L_{\mathrm{e}}\left(\mathbf{p}^{\prime} \rightarrow \mathbf{p}\right)+\int_{A} f\left(\mathbf{p}^{\prime \prime} \rightarrow \mathbf{p}^{\prime} \rightarrow \mathbf{p}\right) L\left(\mathbf{p}^{\prime \prime} \rightarrow \mathbf{p}^{\prime}\right) G\left(\mathbf{p}^{\prime \prime} \leftrightarrow \mathbf{p}^{\prime}\right)\: d A\left(\mathbf{p}^{\prime \prime}\right). \]

这就是表面形式的光线传输方程．

路径积分形式

表面形式的光线传输方程给出了 \(L\left(\mathbf{p}_{i+1} \rightarrow \mathbf{p}_i\right)\) 关于 \(i\) 的递推关系

\[ \begin{aligned} &L\left(\mathbf{p}_{i+1} \rightarrow \mathbf{p}_i\right)\\ =\;&L_{\mathrm{e}}\left(\mathbf{p}_{i+1} \rightarrow \mathbf{p}_i\right)+\int_{A} f\left(\mathbf{p}_{i+2} \rightarrow \mathbf{p}_{i+1}\rightarrow \mathbf{p}_i\right) L\left(\mathbf{p}_{i+2}\rightarrow \mathbf{p}_{i+1}\right) G\left(\mathbf{p}_{i+2}\leftrightarrow \mathbf{p}_{i+1}\right)\: d A\left(\mathbf{p}_{i+2}\right). \end{aligned} \]

从 \(i=0\) 出发，先将等式右边的 \(L\left(\mathbf{p}_{2}\rightarrow\mathbf{p}_{1}\right)\) 用 \(L\left(\mathbf{p}_{3}\rightarrow \mathbf{p}_{2}\right)\) 表示，再将先将等式右边的 \(L\left(\mathbf{p}_{3}\rightarrow \mathbf{p}_{2}\right)\) 用 \(L\left(\mathbf{p}_{4}\rightarrow \mathbf{p}_{3}\right)\) 表示，以此类推，最终将得到如下等式

\[ \begin{aligned} L\left(\mathbf{p}_{1} \rightarrow \mathbf{p}_{0}\right)&=L_{\mathrm{e}} \left(\mathbf{p}_{1} \rightarrow \mathbf{p}_{0}\right) \\ &\quad+\int_{A} L_{\mathrm{e}}\left(\mathbf{p}_{2} \rightarrow \mathbf{p}_{1}\right) f\left(\mathbf{p}_{2} \rightarrow \mathbf{p}_{1} \rightarrow \mathbf{p}_{0}\right) G\left(\mathbf{p}_{2} \leftrightarrow \mathbf{p}_{1}\right)\: d A\left(\mathbf{p}_{2}\right) \\ &\quad+\int_{A} \int_{A} L_{\mathrm{e}}\left(\mathbf{p}_{3} \rightarrow \mathbf{p}_{2}\right) f\left(\mathbf{p}_{3} \rightarrow \mathbf{p}_{2} \rightarrow \mathbf{p}_{1}\right) G\left(\mathbf{p}_{3} \leftrightarrow \mathbf{p}_{2}\right)\\ &\quad\quad\times f\left(\mathbf{p}_{2} \rightarrow \mathbf{p}_{1} \rightarrow \mathbf{p}_{0}\right) G\left(\mathbf{p}_{2} \leftrightarrow \mathbf{p}_{1}\right) \:d A\left(\mathbf{p}_{3}\right) d A\left(\mathbf{p}_{2}\right)\\ &\quad+\cdots \end{aligned} \]

记

\[ T\left(\overline{\mathbf{p}}_{n}\right):=\prod_{i=1}^{n-1} f\left(\mathbf{p}_{i+1} \rightarrow \mathbf{p}_{i} \rightarrow \mathbf{p}_{i-1}\right) G\left(\mathbf{p}_{i+1} \leftrightarrow \mathbf{p}_{i}\right) \]

为路径 \(\overline{\mathbf{p}}_{n}=(\mathbf{p}_n,\mathbf{p}_{n-1},\cdots,\mathbf{p}_{0})\) 的吞吐量（throughput），又令

我们可以将 \(L\left(\mathbf{p}_{1}\to\mathbf{p}_{0}\right)\) 写成级数形式

\[ L\left(\mathbf{p}_{1} \rightarrow \mathbf{p}_{0}\right)=\sum_{n=1}^{\infty} P\left(\overline{\mathbf{p}}_{n}\right). \]

这被称作路径积分形式的光线传输方程．注意到等式右端不再包含 \(L\)，因此，我们实际上以级数形式显式地写出了光线传输方程的解．\(P\left(\overline{\mathbf{p}}_{1}\right)\) 表示从光源直接射向照相机的光，\(P\left(\overline{\mathbf{p}}_{2}\right)\) 表示经过一次散射的直接光，\(P\left(\overline{\mathbf{p}}_{3}\right),P\left(\overline{\mathbf{p}}_{4}\right),\cdots\) 表示经过多次散射的间接光．

路径追踪算法

路径追踪算法（path tracing algorithm）的目标是估计

\[ L\left(\mathbf{p}_{1} \rightarrow \mathbf{p}_{0}\right)=\sum_{n=1}^{\infty} P\left(\overline{\mathbf{p}}_{n}\right). \]

的值．需要解决的问题有两个：第一个是如何估计无穷求和，第二个是如何估计积分项 \(P\left(\overline{\mathbf{p}}_{n}\right)\)．

俄罗斯轮盘赌

无穷求和的困难可以由所谓的俄罗斯轮盘赌（Russian roulette）解决．

假设我们有办法得到任意一个 \(P\left(\overline{\mathbf{p}}_{n}\right)\) 的无偏估计量 \(\widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\)，那么我们可以立刻写出

\[ \sum\limits_{n=1}^{\infty} {P} \left(\overline{\mathbf{p}}_{n}\right) \]

的一个无偏估计量

\[ \sum_{n=1}^{\infty} \widehat{P} \left(\overline{\mathbf{p}}_{n}\right). \]

但问题是在现实中，我们只能获得有限项的 \(\widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\)．俄罗斯轮盘赌的想法，就是以一定概率 \(q\) 丢弃第 \(N\) 项之后的所有项之和 \(\sum\limits_{n=N}^{\infty} \widehat{P} \left(\overline{\mathbf{p}}_{n}\right)\)．而在没有被丢弃的情况下，放大 \(\sum\limits_{n=N}^{\infty} \widehat{P} \left(\overline{\mathbf{p}}_{n}\right)\) 作为补偿．

对于前 \(2\) 项 \(P\left(\overline{\mathbf{p}}_{1}\right),P\left(\overline{\mathbf{p}}_{2}\right)\)，我们并不采用俄罗斯轮盘赌，而是将它们直接计算出来．从第 \(3\) 项开始，我们执行俄罗斯轮盘赌．设 \(\eta_1\) 服从 0-1 分布，

\[ \mathrm{P}(\eta_1=k) = q^{k}(1-q)^{1-k},\quad k=0,1, \]

当 \(\eta_1\) 取 \(1\) 时，我们丢弃之后的项，只取前两项之和

\[ \sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right), \]

而当 \(\eta_1\) 取 \(0\) 时，我们将后面的项乘上一个缩放因子 \(\frac{1}{1-q}\)，即得

\[ \sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\frac{1}{1-q}\sum_{n=3}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right). \]

结合这两种情况，我们的统计量应该写成

\[ \widehat{Q}_3=\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\frac{1_{\eta_1=0}}{1-q}\sum_{n=3}^{\infty}\widehat{P}\left(\overline{\mathbf{p}}_{n}\right) \]

注意到

\[ \mathrm{E}\left[\frac{1_{\eta_1=0}}{1-q}\right]=\frac{\mathrm{P}(\eta_1=0)}{1-q}=\frac{1-q}{1-q}=1, \]

我们有

\[\mathrm{E}\left[ \sum_{n=1}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\right]=\mathrm{E}\left[\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\frac{1_{\eta_1=0}}{1-q}\sum_{n=3}^{\infty}\widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\right], \]

即 \(\widehat{Q}_3\) 也是 \(\sum\limits_{n=1}^{\infty}{P}\left(\overline{\mathbf{p}}_{n}\right)\) 的无偏估计量．

设 \(\eta_2\) 与 \(\eta_1\) 独立同分布，我们还可以用同样的方法继续构造出

\[ \widehat{Q}_4=\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\frac{1_{\eta_1=0}}{1-q}\left(\widehat{P}\left(\overline{\mathbf{p}}_{3}\right)+\frac{1_{\eta_2=0}}{1-q}\sum_{n=4}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\right). \]

显然，\(\mathrm{E}[\widehat{Q}_4]=\mathrm{E}[\widehat{Q}_3]\)，故 \(\widehat{Q}_4\) 也是一个无偏估计量．

一般地，若 \(\eta_1,\eta_2,\cdots\) 独立且都服从成功率为 \(q\) 的 0-1 分布，则

\[ \begin{aligned} \widehat{Q}_{N+2}&=\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\frac{1_{\eta_1=0}}{1-q}\left(\widehat{P}\left(\overline{\mathbf{p}}_{3}\right)+\frac{1_{\eta_2=0}}{1-q}\left(\widehat{P}\left(\overline{\mathbf{p}}_{4}\right)+\cdots+\frac{1_{\eta_N=0}}{1-q}\sum_{n=N+2}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\right)\right)\\ &=\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\frac{1_{\eta_1=0}}{1-q}\widehat{P}\left(\overline{\mathbf{p}}_{3}\right)+\frac{1_{\eta_1,\eta_2=0}}{(1-q)^2}\widehat{P}\left(\overline{\mathbf{p}}_{4}\right)+\cdots+\frac{1_{\eta_1,\eta_2,\cdots,\eta_N=0}}{(1-q)^N}\sum_{n=N+2}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right) \end{aligned} \]

是 \(\sum\limits_{n=1}^{\infty}{P}\left(\overline{\mathbf{p}}_{n}\right)\) 的无偏估计量．令 \(\widehat{Q}=\lim\limits_{N\to\infty}\widehat{Q}_{N}\)，则 \(\widehat{Q}\) 也是无偏估计量．

设

\[ \tau=\inf\{\:i\mid\eta_i=1\:\}, \]

即第 \(\tau+2\) 项及其后的所有项都被丢弃，我们可以将 \(\widehat{Q}\) 写作

\[ \widehat{Q}=\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\tau+1}\frac{1}{(1-q)^{n-2}} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)=\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\infty}\frac{1}{(1-q)^{n-2}} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)1_{\tau>n-2}. \]

不难求出，\(\tau\) 服从几何分布，其概率质量函数为

\[ \mathrm{P}(\tau=k)=(1-q)^{k-1} q,\quad k=1,2, \cdots. \]

于是，我们可以用计算验证 \(\widehat{Q}\) 的无偏性

\[ \begin{aligned} \mathrm{E}\left[\widehat{Q}\right]&=\sum_{n=1}^{2} P\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\infty}\frac{1}{(1-q)^{n-2}} P\left(\overline{\mathbf{p}}_{n}\right)\mathrm{P}(\tau>n-2)\\ &=\sum_{n=1}^{2} P\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\infty}\frac{1}{(1-q)^{n-2}} P\left(\overline{\mathbf{p}}_{n}\right)\sum_{k=n-1}^{\infty}(1-q)^{k-1}q\\ &=\sum_{n=1}^{2} P\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\infty}\frac{1}{(1-q)^{n-2}} P\left(\overline{\mathbf{p}}_{n}\right)(1-q)^{n-2}\\ &=\sum_{n=1}^{\infty} P\left(\overline{\mathbf{p}}_{n}\right). \end{aligned} \]

注意到 \(\mathrm{P}(\tau<\infty)=1\)，即 \(\widehat{Q}\) 为有限项的概率为 \(1\)，我们完全可以在现实中求出 \(\widehat{Q}\) 的估计值．

伪代码

路径追踪算法的伪代码如下：

shade(p, wo):
    # Contribution from the light source.
    uniformly sample the light at x’ (pdf_light = 1 / A);
    shoot a ray from p to x’;
    L_dir = 0.0;
    if the ray is not blocked in the middle:
        L_dir = L_i * f_r * cos θ * cos θ’ / |x’ - p|^2 / pdf_light;

    # Contribution from other reflectors.
    L_indir = 0.0;
    test Russian Roulette with probability P_RR;
    if passing Russian Roulette test:
        uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi);
        trace a ray r(p, wi);
        if ray r hit a non-emitting object at q:
            L_indir = shade(q, -wi) * f_r * cos θ / pdf_hemi / P_RR;
    
    return L_dir + L_indir;