PBR论文简读:A Practical Model for Subsurface Light Transport

PBR论文简读:A Practical Model for Subsurface Light Transport

导语

离线环境下进行真实感渲染的时候,会有各种散射材质,比如人类皮肤、硅胶等。从外观上看,这类材质的通性在于都有一种“透光”的感觉,简单的使用brdf创造出来的材质往往显得很生硬,因此Jensen'01的这篇论文给出了一个实际的次表面散射(subsurface scattering)bxdf来实现这种材质。

理论

BRDF理论假设材质光的进入和离开的起、终点一致,本质上是对BSSRDF的一种近似\(x_o = x_i\),并且给出的积分公式也是简单地在半球面上进行积分,BSSRDF则考虑的是出、入射点不同时,对于区域\(A\)内所有入射光线的radiance的半球面积分:

\[L_o = \int_A\int_{2\pi}S\cdot L_i\cdot (\vec n\cdot \vec {\omega_i})d\omega_i dA(x_i)\]

我们需要的就是文章给出的BSSRDF:\(S\)项,这里Jensen'01的论文将次表面散射的BSSRDF拆解成了Diffusion ApproximationSingle Scattering两个term求和。前者基于的是:通过观察,发现高度散射的介质的光线分布有着各项同性的特征(即使入射光和相位函数是高度各向异性的),后者论文中并没有给出一个拆分的理由,而是从前人的论文中得到的BRDF扩展成为了BBSSRDF,并且后者当且仅当出、入射光的折射发生相交时才会出现,个人理解是对于diffuse的一个specular补偿。

推导

在介质中,光线的传播被我们使用辐射传播方程所描述:

\[(\vec{\omega} \cdot \vec{\nabla}) L(x, \vec{\omega})=-\sigma_{t} L(x, \vec{\omega})+\sigma_{s} \int_{4 \pi} p\left(\vec{\omega}, \vec{\omega}^{\prime}\right) L\left(x, \vec{\omega}^{\prime}\right) d \omega^{\prime}+Q(x, \vec{\omega})\]

我们分别使用三个参数:吸收系数\(\sigma_\alpha\), 散射系数\(\sigma_s\)(这两者的和定义为消散系数\(\sigma_t\))和相位函数\(p(\vec w, \vec w')\)来描述介质的属性,其中相位函数可以看作是描述在介质中一点的光线散射的角度分布[1]的函数,同时文章假设这个函数已经normalize了(对于整个球面积分为1)。紧接着定义了平均余弦函数\(g\)用来描述相位函数的性质\(g=\int_{4 \pi}\left(\vec{\omega} \cdot \vec{\omega}^{\prime}\right) p\left(\vec{\omega} \cdot \vec{\omega}^{\prime}\right) d \omega^{\prime}\):大于0说明散射大部分是前向的,小于0则是反向的,等于零说明在球面上均匀分布(各向同性)。对于一个理想的无穷小光束,radiance根据传播距离\(s\)的衰减定义为指数级的:\(L_{r i}\left(x_{i}+s \vec{\omega}_{i}, \vec{\omega}_{i}\right)=e^{-\sigma_{t} s} L_{i}\left(x_{i}, \vec{\omega}_{i}\right)\)。这样,一阶的散射就可以被作为次级的体积光源使用,可以得出一个点在\(x\)处接收到周围某个点发射的irradiance:\(Q(x, \vec{\omega})=\sigma_{s} \int_{4 \pi} p\left(\vec{\omega}^{\prime}, \vec{\omega}\right) L_{r i}\left(x, \vec{\omega}^{\prime}\right) d \omega^{\prime}\)。接下来就可以对\(x\)点的辐射传播方程的散射参数\(\omega'\)进行积分,获得周围所有光对\(x\)点的影响,将两个irradiance,\(\phi(x)\)\(\vec E(x)\)联系了起来:

\[\vec{\nabla} \cdot \vec{E}(x)=-\sigma_{a} \phi(x)+Q_{0}(x)\]

这里引入的零阶源\(Q_{0}(x)=\int_{A_{\pi}} Q(x, \vec{\omega}) d \omega\)就是前文提到的次级光源,同时将更高阶的源1定义为\(\vec{Q}_{1}(x)=\int_{4 \pi} Q(x, \vec{\omega}) \vec{\omega} d \omega\)

Diffusion Approximation

上一节说到,diffusion被看作是各向同性的,论文中就直接利用这个特性把radiance展开成了radiant fluence和vector irrandiance和(这两个本质上都是irradiance),除以了球体单位角的总和\(4\pi\)即是\(Le\)(原文里面说这两个的常数系数是由fluence和vector irradiance的定义决定的):

\[L(x, \vec{\omega})=\frac{1}{4 \pi} \phi(x)+\frac{3}{4 \pi} \vec{\omega} \cdot \vec{E}(x) \tag{1}\]

将这个式子代入辐射传播方程,类似的,对\(\omega\)进行积分,得到的结果式:\(\vec{\nabla} \phi(x)=-3 \sigma_{t}^{\prime} \vec{E}(x)+\vec{Q}_{1}(x)\),其中\(\sigma_{t}^{\prime}=\sigma_{s}^{\prime}+\sigma_{a} \text { where } \sigma_{s}^{\prime}=\sigma_{s}(1-g)\),然后将上式带入上上个等式,得到经典的diffusion方程:

\[D \nabla^{2} \phi(x)=\sigma_{a} \phi(x)-Q_{0}(x)+3 D \vec{\nabla} \cdot \vec{Q}_{1}(x)\]

在一个无限大的介质里的各项同性点光源,上式\(x\)点的irradiance有一个简单的解如下:

\[\phi(x)=\frac{\Phi}{4 \pi D} \frac{e^{-\sigma_{t r} r(x)}}{r(x)}\]

当考虑传输介质是处于有限区域的时候,就需要考虑边界区域(出射)的情况,这里论文考虑边界只有出射的情况,也就是入射的flux为0:\(\int_{2 \pi_{-}} L\left(x_{s}, \vec{\omega}\right)\left(\vec{\omega} \cdot \vec{n}\left(x_{s}\right)\right) d \omega=0\)(法线向外,但是在内部积分),使用公式(1)展开有:

\[\phi\left(x_{s}\right)-2 D(\vec{n} \cdot \vec{\nabla}) \phi\left(x_{s}\right)=0 \tag{3}\]

但是上面的公式只考虑了surface两边IOR相等的情况,对于两边不等的情况,在表面处会发生反射,根据菲涅尔公式可以计算出平均的菲涅尔漫反射系数\(F_{dr}\),具体方法就是对菲涅尔系数与视线、法线的夹角余弦做积分。可以用菲涅尔公式做一个解析的解,但是论文里面选择了一个基于多项式的近似:

\[F_{d r}=-\frac{1.440}{\eta^{2}}+\frac{0.710}{\eta}+0.668+0.0636 \eta\]

根据前文的净通量为0,同样的有:

\[\phi\left(x_{s}\right)-2 D(\vec{n} \cdot \vec{\nabla}) \phi\left(x_{s}\right)=F_{d r}\left[\phi\left(x_{s}\right)+2 D(\vec{n} \cdot \vec{\nabla}) \phi\left(x_{s}\right)\right]\]

化简:

\[\phi\left(x_{s}\right)-2 A D(\vec{n} \cdot \vec{\nabla}) \phi\left(x_{s}\right)=0\tag{4}\]

实际上观察这个公式和(3)的区别也仅仅在于多了一个平均漫射菲涅尔系数相关的\(A\)。这样我们就可以根据边界的值,通过出射的radiance(\(\vec n \cdot \vec E(x_s)\)等于表面的通量梯度)除以入射的flux计算出BSSRDF diffusion approximation了:

\[R_{d}(r)=-D \frac{(\vec{n} \cdot \vec{\nabla} \phi)\left(x_{s}\right)}{d \Phi_{i}\left(x_{i}\right)}\]

但是有限介质的辐射方程没有精确的解析解,为了简化模型,论文引用了“双极子”的方法模拟体积光源分布,对于次级体积光源,引入了两个点光源:

其中正光源在surface下,负光源在surface上,用来模拟光的衰减效应,两者总体的fluence满足:

\[\phi(x)=\frac{\Phi}{4 \pi D}\left(\frac{e^{-\sigma_{t r} d_{r}}}{d_{r}}-\frac{e^{-\sigma_{t r} d_{v}}}{d_{v}}\right)\]

这样,在有限介质下的边界条件也得到了,紧接着就可以带入之前的公式求出\(R_d(r)\)

\[\begin{aligned} R_{d}(r) &=-D \frac{\left(\vec{n} \cdot \vec{\nabla} \phi\left(x_{s}\right)\right)}{d \Phi_{i}} \\ &=\frac{\alpha^{\prime}}{4 \pi}\left[\left(\sigma_{t r} d_{r}+1\right) \frac{e^{-\sigma_{t r} d_{r}}}{\sigma_{t}^{\prime} d_{r}^{3}}+z_{v}\left(\sigma_{t r} d_{v}+1\right) \frac{e^{-\sigma_{t r} d_{v}}}{\sigma_{t}^{\prime} d_{v}^{3}}\right] \end{aligned}\]

Update 2022-5-12 在实现本论文算法的时候发现这个\(R_d(r)\)的式子是有问题的,第二行中括号前后两个式子的量纲是不一样的(后面的式子多了一个\(z_v\),会导致结果与实际产生10的指数级别的误差,正确的公式应该参照d'Eon的A better dipole中给出的公式(看了这么多实现竟然没有人指出这个问题,有点匪夷所思了))。

可以看出来,\(R_d\)项实际上是只与散射材料本身和双极子与表面距离相关的函数(非严格意义上的反比例,这样就给出了一个很好的性质,让我们在做BSSRDF积分的时候不需要完整地考虑整个surface,而是主要对临近的表面积分即可),进而可以求出BSSRDF中的diffusion term:

\[S_{d}\left(x_{i}, \vec{\omega}_{i} ; x_{o}, \vec{\omega}_{o}\right)=\frac{1}{\pi} F_{t}\left(\eta, \vec{\omega}_{i}\right) R_{d}\left(\left\|x_{i}-x_{o}\right\|\right) F_{t}\left(\eta, \vec{\omega}_{o}\right)\]

这一项主要是代表了光线在介质中多次散射产生的衰减。

Single Scattering

论文中还专门将一次散射拆成了一个单独的项\(S^{(1)}\)提了出来,这里个人暂时无法理解这么做的理论依据,只能理解为单次散射无法看作各向同性的新光源,而是作为一个specular项存在:

\[\begin{aligned} L_{o}^{(1)}\left(x_{o}, \vec{\omega}_{o}\right) &=\sigma_{s}\left(x_{o}\right) \int_{2 \pi} F p\left(\vec{\omega}_{i}^{\prime} \cdot \vec{\omega}_{o}^{\prime}\right) \int_{0}^{\infty} e^{-\sigma_{t c^{s}}} L_{i}\left(x_{i}, \vec{\omega}_{i}\right) d s d \vec{\omega}_{i} \\ &=\int_{A} \int_{2 \pi} S^{(1)}\left(x_{i}, \vec{\omega}_{i} ; x_{o}, \vec{\omega}_{o}\right) L_{i}\left(x_{i}, \vec{\omega}_{i}\right)\left(\vec{n} \cdot \vec{\omega}_{i}\right) d \omega_{i} d A\left(x_{i}\right) \end{aligned}\]

论文认为单次散射的出射radiance \(L_o\)等于对折射路径\(s\)上的入射radiance进行积分(第一行的内部积分),再对入射的半球面立体角进行积分,总体上也就等价于固定了\(x_o\),再对入射的立体角和单位入射面积进行二重积分,便可以得到\(S^{(1)}\)的隐式表达。

论文里还指出了diffusion term和single scattering有着不同的平均出射距离,并且两者的距离尺度有着很大的差别,对于有效消散系数\(\sigma_{tr}\)远远小于消散系数\(\sigma_t\),single scattering随着\(x_o\)增加而衰减的情况会远明显于diffusion。

BRDF近似

完成上面两项的推导之后,论文提出了使用BRDF模拟BSSRDF的做法:假设surface上受到的光照是均匀的,这样就可以通过对于整个surface积分来实现计算(\(r\)从0积分到正无穷),对于single scattering也是同样的,就可以得到近似的BRDF:

\[ f_{r}\left(x, \vec{\omega}_{i}, \vec{\omega}_{o}\right)=f_{r}^{(1)}\left(x, \vec{\omega}_{i}, \vec{\omega}_{o}\right)+F \frac{R_{d}}{\pi}\]

对于不透明,消散系数很大的材质(使得平均自由路径很短)来说,这个方法可以较好地模拟出BSSRDF(实际积分的权重多数集中于入射、反射点)。

BSSRDF实现

相较于BRDF,虽然我们已经简化了BSSRDF原本要求的入射区域A的积分而变成了在入射的位置进行多次采样,论文提出使用类似于分布式光线追踪的思路,在一个shadow ray的两端进行多次采样,并且针对single scattering和diffusion两项进行不同的采样。其中diffusion采样直接使用蒙特卡洛方法根据离出射点的距离,使用\(\sigma_{t r} e^{-\sigma_{t r} d}\)的权重进行采样。对于single scattering,由于入射和出射的折射光必须相交,论文对其进行了重组参数处理。

对于任意的几何表面,为了避免高度曲折情况导致的artifact,论文认为应该始终选择将双极子放置在最小距离\(1/\sigma'_t\)处,这样可以处理在边界处出现的奇点问题。同时,如果物体表面有纹理材质,那么论文将纹理近似为仅仅覆盖在物体的表面,因此内部的体积光散射将不会收到纹理的variance造成的影响。

作者

Carbene Hu

发布于

2022-04-09

更新于

2022-05-12

许可协议

评论