【学习】快速次表面散射（SSS）的原理与实现

本文链接：https://blog.csdn.net/coldbluer/article/details/138170142

一直以来都对次表面散射(Subsurface Scattering,简写为SSS)的具体实现非常感兴趣。所幸最近有时间，研究了一下网上的视频和文章，这里简单作一下记录。自己也实现一下试试。
基本的实现最初参考了B站视频【TA进阶】快速次表面散射（SSS）的原理与实现，但是后面基本上都在查文章资料了，也会补充一些自己的思考以及检索到的一些文章。
参考文章列表：
[图形技能树-渲染]次表面散射
 《GPU Gems》第16章：对次表面散射的实时近似
 实时皮肤渲染综述
 实时皮肤渲染——BTDF
基于物理着色（四）- 次表面散射

SSS材质，一般用于表现玉石、蜡烛、人体皮肤的材质效果。在强光照射下其材质特性极为明显。这种效果
在这里插入图片描述

次表面散射原理

在这里插入图片描述
考虑光线进入SSS物体的三种宏观行为抽象。首先是大家最熟悉的BRDF（上图b），关注模型表面的漫反射、高光、环境光等效果；BSSRDF（Bidirectional Surface Scattering Reflectance Distribution Function，上图a）更多关注光源和相机在物体同一侧时光源打入介质的晕染效果；BTDF（Bidirectional Transmittance Distribution Function，上图c）则更多关注物体间隔在相机和光源中间，光线透射物体的效果。
要想实现良好的SSS效果，最完美的方法肯定是采用光线追踪来实现，但是这无疑需要消耗巨量的性能，很难直接应用于实际开发。接下来只考虑各分量的简易实现方法。
BRDF就不说了，按照传统的光照模型计算即可。这里主要考虑BSSRDF和BTDF。我学习的视频主要针对移动端进行快速SSS实现，主要是实现了BTDF简化方法。

BTDF的快速实现方法

视频中给出一个公式：
SSS = 背光 + 扰动 + 扩散
个人理解背光+扰动即BTDF，而扩散属于BSSRDF的部分
在这里插入图片描述

如果只考虑单纯的背光，则根据Lambert定律，观察方向接收到的背光强度为，V与-L的点乘
背光强度 = DOT（V,-L），由于背光情况下DOT(N,L)是背面不受光照，贡献为，这里暂时不考虑

接下来是扰动项，视频中给出示意图和公式如下图所示：
在这里插入图片描述
其中δ为扰动项。
这里的公式非常简洁，但是具体的实现一开始让我不太能理解，为什么要求L和N的半途矢量H，然后再用-H与V点乘得到结果？一开始以为H是像Blinn-Phong中的半途矢量类似，如图示中用H去近似L在表面镜像位置的方向，模拟此方向有另外一个光源，但是为什么最终计算要用-H来点乘V呢？
这里我找到了这个方法的原始资料：
实时皮肤渲染——BTDF
Fast Subsurface Scattering in Unity (Part 2)
他们引用的方法是GDC2011的文章Approximating Translucency for a Fast, Cheap and Convincing Subsurface Scattering Look
仔细看了一下并没有讲原因，原始文章暂时没翻墙去读。

这里我斗胆提出一个自己的理解，如有谬误欢迎指正
在不考虑扰动的情况这个背光非常好理解，就是光线在均匀介质中的传播。
为什么会有扰动呢，背光的偏折方向为什么和法线方向有关？这里我想到了折射定律，SSS材质这里的的光照虽然更复杂，但是作简化考虑也需要符合折射定律！
考虑到这种BTDF的观察通常是在空气中观察目标物体，也就是说光线是从高IOR的材质出射到低IOR的材质，不考虑折射时有θi = θt，然后目标物质折射率越大则θi越大，相当于实际的出射方向向靠近平面的方向扰动，即法线的反方向！
这里使用δ控制扰动的程度，向下扰动得越明显，δ的值越大，有点像用δ来做不同折射率下出射方向的一个插值
由此近似可以得到Lout = -(L+δN）
在这里插入图片描述
到这里得到V点乘-H可以得到BTDF在观察方向看到的分量
后面就简单了，参考Phong模型的高光计算，我们用幂次来加剧这里方向对结果的变化影响。使用公式
背光强度 = pow(saturate(V·-(L+δN)),p)*s
其中p是控制变化程度的参数，类似Phong模型里高光的gloss参数
到这里就是BTDF的方向计算就完成了，然后还要考虑在传播方向上强度的衰减
真正的光照模型是非常复杂的，这里也是要考虑一个简化的模型，即根据透射距离决定吸收程度。比如在背光时人头上的皮肤以耳朵的透光最为明显。

接下来的问题就是如何获取这个距离呢，目前考虑有两种思路：

1.厚度值贴图

视频教程中的思路，也就是对模型提前烘培一个厚度值贴图，将采样而得的厚度值映射到一个衰减值，将最终的衰减值乘到上面的的背光强度可得最后的表达式。这里样图是厚的地方白色，薄的地方黑色，那么结果直接使用采样到的厚度thickness即可，如果贴图是反过来的1-thickness即可
在这里插入图片描述

最终的背光贡献 = 灯光衰减值*(背光强度+ambient)*thickness
这里的ambient是为了防止最终结果过黑而加上的，毕竟上面点乘的结果很容易截取到0，容易导致大面积的看不到背光效果。
严格的来说，这里的厚度值在每一点是一个统一的值，而真实情况中同一点从不同光照方向穿透的厚度并不一样，但是这里为了快速应该是做了一个近似。如果想进一步准确记录不同方向的厚度，可以考虑使用球谐函数？一张贴图四个通道可以记录到2阶的球谐函数。再高阶的话代价就有点太大了。这样烘培工具应该需要单独写，这里留个坑。
关于厚度贴图的烘培,Blender和SubstancePainter等软件都可以直接烘培，但是我查了一下有些软件中好像白色代表厚黑色代表薄，无所谓直接反转结果即可。
原始论文中也给出了一种生成厚度图的方法：
1.反转表面法线
2.烘培AO
3.将烘培结果颜色翻转存到贴图中
在这里插入图片描述
这个方法说实话不是很理解，回头实际做个测试看看

2使用深度映射模拟吸收

方法来自这里，有点像阴影贴图。
为了模拟这种效果，我们需要测量光在物质中传播的距离。而估算传播距离可以使用深度映射技术，此技术非常类似于阴影映射，而且可用于实时渲染：如下图，在深度遍中，我们以光源为摄像机渲染场景，存储从光源到当前片元的距离。在渲染遍中，把当前片元坐标转换到光源空间，采样深度贴图得到当前片元对应入射点到光源的距离(di)，再计算出当前片元(出射点)到光源的距离(do)，后者减去前者，即光线在物体中传播的距离(s)。
在这里插入图片描述
这个方法优点是不需要额外烘焙厚度贴图，而且考虑到了光的入射方向，能体现不同方向的厚度不同，同时也可用于动画物体的实时计算（物体能变形的话厚度贴图必然失效）。但也存在缺点：只对凸的物体有效，不能正确描绘物体上的洞（问题不大）。
文章中有详细的代码和方法，可自行阅读。

BSSDF的实现

与Blinn-Phong模拟BRDF类似，BSSRDF也有几种经典的拟合模型，人们根据观察认为，BSSRDF对最终像素的主要作用在于一定程度的blur和表面颜色根据光源颜色产生一些偏移。
个人对此的理解是，BSSDF考虑光线照射的部分，其中会有部分散射光线扩散到旁边较暗的部分中去，被直接照射的部分之间光线也会互相影响，其结果就极像对最终图像做一个模糊然后再你混合叠加。
以下是一些blur与颜色的变换方法来模拟BSSRDF：
texture space blur：在原始纹理中计算BSSRDF的等价blur做成新纹理
screen space blur：同样是blur，不过是先光栅化到屏幕空间再realtime计算
pre-integrated预积分方法：部分变换中间结果预先存储在表里，本质还是blur
light warping：用经验表达式变换颜色

texture space blur

在这里插入图片描述
上图是这种做法的流程图。先计算整张皮肤漫反射的结果（比如用lambert diffuse）并展开成一张纹理(也叫转换到texture space)，进行若干次blur（demo用了六次），将每次blur结果加在一起作为diffuse，最终加上specular。
具体的实现在《GPU Gems》第16章：对次表面散射的实时近似中16.4部分有详细的解释和代码，可以自行查阅

screen space blur

还是blur，只不过不在texture space做，而是在screen space做，将需要blur的部分用stencil buffer标记。这样的话开销是固定的，在屏幕几乎没有皮肤时最不值。
在这里插入图片描述

pre-integrated预积分方法

前面的两种方法都需要进行blur,也需要Blit，其开销很大，在移动平台上显然使用起来有困难。预积分方法完全从结果出发，将结果烘焙在贴图上，以光照的变量作为索引（横轴竖轴），制成一个查找表（look-up table, LUT）。
在这里插入图片描述
如上图所示，左侧是LUT，输入横轴纵轴的值，对应像素的颜色就是diffuse color（当然后续还要乘以物体颜色光颜色之类的）。横轴是N dot L的取值，他表示法线与光线夹角的cos值。图中的这些颜色怎么理解呢？我们忽略纵轴，只看NdotL对颜色值的影响，可以看到右侧NdotL>0时，N和L夹角小于90度，表明光照在物体表面，所以结果是白色；反之，NdotL<0，N和L夹角大于90度，意味着这个点在物体背面，所以光照结果是黑色的。
纵轴1/r表示曲率，对应到右侧的人脸上就是额头、脸颊相对平整，曲率低，而鼻沟、眼皮和耳朵曲率高。曲率高的地方LUT的y值更大，此时NdotL在明暗交界处会偏红，暗部会更亮，也就是更透红光。
为了使用该方法，我们需要提前渲染一张角色的curvature map（曲率贴图），就像图中的人脸一样。Substance painter等dcc软件可以生成。进行皮肤渲染的时候，在shader中计算NdotL作为x值（[-1,1]映射到[0,1]），读取curvature map作为y值，读取look-up table的颜色值，这个值在光照计算中替换掉原来的光照模型的NdotL。也就是说新的NdotL不仅是[0,1]的，而且是有色相的。
单独观察这个查找表能看出来，先不考虑曲率高低的区别只看横轴，发生BSSDF的部分主要集中在光照从明到暗变化的中间区域，对其作一个连续带颜色的整体提亮即可近似的模拟出BSSDF，由此可以引出 Light Warping方法

Light Warping方法

这是一种比预积分更省的方法，它不需要look-up table也不需要curvature map，但是也像预积分方法一样用重新映射NdotL的手段模拟次表面散射。如下图所示，我们直接对NdotL的结果做一个数学运算，用一个参数控制暗处的提亮。
在这里插入图片描述
warp的公式，就是图中所写的y=(x+warp)/(1+warp)，是怎么来的呢？这是纯粹的人为设计，可以看到warp=0时，NdotL输入是多少输出也是多少，当warp增加时，暗部会越来越亮，但最大值还是1，不会溢出。想要色相变化，再乘以一个颜色参数即可。
在照明度接近0时，可以表达那种倾向红色的微小颜色漂移，这是模拟皮肤散射的一种廉价方法。
其实看到这里还有一个思路，是不是可以用贴图映射来做？取一个一维贴图，根据NDotL去采样即可，不过读取纹理的性能肯定比不上后面的近似计算。

在这里插入图片描述
这里贴一下知乎上的代码

half4 frag (v2f i) : SV_Target
{
    half3 N = normalize(i.worldNormal);
    half3 L = normalize(UnityWorldSpaceLightDir(i.worldPos));
    half3 V = normalize(UnityWorldSpaceViewDir(i.worldPos));
    half3 H = normalize(L + V);

    half NdotL = saturate(dot(N, L));
    half NdotH = saturate(dot(N, H));

    half diffuse = NdotL;
    half specular = pow(NdotH, _Shininess);
    half scatter = 0;
#ifdef _SUBSURFACESCATTERING
    half NdotL_wrap = (NdotL + _Wrap) / (1 + _Wrap);
    diffuse = max(0, NdotL_wrap);
    half range1 = smoothstep(0, _ScatterWidth, NdotL_wrap);
    half range2 = smoothstep(_ScatterWidth * 2, _ScatterWidth, NdotL_wrap);
    scatter = range1 * range2;
#endif

    half4 col = tex2D(_MainTex, i.uv);
    half4 lit = diffuse * _LightColor0 * col + specular * _SpecularColor + scatter * _ScatterColor;

    return lit;
}

可以看出，核心公式就是这里(NDotL+ wrap)/(1+wrap)，将原本的NDotL替换成(NDotL+ wrap)/(1+wrap)的形式，但是这里的range1 range2我一开始没看明白。
使用扰动后的NDotL做插值，这里的算式简化后scatter=1-(1-NdotL_wrap)²，即y=1-(1-x)²，看起来是将NdotL_wrap结果中间部分整体提亮而且仍然保持值域[0,1]不变，因为SSS的BSSDF主要影响中间部分，这里也是强化影响同时满足能量守恒。
scatterwidth的平方项控制最终的结果数值大小。
在这里插入图片描述
这篇文章还给出了几个其他的计算公式：
f(θ)=0.25*pow((cos(θ)+1),2)
f(θ,a)=pow((cosθ+a)/(1+a),1+a)
可以作为不同的效果进行测试。这里主要参考了《GPU Gems》第16章：对次表面散射的实时近似
下面这个公式的实现代码如下：

inline half4 LightingStandardTranslucent(SurfaceOutputStandard s, half3 viewDir, UnityGI gi)
{
    // PBR
    half4 pbr = LightingStandard(s, viewDir, gi);

    // Subsurface Scattering
#ifdef _SUBSURFACESCATTERING
    float3 L = gi.light.dir;
    float3 N = s.Normal;
    half NdotL = dot(N, L);
    half alpha = _SSS_Radius;
    half theta_m = acos(-alpha); // boundary of the lighting function
    half theta = max(0, NdotL + alpha) - alpha;
    half normalization_jgt = (2 + alpha) / (2 * (1 + alpha));
    half wrapped_jgt = (pow(((theta + alpha) / (1 + alpha)), 1 + alpha)) * normalization_jgt;
    //half wrapped_valve = 0.25 * (NdotL + 1) * (NdotL + 1);
    //half wrapped_simple = (NdotL + alpha) / (1 + alpha);
    half3 subsurface_radiance = _SSS_Color * wrapped_jgt * pow((1 - NdotL), 3);
    pbr.rgb = pbr.rgb * (1 - _SSS_Value * _SSS_Value) + gi.light.color * subsurface_radiance * _SSS_Value;
#endif

    return pbr;
}

这里half theta = max(0, NdotL + alpha) - alpha这句将cosθ限制在(-alpha,1)之间，看起来这里似乎考虑了光线从背面来的散射？
normalization_jgt = (a+2)/(2a+2)不是很好理解，看起来似乎是对wrapped_jgt做一些修正，我将alpha取0.1，0.3，…，0.9然后对wrapped_jgt绘制了一下函数,整体是一个上扬的状态。而normalization_jgt是图示中下降的函数。实际上光照强的部分BSSDF不应该明显？所以相乘做一个修正？这里不是很理解，欢迎大佬进行一下指正。
在这里插入图片描述
上面的代码最后几句还有pow((1 - NdotL), 3)这个参数，我看了在引用的国外文章中并没有，也许是Trick。至于(1 - _SSS_Value * _SSS_Value)这个我倾向于认为作者是写错了，这样就能量不守恒了。