从零开始的预积分次表面散射

HkingAuditore

发布于 2023-10-26 17:35:46

3170

发布于 2023-10-26 17:35:46

没什么技术力，可能会有错，希望各位大佬批评指正Orz

次表面散射（Subsurface scattering，SSS），是光在传播时的一种现象。光一般会穿透物体的表面，并在与材料之间发生交互作用而被散射开来，在物体内部在不同的角度被反射若干次，最终穿出物体[1]。对于大理石，皮肤，树叶，蜡和牛奶等材料，次表面散射对于提升材质的质感而言非常重要。

皮肤的三层渲染模型

因为散射这种效果太受材质内部性质影响了，Donner（2005）的论文《Light Diffusion in Multi-Layered Translucent Materials》[2]就证明了单层模型不足以表现出良好的渲染效果，因此提出了一个三层的皮肤渲染模型。对于皮肤而言，光路表现出如图的散射路线。光射入皮肤后会经过一系列内部反射，最后从另外一处射出来，而光射出时的颜色取决于其散射光路，对于人类皮肤而言，散射出的颜色主要是红色[3]。

看完皮肤结构再来看光。假设散射效果与光的入射角度和其与入射点距离存在某种函数关系，但对于一个质地均匀的材质，光在它上面所有方向上的散射都相同的，而且由于在物体内部散射时，入射光的方向性几乎丢失，所以散射效果与光照角度无关；基于这些推论，光的散射效果只与入射光强度和距离有关。

现在我们对着这个质地均匀物体射出一束光，左图中光源外围红色的一圈就是散射出的效果[4]。且不同颜色光具有不同的散射轮廓，如右图。从图上看，皮肤对红光的散射比绿光和蓝光都更强烈。在这里，我们用r表示观测点与光的入射点之间的距离。

左边这张图其实我们也能做。它是怎么来的呢，GPU Gems 3给出来的思路是，用多个高斯函数的和去拟合扩散曲线，即对扩散曲线

，有

R\left(r\right)\approx\sum_{i=1}^{k}{w_iG\left(v_i,r\right)}

。（曲线来源是Jensen的论文《A Practical Model for Subsurface Light Transport》[5]，暂时还没复现，先挖个坑）。这个找近似高斯函数和的方法就是用matlab之类的查找，论文里已经给出了一个

，只要找到

k

个高斯函数使得

\int_{0}^{\infty}{r\left[R\left(r\right)-\sum_{i=1}^{k}{\omega_iG\left(v_i,r\right)}\right]^2dr}

这个积分取到最小值就好了。英伟达已经把皮肤的参数给找出来了[4]。

GPU Gems 3里给出了这样的拟合高斯和图:

不过我这边按照公式复现出来图是长这样的，我怀疑是因为这里要对RGB的各自的权值做归一化，把绿光和蓝光原本的散射权重给隐藏掉了，所以GPU gems 3里作的这张图应该是用未进行归一化的权值算的。

我这里顺便验证了一下GPU Gems 3中给出的各个高斯函数，分别做了几张Diffuse Profile，效果是对得上的。值得注意的是，这里用的扩散函数不是固定的，可以自己按需求选择[6]。我这里还拿了迪士尼的公式[7]做了个实验：

其实，根据上面的这几张Diffuse Profile，我们已经可以根据任意的输入

r

来查找对应的散射效果了。这张Diffuse Profile其实只存储了

r

向

L

映射的一维信息。然而，在实际做渲染时，我们面对的一般是平行光，目标皮肤模型也通常是个曲面，这就意味着我们一般面临的不是“观测点与入射点”的问题。

根据之前的文献和推理，既然光在射入材质时我们已经忽略了它的方向关系，那我们不妨建立一个“辐照度-距离-散射”的对应关系，在渲染中使用辐照度和距离查找某张LUT，得到散射效果。

在这个对应关系里，辐照度是很好求的，距离该如何求呢？对于平行光，我们已经无法像找激光入射点一样解决这个问题了，那不妨设整个材质是通透的，平行光给这个材质均匀地施加了光照，材质在任何一点都具有相同的散射性质，这样对于某个点接收的散射量，我们都可以根据这个点自身的信息求出，而不需要考虑模型整体的光照情况，这样点就满足

。这样操作虽然不够物理，但是计算量得到了减小，效果也能满足不错的近似。

数学推导

积分基本形式

首先，我们假设存在一个半径为

r

的圆，圆上所有点都具有相同的散射性质。那么对于圆上任意一点

P

，它将受到圆上所有点散射而来的光。现在，给这个圆施加平行光

L

，圆内将发生散射。

现在，设有一点

Q

,其法线为

，

与

P

的法线

N

夹角为

x

，平行光线

L

与

N

的夹角是

。设散射系数为

。根据之前之前的假设，

P

点将受到圆上所有点，包括

Q

点传递过去的散射。可以求得点

Q

受到的直接光照就是

L_Q=L_{intensity}\cdot\left(N_Q\cdot L\right)

，由于法向量、光线向量为归一化向量，可得

L_Q=L_{intensity}\cdot\cos{\left\langle N_Q,L\right\rangle}=L_{intensity}\cdot\cos{\left(\theta-x\right)}

。现在设平行光

L

的光强为1，则有

。那么根据之前的定义与推论，点

Q

将给点

P

提供

\cos{\left(\theta-x\right)}\cdot s\left(x\right)

的散射光，意思就是点将其受到的光照根据材质的散射系数散射出去。

如果不考虑二次散射，那么只有受到光照的点才能向四周发出散射，则

P

点受到的散射量就是圆上所有点的散射量之和，积分可得

\int_{\theta-\frac{\mathrm{\Pi}}{2}}^{\theta+\frac{\mathrm{\Pi}}{2}}{\cos{\left(\theta+x\right)}\cdot s\left(x\right)dx}

。不过这样的形式显然太复杂了，实际上，对于未受到光照的点，其直接光照量为0，它乘上散射系数后也是0，积分的时候其实不必把这段单独处理，只需要把

限制到

即可，即：

\int_{-\pi}^{\pi}{saturate\left(\cos{\left(\theta+x\right)}\right)\cdot s\left(x\right)dx}

。

好了，现在已经完成了对圆环散射的积分；不过毕竟是渲染三维物体，所以还需要对球面进行积分。虽然Panner在他的PPT[8]中认为在环上积分和在球面上积分效果区别不大，不过这里还是打算做到底把球面积分的路也走一遍。

如图，我们在对

PQ

环进行积分的同时对

Q

所在的垂直于

PQ

环的

环积分。先前在积

PQ

环时我们设

为

x

进行积分；现在为了方便进行二重积分，我们设

为

，设

为

。

首先，对应原积分里的

，二重积分中我们需要用

\cos\left\langle L,N_{Q^{'}} \right\rangle

表示光照量。由于这是一个积分，

L

又是平行光，因此我们可以视

与

L

构成的平面垂直于

与

构成的平面，并由此抽象出一个过同顶点

Q

的三条棱两两垂直的三棱锥，如图。

在这个三棱锥中，我们不难发现

L

在

上的投影就是

L

在

上的投影再向

上的投影，并且由于

L

、

与都是归一化的向量，因此可以得出

\cos{\left\langle L,N_{Q^\prime}\right\rangle}=\cos{\left(\theta+\alpha\right)}\cdot\cos{\beta}

。

到此，我们可以得出二重积分的形式：

\int_{0}^{\pi}{\int_{-\pi}^{\pi}{saturate\left(\cos{\left(\theta+\alpha\right)}\cdot\cos{\beta}\right)\cdot s\left(\alpha,\beta\right)d\alpha}d\beta}

散射系数

好的，那现在就剩下散射系数

没有求解了。

在上文已经介绍了以距离

d

（上文写的是

r

，这里为了与理想圆球的半径

r

区分改成

d

）为自变量的扩散函数

，我们设材质的散射系数

s\left(x\right)=kR\left(d\right)+\lambda

，其中

为散射时的能量损耗。根据能量守恒定理，对一重积分而言，对于辐射通量为1的平行光，环上散射及耗散的能量总和应当也为1，即有

\int_{-\pi}^{\pi}\left(kR\left(d\right)+\lambda\right)dx=1

，化简得

k=\frac{1-2\pi\lambda}{\int_{-\pi}^{\pi}R\left(d\right)dx}

。对于抽象出的半径为

r

的理想环形，可以得出

P

点与

Q

点之间的距离

，代入即得

k=\frac{1-2\pi\lambda}{\int_{-\pi}^{\pi}R\left(2r\sin{\frac{x}{2}}\right)dx}

。

再看对球面积分的情况。对二重积分而言，能量守恒的积分式子应为

\int_{0}^{\pi}{\int_{-\pi}^{\pi}\left(kR\left(d\right)+\lambda\right)d\alpha d\beta}=1

，化简得

k=\frac{1-2\pi^2\lambda}{\int_{0}^{\pi}{\int_{-\pi}^{\pi}R\left(d\right)d\alpha d\beta}}

。参考下图，由于

，可证

，

Q^\prime K=\sqrt{r^2+r^2\cos^2{\alpha}-2r^2\cos{\alpha\cdot\cos{\beta}}}

，得距离

d

即

Q^\prime P=r\sqrt{2-2\cos{\alpha\cdot\cos{\beta}}}

，代入原积分得

k=\frac{1-2\pi^2\lambda}{\int_{0}^{\pi}{\int_{-\pi}^{\pi}R\left(r\sqrt{2-2\cos{\alpha\cdot\cos{\beta}}}\right)d\alpha d\beta}}

。

把上面的式子整理一下，可以得到对环积分的最终形式为：

\frac{\left(1-2\pi\lambda\right)\int_{-\pi}^{\pi}{saturate\left(\cos{\left(\theta+x\right)}\right)\cdot R\left(2r\sin{\frac{x}{2}}\right)dx}}{\int_{-\pi}^{\pi}R\left(2r\sin{\frac{x}{2}}\right)dx}+\lambda\int_{-\pi}^{\pi}saturate\left(\cos{\left(\theta+x\right)}\right)dx

以及对球面积分的最终形式：

\frac{\left(1-2\pi^2\lambda\right)\int_{0}^{\pi}{\int_{-\pi}^{\pi}{saturate\left(\cos{\left(\theta+\alpha\right)}\cdot\cos{\beta}\right)\cdot R\left(r\sqrt{2-2\cos{\alpha\cdot\cos{\beta}}}\right)d\alpha}d\beta}}{\int_{0}^{\pi}{\int_{-\pi}^{\pi}R\left(r\sqrt{2-2\cos{\alpha\cdot\cos{\beta}}}\right)d\alpha d\beta}}+\lambda\int_{0}^{\pi}{\int_{-\pi}^{\pi}saturate\left(\cos{\left(\theta+\alpha\right)}\cdot\cos{\beta}\right)d\alpha d\beta}

终于......到此为止，我们在积分上的工作就完成了。它看起来很复杂，但是对代码来说已经足够简单了。现在只要根据GPU Gems 3中给的

R\left(r\right)=\sum_{i=1}^{k}{w_iG\left(v_i,r\right)}

、高斯系数、以及RGB三色光的散射系数，写一个简单的程序就能把这张LUT图弄出来了。

写个Python出图

为了方便我们在渲染时进行采样，这张LUT在竖直方向上以

采样，水平方向上以受光强度

采样[6]。

我用Python做了一个出图的工具，不得不说Python原生速度是真的慢...后来上了多线程和numba速度才算起飞。核心部分是这么写的：

@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有些不一样。

核心部分：

@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版本也会带有青色，而且函数曲线画出来在

r

接近0的地方就是绿色蓝色光的散射明显大于红色，但感觉都没有GPU里算的那么明显Orz。如果不想慢慢调教颜色参数的话就在PS里调吧（）

渲染

好了，现在可以把图放进渲染里用了。

首先是对LUT如何采样的问题。UV.X方向上直接用

采应该没什么问题，UV.Y方向上需要处理一下。因为在生成LUT图时，竖直方向上我们是以理想小球半径的倒数，即用

来排布的，在采样时，我们就应当反求出每个像素点对应的小球半径。

这里用一个相似三角形就可以解出来了：

\frac{\left|\vec{N}\right|}{r}=\frac{\left|\Delta\vec{N}\right|}{\left|\Delta\vec{P}\right|}

。这玩意其实就是曲率，详细的可以看Jeffrey Zhuang大佬的证明[9]。Shader代码就很简单了，用上面推出来的算式直接采样：

    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);
    }

效果对比：