前言

前文提到了如何实现适合于移动端的预积分次表面散射效果，但这种效果并不是特别好，对于PC端还有其他效果更好的实时渲染方案，且两大的引擎Unity、UE都用到了此方案——即Screen Space Subsurface Scattering(SSSSS, 5S)。但光是5S需要消耗不少性能，基于此基础上，两大引擎又做了优化，使用高斯混合模型，以及高斯核可分离性，让 SSS 可以以比较低的代价实时的运行起来——即Separable Subsurface Scattering)(4S)

UE的方案

UE对皮肤、玉石等具有次表面散射效应的物体提供了三种渲染方案：

Subsurface
Preintegrated Skin
Subsurface Profile

三者性能消耗逐渐递增。其中，Preintegrated Skin在前文已经介绍过。而4s相对传统3s主要区别在于：

适用于延迟渲染，否则许多参数获取十分麻烦
基于屏幕空间
要求输入分离的Diffuse和Specular
算法核心的卷积（滤波）运算，工作在后处理之后，tonemapping之前

次表面散射算法

关于次表面散射的成因及皮肤组织建模在上篇已经提过，这里不再提及

高光

高光部分依旧采用Cook-Torrance的BRDF，对应公式如下：

漫反射

BSSRDF

漫反射部分使用到的是BSSRDF，相较于BRDF，BSSRDF每一次反射在物体表面上每一个位置都要做一次半球面积分，是一个嵌套积分：

它的含义是：当一束光以任意角度入射到某个确定的微表面p上时，有多少Radiance能够从任意一个给定的微表面q上以某个给定的出射角度反射出去

对BSSRDF渲染方程的解析在上篇也有提及，这里不再解释

对于BSSRDF渲染方程的S函数部分可以笼统的归纳为以下形式：

表示以 $ω_i$ 角度向 $x_i$ 位置入射的radiance，在 $x_o$ 处以 $ω_o$ 角度出射的irradiance为多少。这部分看似比较简单，但是它是次表面散射中最为重要的一部分，也是最复杂的一部分

对于S函数，有如下经验公式：

其中：

$F_t$ ：Fresnel透射效应，用于模拟入射和出射过程的消耗，一般只和材质属性和角度有关
R_d(|p_i - p_o|)：散射函数，主要与入射点和出射点的距离有关。它的真身如下：
- D：漫反射常量

可以看到BSSRDF渲染方程是十分复杂的，想用此做到实时是不可能的。虽然完全物理的做不到，但可以用近似法来逼近以降低消耗

其中：

$F_r(cos_{θ_o})$ ：Fresnel反射项
$S_p(p_o, p_i)$ ：点p处的散射函数。可以进一步表示为 $S_r(|p_0 - p_i|)$

而 $S_r$ 与介质的许多属性相关，可用公式简化并表达：

简化后的Sr只跟ρ和r有关，每种材料的ρ和r可组成一个BSSRDF表：
$S_ω(ω_i)$ ：含有缩放因子的Fresnel项。公式：

其中c是一个嵌套的半球面积分：

SeparableSSS

Diffusion Profile

实时渲染中，次表面散射的本质是利用附近pixel进行加权平均，具有以下特点：

先对每个像素进行常规的漫反射计算
再根据某种特殊的函数Rd(r)和上述漫反射结果，加权计算周围若干个像素对当前像素的次表面散射贡献

可以看到，次表面散射和Blur很像，都利用卷积核进行加权。其实，这里的Rd(r)函数和卷积核是同一概念，该函数可以绘制得到所谓的扩散剖面（Diffusion Profile）。如下图所示：

那Rd(r)函数究竟有何物理含义呢？对于一块面积无限，介质均匀，厚度无限的理想表面，当一速激光照射上去时会引发光线向周围的扩散，形成以照射点为中心的稳定光晕，建模出来可以得到以下结果：

所谓的Diffusion Profile，即上图中投影到x-z或y-z平面上的一块投影。它反映了散射光强随距离的分布趋势，实验表明这种趋势只受到材质属性、扩散距离的影响，与光线的入射/出射位置或者光线的入射/出射方向无关，因而是各向同性（isotropy）的分布函数

Separable Convolution

虽然BSSRDF用到了近似逼近来优化性能，但在后续还需进行2D卷积，计算量还是很大，有没有什么办法对2D卷积进行优化呢？还真有，也就是本文着重介绍的SeparableSSS

SeparableSSS是Jimenez和Gutierrez在2015年的论文中提出的，主要贡献是将2D卷积拆分成2个相关的1D卷积，该方法称为卷积分离（Separable Convolution）。主要形式如下图所示：

证明过程如下图所示：

其中：

$E$ ：radiance
$R_d$ ：散射函数
$x、y$ ：出射点位置
$x'、y'$ ：积分变量
$dx'、dy'$ ： $R^2$ 的一块微元面
$A$ ：对 $R_d$ 函数的近似
第三行：将 $A$ 拆分为多个低维函数，并求和
最后一行：也是关键所在，因为卷积核a(x)、a(y)没有关联，所以可以分开积分

其实，Separable Convolution就是在模糊中经常用到的高斯模糊，水平加权一次，竖直再加权一次。只是并不是所有卷积核都可以Separable

$R_d$的近似函数A

为了得到近似函数A，需要先确定 $R_d$ 函数。数学上定义的散射函数是十分复杂的，它是由Jensen在2001年提出的散射模型，存在大量的公式推导，这里就省去了，直接下结论。如下图所示：

这里提到了偶极子（Dipole）的概念，所谓Dipole即一对互为正负的点光源。其中：

将正光源放于表皮层的下方 $Z_r$ 深处。理想状态下，正光源向外均匀辐射光强，同时受到介质的影响，衰减散射
将负光源放于表皮层与正光源对立的一侧，高度为 $Z_v$ 。负光源用于吸收从表面散射的光强，以调节强度分布

说了这么多，还是不太清楚Dipole到底是个啥东西，本质上讲Dipole表示的是皮肤表面x处的光通量，它可以用以下公式表达：

其中：

$D$ ：散射常量
$σtr$ ：透射吸收因子
$d_r$ ：点x到正光源的直线距离
$d_v$ ：点x到负光源的直线距离

进一步可以得到 $R_d$ 函数的公式：

对于牛奶、大理石这类材质，一个Dipole足以描述散射函数，消耗貌似不大。但对于皮肤这样多层结构的材质，需要使用3个Dipole才能达到理想效果，消耗是无法接受的，综合来说不适合用于实时

优化

但Diffusion Profile的分布和高斯函数类似，如下图所示：

这意味着，可以通过结合不同参数的高斯函数来近似Diffusion Profile。高斯函数还有个优点：可分离性和径向对称性，这使得高斯函数满足卷积分离的需求

使用6个高斯函数即可拟合皮肤等带有3层Dipole的Diffusion Profile

其中：

$w_i$ ：控制参与求和的高斯函数形状
$v_i$ ：方差
$r$ ：散射距离

在论文中对皮肤的拟合结果最终给出如下的参数，可见R、G、B通道拟合出的曲线有所不同，而R通道曲线的扩散范围最远，这也是皮肤显示出红色的原因：

需要注意的是，在实际计算时，只需计算R通道的散射效果。因为即便如此，对于蓝绿色的散射效果也是可以接受的

卷积核大小

在实现Blur效果时，所用到的卷积核是固定大小的。但在实现SSSS时，情况不同了，因为这里卷积核代表基于物理的Diffusion Profile，若是固定的，得到的模糊效果整体都是一样的，且因为透视投影的原因，不可能每个pixel所占面积都是相同的，所以需要根据不同的区域属性来调整卷积核大小

决定卷积核大小的因素如下：

缩放系数

物理意义是卷积核大小应随pixel深度的增加而减小，且随视场角的变小而变大
- $f_y$ ：相机的视场角
- $p_d$ ：当前pixel的深度（非[0,1]）
SSS宽度：限定在什么样的尺度上能够看到次表面效果
SSS强度：通常以贴图的形式存在，用以表示不同区域的次表面强度

几何过滤

在屏幕空间几个相邻的点，次表面散射材质是连续的；但在世界空间上不一定连续。为了避免混淆几何上迥异的2块次表面材质，造成过渡不自然。解决方案是比较当前采样点和周围采样点的像素颜色和深度值

其中：

$i_c$ ：当前pixel的颜色
$o_c$ ：周边pixel的颜色
$i_d$ ：当前pixel的深度
$o_d$ ：周边pixel的深度

工业界的流程

大名鼎鼎的UE使用到的流程如下：

大致含义是先对高光和漫反射光分开计算，对漫反射光再次分离，分成水平竖直两次计算，最后叠加

实现

这里先不分离漫反射和高光

高斯函数

/// <summary>
/// 高斯函数. 增加了FalloffColor颜色，对应不同颜色通道的值
/// </summary>
/// <param name="variance"> 方差 </param>
/// <param name="radius"> 次表面散射的最大影响距离。单位mm </param>
/// <param name="falloffColor"> 光线随距离增加而衰减的程度, 数值越小表示对应方向上光线衰减得越快，数值越大表示衰减得越慢 </param>
/// <returns></returns>
static Vector3 SeparableSSS_Gaussian(float variance, float radius, Vector3 falloffColor)
{
    Vector3 result = Vector3.zero;

    for (int i = 0; i < 3; ++i)
    {
        float rr = radius / (0.001f + falloffColor[i]);
        result[i] = Mathf.Exp(-(rr * rr) / (2f * variance)) / (2f * Mathf.PI * variance);
    }

    return result;
}

近似Diffusion Profile

/// <summary>
/// 6个高斯函数拟合3个dipole曲线
/// </summary>
/// <param name="radius"> 次表面散射的最大影响距离。单位mm </param>
/// <param name="falloffColor"> 光线随距离增加而衰减的程度, 数值越小表示对应方向上光线衰减得越快，数值越大表示衰减得越慢 </param>
/// <returns></returns>
static Vector3 SeparableSSS_Profile(float radius, Vector3 falloffColor)
{
    // 去掉0.233f * SeparableSSS_Gaussian(0.0064f, radius, falloffColor)。理由是这个是直接反射光，且考虑了strength参数
    return 0.1f * SeparableSSS_Gaussian(0.0484f, radius, falloffColor) +
           0.118f * SeparableSSS_Gaussian(0.187f, radius, falloffColor) +
           0.113f * SeparableSSS_Gaussian(0.567f, radius, falloffColor) +
           0.358f * SeparableSSS_Gaussian(1.99f, radius, falloffColor) +
           0.078f * SeparableSSS_Gaussian(7.41f, radius, falloffColor);
}

计算卷积核数据

确定卷积范围：kernelList[i].w

// 卷积核先给定一个默认的半径范围，不能太大也不能太小，根据nTotalSamples调整Range(单位是毫米mm)
// 得到一个前半为负，后半为正的步长数组，并通过一个步长函数来控制步长值
float range = nTotalSamples > 20.0f ? 3.0f : 2.0f;
const float exponent = 2.0f;    // 指定步长函数的形状.高斯分布的简化版，距离原点越远的样本步长越大
float step = 2.0f * range / (nTotalSamples - 1);    // 每次卷积的偏移值
for (int i = 0; i < nTotalSamples; ++i)
{
   float o = -range + (float)i * step; // 第i步卷积的总偏移值
   float sign = o < 0.0f ? -1.0f : 1.0f;
   // 将当前的range和最大的Range的比值存入alpha通道
   kernelList.Add(new Vector4(0.0f, 0.0f, 0.0f, range * sign * Mathf.Abs(Mathf.Pow(o, exponent))));
}

采用Mathf.Pow(o, exponent)，正是利用到高分分布的特性，中心周围点的贡献度最大，越往两边贡献度越低（类似于重要性采样）

计算权重：kernelList[i].xyz

计算每卷积范围后，需要知道卷积范围内，每个点的比重

// 计算Kernel权重
for (int i = 0; i < nTotalSamples; ++i)
{
   float w0 = i > 0 ? Mathf.Abs(kernelList[i].w - kernelList[i - 1].w) : 0.0f;     //左邻居的步长差
   float w1 = i < nTotalSamples - 1 ? Mathf.Abs(kernelList[i].w - kernelList[i + 1].w) : 0.0f; //右邻居的步长差
   float area = (w0 + w1) / 2.0f;  // 该采样点的所占长度

   Vector3 t = area * SeparableSSS_Profile(kernelList[i].w, falloffColor); // 用卷积步长计算得到颜色
   kernelList[i] = new Vector4(t.x, t.y, t.z, kernelList[i].w);
}

将最中间的元素移到数组的第一位，在后续计算Blur时，优化性能

// 使用nSamples / 2表示中间的位置，再从中间元素的位置开始，往前遍历半个数组，将所有元素都向后移动一个位置
// 确保kernel数组的中间元素处于最容易访问的位置，后续对kernel数组的操作更高效
Vector4 centerKernel = kernelList[nTotalSamples / 2];
for (int i = nTotalSamples / 2; i > 0; --i)
{
   kernelList[i] = kernelList[i - 1];
}
kernelList[0] = centerKernel;

将权重归一化，即所有点的权重和为1

// kernelList[i].xyz = subsurfaceColor * kernelList[i].xyz / ∑kernelList[i].xyz
Vector4 sum = Vector4.zero;
for (int i = 0; i < nTotalSamples; ++i)
{
   sum.x += kernelList[i].x;
   sum.y += kernelList[i].y;
   sum.z += kernelList[i].z;
}
for (int i = 0; i < nTotalSamples; ++i)
{
   kernelList[i] = new Vector4(kernelList[i].x / sum.x, kernelList[i].y / sum.y, kernelList[i].z / sum.z, kernelList[i].w);
}

将权重比和一个可开放且可修改的参数lerp，从而手动控制漫反射效果

var tempKernel = kernelList[0];
tempKernel.x = Mathf.Lerp(1.0f, kernelList[0].x, subsurfaceColor.x);
tempKernel.y = Mathf.Lerp(1.0f, kernelList[0].y, subsurfaceColor.y);
tempKernel.z = Mathf.Lerp(1.0f, kernelList[0].z, subsurfaceColor.z);
kernelList[0] = tempKernel;

for (int i = 1; i < nTotalSamples; ++i)
{
   tempKernel = kernelList[i];
   tempKernel.x *= subsurfaceColor.x;
   tempKernel.y *= subsurfaceColor.y;
   tempKernel.z *= subsurfaceColor.z;
   kernelList[i] = tempKernel;
}

分离的高斯模糊

分为水平和竖直的模糊

缩放系数：

_SSSSMaterial.SetFloat("_DistanceToProjectionWindow", 1.0f / (0.5f * Mathf.Tan(math.radians(renderingData.cameraData.camera.fieldOfView))));

float  SSSSIntensity = _SSSSScale * _CameraDepthTexture_TexelSize.x;    // 水平方向为_CameraDepthTexture_TexelSize.x，竖直方向为_CameraDepthTexture_TexelSize.y
float4 sceneColor = _MainTex.SampleLevel(Smp_ClampU_ClampV_Linear, uv, 0);

float3 SSSSBlur(float3 sceneColor, float2 uv, float2 SSSSIntensity)
{
    float eyeDepth = GetLinearEyeDepth(GetDeviceDepth(uv));
    float blurLength = _DistanceToProjectionWindow / eyeDepth;  // 缩放系数

    float2 uvOffset = blurLength * SSSSIntensity;
    float3 blurColor = sceneColor * _KernelArray[0].rgb;

    [unroll(64)]
    for(int i = 1; i < _Sample_Nums; ++i)
    {
        float2 SSSSUV = uv + uvOffset * _KernelArray[i].a;

        float3 SSSSSceneColor = _MainTex.SampleLevel(Smp_ClampU_ClampV_Linear, SSSSUV, 0).rgb;
        float  SSSSEyeDepth = _CameraDepthTexture.SampleLevel(Smp_ClampU_ClampV_Linear, SSSSUV, 0);

        float DPTimes300 = _DistanceToProjectionWindow * 300;    // 经验参数
        float SSSSSCale = saturate(DPTimes300 * SSSSIntensity * abs(eyeDepth - SSSSEyeDepth));
        SSSSSceneColor = lerp(SSSSSceneColor, sceneColor, SSSSSCale);
        blurColor += _KernelArray[i].rgb * SSSSSceneColor;
    }

    return blurColor;
}

Diffuse光照选择

UE选择的是Disney提出的Burley模型

// [Burley 2012, "Physically-Based Shading at Disney"]
// Lambert漫反射模型在边缘上通常太暗，而通过尝试添加菲涅尔因子以使其在物理上更合理，但会导致其更暗
// 根据对Merl 100材质库的观察，Disney开发了一种用于漫反射的新的经验模型，以在光滑表面的漫反射菲涅尔阴影和粗糙表面之间进行平滑过渡
// Disney使用了Schlick Fresnel近似，并修改掠射逆反射（grazing retroreflection response）以达到其特定值由粗糙度值确定，而不是简单为0
float3 Diffuse_Burley( float3 DiffuseColor, float Roughness, float NoV, float NoL, float VoH )
{
    float FD90 = 0.5 + 2 * VoH * VoH * Roughness;
    float FdV = 1 + (FD90 - 1) * pow5( 1 - NoV );
    float FdL = 1 + (FD90 - 1) * pow5( 1 - NoL );
    return DiffuseColor * ( (1 / PI) * FdV * FdL );
}

效果

仅漫反射无高光，再SSSS效果如下：

高光

高光部分有些繁琐，不能走正常的延迟渲染流程，需要修改延迟光照算法，先计算diffuse，然后SSSS，最后再进行延迟光照来计算specular。重点在于specular算法如何抉择

NDF

UE的皮肤渲染采用双镜叶高光（Dual Lobe Specular），由两个独立的高光镜叶组成，各自使用不同的roughness——即使用两个独立的roughness，计算两个相同的NDF，再对两个roughness混合，UE的混合比例是固定的：主roughness为0.85，次roughness为0.15

half  _LodeA;
half  _LodeB;

// GGX / Trowbridge-Reitz
// [Walter et al. 2007, "Microfacet models for refraction through rough surfaces"]
// 在流行的模型中，GGX拥有最长的尾部。而GGX其实与Blinn (1977)推崇的Trowbridge-Reitz（TR）（1975）分布等同。然而，对于许多材质而言，即便是GGX分布，仍然没有足够长的尾部
float NDF_GGX( float roughness2, float NoH )
{
    const float a2 = pow2(roughness2);
    const float NoH2 = pow2(NoH);
    const float d = PI * pow2(NoH2 * (a2 - 1.f) + 1.f);

    if(d < FLT_EPS) return 1.f;

    return a2 / d;
}

float lobeARoughness = brdfData.roughness * brdfData.LobeA;
lobeARoughness = lerp(1.f, lobeARoughness, saturate(brdfData.opacity * 10.0f));
float lobeBRoughness = brdfData.roughness * brdfData.LobeB;
lobeBRoughness = lerp(1.f, lobeBRoughness, saturate(brdfData.opacity * 10.0f));
float lobeMix = 0.15f;

// AverageRoughness: roughness
float AverageAlpha2 = Pow4(AverageRoughness);
float Lobe0Alpha2 = Pow4(Lobe0Roughness);
float Lobe1Alpha2 = Pow4(Lobe1Roughness);
float D = lerp(NDF_GGX(Lobe0Alpha2, brdfData.NoH), NDF_GGX(Lobe1Alpha2, brdfData.NoH), LobeMix);

几何项

使用 Joint-Smith项（联合史密斯）的一种近似高效方案：

float Vis_SmithJointApprox( float a2, float NoV, float NoL )
{
    float a = sqrt(a2);
    float Vis_SmithV = NoL * ( NoV * ( 1 - a ) + a );
    float Vis_SmithL = NoV * ( NoL * ( 1 - a ) + a );
    return 0.5 * rcp( Vis_SmithV + Vis_SmithL );
}

float G = Vis_SmithJointApprox(AverageAlpha2, brdfData.NoV, brdfData.NoL);

Fresnel

用普通的Schlick-Fresnel方案即可：

float3 SchlickFresnel(float HdotV, float3 F0)
{
    float m = clamp(1-HdotV, 0, 1);
    float m2 = m * m;
    float m5 = m2 * m2 * m; // pow(m,5)
    return F0 + (1.0 - F0) * m5;
}

float3 F = SchlickFresnel(brdfData.HoV, brdfData.F0);

效果

得到的效果如下：

项目链接

https://github.com/chenglixue/Separable-Subsurface-Scatter/tree/main

Reference

剖析Unreal Engine超真实人类的渲染技术Part 1 - 概述和皮肤渲染

简谈实时皮肤渲染之SSSS(可分离次表面散射模型-SeparableSSS)

【详解4S皮肤】可分离次表面散射Separable SSS（含透射项）

Azusa

更好的次表面渲染 Separable Subsurface Scatter

前言