没什么技术力，可能会有错，希望各位大佬批评指正Orz
次表面散射（Subsurface scattering，SSS），是光在传播时的一种现象。光一般会穿透物体的表面，并在与材料之间发生交互作用而被散射开来，在物体内部在不同的角度被反射若干次，最终穿出物体[1]。对于大理石，皮肤，树叶，蜡和牛奶等材料，次表面散射对于提升材质的质感而言非常重要。
皮肤的三层渲染模型
因为散射这种效果太受材质内部性质影响了，Donner（2005）的论文《Light Diffusion in Multi-Layered Translucent Materials》[2]就证明了单层模型不足以表现出良好的渲染效果，因此提出了一个三层的皮肤渲染模型。对于皮肤而言，光路表现出如图的散射路线。光射入皮肤后会经过一系列内部反射，最后从另外一处射出来，而光射出时的颜色取决于其散射光路，对于人类皮肤而言，散射出的颜色主要是红色[3]。
看完皮肤结构再来看光。假设散射效果与光的入射角度和其与入射点距离存在某种函数关系，但对于一个质地均匀的材质，光在它上面所有方向上的散射都相同的，而且由于在物体内部散射时，入射光的方向性几乎丢失，所以散射效果与光照角度无关；基于这些推论，光的散射效果只与入射光强度和距离有关。
现在我们对着这个质地均匀物体射出一束光，左图中光源外围红色的一圈就是散射出的效果[4]。且不同颜色光具有不同的散射轮廓，如右图。从图上看，皮肤对红光的散射比绿光和蓝光都更强烈。在这里，我们用r表示观测点与光的入射点之间的距离。
左边这张图其实我们也能做。它是怎么来的呢，GPU Gems 3给出来的思路是，用多个高斯函数的和去拟合扩散曲线，即对扩散曲线
，有
。（曲线来源是Jensen的论文《A Practical Model for Subsurface Light Transport》[5]，暂时还没复现，先挖个坑）。这个找近似高斯函数和的方法就是用matlab之类的查找，论文里已经给出了一个
，只要找到
个高斯函数使得
这个积分取到最小值就好了。英伟达已经把皮肤的参数给找出来了[4]。
GPU Gems 3里给出了这样的拟合高斯和图:
不过我这边按照公式复现出来图是长这样的，我怀疑是因为这里要对RGB的各自的权值做归一化，把绿光和蓝光原本的散射权重给隐藏掉了，所以GPU gems 3里作的这张图应该是用未进行归一化的权值算的。
我这里顺便验证了一下GPU Gems 3中给出的各个高斯函数，分别做了几张Diffuse Profile，效果是对得上的。值得注意的是，这里用的扩散函数不是固定的，可以自己按需求选择[6]。我这里还拿了迪士尼的公式[7]做了个实验：
其实，根据上面的这几张Diffuse Profile，我们已经可以根据任意的输入
来查找对应的散射效果了。这张Diffuse Profile其实只存储了
向
映射的一维信息。然而，在实际做渲染时，我们面对的一般是平行光，目标皮肤模型也通常是个曲面，这就意味着我们一般面临的不是“观测点与入射点”的问题。
根据之前的文献和推理，既然光在射入材质时我们已经忽略了它的方向关系，那我们不妨建立一个“辐照度-距离-散射”的对应关系，在渲染中使用辐照度和距离查找某张LUT，得到散射效果。
在这个对应关系里，辐照度是很好求的，距离该如何求呢？对于平行光，我们已经无法像找激光入射点一样解决这个问题了，那不妨设整个材质是通透的，平行光给这个材质均匀地施加了光照，材质在任何一点都具有相同的散射性质，这样对于某个点接收的散射量，我们都可以根据这个点自身的信息求出，而不需要考虑模型整体的光照情况，这样点就满足
。这样操作虽然不够物理，但是计算量得到了减小，效果也能满足不错的近似。
数学推导
积分基本形式
首先，我们假设存在一个半径为
的圆，圆上所有点都具有相同的散射性质。那么对于圆上任意一点
，它将受到圆上所有点散射而来的光。现在，给这个圆施加平行光
，圆内将发生散射。
现在，设有一点
,其法线为
，
与
的法线
夹角为
，平行光线
与
的夹角是
。设散射系数为
。根据之前之前的假设，
点将受到圆上所有点，包括
点传递过去的散射。可以求得点
受到的直接光照就是
，由于法向量、光线向量为归一化向量，可得
。现在设平行光
的光强为1，则有
。那么根据之前的定义与推论，点
将给点
提供
的散射光，意思就是点将其受到的光照根据材质的散射系数散射出去。
如果不考虑二次散射，那么只有受到光照的点才能向四周发出散射，则
点受到的散射量就是圆上所有点的散射量之和，积分可得
。不过这样的形式显然太复杂了，实际上，对于未受到光照的点，其直接光照量为0，它乘上散射系数后也是0，积分的时候其实不必把这段单独处理，只需要把
限制到
即可，即：
。
好了，现在已经完成了对圆环散射的积分；不过毕竟是渲染三维物体，所以还需要对球面进行积分。虽然Panner在他的PPT[8]中认为在环上积分和在球面上积分效果区别不大，不过这里还是打算做到底把球面积分的路也走一遍。
如图，我们在对
环进行积分的同时对
所在的垂直于
环的
环积分。先前在积
环时我们设
为
进行积分；现在为了方便进行二重积分，我们设
为
，设
为
。
首先，对应原积分里的
，二重积分中我们需要用
表示光照量。由于这是一个积分，
又是平行光，因此我们可以视
与
构成的平面垂直于
与
构成的平面，并由此抽象出一个过同顶点
的三条棱两两垂直的三棱锥，如图。
在这个三棱锥中，我们不难发现
在
上的投影就是
在
上的投影再向
上的投影，并且由于
、
与 都是归一化的向量，因此可以得出
。
到此，我们可以得出二重积分的形式：
散射系数
好的，那现在就剩下散射系数
没有求解了。
在上文已经介绍了以距离
（上文写的是
，这里为了与理想圆球的半径
区分改成
）为自变量的扩散函数
，我们设材质的散射系数
，其中
为散射时的能量损耗。根据能量守恒定理，对一重积分而言，对于辐射通量为1的平行光，环上散射及耗散的能量总和应当也为1，即有
，化简得
。对于抽象出的半径为
的理想环形，可以得出
点与
点之间的距离
，代入即得
。
再看对球面积分的情况。对二重积分而言，能量守恒的积分式子应为
，化简得
。参考下图，由于
，可证
，
，得距离
即
，代入原积分得
。
把上面的式子整理一下，可以得到对环积分的最终形式为：
以及对球面积分的最终形式：
终于......到此为止，我们在积分上的工作就完成了。它看起来很复杂，但是对代码来说已经足够简单了。现在只要根据GPU Gems 3中给的
、高斯系数、以及RGB三色光的散射系数，写一个简单的程序就能把这张LUT图弄出来了。
写个Python出图
为了方便我们在渲染时进行采样，这张LUT在竖直方向上以
采样，水平方向上以受光强度
采样[6]。
我用Python做了一个出图的工具，不得不说Python原生速度是真的慢...后来上了多线程和numba速度才算起飞。核心部分是这么写的：
代码语言：javascript复制@jit(nopython=True)
def sp_integrate(accuracy=.1, thickness=1.0, k=1.0, theta=0.0, cost=.0):
delta_b = np.pi * accuracy
delta_a = 2 * np.pi * accuracy
b = 0
total_weights = np.array([.0, .0, .0])
total_light = np.array([.0, .0, .0])
total_cost = np.array([.0, .0, .0])
# 二重积分
while b <= np.pi:
a = -np.pi
sub_weights = np.array([.0, .0, .0])
sub_light = np.array([.0, .0, .0])
sub_cost = np.array([.0, .0, .0])
while a <= np.pi:
weight = R(b, a, thickness) * k
sub_weights = weight * delta_a
light = saturate(sample_light(theta, b, a))
sub_cost = light * delta_a
sub_light = light * weight * delta_a
a = delta_a
total_weights = sub_weights * delta_b
total_light = sub_light * delta_b
total_cost = sub_cost * delta_b
b = delta_b
return (total_light, total_weights, cost * total_cost)
@jit(nopython=True)
def ring_integrate(accuracy=.1, thickness=1.0, k=1.0, theta=0.0, cost=.0):
delta_x = 2 * np.pi * accuracy
total_weights = np.array([.0, .0, .0])
total_light = np.array([.0, .0, .0])
total_cost = np.array([.0, .0, .0])
x = -np.pi
while x <= np.pi:
weight = R(a=x, r=thickness) * k
total_weights = weight * delta_x
light = saturate(sample_light(theta, x))
total_light = light * weight * delta_x
total_cost = light * delta_x
x = delta_x
return (total_light, total_weights, cost * total_cost)
@jit(nopython=True)
def integrate(theta=0.0, thickness=1.0, accuracy=.1, use_sphere=False, k=1.0, cost=.0):
if use_sphere:
total_light, total_weights, total_cost = sp_integrate(accuracy, thickness, k, theta, cost=cost)
rgb = ((1 - 2 * np.pi * np.pi * cost) * total_light) / total_weights total_cost
rgb = np.power(rgb, 1 / 2.2)
return rgb
else:
total_light, total_weights, total_cost = ring_integrate(accuracy, thickness, k, theta, cost=cost)
rgb = ((1 - 2 * np.pi * cost) * total_light) / total_weights total_cost
rgb = np.power(rgb, 1 / 2.2)
return rgb完整的代码在：
预积分LUT工具_Python
出来的效果勉勉强强吧，觉得不合适也可以在代码里手动调RGB三色光散射的参数：
虽然用了numba之后已经快多了，不过渲球面积分还是非常……龟速Orz，这里还写了一个用GPU加速的版本，不过可能是精度的原因，导致GPU出的图看起来和CPU有些不一样。
核心部分：
代码语言：javascript复制@nb.jit(nopython=True)
def sp_integrate(idx, lw, vli, accuracy=.1, thickness=1.0, k=1.0, theta=0.0):
delta_b = np.math.pi * accuracy
delta_a = 2 * np.math.pi * accuracy
b = 0.0
total_weights = 0
total_light = 0
# 对球面积分
# 二重积分
while b <= np.math.pi:
a = -np.math.pi
while a <= np.math.pi:
d = abs(thickness * np.math.sqrt(2 - 2 * np.math.cos(a) * np.math.cos(b)))
sub_weights = 0.0
sub_light = 0.0
weight = R_origin(idx, lw, vli, d) * k
sub_weights = weight * delta_a
light = saturate(sample_light(theta, b, a)) * weight
sub_light = light * delta_a
a = delta_a
total_weights = sub_weights * delta_b
total_light = sub_light * delta_b
b = delta_b
return total_light / total_weights
@nb.jit(nopython=True)
def ring_integrate(idx, lw, vli, accuracy=.1, thickness=1.0, k=1.0, theta=0.0):
delta_a = 2 * np.math.pi * accuracy
total_weights = 0.0
total_light = 0.0
a = -np.math.pi
while a <= np.math.pi:
d = abs(thickness * np.math.sqrt(2 - 2 * np.math.cos(a)))
weight = R_origin(idx, lw, vli, d) * k
total_weights = weight * delta_a
light = saturate(sample_light(theta, a, 0)) * weight
total_light = light * delta_a
a = delta_a
return total_light / total_weights
# @nb.jit()
@nb.jit(nopython=True)
def integrate(idx, lw, vli, theta=0.0, thickness=1.0, accuracy=.1, use_sphere=False, k=1.0):
result = 0
if use_sphere:
result = sp_integrate(idx, lw, vli, accuracy, thickness, k, theta)
else:
result = ring_integrate(idx, lw, vli, accuracy, thickness, k, theta)
rgb = np.math.pow(result, 1 / 2.2)
return rgb
@nb.jit(nopython=True)
def integrate_cost(theta=0.0, accuracy=.1, use_sphere=False,cost = 0.0):
if use_sphere:
delta_b = np.math.pi * accuracy
delta_a = 2 * np.math.pi * accuracy
b = 0.0
total_cost = 0.00000
# 对球面积分
# 二重积分
while b <= np.math.pi:
a = -np.math.pi
sub_cost = 0.0000
while a <= np.math.pi:
light = saturate(sample_light(theta, b, a))
sub_cost = light * delta_a
a = delta_a
total_cost = sub_cost * delta_b
b = delta_b
return cost * total_cost
else:
delta_a = 2 * np.math.pi * accuracy
total_cost = 0.0
a = -np.math.pi
while a <= np.math.pi:
light = saturate(sample_light(theta, a, 0))
total_cost = light * delta_a
a = delta_a
return cost * total_cost
@cuda.jit
def generate_lut(result, accuracy, use_sphere, k, cost, max_r):
row, col = cuda.grid(2)
width = result.shape[0]
height = result.shape[1]
if row < width and col < height:
uv_x =((col / height) - .5) * 2.0
uv_y = 1.0 - (row / width)
theta = np.math.acos(uv_x)
rr = 1.0 / (uv_y * max_r .0000001)
c = ((1.0 - 2.0 * np.math.pi * np.math.pi * cost* accuracy) if use_sphere else (1.0 - 2.0 * np.math.pi * cost* accuracy))
r = integrate(0, light_weights_tuples, vl, theta, rr, accuracy, use_sphere, k) * c
g = integrate(1, light_weights_tuples, vl, theta, rr, accuracy, use_sphere, k) * c
b = integrate(2, light_weights_tuples, vl, theta, rr, accuracy, use_sphere, k) * c
co = integrate_cost(theta,accuracy,use_sphere,cost)
r = saturate(r co)
g = saturate(g co)
b = saturate(b co)
result[row, col][0] = r * 255
result[row, col][1] = g * 255
result[row, col][2] = b * 255完整代码：
GPU加速版本
青色会更明显一点，虽然CPU版本也会带有青色，而且函数曲线画出来在
接近0的地方就是绿色蓝色光的散射明显大于红色，但感觉都没有GPU里算的那么明显Orz。如果不想慢慢调教颜色参数的话就在PS里调吧（）
渲染
好了，现在可以把图放进渲染里用了。
首先是对LUT如何采样的问题。UV.X方向上直接用
采应该没什么问题，UV.Y方向上需要处理一下。因为在生成LUT图时，竖直方向上我们是以理想小球半径的倒数，即用
来排布的，在采样时，我们就应当反求出每个像素点对应的小球半径。
这里用一个相似三角形就可以解出来了：
。这玩意其实就是曲率，详细的可以看Jeffrey Zhuang大佬的证明[9]。Shader代码就很简单了，用上面推出来的算式直接采样：
代码语言：javascript复制 float3 GetFakeSSS(float3 normal,float3 lightDir,float3 worldPos,half intensity,half thickness)
{
half NdotL = dot(normal,lightDir) ;
half t = (1.0-thickness) * saturate(dot(-normal,lightDir));
half cuv = length(fwidth(normal)) / length(fwidth(worldPos));
float3 dif = SAMPLE_TEXTURE2D(_SSSMap, sampler_SSSMap, float2(NdotL t, intensity * cuv));
return dif * (1 t);
}效果对比：
