次表面散射

继续写之前的基于物理着色系列，这一篇谈谈次表面散射（Subsurface Scattering）的理论以及渲染的算法。次表面散射是现实中一种非常常见的材质外观，所有有着半透明外观的物体，例如玉石，大理石，蜡烛，可乐，苹果，牛奶（见下图）包括人类的皮肤，等等。现实中这种半透明的材质其实非常普遍，所以如何为它们的外表建立正确的数学模型是照片级别真实感渲染的重要一部分。

下图则是电影《功夫熊猫3》中的一幕，画面中的玉石也是次表面散射材质的典型例子之一。

与此同时，次表面散射现象的模拟也比在前三篇文章里介绍过的一般的表面反射复杂很多，因为要正确的模拟这种现象，光线不止再物体的表面发生散射，而是会先折射到物体内部，然后再物体内部发生若干次散射，直到从物体表面的某一点射出，例如下图中的（b）。所以对于次表面散射性质的材质来说，光线出射的位置和入射的位置是不一样的，而且每一点的亮度取决于物体表面所有其他位置的亮度，物体的形状，厚度等。之前文章里介绍的表面散射都是基于BSDF模型（双向散射分布函数），而BSDF只能用于描述物体表面某一点的散射性质，所以它无法描述像次表面散射这种现象。

Participating Media

要模拟这样的现象，最简单，最精确，也是计算量最大，最慢的方法，就是直接在物体内部的空间求解带有Participating Media的渲染方程。为了文章的完整性，这里简单介绍一下Participating Media的数学定义，想详细了解可以参考pbrt第11章[1]。

Participlating Media的存在，意味着传统的渲染方程不再是传统的以场景中所有表面为定义域的积分，因为光线在有介质的空间中传播的时候，不用和表面接触，也会被空间中的介质所吸收，散射。通常我们用一下几个参数来描述空间中介质的性质。吸收系数 ，和散射系数 。前者描述了空间中每单位长度中光子被吸收的概率，后者则描述了每单位长度中光子被散射、既传播方向被改变的概率。散射系数包括了其他方向的光线入射到当然方向的概率，也包括当前方向出射到其他方向的概率。另一个参数是 ，就是衰减系数，它是吸收系数和散射系数的和，代表了单位距离内光子出射到其他方向或者被吸收掉的总概率。所以单位距离内辐射度的变化率可以写成：

第一眼看去貌似有些吓人，实际上说的就是，在介质中，单位长度上Radiance强度的变化会因为被介质吸收而减少，同时一部分能量会随机向外散射出去，最后来自各个其他方向的光线也有一定概率散射到当前方向。其中方程最右边积分中的 是phase function，可以理解成空间中的BSDF，定义了空间中一点两个方向间发生散射的概率。另一个值得一提的是，每个频率的光线都有对应的系数。也就是说在用RGB渲染时，每个颜色都要有各自的吸收和散射系数。

这篇文章封面的图片以及下面这张图都是用直接用模拟Participating Media的方法渲染出的玉石。玉石材质是由一个无色的，折射率为1.6的电介材质表面包裹着绿色散射系数较大的均匀介质。由于下图中介质的散射系数要远大于封面中的，所以光线散射的更加频繁，散射间传递的距离更短，导致外观比封面中的玉石更加厚密。

渲染Participating Media的方法有很多，例如最简单的Volumetric Path Tracing，Volumetric Photon Mapping等等。两个参考，[2]是最简单的Volumetric Path Tracing的tutorial，[3]是State of the art的Participating media estimator。因为话题很大，这里不再展开，以后有时间会再专门写写相关的内容。

BSSRDF的定义

模拟次表面散射的另一大类方法就是通过BSSRDF（Bidirectional Surface Scattering Reflectance Distribution Function）对材质建模。和BRDF的不同之处在于，BSSRDF可以指定不同的光线入射位置和出射的位置。Jensen在2001年的论文[4]可以说那是次表面材质建模最重要的一篇论文，推导了许多重要的物理公式，计算模型，渲染时的参数转换，以及测量了许多生活中常见材质的散射系数等等。大部分后来的论文都是在基于这篇文章中的理论的一些提升。不过这篇论文中有大量的数学以及物理公式，对没有一定数学基础和没有接触过材质建模，体渲染的人来说读起来很快就会摸不着头脑。好在大部分公式的意义都在于推导和证明BSSRDF的计算，所以这篇文章会挑出几个和实现次表面散射直接相关的重要公式简单介绍，如果想了解全部细节，还是请参考[4]原文。

BSSRDF的意义在于快速的近似。对于许多吸收系数特别低，散射系数特别高的材质来说（大部分浑浊的半透明物体例如牛奶，大理石，苹果，肉等等），一束光线会在物体内部散射成千上万次散射，才会彻底被吸收。直接用最暴力的Path Tracing去模拟则意味着每个Sample都需要散射上千次，收敛速度很慢。

对于之前的BRDF来说，一次反射光照的计算是在光线交点的法线半球上的球面积分。而对于BSSRDF来说，每一次反射在物体表面上每一个位置都要做一次半球面积分，是一个嵌套积分。

其中BSSRDF的定义 是：

其中包括了光纤入射和出射材质表面时候的折射Fresnel项 ，一个归一化的 项，和一个diffuse项 。除了diffuse项以外，其余的项都非常直观。下面就着重说说这个diffuse项。首先 只接受一个标量参数，这个参数的意义是光线入射位置和初设位置的曼哈顿距离。直观的理解就是，BSSRDF尝试将光线在物体表面内部中数千次的散射后所剩余的能量用一个基于入射点和出射点之间距离的函数去近似。

这个近似则是基于几个假设。第一个假设是，次表面散射的物体是一个曲率为零的平面。第二个假设是，这个平面的厚度，大小都是无限。第三个假设是，平面内部的介质参数是均匀的，最后一个假设是光线永远是从垂直的方向入射表面。正因为有这些假设，所以很容易把出射光的强度与出射点和入射点之间的距离用一个函数去近似。下图是用unbiased monte carlo方法模拟的一束光照在平面均匀介质

当然实际渲染中这种平面不存在，所有模型都不是无限厚度的，也都是有各种形状的，现实中所有物体内部的散射系数也都肯定不会是均匀的。这也是BSSRDF近似渲染的最主要的错误来源。下图是我在EDXRay里用两种完全不同的方法渲染的颜色，材质以及介质参数完全一样的两个佛像。左边的佛像是基于BSSRDF，右边的则是无偏的Participating Media，两者还是有细微的差别的。BSSRDF因为假设所有表面都的厚度都是无限的，在实际模型中如果有厚度比较薄的地方，误差会更为明显。这也是为什么BSSRDF渲染的结果相比Participating Media要暗淡一些。除了这些误差以外，BSSRDF也没办法渲染出焦散，无法正确的渲染表面以下的几何模型（注意封面shader ball左侧在球内部的写有sub-surface的黑色模型）。

那么究竟 的定义是什么呢。理论上来说它需要尽可能的只通过离光线入射点的距离就算出出射点的能量。首先在场景及其简化的情况下，例如一个点光源完全被包围在散射系数均匀的介质中，散射的能量和光源的距离的关系可以用下面的公式算出来。

但是这个积分并没有办法解析的求解，数值方法求得的解也会发散。再加上没有任何公式能直接根据距离准确的计算任何在有限空间内，表面形状任意时候的次表面散射的能量。所以寻求某种近似是必须的。其中一个这些年非常流行的近似算法就是用求解Diffusion方程的方法。

,

相比上面的精确解，Diffusion的近似简单了许多，只是一个距离自然指数除以距离。这个方程也是Jensen 2001[4]中Dipole近似的核心。实际上过去的所有次表面散射相关的论文都在不断的提出新的方法来更准确的理解这个方程的性质和提高它的精确度。

然而这篇文章并不想着重介绍Dipole[4]以及其延伸算法例如Multipole[5]，Directional Dipole[6]。所以读到这也不必太在意上面两个眼花缭乱的公式。

Normalized Diffusion

Pixar的Christensen（曾经差点和他成为同事）在去年发表的一篇论文中[7]，试验出了一个函数，它没有基于任何Diffusion或者光线传递的理论，就直接可以比几乎以往所有的算法都能更精确的直接逼近介质渲染方程的正确解。不仅如此，这个方程本身只含有两个指数函数，实现和计算都非常简单。而且它还可以的解析的求Cumulative Distribution Function(CDF)，渲染时要Importance Sampling也非常方便。实际上上面红色的佛像就是用这个方法渲染的。当然它同样也要满足以上提倒的那几个假设。因为这个函数从0到 的积分永远是1，所以这个方法的名字叫做Normalized Diffusion。它的定义如下：

，

就是这么简单。其中 是一个控制这个函数形状的参数。我们只需要根据不同的介质散射系数计算出适合的 ，就可以非常准确的逼近Participating Media的无偏解。下图截自[7]，意思就是说上面那个如此简单的公式可以做到比传统的数个复杂许多的算法更精确度逼近基于蒙特卡洛暴力积分的无偏解。A是材质表面的反射率（既RGB值）。注意，原来的Diffusion方法所计算的都只包括了材质中的多重散射（Multiple Scattering），直接散射（Single Scattering）则要单独用Ray Marching的方法求解。而Normalized Diffusion在函数里就默认包含了直接散射，实现和计算上都更加简单。

Normalized Diffusion所接受的介质参数和传统的有所不同。通常描述介质用的是吸收系数 ，和散射系数 。但是由于它们和最终材质外观的联系并不是一个线性关系，下面的表格可以看出 和 与材质的颜色的关系非常摸不着头脑。所以Artists并不能很方便的直接调整这些系数。所以其实大部分渲染器都会把这些系数映射成材质表面的Albedo 和平均散射距离 （diffuse mean free path）供Artists调整。不仅如此，也只有这也，才可以将纹理上的颜色应用到次表面散射材质上去。

从Albedo和平均散射距离到介质系数的映射在[4]中就有介绍。 其中 是reduced albedo， 。 和 则是散射系数受到phase function影响之后的reduced scattering coefficient。 , 。 则和表面的折射率有关， ， 。

平均散射距离 ，而 。

上面那堆公式就足够能解出材质的散射系数。这里可以看出Normalized Diffusion另一个优势是，方程的输入直接就是Diffuse Albedo和平均散射距离，不需要求解上方那些复杂的公式。

值本身的计算，也是通过函数拟合。给定一个材质的反射率 ，以及介质的平均散射距离 ，再假定光线是以光束的形式垂直射入物体的话， 。想了解其他拟合的计算方法，例如光线是以漫射的形式入射表面等，请参考[7]。

因为这里的 积分时永远等于1，所以渲染的时候只要直接乘以材质的Diffuse Albedo就可以了。当平均散射距离为0的时候，会收敛到和漫反射材质一样的结果。下图是五只用Normalized Diffusion渲染的兔子，平均散射距离从左到右由小变大。

终于，文章写到这里，BSSRDF的公式的定义终于介绍完整了！最后简单的谈谈具体如何用BSSRDF进行采样和渲染。

BSSRDF重要性采样

利用BSSRDF渲染的方法大致有两种。Jensen在[4]中提出了一种直接在物体的入射点周围根据 重要性采样出射点的方法。接着Jensen又在[8]中提出了一种基于点云的方法，先将许多点随机分布在次表面的模型上，然后在每个点预计算irradiance，最后利用这些点云信息快速查找积分。这个算法的优势是速度比前者快两三个数量级，缺点则有额外的内存开销，需要多一个预计算的Pass，实现起来较为麻烦，不适合progressive渲染，同时点云的分布影响渲染质量，容易在animation中出现flickering等等。

所以这里介绍的方法还是基于直接对BSSRDF进行重要性采样。Solid Angle在2013年提出一个多轴投射的采样方法[9]。首先根据 在入射点周围的圆盘上根据 的分布采样一个距离 ，接着将这个点垂直投影到模型的表面上。因为BSSRDF是在物体的表面积上积分，所以在遇到和入射点法线方向不同的表面时，样本的pdf需要乘以入射和出射点法线的点积。当入射与出射表面接近垂直的时候，pdf就会非常小。所以[9]中的做法是用一定的概率在binormal和tangent方向进行投影，然后用多重重要性采样计算样本的pdf，就可以在凹凸不平的表面达到更快的收敛速度。

至于对散射距离的采样，Normalized Diffusion的CDF可以直接用下面的公式计算: .

只不过采样需要的 无法解析的表达，这里可以直接用预计算查表的方法解决。计算的时候让 等于1，查表时乘以实际的 就能得到散射距离。有了以上采样的方法，就可以把次表面散射采样的路径作为一种特殊路径集成到Path tracing或者Bidirectional path tracing的积分中去。

最后放多几张次表面散射材质的美图。本文中所有渲染结果均使用本人自己开发的渲染器EDXRay。

参考文献