渲染方程中的路径积分

路径追踪是第一个光线传输算法,并广泛应用于Realistic Rendering中。路径追踪将渲染方程的积分式转换成一个路径积分的形式,并且通过蒙特卡洛积分方法对每一项进行积分,从而得到一个统计上无偏的渲染结果。

渲染方程中的路径积分

路径追踪是第一个光线传输算法,并广泛应用于Realistic Rendering中。路径追踪将渲染方程的积分式转换成一个路径积分的形式,并且通过蒙特卡洛积分方法对每一项进行积分,从而得到一个统计上无偏的渲染结果。因为路径积分的无偏属性,路径积分也常常被用作为渲染的Ground Truth,用于对比其他的渲染方法。

渲染方程

要完全推导路径积分并且用代码表示路径积分的计算过程,我们需要从渲染方程入手。

首先需要明确的一点是,渲染的目标是计算最终进入眼睛(摄像机)的辐射亮度(Radiance),即单位面积,单位立体角上的光通量(辐射功率)。另外一个重要的概念是辐射照度,辐照度(Irradiance),即单位面积上光通量(辐射功率)。根据几何光学的线性假设,当光射到物体表面产生反射的时候,其反射方向的辐射亮度的微分,和入射方向上辐射照度的微分成线性关系。这两者比例常数由入射角$\omega_i$和反射角$\omega_o$和物体表面属性直接给出,我们把这个函数称之为BRDF(双向反射分布函数),$f(p, \omega_i, \omega_o)$。双向反射分布函数理解起来比较具有挑战性,一个简单的方法是将其辐照度和辐亮度类比为连续型概率的中的随机变量和累积概率函数,那么累积概率函数和随机变量的微分比值,即BRDF,就是概率密度函数。

$$dL_o(p, \omega_o) = dE(p, \omega_i) = f(p, \omega_i, \omega_o)L_i(p, \omega_i)cos\theta_id\omega_i$$

有了BRDF函数,我们就建立起了,对于物体表面上的一个点,点$p$,其反射的辐射亮度$L_o(p, \omega_o)$和入射的辐射亮度之间的关系,进一步的,如果我们考虑$p$上的自发光的辐射亮度$L_e(p, \omega_o)$,那我我们就可以建立起渲染方程

$$L_o(p, \omega_o) = L_e(p, \omega_o) + \int f(p, \omega_i, \omega_o)L_i(p, \omega_i)cos\theta_id\omega_i$$

如果我们只考虑不透明反射模型,而不考虑更复杂的折射和次表面反射,那么我们的积分范围即是围绕表面法线的半球积分。

方向积分形式和面积积分形式

从渲染方程中,对于$L_e$项是相对简单的,我们可以从物体表面材质定义获得。对于积分项,我们则需要知道每个方向的$L_i$。首先需要明确一点,在真空环境中,光在传播过程中,辐射亮度并不发生变化。对于$L_i$项来说,其就等于它来自方向的点$p'$的反射的辐射亮度。为了更好的表示这一点,我们用函数$t(p, w_i) = p'$表示。函数$t$指的是,从$p$点出发,沿着$w_i$方向,第一个相交的点,即为$p'$。由此,渲染方程可以改写为

$$L_o(p, \omega_o) = L_e(p, \omega_o) + \int f(p, \omega_i, \omega_o)L_o(t(p, \omega_i), \omega_i')cos\theta_id\omega_i$$

然而我们无法显式定义函数$t$,于是上面这个式子就无法显式定义,对于上面函数的展开,则必然包含着一个关于函数$t$的递归项无法消除。这使得我们仅仅可能通过递归黎曼求和的形式得到积分结果,而无法使用诸如蒙特卡洛近似的方式进行数值上的快速拟合。

为了显式改写该函数,我们重新定义一些记号,用于方更方便阅读。

首先是反射的辐射量度和入射的辐射亮度。在光线传输过程中提到,对于真空来说,辐射亮度在传输过程中不发生变化。那么对于一个点的入射辐射亮度,即是其来自方向反射(或者自发光)的辐射亮度。于是我们是用$L(p_1, p_0)$记号替代$L_i$和$L_o$,其表示点$p_1$到点$p_0$的辐射亮度。

同样的,BRDF方程记号我们也可以做相应的变化,我们使用$f(p_2, p_1, p_0)$表示从点$p_2$出发,击中表面点$p_1$,发生反射到达$p_0$的光路的BRDF。显然,该表达方式蕴涵了之前的立体角,从而在数学上是等价的。

经过上面两次记号的变化,在表达式上我们已经消除了立体角$\omega$的存在,然而微分项$d\omega_i$依然存在,我们需要将该微分项同样转换为关于点的表达式。由于

$$\omega = \frac{Acos\theta'}{d^2}$$

我们便将立体角的微分转换为对应顶点处的面积微分,为此,我们引入一个函数

$$G(p_0, p_1) = V(p_0, p_1)\frac{|cos\theta_0||cos\theta_1|}{||p_0-p_1||^2}$$

其中函数$V(p_0, p_1)$为可见性函数,如果两个点相互可见,则$V=1$,否则$V=0$。

于是渲染方程可以改写为

$$L(p_1, p_0) = L_e(p_1, p_0) + \int_A f(p_2, p_1, p_0)L(p_2, p_1)G(p_1, p_2)dA(p_2)$$

$A$表示场景中的所有表面。

在该方程中,首先我们已经完全消除了立体角变量,建立起关于空间中的点的表达式。同时,我们从方程中移除了隐式函数$t$,我们不再依赖$t$去寻找下一个顶点,而是通过一个可见行函数$V$对空间中的两点做可见性测试。由此,在数学意义上完全等效的同时,我们可以显式对该式展开。

路径积分

我们现在对面积积分方程进行递归展开

$$\begin{eqnarray}L(p_1, p_0) &=& L_e(p_1, p_0) + \int_A f(p_2, p_1, p_0)L(p_2, p_1)G(p_1, p_2)dA(p_2) \\&=& L_e(p_1, p_0) + \int_A f(p_2, p_1, p_0)\Bigg[L_e(p_2, p_1) + \int_A f(p_3, p_2, p_1)L(p_3, p_2)G(p_2, p_3)dA(p_3)\Bigg]G(p_1, p_2)dA(p_2) \\&=& L_e(p_1, p_0) + \int_A f(p_2, p_1, p_0)L_e(p_2, p_1)G(p_1, p_2)dA(p_2) \\&+& \int_A \int_A L(p_3, p_2)f(p_3, p_2, p_1)G(p_2, p_3)f(p_2, p_1, p_0)G(p_1, p_2)dA(p_3)dA(p_2) + ...\end{eqnarray}$$

不难察觉到,我们还可以继续展开$L(p_3, p_2)$,展开项为一个二重积分的自发光项和三重积分的反射项。对于多重积分项,我们不难发现,其具有一定的重复结构,我们将积分项定义为函数$P(\bar{p}_n)$,其表示具有$n+1$个顶点的光线路径的路径积分,其路径顶点顺序为$p_0,p_1,...,p_n$。其中,$p_0$表示摄像机。于是

$$P(\bar{p}_n) = \int_A \int_A ... \int_A L_e(p_n, p_{n-1}) \times \Bigg(\prod_{i=1}^{n-1}f(p_{i+1}, p_i, p_{i-1})G(p_i, p_{i+1})\Bigg)dA(p_2)...dA(p_n)$$

其中多重积分项为$n-1$重积分,BRDF函数的乘积表示光路的吞吐量,我们可以用$T(\bar{p}_n)$表示。于是

$$P(\bar{p}_n) = \int_A \int_A ... \int_A L_e(p_n, p_{n-1})T(\bar{p}_n)dA(p_2)...dA(p_n)$$

那么对于渲染方程来说,其可以直接表示为

$$L(p_1, p_0) = \sum_{n=1}^{\infty}P(\bar{p}_n)$$

路径追踪

我们已经知晓,最后进入眼睛(摄像机)的辐射亮度的路径追踪形式。然而,该形式是一个无穷级数,好在该级数是一个单调收敛级数,另外,一个直观的事实,对于路径较长的项来说,其吞吐由于BRDF的叠乘变得十分渺小,其对最后的数值贡献也显得不那么重要。因此我们可以在一定计算之后,将其截断,从而减少无穷项的计算量。为了保证路径积分的无偏属性,我们使用渐进式俄罗斯轮盘进行阶段,即:

$$L(p_1, p_0) = P(\bar{p}_1) + P(\bar{p}_2) + ... + P(\bar{p}_n) + \frac{1}{1-q_{n+1}}\bigg[P(\bar{p}_{n+1}) + \frac{1}{1-q_{n+2}}(P(\bar{p}_{n+2})+...)\bigg]$$

对于$n$之后的每一项,我们有$q_i$的概率停止计算,其中$q_i$是一个递增数列,从而保证越往后的项,截断概率越大。一般来说计算前三项,我们便可以得到一个足够准确的结果。

另外一点,路径积分和渲染方程是完全等价的,不过在路径积分形式中,我们只保留了$L_e$项,而所有的$L_i$都被展开。由此可知,如果任何一项$P(\bar{p}_n)$的最后一个顶点$\bar{p}_n$不是光源的话,那么该项则为零。所以在路径积分中,一个隐藏的条件是,所有的项最后一个顶点都应该是光源。

对于$P(\bar{p}_{n})$来说,我们已经知道这是一个$n-1$重积分,我们自然也可以使用蒙特卡洛随机积分方法进行积分。对于这一个$n-1$重积分来说,我们需要对$n-1$个点进行采样,并且计算每一个点的PDF,从而得到积分的近似。

一个最简单的采样方法是全局随机采样。很显然,对于整个模型来说,全局随机采样的PDF为

$$PDF_{p_i}=\frac{1}{\sum_i A_i}$$

值得注意的是,这个策略计算结果是无偏的。然而这个策略却容易导致两个严重的问题。第一个问题是,对于一个复杂场景来说,其可能导致非常大的方差。这是因为对于随机的两个点来说,他们有很大可能是互相不可见的,那么这会直接导致结果为0。另外,就算他们是可见的,但由于距离或者角度影响,可能导致他们对光路的实际贡献是非常小的,同样也使得结果为0。另外一个问题则可能会导致严重的错误。这是因为在BRDF函数中有一些特殊的Delta分布的函数,这一部分函数不能使用随机采样来拟合,因此该方法直接导致了这部分函数无法正常工作。

渐进式光路重建

既然随机选择的顶点序列会造成方差过大的问题,我们可以通过渐进式光路重建的方式来寻找一条更加稳定,方差更小的路径,即让每一条光路都尽可能的贡献最大辐射亮度。一个简单可行的方法是,我们从眼睛(摄像机)$p_0$顶点出发,然后和场景相交得到第一个顶点$p_1$以及其所对应的BRDF方程。由于BRDF方程的可逆性,我们可以对BRDF方程采样,得到出射光线,然后根据出射光线,继续和场景相交得到$p_2$顶点。持续这个步骤直到面积积分被俄罗斯轮盘截断。

因为在渐进式光路重建方法中,每一个点的光路的出射方向是通过BRDF采样得到的,然而BRDF采样是一个根据立体角$\omega_i$的采样方程,路径积分则是一个场景面积的积分形式,因此我们需要将立体角$\omega_i$的概率密度转换到面积的概率密度。根据立体角和面积的微分关系,我们有

$$p_A(p_{i+1}) = p_{\omega}(p_{i+1}-p_{i})\frac{|cos\theta_i|}{||p_{i+1}-p_{i}||^2}$$

蒙特卡洛近似

对于任意一个积分,我们可以通过蒙特卡洛方法近似。蒙特卡洛是一种随机算法,对于积分

$$F = \int_a^bf(x)dx$$

其有蒙特卡洛近似

$$F^N = \frac{1}{N}\sum_{i=1}^N\frac{f(X_i)}{pdf(X_i)}$$

蒙特卡洛近似是一个无偏近似。当随机变量$X_i$的概率密度函数和函数$f$完全相同的时候,蒙特卡洛近似就等于积分值。然而,选择一个概率密度相同的随机变量是不可能的,但是我们可以让其尽可能接近原函数,使得蒙特卡洛近似收敛的速度更快。

在路径积分中,我们有

$$P(\bar{p}_n) = \int_A \int_A ... \int_A L_e(p_n, p_{n-1}) \times \Bigg(\prod_{i=1}^{n-1}f(p_{i+1}, p_i, p_{i-1})G(p_i, p_{i+1})\Bigg)dA(p_2)...dA(p_n)$$

首先我们知道,根据渐进式光路重建方法,任意相邻的两点都是互相可见的,这意味着原函数中的$V$项永远等于1,可以直接忽略。

所以我们得到

$$\begin{eqnarray} P(\bar{p}_n) &\approx& \frac{L_e(p_n, p_{n-1}) \times \Bigg(\prod_{i=1}^{n-1}f(p_{i+1}, p_i, p_{i-1})G(p_i, p_{i+1})\Bigg)}{p_A(p_2)p_A(p_3)p_A(p_4)...p_A(p_n)} \\ &\approx& L_e(p_n, p_{n-1})\frac{G(p_{n-1}, p_n)\prod_{i=1}^{n-1}f(p_{i+1},p_i,p_{i-1})\prod_{i=1}^{n-2}\frac{|cos\theta_i||cos\theta_{i+1}|}{||p_i-p_{i+1}||^2}}{p_A(p_n)\prod_{i=2}^{n-1}\frac{p_\omega(p_i-p_{i-1})|cos\theta_i|}{||p_i-p_{i-1}||^2}} \\ &\approx& \frac{L_e(p_n, p_{n-1})G(p_{n-1}, p_n)f(p_{n},p_{n-1},p_{n-2})}{p_A(p_n)} \times \bigg(\prod_{i=1}^{n-2}\frac{f(p_{i+1},p_i,p_{i-1})|cos\theta_i|}{p_\omega(p_{i+1}-p_i)}\bigg) \end{eqnarray}$$

我们令

$$\beta_n=\prod_{i=1}^{n-2}\frac{f(p_{i+1},p_i,p_{i-1})|cos\theta_i|}{p_\omega(p_{i+1}-p_i)}$$

$\beta_n$就是在蒙特卡洛近似中的每一条光路的吞吐量。而每条光路的蒙特卡洛近似就等于

$$P(\bar{p}_n) = \frac{L_e(p_n, p_{n-1})G(p_{n-1}, p_n)f(p_{n},p_{n-1},p_{n-2})\beta_n}{p_A(p_n)}$$