什么是次表面散射？

在最为常见的BRDF中，描述了光线打到物体表面，如何反射；但针对不同物体，光线并非全部都反射或被吸收。存在一些物体，光线会穿过它们，发生折射从另一面射出，这种现象被称之为透射，即为BTDF；BSDF对BRDF 和 BTDF做了融合，更有BSSDF基于BSDF进行了升级

其中，次表面散射用BSSDF来表示。它主要用于模拟皮肤、玉石、牛奶等半透光性物质。对于这些物质，光线打中物体后的折射光不太可能直接从另一侧透射出人体，更可能在某一层介质中经过了多次光路改变，最后出射光以一个相对随机的角度再次回到入射光一侧的空气中

在次表面散射中，漫反射来自于这些重新回到空气中的折射光。当观察尺度足够小，便不能近似的认为这些漫反射光的出射点恰好位于入射光所在的入射点处，而是如上图所示，出射光被认为分布在入射光周围一定范围内（绿色圆圈）。如果上图中绿色圆圈的区域恰好落在某个像素内，得到的依然是普通的漫反射效果

对于皮肤这类特殊的次表面散射物质，当光线穿透油脂和上皮层，遇到人体皮肤中富含血管和组织液的真皮部分，光线很可能会沿着局部的均质(血液，组织液)做较长范围的折射迁移，从而在绿色圆圈外射出

皮肤组成

为了还原皮肤需要了解皮肤如何组成的，以及它们对光线如何作用。皮肤的组成结构大致如下：

表面油脂层（Thin Oily Layer，蓝色部分）：模拟皮肤的反射
表皮层（Epidermis，浅橙色部分）：模拟次表面散射的贡献层
真皮层（Dermis，深橙色部分）：模拟次表面散射的贡献层

皮肤存在两种特性：

随机性：由于材质内部介质通常是不连续的，以及大量毛细血管和里面的血液构成，光路在内部的折射/反射往往是无规律可循的，导致最终出射时方向的随机化
电介质的着色性(漫反射光会被染上颜色)：光线在物体内部传播时，一部分波长的光能被吸收

次表面散射算法

SSS

类似于 BRDF，人们专门为SSS引入了新的方程BSSDF，公式如下： $L_o(x_o, o) = \int_A R(| x_i - x_o |) \int_Ω L_i(x_i, i)(i \cdot n_i) di dA$

$L_o(x_o, o)$ ：最终得到的出射位置与出射角度下的亮度 $L_o$
在给定一个计算面积 $A$ 的前提下，对半球区域 $Ω$ 所有入射光进行积分
- 在积分时需要考虑光线对表面投影造成的衰减—— $dot(i, n_i)$ ，其中 $L_i(x_i, i)$ 表示入射光线的入射位置和入射角度下的亮度 $L_i$
- 当前积分最终得到的是单位面积的irradiance
随后是对面积A积分，函数R和半球区域 $Ω$ 的irradiance的乘积
- 函数R的实参是 $| x_i - x_o |$ ，意为 $x_i$ 与 $x_o$ 的距离
- R为模糊函数（各种模糊都可以，但需保证能量守恒），用于后续近似散射函数，但该函数有个限制，由于该函数只关心入射点和出射点，不关心出射方向和入射方向，因此只能用于各向同性物体
- 当前积分最终得到的是面积A的irradiance

预积分(diffuse)

条件

因为需要根据某一方向得到某个方向出射的irradiance，所以可以看到BSSDF的计算方式是比较昂贵的。有没有什么方法可以优化呢？对于SSS常用到的方法是预积分——将BSSDF的部分或整体预先计算出来，并以一定方式保存进纹理

但为了确保预计算的内容不被模型和表面参数所束缚，就必须限定积分对象的材质，并且仅预计算那些不变的组成部分。要满足以上条件，需要做到以下三点：

弧度变化的模型
法线扰动：表面法线完全没有扰动的话，那么微表面区域完全镜面和平坦，SSS失去了形成细节特征的条件；如果表面法线扰动过于锐利，那么看起来效果会不好
光线遮挡：必须包含阴影，从而形成局部由暗至明的转向，避免全方位的亮光造成SSS效果无法被察觉

过程

预积分过程分为三步：

处理光照：确定入射函数 $L_i$
处理微表面：确定散射分布函数 $R$
处理遮蔽：计算光照的投影 $dot(i, n_i)$

对象

预积分积分的是皮肤表面漫反射分量，与BRDF注重高光不同

推导

由于预积分公式有点长，因此这里主要介绍最重要的部分

先从物理模型开始：

上图(来自GPU Pro2)中，几个点需要注意：

左图表示BSSDF中的R函数（近似的散射函数），其中x意为出射点和入射点的距离，返回值为贡献强度
- 红蓝绿三色的线条代表不同的R函数
右图为抽象的几何模型（可以想象为人的头），用于计算N方向的漫反射量
- 法线N和入射光L的夹角设为 $θ$
- 圆形半径设为 $r$
- N+θ意为从N开始逆时针旋转θ角对应的位置

这里先直接给出预积分公式： $D(θ, r) = \frac { ∫^{\frac{Π}{2}} _{-\frac{Π}{2} } cos(θ+x) R(2r sin(\frac{x}{2})) dx }{ ∫^{\frac{Π}{2}} _{-\frac{Π}{2} } R(2r sin(\frac{x}{2})) dx }$

与上图不同的是，该公式将积分区域缩小到了半球区域，多了参数r。这是因为上图左侧主要描述x对散射的贡献强度

如下图所示，对于圆形上任意一点Q，它的入射光强度如下：假设OQ和ON的夹角为x(区域 $[-\frac{Π}{2}, \frac{Π}{2}]$ )，可以求得Q点的入射光强度为 $L_Q = L_i(dot(OQ, OL)) = L_i cos(θ+x)$

但主要目标并不是求半圆上任意一个点的入射光强度，而是求半圆上所有入射光散射后对N的贡献度。因此还需要一个散射率 $q(x)$ ，它表示Q点到N点的散射概率；但这还不够，因为这里点Q的弧长并不为1(因为连续概率密度函数的原因，不能直接说某个点的概率为多少，但可以用一小段区间来表示)，所以需要利用微分，可得点Q的弧长为 $Δx$ ，对N的贡献度为 $L_i cos(θ + x) q(x) dx$

因为半球上很可能存在许多像Q点类似的点，因此需要对半球做个积分，可得： $D(θ) = ∫ ^{\frac{Π}{2}} _{- \frac{Π}{2}} L_i cos(θ + x) q(x) dx$

但在前面提到过，近似的散射函数R，它的参数是入射点和出射点的距离，而目前x表示的是OQ和ON的角度差，这也并不难，将角度差转换为点Q和点P的距离即可—— $length = 2r sin(\frac{θ}{2})$ 。可以得到 $D(θ) = ∫ ^{\frac{Π}{2}} _{- \frac{Π}{2}} L_i cos(θ + x) q(2r sin(\frac{x}{2})) dx$

公式剩下的部分不再推导，主要围绕两个点展开：

近似的散射函数R需要让结果看着不错
近似的散射函数R需要满足能量守恒

生成预积分图

这里就懒得写成工具了，直接用脚本配合compute shader实现

using System;
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using UnityEngine.Rendering;

public class PreIntegrateSSSLUT : MonoBehaviour
{
    public int resolution = 1024;
    public RenderTexture SSSLUTRT;
    public Texture2D SSSLUT;

    void Start()
    {
        if (SSSLUTRT != null) RenderTexture.ReleaseTemporary(SSSLUTRT);
        if (SSSLUT != null) Texture2D.Destroy(SSSLUT);
        ComputeShader computeShader = Resources.Load<ComputeShader>("SSS/CS_PreIntegrateSSSLUT");
        if (computeShader == null)
        {
            Debug.LogError("PreIntegrateSSSLUT compute shader is missing");
        }

        SSSLUTRT = RenderTexture.GetTemporary(resolution, resolution, 0, RenderTextureFormat.ARGBHalf, RenderTextureReadWrite.Default);
        SSSLUTRT.enableRandomWrite = true;
        SSSLUT = new Texture2D(resolution, resolution, TextureFormat.RGB24, true, true);

        var kernelIndex = computeShader.FindKernel("PreIntegrateSSSLUT");
        computeShader.SetTexture(kernelIndex, "_RW_OutputTex", SSSLUTRT);
        computeShader.SetVector("_SSSLUTSize", new Vector4(resolution, resolution, 1f / resolution, 1f / resolution));
        computeShader.Dispatch(kernelIndex, resolution / 8, resolution / 8, 1);

        RenderTexture.active = SSSLUTRT;
        SSSLUT.ReadPixels(new Rect(0, 0, resolution, resolution), 0, 0);
        RenderTexture.active = null;

        System.IO.File.WriteAllBytes("Assets/Resources/Tex/SSSLut.jpg", SSSLUT.EncodeToJPG());
    }

    private void OnDestroy()
    {
        RenderTexture.ReleaseTemporary(SSSLUTRT);
        Texture2D.Destroy(SSSLUT);
    }
}

#pragma once
#include "Assets/Resources/Library/Common.hlsl"

#pragma kernel PreIntegrateSSSLUT

RWTexture2D<float4> _RW_OutputTex;

float4 _SSSLUTSize;

float3 Gaussian(float v,float r)
{
    return 1.0/sqrt(2.0 * PI * v) * exp(-(r*r)/(2.0*v));
}
float3 ApproximateScatter(float r)
{
    return float3(0.0, 0.0, 0.0)
        + Gaussian(0.0064, r) * float3(0.233, 0.455, 0.649)
        + Gaussian(0.0484, r) * float3(0.100, 0.336, 0.344)
        + Gaussian(0.187, r) * float3(0.118, 0.198, 0.0)
        + Gaussian(0.567, r) * float3(0.113, 0.007, 0.007)
        + Gaussian(1.99, r) * float3(0.358, 0.004, 0.0)
        + Gaussian(7.41, r) * float3(0.233, 0.0, 0.0);
}

float3 GetLUT(float2 uv)
{
    float NoL = uv.x;
    float oneOverR = uv.y;

    float theta = acos(NoL * 2.f - 1.f);
    float radius = 1.f / oneOverR;

    float3 totalWeight = 0.f, totalLight = 0.f;
    for(float x = -PI / 2; x <= PI / 2; x+= PI * 0.001f)
    {
        float sampleDistance = abs(2.f * radius * sin(x / 2));
        float3 weight = ApproximateScatter(sampleDistance);

        totalWeight += weight;
        totalLight += saturate(cos(theta + x)) * weight;
    }

    float3 result = totalLight / totalWeight;

    return result;
}

[numthreads(8,8,1)]
void PreIntegrateSSSLUT (uint3 id : SV_DispatchThreadID)
{
    float2 uv = ((float2)id.xy + 0.5f) / _SSSLUTSize;
    _RW_OutputTex[id.xy] = float4(GetLUT(uv), 1.f);
}

可以得到如下结果：

横坐标是NoL，纵坐标是 $\frac{1}{r}$

如何使用

如何采样该预积分图呢？其实就是需要知道NoL、 $\frac{1}{r}$

NoL：这部分很简单，只需要获得光源向量、物体的normal。但需注意，normal尽量使用法线贴图，这样可以采样法线贴图的mipmap，防止出现生硬的边缘
$\frac{1}{r}$ ：因为 $\frac{1}{r}$ 范围在[0,1]，所以这意味着r的范围在 $[1, \infty]$ 。当r趋近于最小单位1时， $\frac{1}{r}$ 趋近于1，从上图可以看出，得到的散射效果都提高了，且带有红色；当r越来越大， $\frac{1}{r}$ 趋近于0，从上图可以看出，得到的散射效果更接近于普通的明暗光照，并无明显的SSS效果

但问题在于如何得到 $\frac{1}{r}$ ?这里用到的是初中的数学知识——映射
如下图所示，先计算当前pixel和它附近3个pixel(因为GPU在计算时，是2x2个区域计算的)的normal和position的差值( $ΔN、ΔP$ )，因为normal在实际计算中会被归一化，所以这里可以得到一个映射关系—— $\frac{ΔN}{ΔP} = \frac{1}{r}$

现在剩下的问题是如何计算 $ΔN、ΔP$ ?

这里就需要用到shader中的一个函数fwidth()，它可以计算当前2x2区域中，pixel的特定值在水平和竖直方向上的差值和。其中，fwidth()又由ddx()、ddy()计算而来，因为ddx函数是计算水平方向的两个相邻像素的差值，ddy是计算垂直方向的两个相邻像素的差值，这意味着在2x2区域中，相邻的pixel的ddx()、ddy得到的结果是相同的，而对于fwidth()相同的pixel在对角线上
```
fwidth(v) = abs(ddx(v))+ abs(ddy(v))
```
对于pixel的数据，如这里的normal和position，如果想基于它们的变化来实现特定的功能，不能仅仅采用ddx或ddy，因为它只代表某一数据在x或y方向上的变化，为了更加全面需要用到fwidth()

那么计算 $ΔN、ΔP$ 的代码就比较明了了，大致实现如下：
```
float curve = saturate( _CurveFactor * (length(fwidth(normalWS)) / length(fwidth(positionWS))) );
float NoL = dot(blurNormalWS, lightDirection);
float3 diffuse = _SSSLUTTex.SampleLevel(Smp_ClampU_ClampV_Linear, float2(NoL * 0.5f + 0.5f, curve), 0);
```
$\frac{1}{r}$ 也可以不用实时计算，可以预先将模型表面各点的曲率烘焙到uv贴图上

Specular

皮肤的高光和漫反射光不同，只有很少一部分的光在接触到皮肤表面后会被镜面反射，这种反射光具有若干性质，可以归纳如下：

反射光不带有任何颜色信息，这是电介质作为反射层的固有特性
反射光中的大部分来自于皮肤最表面的一层油脂层(Thin Oily Layer)
反射光遵循Fresnel效应，掠射角处反射量增大

高光的质感在于选择合适的NDF、Geometry、Fresnel函数，皮肤的渲染也不例外。这里主要谈及GPU Gem 3中的Kelemen/Szirmay-Kalos BRDF模型

在PBR中，对于渲染方程的高光项，主要由三个部分组成，以BRDF为例：

specularLight += lightColor[i] * lightShadow[i] * rho_s * specBRDF( N, V, L[i], eta, m) * saturate( dot( N, L[i] ) );

lightColor[i] * lightShadow[i]：直接光强度
rho_s：非物理因子用于控制整体强度
specBRDF：由NDF、Geometry、Fresnel三个函数组成
dot( N, L[i] )：控制入射光能的衰减强度

整体公式

公式采用一种更为简版的Cook-Torrance模型BRDF

P项：NDF项
F项：Fresnel项
余下部分：Geometry项

Fresnel

Kelemen/Szirmay-Kalos BRDF模型中的Fresnel函数依旧采用Schlick菲尼尔近似，当视方向V越与微表面法线H垂直所形成的反射光照强度越高

float fresnelReflectance( float3 H, float3 V, float F0 ) 
{   
    float base = 1.0 - dot( V, H );   
    float exponential = pow( base, 5.0 );   
    return exponential + F0 * ( 1.0 - exponential ); 
}

其中，F0表示当皮肤遇到垂直入射光照时的反射率，GPU Gem 3中用的值为0.028

Geometry

将Cook-Torrance模型的Geometry和 $\frac{1}{4NoLNoV}$ 连同一起简化为 $\frac{1}{dot(h,h)}$ ，其中 $h = v + l$ （未被归一化的半程向量）

NDF

对于Fresnel项和Geometry项而言，由于它们需要用到大量PS阶段的数据，因此它们只能放在PS阶段进行计算。但NDF项就不一样了，它只涉及到一个微表面法线H，可以对NDF进行预计算从而减少性能消耗

在Kelemen/Szirmay-Kalos中的NDF项本质是Beckmann：

其中， $α$ 为宏观法线N和微观法线H的夹角；m为roughness

预计算NDF的实现如下：

public enum Quality
{
    low = 256,
    middle = 1024,
    high = 2048,
    veryHigh = 4096
}
public Quality quality = Quality.middle;
private int resolution;
public RenderTexture NDFLUTRT;
public Texture2D NDFLUTTex;

void Start()
{
    if (NDFLUTRT != null) RenderTexture.ReleaseTemporary(NDFLUTRT);
    if (NDFLUTTex != null) Texture2D.Destroy(NDFLUTTex);
    ComputeShader computeShader = Resources.Load<ComputeShader>("SSS/CS_SSSNDFLUT");
    if (computeShader == null)
    {
        Debug.LogError("SSS NDF LUT compute shader is missing");
    }

    resolution = (int)quality;

    NDFLUTRT = RenderTexture.GetTemporary(resolution, resolution, 0, RenderTextureFormat.RFloat, RenderTextureReadWrite.Linear);
    NDFLUTRT.enableRandomWrite = true;
    NDFLUTTex = new Texture2D(resolution, resolution, TextureFormat.RFloat, true, true);

    var kernelIndex = computeShader.FindKernel("SSSNDFLUT");
    computeShader.SetTexture(kernelIndex, "_RW_OutputTex", NDFLUTRT);
    computeShader.SetVector("_NDFLUTSize", new Vector4(resolution, resolution, 1f / resolution, 1f / resolution));
    computeShader.Dispatch(kernelIndex, resolution / 8, resolution / 8, 1);

    RenderTexture.active = NDFLUTRT;
    NDFLUTTex.ReadPixels(new Rect(0, 0, resolution, resolution), 0, 0);
    RenderTexture.active = null;

    System.IO.File.WriteAllBytes("Assets/Resources/SSS/SSSNDFLut.jpg", NDFLUTTex.EncodeToJPG());
}

private void OnDestroy()
{
    RenderTexture.ReleaseTemporary(NDFLUTRT);
    Texture2D.Destroy(NDFLUTTex);
}

#pragma kernel SSSNDFLUT

RWTexture2D<float4> _RW_OutputTex;

float4 _NDFLUTSize;

float PHBeckmann( float NoH, float m ) 
{   
    float alpha = acos( NoH );
    float ta = tan( alpha );
    float val = 1.0 / (pow2(m) * pow4(NoH)) * exp(-(ta*ta) / pow2(m));   
    return val; 
} 

float GetNDFLUT(float2 uv)
{
    float NoH = uv.x;
    float roughness = uv.y;

    // Scale the value to fit within [0,1] – invert upon lookup. 
    return 0.5 * pow( PHBeckmann( NoH, roughness ), 0.1 ); 
}

[numthreads(8,8,1)]
void SSSNDFLUT (uint3 id : SV_DispatchThreadID)
{
    float2 uv = ((float2)id.xy + 0.5f) / _NDFLUTSize;
    _RW_OutputTex[id.xy] = GetNDFLUT(uv);
}

最终可以得到如下结果：

其中，横坐标代表NoH，纵坐标代表roughness

如何使用

照着一开始提到的公式，往里套即可

float3 halfVectorUnNor = light.direction + (GetCameraPositionWS() - posWS);
float NoH = dot(lightData.normalWS, halfVectorUnNor);

float NDF = pow(2.f * _SSSNDFLutTex.SampleLevel(Smp_ClampU_ClampV_Linear, float2(NoH, brdfData.roughness), 0), 10);
float F = SchlickFresnel(brdfData.HoV, brdfData.F0);
float G = dot(halfVectorUnNor, halfVectorUnNor);
float3 specular = max(NDF * F * rcp(G), FLT_EPS);

result += specular * brdfData.NoL * _SubsurfaceSpecularIntensity;