【实时渲染】皮肤渲染 理论及UE实现(SeparableSSS和Burley Normalize SSS)

UE中皮肤渲染采用的是Separable Subsurface Scattering和Burley Normalize SSS，本节我们先从理论上了解一下BSSRDF，接着分析以下Separable Subsurface Scattering和Burley Normalize SSS算法，并对UE中的实现进行分析。

渲染方程

BRDF

先从最熟悉的BRDF开始，假设表面上有一个点$p$，来自$w_i$方向的radiance是$L_i$，那么
$dL=\frac{d\Phi }{d{A}^{T}dw}=\frac{dE}{cos\theta dA}$
=>$dE(p,w_i)=L_i(p,w_i)cos{\theta }_{i}d{w}_i$

假设出射的radiacne $L_o$符合下面的式子：$dL(p,w_o)\propto dE(p,w_i)$。这样就得到了$f_r=\frac{dL(p,w_o)}{dE(p,w_i)}=\frac{dL(p,w_o)}{L_i(p,w_i)cos{\theta }_{i}d{w}_i}$,从这几个式子可以看出$f_r$位置上只与p点有关，p点接受入射光并反射到对应方向。

$L_o=={\int }_{S^2}{f}_r(p,w_i,w_o)L_i(p,w_i)|cos{\theta }_i|d{w}_i$
它的积分区域在整个球面上，是一个二维的积分

BSSRDF

BSSRDF与BRDF不同的是，在BSSRDF，光线在材质内部进行多次散射并最终离开材质内部，所以它的入射点和初始点不同。光线在内部的过程我们先不考虑，只考虑进出这两个关键点。
$S(p_i,p_o,w_i,w_o)=\frac{dL(p_o,w_o)}{dE(p_i,w_i)}$

先朴素地思考如何计算一个点$p_o$的次表面散射，首先在整个表面上采样一点$p_i$，对$p_i$上半球方向采样，根据出入点和出入方向计算S，并累加到最终结果上。之后再将整个表面的的样本累加。也就是对整个表面上的点采样，在对采样点的整个半球采样，这样维度就从原来整个半球的两维变成四维。
$L_o={\int }_A{\int }_{H^2(n)}S(p_i,p_o,w_i,w_o)L_i(p,w_i)|cos{\theta }_i|d{w}_{i}dA$

Separable BSSRDFs

光线在材质内部的传播是一个很复杂的过程，我们把这个过程分为三个部分，第一部分是光线从p进入此表面内部，进入到表面内部的能量是$1-F$，这里的F是菲涅尔项，在加一个缩放系数就得到：
$S_w(w_i)=\frac{1-F_r(cos{\theta }_i)}{c\pi }$

归一化系数逆推得到：
${\int }_{H^2}{S}_w(w)cos\theta dw=1$
=>$c={\int }_0^{2\pi }{\int }_0^{\frac{\pi }{2}}\frac{1-F_r(cos\theta )}{\pi }sin\theta cos\theta d\theta d\varphi$
=>$1-2{\int }_0^{\frac{\pi }{2}}1-F_r(cos\theta )sin\theta cos\theta d\theta d\varphi$

光线进入到材质内部之后会进行多次散射，这一部分的过程非常难建模，其实在后面会提到的几种实时渲染算法都是在对这部分过程进行建模，我们先忽略掉这一过程，把这一过程剩下能量的比率记为：
$S_p(p_o,p_i)$
接着进一步假设这个表面是个平面，并且这个比率只与两者的距离有关。
$S_p(p_o,p_i)\approx S_r(||p_o-p_i||)$

第三个过程是光线从材质内部出来，这一部分和第一个过程类似，即：
$1-F_r(cos{\theta }_o)$

这样整个方程就变成下面的样子
$L_o={\int }_A{\int }_{H^2(n)}S(p_i,p_o,w_i,w_o)L_i(p,w_i)|cos{\theta }_i|d{w}_{i}dA$
$={\int }_A{\int }_{H^2(n)}(1-F_r(cos{\theta }_o))(S_r(||p_o-p_i||))(\frac{1-F_r(cos\theta )}{\pi })L_i(p,w_i)|cos{\theta }_i|d{w}_{i}dA$
$=(1-F_r(cos{\theta }_o)){\int }_A(S_r(||p_o-p_i||)){\int }_{H^2(n)}(\frac{1-F_r(cos\theta )}{\pi })L_i(p,w_i)|cos{\theta }_i|d{w}_{i}dA$

可以发现这里面最难处理的就是$S_r$，它的参数很多$S_r(\eta ,g,\rho ,{\sigma }_t,r)$,这是一个高维函数，我们需要移除或固定某些参数来降低维度。首先是当给定${\sigma }_t$后，r和${\sigma }_t$存在一个对应关系，unitless optical radius即$r_{optical}={\sigma }_{t}r$。由于$S_r$是在一个2D的极坐标系上定义的，$S_r$可以简化为$S_r(\eta ,g,\rho ,{\sigma }_t,r)={\sigma }_t^2{S}_r(\eta ,g,\rho ,1,r_{optical})$

接着在固定$\eta ,g$，那么参数只剩下$\rho ,r_{optical}$，也就是给定$\rho ,r_{optical}$，我们要求出一个$S_r$。

注意这里的$\rho$描述的是单次散射之后的能量减少，而材质整体的$\rho$则需要考虑任意次数的散射，这里我们定义单次散射的albedo为$\rho$，整体的albedo为${\rho }_{eff}$,它们符合如下公式：
${\rho }_{eff}={\int }_0^{2\pi }{\int }_0^{\infty }r{S}_r(r)drd\varphi =2\pi {\int }_0^{\infty }r{S}_r(r)dr$

如何求$S_r$呢？Tabulated BSSRDF采用查表的方式，下面是如何生成这个表d 的步骤。

首先生成一系列半径r，生成方式是首先指定前两个样本的r，后续的样本成指数递增。接着生成一系列的$\rho$

$r_0=0$
$r_1=2.5e-3f$
$r_i=1.2r_{i-1}$

$\rho =\frac{1-8e^{-8i/(N-1)}}{(1-e^{-8})}$

然后对每对$\rho$和r，计算对应的$S_r$,这个表的大小是$N_{\rho }\times N_r$。这个$S_r$又如何计算呢？将$S_r$看作单次散射和多次散射的相加，即$S_r(r)=S_{ss}(r)+S_{ms}(r)$。接下来就是计算$S_{ss}(r)$和$S_{ms}(r)$，计算的方法来自这篇论文，下面是计算流程，只放MS的部分

输入${\sigma }_a{\sigma }_{s}g\eta r$

计算所需参数（步骤一）

(1)
${\sigma }_s^{{}^{\prime }}={\sigma }_s(1-g)$
${\sigma }_t^{{}^{\prime }}={\sigma }_a+{\sigma }_s^{{}^{\prime }}$
${\rho }^{{}^{\prime }}=\frac{{\sigma }_s^{{}^{\prime }}}{{\sigma }_t^{{}^{\prime }}}$
(2)
$D_g=\frac{2{\sigma }_a+{\sigma }_s^{{}^{\prime }}}{3({\sigma }_a+{\sigma }_s^{{}^{\prime }}{)}^2}$
（3）
${\sigma }_{tr}^{{}^{\prime }}=\sqrt{\frac{{\sigma }_a}{{\sigma }_g}}$
(4)
$Z_e=-2D_g\frac{1+3F_{r,2}(\eta )}{1-2F_{r,1}(\eta )}$
(5)
$C_{\varphi }=\frac{1}{4}(1-2F_{r,1}(\eta ))$
$C_d=\frac{1}{2}(1-3F_{r,2}(\eta ))$

采样nSample和点（步骤二）

for(i=1->nSample):
(1)
$Z_r=-\frac{log(1-\frac{i+0.5}{nSamples})}{{\sigma }_t^{{}^{\prime }}}$
$Z_v=-Z_r+2Z_e$
(2)
$d_r=\sqrt{r^2+Z_r^2}$
$d_v\sqrt{v^2+Z_v^2}$
(3)
${\Phi }_D=\frac{1}{4\pi D_g}(e^{-\frac{{\sigma }_{tr}}{d_r}}-e^{-\frac{{\sigma }_{tv}}{d_v}})$
$n{E}_D=\frac{1}{4\pi }(\frac{Z_r(1+d_r{\sigma }_{tr})}{d_r^3}{e}^{{\sigma }_{tr}{d}_r}-\frac{Z_v(1+d_v{\sigma }_{tv})}{d_v^3}{e}^{{\sigma }_{tv}{d}_v})$
（4）
$Ed+={\Phi }_D{C}_{\varphi }+n{E}_D{C}_d$

返回结果

$Ed/nSamples$

上面的方法是离线渲染中采用的，下面的两篇paper是用在实时渲染的

Separable Subsurface Scattering

$M(x,y)=\int \int E(x,y)R(r)dxdy$这里的M(x,y)是p点出射的radiance，E（x，y）是p点周围的irradiance。这个方程可以写成2D卷积的形式。
$M(x,y)=E(x,y)\ast R(r)$
2维的卷积对于实时来说性能消耗太大，我们使用高斯卷积来将它分成水平垂直两个分离的部分。
$R(r)={\Sigma }_{i=1}^k{w}_{i}G(v_i,r)$

单个高斯很难模拟diffuse profile，我们采用多个高斯来拟合。但是多个高斯合起来的函数是不可分离的，所以必须用多个Pass对不同的高斯函数进行处理。多个Pass就带来了性能上的消耗。

Separable Subsurface Scattering是动视在2015年发布的一个Paper，这个Paper的主页有点迷惑，这个主页上放了2012年一个Github实现代码，我看了一下代码和这篇论文并没有很大的关系，但是UE种是按Github实现代码做的。

    return  // 0.233f * gaussian(0.0064f, r) + /* We consider this one to be directly bounced light, accounted by the strength parameter (see @STRENGTH) */
               0.100f * gaussian(0.0484f, r) +
               0.118f * gaussian( 0.187f, r) +
               0.113f * gaussian( 0.567f, r) +
               0.358f * gaussian(  1.99f, r) +
               0.078f * gaussian(  7.41f, r);

Separable Subsurface Scattering将多个高斯函数的累积作为一个新的函数，在水平和竖直方向进行累积，整个算法就是这些内容，之前是这6个高斯函数分到6*2个Pass中或者是在一个Pass上进行2D的Sample，现在是直接分开，就是水平和竖直两个部分。这里我比较困惑的一点就是这个新的函数为什么是可分离的。

2015年的论文就比较清晰，直接将$R(r)$做SVD分解，只取1-rank，那么不论是不是可分离的都可以用两个水平或竖直的Pass处理。

Burley Normalize SSS

选择高斯函数去拟合diffusion profile有以下几个好处：
1.高斯的卷积可以分离，降低了复杂度
2.不断地用一个”小“的高斯进行卷积相当于用一个大的高斯进行卷积。
使用高斯的Separable Subsurface Scattering也有几个问题：
1.不同高斯核的混合不能分离
2.SSSS的参数过多，并且其物理意义不明，难以调整。
3.至少需要4个Pass，消耗较大。

Burley Normalize SSS使用下面的diffusion profile：
$R(r)=As\frac{e^{-sr}+e^{-sr/3}}{8\pi r}$

完整的公式为：
$R(r)=As\frac{e^{-sr/l}+e^{-sr/3l}}{8\pi r}$
l是mean free pass，作者认为要去用s拟合曲线将l取为1就足够了

Burley Normalize SSS的累积分布函数CDF很容易计算：
$cdf(r)=\frac{{\int }_0^{r}R(t)2\pi tdt}{{\int }_0^{\infty }R(t)2\pi tdt}$
$=\frac{A/4(4-e^{-r/d}-3e^{-r/(3d)})}{A}$
$=(1-\frac{1}{4}{e}^{-r/d}-\frac{3}{4}{e}^{-r/(3d)})$

r采样

cdf属于$[0,1]$，给定一个0-1的随机数，可以得到其对应的半径为
$cdf(r,d)=random$
$r=cd{f}^{-1}(random,d)$

UE中给出了3中方法：
第一种是近似方法，$cd{f}^{-1}$近似为以下函数：
$cd{f}^{-1}(d,random)=d((2-2.6)random-2)log(1-random)$

第二种采用数值计算方法（Halley）进行求解：
$x_{n+1}=x_n-\frac{2f(x)f^{{}^{\prime }}(x)}{2f^{{}^{\prime }}(x)f^{{}^{\prime }}(x)-f(x)f^{{}^{\prime \prime }}(x)}$

第三种计算解析解：
令$F(X)=\frac{1}{4}{e}^{-X}+\frac{3}{4}{e}^{-X/3}$
$cdf(X)=1-F(X/D)$
也就是计算$cdf(r,D)=random$
=>$F(X/D)=1-random$
=>$X=D\ast F^{-1}(1-random)$
$F^{-1}$可直接计算
令$U=1-random$
$G=1+4U(2U+\sqrt{1+4U^2})$
$F^{-1}(1-random)=F^{-1}(U)=3log(\frac{1+G^{-1/3}+G^{1/3}}{4U})$

原始方法的样本点大多集中在中心，UE中做了改进，首先计算出center sample在世界空间中的半径R，然后计算出R对应的CDF®=T，之后把中心的Smple和其他区域分开。
不在中心的样本在$[T,1]$中取，并将计算结果呈上缩放系数1-T。

s拟合

zai知道r后接下来就是计算$R(r)$了，在R( r)中还有一个参数没有计算，原始论文中给出了三种拟合s的公式：

第一种假设的情况是光线垂直进入表面，它对MC的平均相对误差为5.5%
$s=1.85-A+7|A-0.8{|}^3$

第二种情况假设光线是任意方向的，它对MC的平均相对误差为3.9%
$s=1.9-A+3.5(A-0.8{)}^2$

第二种方法用diffuse mean free path (dmfp)$l_d$代替mean free path (dmfp)$l$拟合s，它对MC的平均相对误差为7.7%
$s=3.5+100(A-0.33{)}^4$

Bilateral Filtering

基本的工作已经完成了，不过这个模型还有一个很大的缺陷就是之前假设表面是平面，并且在平面上进行采样，这显然错误，现在需要把深度信息考虑进去。

$I\approx \frac{2\pi }{n}{\Sigma }_{i=1}^n\frac{R(r_i,d_i)}{p(r_i,d_i)}L=\frac{2\pi }{n}{\Sigma }_{i=1}^n\frac{R(\sqrt{r_i^2+d_i^2})}{p(r_i,d_i)}L$

接下来就碰到困难了，KaTeX parse error: Expected 'EOF', got '}' at position 11: p(r_i,d_i)}̲怎么算，之前的pdf是按平面计算的，它们的对应关系很难求。Unity在报告中采用强制能量守恒，添加一个归一化常数。

$I=\frac{{\Sigma }_{i=1}^n\frac{R(\sqrt{r_i^2+d_i^2})}{p(r_i,d_i)}L}{{\Sigma }_{i=1}^n\frac{R(\sqrt{r_i^2+d_i^2})}{p(r_i,d_i)}}$

UE中Bilateral Filtering相关代码如下：

float3 WeightingFactor = 0.0f;
LOOP for (int i = 0; i < NumOfSamples; ++i)
{
	``````
	float DeltaDepth = (SampledRadianceAndDepth.w - DepthAtDiscCenter)*BURLEY_CM_2_MM;
	float RadiusSampledInMM = sqrt(BurleySampleInfo.RadiusInMM * BurleySampleInfo.RadiusInMM + DeltaDepth * DeltaDepth);
	``````
	float3 DiffusionProfile = GetDiffuseReflectProfileWithDiffuseMeanFreePath(BurleyParameter.DiffuseMeanFreePath.xyz, S3D.xyz, RadiusSampledInMM);
	float3 SampleWeight = (DiffusionProfile / BurleySampleInfo.Pdf) * Mask * NormalWeight;
	RadianceAccumulated.xyz += SampleWeight * (SampledRadianceAndDepth.xyz*TintColor);


	WeightingFactor += SampleWeight;
	``````
}
// 0.99995f is a compensitation to make it energe conservation.
const float EnergyNormalization = 1.0f / 0.99995f;

RadianceAccumulated.xyz *= (WeightingFactor == 0)? 0:(1.0f /WeightingFactor*EnergyNormalization);

Reference：
PBRT11-4BSSRDF
PBRT15-5Subsurface Scattering Using the Diffusion Equation
Separable Subsurface Scattering
Approximate Reflectance Profiles for Efficient Subsurface Scattering
Efficient screen space subsurface scattering Siggraph 2018