19 min to read
多重散射PBR
介绍
光的散射包含了反射和折射，我们在实时渲染引擎中所计算的PBR通常只包含了一次散射，也就是一次反射。而在真实世界中物体表面通常有多次散射，常规的一次散射会出现能量不守恒的现象，一般来说物体越粗糙散射的次数越多，但如果只有一次散射，会导致最后的渲染结果能量丢失，变得粗糙度越高物体越黑。
在我的渲染引擎中参考1实现了一个多重散射的PBR（仅在IBL）来达到能量守恒。
背景
反射方程：
下面是一个我们通常定义的某点p的出射辐率：\(L_o(p, \omega_o) = L_e(p, \omega_o) + \int_{\Omega} f_r(p, \omega_i, \omega_o) L_i(p, \omega_i) (\omega_i \cdot n) d\omega_i\) 其中：
- $L_o(p, \omega_o)$ 是从表面点 $p$ 在方向 $\omega_o$ 离开的辐射亮度。
- $L_e(p, \omega_o)$ 是点 $p$ 在方向 $\omega_o$ 的自发光亮度。
- $\Omega$代表一个和
n对齐的半球。 - $\int_{\Omega}$ 表示对半球 $\Omega$ 上所有方向的积分，积分代表这个点要接受从整个半球所有射出光，可以用黎曼和(Riemann sum)的方式来离散求解。
- $f_r(p, \omega_i, \omega_o)$ 是点 $p$ 处的BRD），描述了入射光 $\omega_i$ 如何被反射到方向 $\omega_o$。
- $L_i(p, \omega_i)$ 是入射辐射亮度，即从方向 $\omega_i$ 到点 $p$ 的辐射亮度。
- $(\omega_i \cdot n)$ 是入射光方向（在半球中是一个立体角）和表面法线 $n$ 的点积，代表了光线入射角度的余弦值，用于计算光线与表面的相互作用。
- $d\omega_i$ 是入射方向的微小立体角元素。
h是一个介于l（光照方向）和v（视角方向）之间的半程向量。BRDF:
BRDF定义了入射光线与反射光线之间的关系，也就是物体表面最后发出来的光是什么样的。有非常多BRDF，如下举例几种常见的：
-
Lambertian:\(f_r(\omega_i, \omega_o) = \frac{\rho}{\pi}\)其中，$\rho$ 是反射率，也就是光照颜色。
-
Phong:\(f_r(\omega_i, \omega_o) = k_s \cdot (\max(0, n \cdot h))^{\alpha}\)其中，$k_s$ 是镜面反射系数，$\alpha$ 是高光锐度参数。
- Cook-Torrance:\(f_r(\omega_i, \omega_o) = \frac{D \cdot F \cdot G}{4(\omega_i \cdot n)(\omega_o \cdot n)}\)这个部分只看镜面反射，这里的(D)、(F)和(G)分别表示微面分布、菲涅尔项和几何遮蔽项。
IBL:
通常使用立方体贴图，将整个环境看作一个大光源。也使用上述的反射方程。
直接光
单次散射：
对于直接光，我们使用的BRDF包含漫反射和镜面反射两个部分，Lambertian BRDF作为漫反射部分，Cook-Torrance作为镜面反射部分，主要参考了UE4中的做法2。我们的BRDF就变成了：\(f_r = k_d f_{\text{lambert}} + k_s f_{\text{cook-torrance}}\)
其中$k_d$和$k_s$分别代表直接光漫反射和镜面反射的占比。我们的直接光包含定向光，点光源和聚光灯（都是从一个点发出的光），因为这些光源都带有固定的方向，所以我们在求反射方程的时候就不需要计算积分(只有一个方向有光，其他方向为0)，只需把每个光源的单方向反射方程结果加在一起就可以了。将l光照方向作为如何方向，v视角方向作为出射方向，于是我们的反射方程就变成了（暂时忽略自发光）：\(L_o(p, v) = f_r * L_i(p, l) * (l \cdot n)\)
Cook-Torrance BRDF的镜面反射部分包含三个函数，此外分母部分还有一个标准化因子 。字母D，F与G分别代表着一种类型的函数，各个函数分别用来近似的计算出表面反射特性的一个特定部分。三个函数分别为法线分布函数(Normal Distribution Function)，菲涅尔方程(Fresnel Rquation)和几何函数(Geometry Function)：
- D 法线分布函数：估算在受到表面粗糙度的影响下，朝向方向与半程向量一致的微平面的数量。这是用来估算微平面的主要函数。
- G 几何遮挡函数：描述了微平面自成阴影的属性。当一个平面相对比较粗糙的时候，平面表面上的微平面有可能挡住其他的微平面从而减少表面所反射的光线。
- F 菲涅尔方程：菲涅尔方程描述的是在不同的表面角下表面所反射的光线所占的比率。
D,F,G同样有不同的选择3，这里我们使用和UE4相同的DFG: D使用Trowbridge-Reitz GGX，F使用Fresnel-Schlick近似(Fresnel-Schlick Approximation)，而G使用Smith’s Schlick-GGX。
D 法线分布函数：
D近似地表示了与h取向一致的微平面比率（0-1）。我们使用UE4中的Trowbridge-Reitz GGX：\(NDF_{GGXTR}(n, h, \alpha) = \frac{\alpha^2}{\pi \left( (n \cdot h)^2 (\alpha^2 - 1) + 1 \right)^2}\)
其中$\alpha$=粗糙度，当粗糙度很低的时候反射就会非常集中，粗糙度高的的时候分社就会非常分散，颜色会更加灰暗。
G 几何函数：
G近似地求得了微表面间相互遮蔽的比率 （0-1），遮挡会损耗光线的能量。这里使用UE4中的Schlick-GGX： \(G_1(v) = \frac{n \cdot v}{(n \cdot v)(1 - k) + k}\)k 是粗糙度的映射，这个映射只在直接光中使用。 \(k = \frac{(Roughness + 1)^2}{8}\)最后要把光照方向和几何方向的遮挡同时考虑进去：\(G(l, v, h) = G_1(l) \cdot G_1(v)\)
F 菲涅尔函数：
菲涅尔函数描述的是被反射的光线对比光线被折射的部分所占的比率，这个比率会随着我们观察的角度不同而不同。当与法线的夹角越接近90°时反光越多。当垂直法线观察的时候，任何物体或者材质表面都有一个基础反射率，非导电体通常使用0.04，金属使用他的albedo。这里使用Schlick近似模型，再加上UE4中的Spherical Gaussian approximation2。\(F(v, h) = F_0 + (1 - F_0) \cdot 2^{(-5.55473(v \cdot h) - 6.98316)(v \cdot h)}\) $F_0$就是基础反射率。我们像如下代码，用输入的金属度来决定F0：
vec3 F0 = vec3(0.04);
F0 = mix(F0, albedo, metalness);
在最后加上漫反射和镜面反射的结果时，我们要记住金属只有镜面反射，F已经代表了反射的比率（最后不需要再乘$k_s$）
vec3 kS = F;
vec3 kD = vec3(1.0) - kS;
kD *= 1.0 - metalness;
最后我们计算了$f_r$, $L_i(p, l)$ 是光照颜色衰减, l与n的点乘是夹角衰减。所有项相乘我们就得到了最后的结果。\(L_o(p, v) = f_r * L_i(p, l) * (l \cdot n)\)
间接光
单次散射：
间接光也分为漫反射和镜面反射，先看漫反射。间接光漫反射也有很多方法，包括球协光照函数和IBL。这里我们看IBL。和上面的单点光源不同的是，这里对于一个点p我们需要考虑整个半球方向的光源，所以需要解决积分。我们可以对环境贴图进行卷积，每个方向的环境贴图卷积就是场景的辐照度，然后再采样卷积后的环境贴图（辐照度贴图），相当于提前对环境贴图进行了积分。我们参考learn opengl中提到的立方体卷积方式用compute shader对立方体贴图进行卷积处理。
对于立方体贴图的每个纹素，在纹素所代表的方向的半球 Ω 内生成固定数量的采样向量，并对采样结果取平均值。数量固定的采样向量将均匀地分布在半球内部。注意，积分是连续函数，在采样向量数量固定的情况下离散地采样只是一种近似计算方法，我们采样的向量越多，就越接近正确的结果。 反射方程的积分 ∫ 是围绕立体角 dw 旋转，而这个立体角相当难以处理。为了避免对难处理的立体角求积分，我们使用球坐标 θ 和 ϕ 来代替立体角。
于是我们的漫反射方程就变成了这样，对于球坐标的两个轴进行求积分。
\(L_o(p,\phi_o,\theta_o) = \frac{k_dc}{\pi} \int_{0}^{2\pi} \int_{0}^{\frac{\pi}{2}} L_i(p,\phi_i,\theta_i) \cos(\theta) \sin(\theta) d\phi d\theta\)
$\sin(\theta)$用于权衡（环境贴图）区域的贡献度。用黎曼和：\(L_o(p,\phi_o,\theta_o) = {k_d} \frac{c{\pi}}{n_1 n_2} \sum_{\phi=0}^{n_1} \sum_{\theta=0}^{n_2} L_i(p,\phi_i,\theta_i) \cos(\theta) \sin(\theta) d\phi d\theta\)我们通过 cos(theta) 缩放采样颜色值，因为角度较大时光线较弱，并通过 sin(theta) 缩放采样颜色值，以考虑较高半球区域中较小的样本区域,确保了在进行角度积分时，正确地考虑了不同纬度上的表面区域大小(越靠近极点面积越小)。以下是compute shader代码，要注意设置opengl cubmap时的格式。
#version 450 core
layout(local_size_x = 16, local_size_y = 16) in;
layout(binding = 0) uniform samplerCube sourceCubemap;
layout(rgba8, binding = 1) writeonly uniform imageCube destinationCubemap;
#define PI 3.14159265359
#define TWO_PI 6.2831853071795864769252867665590
#define HALF_PI 1.5707963267948966192313216916398
void main()
{
// 决定当前shader处理哪一个面
ivec3 globalId = ivec3(gl_GlobalInvocationID.xyz);
int face = globalId.z;
// 得到当前面UV
vec2 uv = (vec2(globalId.xy) + vec2(0.5)) / vec2(imageSize(destinationCubemap).xy);
uv = uv * 2.0 - 1.0; // Transform uv to [-1, 1]
// 根据面得到法线
vec3 normal;
if (face == 0) normal = vec3(1, -uv.y, -uv.x); // Positive X
else if (face == 1) normal = vec3(-1, -uv.y, uv.x); // Negative X
else if (face == 2) normal = vec3(uv.x, 1, uv.y); // Positive Y
else if (face == 3) normal = vec3(uv.x, -1, -uv.y); // Negative Y
else if (face == 4) normal = vec3(uv.x, -uv.y, 1); // Positive Z
else if (face == 5) normal = vec3(-uv.x, -uv.y, -1); // Negative Z
// 计算tangent 和 bitangent
vec3 up = abs(normal.z) < 0.999 ? vec3(0.0, 0.0, 1.0) : vec3(1.0, 0.0, 0.0);
vec3 right = normalize(cross(up, normal));
up = normalize(cross(normal, right));
vec3 irradiance = vec3(0.0);
float sampleDelta = 0.025;
float nrSamples = 0.0;
for(float phi = 0.0; phi < 2.0 * PI; phi += sampleDelta)
{
for(float theta = 0.0; theta < 0.5 * PI; theta += sampleDelta)
{
// 球面坐标到cartesian (tangent space)
vec3 tangentSample = vec3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
// tangent space to world space
vec3 sampleVec = tangentSample.x * right + tangentSample.y * up + tangentSample.z * normal;
// 采样
irradiance += texture(sourceCubemap, sampleVec).rgb * cos(theta) * sin(theta);
nrSamples++;
}
}
irradiance = PI * irradiance * (1.0 / float(nrSamples));
imageStore(destinationCubemap, globalId, vec4(irradiance, 1));
}
globalID.z代表目前在cubmap的哪一个面，每一个的normal都朝外，垂直与面。
镜面反射和间接光漫反射一样使用IBL，同样需要计算整个半球的积分，但因为镜面反射和视角方向相关，采样时需要用到光照方向和视角方向，我们没法简单地采样一张贴图。同时我们也需要快速求得半球积分。
Split sum:
这里我们使用重要性采样和UE4的split sum2方法。split sum是把反射方程拆分成了两个部分：\(L_o(p,\omega_o) = \int_{\Omega} Li(p,\omega_i) d\omega_i * \int_{\Omega} f_r(p,\omega_i,\omega_o) (n \cdot \omega_i) d\omega_i\)
重要性采样：
按照与被积函数相关的某种分布来抽取样本，而不是简单地在整个积分区间内均匀抽样。通过这样做，可以让采样更加集中于对积分贡献较大的区域，从而提高估计的效率和准确度，叫做重要性采样。
对于间接光漫反射，我们通过生成均匀分布在半球 Ω 上的样本向量来对环境图使用球坐标进行卷积 ，这种方法对于镜面反射来说效率较低。对于镜面反射，根据表面的粗糙度，光线会紧紧或大致围绕反射r法线 n:
由于大多数光线最终会在h周围反射的镜面波瓣中（微表面模型理论），所以我们不需要采样这之外的样本。
蒙特卡洛积分：
蒙特卡洛积分是一种利用随机抽样来近似计算积分的方法。于是我们可以把我们的反射方程变成如下：\(\int_{\Omega} f_r(p,\omega_i,\omega_o) (n \cdot \omega_i) d\omega_i\approx \frac{1}{N} \sum_{k=1}^{N} \frac{f_{\text{brdf}}(p,\omega_i,\omega_o) L_i(p,\omega_i) (n \cdot \omega_i)}{\text{pdf}(\omega_i,\omega_o)}\)pdf代表概率密度函数，它告诉我们特定样本在整个样本集中出现的概率。我们可以通过选择pdf来进行重要性采样。split sum也就变成了：\(\int_{\Omega} f_r(p,\omega_i,\omega_o) (n \cdot \omega_i) d\omega_i \approx \frac{1}{N} \sum_{k=1}^{N} {L_i(p,\omega_i)} + \frac{1}{N} \sum_{k=1}^{N} \frac{f_{\text{brdf}}(p,\omega_i,\omega_o) L_i(p,\omega_i) (n \cdot \omega_i)}{\text{pdf}(\omega_i,\omega_o)}\)
预过滤环境贴图：
split sum中的左边部分也叫做预过滤环境贴图（环境光辐照率），预先计算的环境卷积图，但这次考虑了粗糙度，粗糙度越高环境贴图越模糊。由于我们事先不知道卷积环境贴图时的视图方向，因此 Epic Games 采用了进一步的近似方法，即假设视图方向（以及镜面反射方向）等于输出样本方向 ωo。
vec3 N = normalize(w_o); vec3 R = N; vec3 V = R;
相当于在卷积时把粗糙度看作0，这个假设可以带来合理的结果。
我们对$\int_{\Omega} Li(p,\omega_i) d\omega_i$进行蒙特卡洛积分和重要性采样。为了让蒙特卡洛积分采样更加平均我们使用Hammersley Sequence 4，同时我们将根据表面的粗糙度生成偏向微表面h的一般反射方向的样本矢量(重要性采样)。所以我们的整个过程就是：取Hammersley Sequence值(一个2D的点)，映射到球面坐标，然后在切线空间中生成样本向量h，变换到世界空间，并对环境贴图的辐射率进行采样。代码如下：
vec3 ImportanceSampleGGX(vec2 Xi, vec3 N, float roughness)
{
float a = roughness*roughness;
// 计算采样点的球面坐标
float phi = 2.0 * PI * Xi.x;
float cosTheta = sqrt((1.0 - Xi.y) / (1.0 + (a*a - 1.0) * Xi.y));
float sinTheta = sqrt(1.0 - cosTheta*cosTheta);
// 球面坐标到cartesian
vec3 H;
H.x = cos(phi) * sinTheta;
H.y = sin(phi) * sinTheta;
H.z = cosTheta;
// tangent space
vec3 up = abs(N.z) < 0.999 ? vec3(0.0, 0.0, 1.0) : vec3(1.0, 0.0, 0.0);
vec3 tangent = normalize(cross(up, N));
vec3 bitangent = cross(N, tangent);
// to world-space
vec3 sampleVec = tangent * H.x + bitangent * H.y + N * H.z;
return normalize(sampleVec);
}
// 收敛速度更快
vec2 Hammersley(uint i, uint N)
{
uint bits = (i << 16u) | (i >> 16u);
bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
float rdi = float(bits) * 2.3283064365386963e-10;
return vec2(float(i) /float(N), rdi);
}
void main()
{
// 决定当前shader处理哪一个面
ivec3 globalId = ivec3(gl_GlobalInvocationID.xyz);
int face = globalId.z;
// 映射到对应面的UV范围
vec2 uv = (vec2(globalId.xy) + vec2(0.5)) / vec2(imageSize(prefilteredMap).xy);
uv = uv * 2.0 - 1.0; // Transform uv to [-1, 1]
// 法线垂直于当前面
vec3 N;
if (face == 0) N = vec3(1, -uv.y, -uv.x); // Positive X
else if (face == 1) N = vec3(-1, -uv.y, uv.x); // Negative X
else if (face == 2) N = vec3(uv.x, 1, uv.y); // Positive Y
else if (face == 3) N = vec3(uv.x, -1, -uv.y); // Negative Y
else if (face == 4) N = vec3(uv.x, -uv.y, 1); // Positive Z
else if (face == 5) N = vec3(-uv.x, -uv.y, -1); // Negative Z
vec3 R = N;
vec3 V = R;
float totalWeight = 0.0;
vec3 prefilteredColor = vec3(0.0);
uint numSamples = 1024; // Number of samples for the prefiltering
for (uint i = 0; i < numSamples; ++i)
{
// 采样一个点
vec2 Xi = Hammersley(i, numSamples);
// 根据粗糙度重要性采样方向
vec3 H = ImportanceSampleGGX(Xi, N, roughness);
vec3 L = 2.0 * dot(V, H) * H - V;
float NdotL = max(dot(N, L), 0.0);
if (NdotL > 0.0)
{
prefilteredColor += texture(environmentMap, L).rgb * NdotL;
totalWeight += NdotL;
}
}
prefilteredColor /= totalWeight;
imageStore(prefilteredMap, ivec3(gl_GlobalInvocationID.xyz), vec4(prefilteredColor, 1.0));
}
注意因为我使用了compute shader，所以需要手动决定法线的朝向，根据立方体贴图不同面来决定法线方向。当我们得到样本h后就可以通过计算反射向量来得到光照方向l，也就是我们采样环境贴图的方向。我们要针对粗糙度生成不同模糊程度的预过滤贴图，使用预过滤贴图时也要用粗糙度来选择不同模糊程度：
vec3 prefilteredColor = textureLod(prefilterMap, reflect(-V, N), roughness * 4.0).rgb * ks;
预计算BRDF：
split sum右部分，也就是我们也可以用同样的蒙特卡洛积分和和重要性采样的方法来进行预计算。预计算需要用到BRDF的输入视角与法线的点乘（G），粗糙度和F0。根据learnopengl中提到的方法 5，我们可以把F0移出积分，右半部分就变成了这样（没有写成蒙塔卡罗积分形式）：\(F_0 \int_{\Omega} f_r(p, \omega_i, \omega_o) \left(1 - (1 - \omega_o \cdot h)^5\right) (n \cdot \omega_i) d\omega_i + \int_{\Omega} f_r(p, \omega_i, \omega_o) \left(1 - \omega_o \cdot h\right)^5 (n \cdot \omega_i) d\omega_i\)我们这里的F用Schlick’s Fresnel，G使用Schlick-GGX。
// Karis 2014
vec2 IntegrateBRDF(float roughness, float NdotV)
{
vec3 V;
V.x = sqrt(1.0 - NdotV * NdotV); // sin
V.y = 0.0;
V.z = NdotV; // cos
// N points straight upwards
const vec3 N = vec3(0.0, 0.0, 1.0);
// 左和右
float A = 0.0;
float B = 0.0;
const uint numSamples = 1024u;
for (uint i = 0u; i < numSamples; i++) {
vec2 Xi = Hammersley(i, numSamples);
vec3 H = ImportanceSampleGGX(Xi, N, roughness);
vec3 L = 2.0 * dot(V, H) * H - V; // 反射
float NdotL = max(dot(N, L), 0);
float NdotH = max(dot(N, H), 0);
float VdotH = max(dot(V, H), 0);
if (NdotL > 0.0)
{
float G = GeometrySmith(N, V, L, roughness);
float G_Vis = (G * VdotH) / (NdotH * NdotV);
float Fc = pow(1.0 - VdotH, 5.0);
A += (1.0 - Fc) * G_Vis;
B += Fc * G_Vis;
}
}
// return 4.0 * vec2(A, B) / float(numSamples);
return vec2(A, B) / float(numSamples);
}
void main()
{
// Normalized pixel coordinates (from 0 to 1)
vec2 uv = vec2(gl_GlobalInvocationID.xy + 1) / vec2(imageSize(brdfLUT).xy);
vec2 integratedBRDF = IntegrateBRDF(uv.x, uv.y);
imageStore(brdfLUT, ivec2(gl_GlobalInvocationID.xy), vec4(integratedBRDF, 0.0, 0.0));
}
BRDF的预计算结果是一个2D贴图，x轴是粗糙度，y轴是N点乘V。IntergrateBRDF中计算的V相当于一个平面上的V，利用点乘结果（cos）求得在2D单位圆的坐标。N一直朝上，可以视为我们在一个平面上积分。A,B分别是上述两个部分，代码就是在复现公式。
最后我们使用这个预计算贴图时，将公式中的F0视为考虑粗糙度的F。
vec2 envBRDF = texture(brdfLUT, vec2(max(dot(N, V), 0.0), roughness)).rg;
vec3 specular = (F * envBRDF.x + envBRDF.y);
因为这里已经乘了F，最后加上间接光镜面反射时就不用乘Ks了。需要注意的是在IBL中计算F时我们要考虑粗糙度,防止出现粗糙度越大菲涅尔现象反而更明显。把上述的部分都结合起来就是我们的单次散射PBR。
多次散射：
在真实世界中物体表面通常有多次散射，常规的一次散射会出现能量不守恒的现象，一般来说物体越粗糙散射的次数越多，但如果只有一次散射，会导致最后的渲染结果能量丢失，变得粗糙度越高物体越黑。
金属材质
首先，让我们考虑一个完美反射体，没有漫反射叶片，其中 F=1。在这种情况下，反射的总能量不论反弹次数，都必须等于入射能量（由能量守恒定律）： \(1 = E_{ss} + E_{ms} \Rightarrow E_{ms} = 1 - E_{ss}\) 其中 $E_{ss}$ 是我们的定向反照率（directional albedo），但是当 F=1 时： \(E_{ss} = \int_{\Omega} \frac{D(h)G(v, l, h)}{4\langle n \cdot l \rangle \langle n \cdot v \rangle} dl\) 为了定义 $E_{ms}$，我们可以以相同的方式表达它，但是使用一个未知的BRDF： \(E_{ms} = \int_{\Omega} f_{ms} \langle n \cdot l \rangle dl = 1 - E_{ss}\) 所以我们可以有效地将这个额外的叶片加到我们的反射方程中： \(L_o(v) = \int_{\Omega} (f_{ss} + f_{ms})L_i(l) \langle n \cdot l \rangle dl = \int_{\Omega} f_{ss} L_i(l) \langle n \cdot l \rangle dl + \int_{\Omega} f_{ms} L_i(l) \langle n \cdot l \rangle dl\) 我们已经讨论了左边的积分，所以让我们专注于第二个积分。我们将再次使用Karis引入的分裂求和方法，但这里Fdez-Agüera假设我们可以认为次级散射事件主要是漫反射，因此使用辐照度作为 $L_i$ 积分的解决方案： \(\int_{\Omega} f_{ms} L_i(l) \langle n \cdot l \rangle dl = (1 - E_{ss}) \int_{\Omega} \frac{L_i(l)} \pi \langle n \cdot l \rangle dl\) Fdez-Agüera进一步指出，这个近似对于狭窄的光源，如分析点光源或定向光源是失败的。有趣的是，Stephen McAuley聪明地观察到这个多次散射项只会在较高的粗糙度上占主导地位，其中我们的辐射图将会变得高度模糊和相当漫反射。他的比较显示了很小的差异，所以如果你之前没有使用辐照度，你可以在这里跳过它。然而，我们已经创建了我们的辐照度图，所以我们将在这里再次使用它。 如果我们使用之前对单次散射事件的近似，并引入这个多次散射事件的近似，我们最终得到以下结果： \(L_o(v) = (f_a + f_b) \text{radiance} + (1 - E_{ss}) \text{irradiance}\) 其中 $f_a, f_b$ 是我们在LUT中存储的Karis的缩放和偏置项。 然而，到目前为止，我们限制自己于一个完美反射的金属。为了将其扩展到通用金属，我们需要重新引入 F。与其他人一样，Fdez-Agüera将 F 分为两项，$F_{ss}$ 和 $F_{ms}$，使得： \(E = F_{ss} E_{ss} + F_{ms} E_{ms}\) 然而，与之前不同，我们不能简单地设置 E=1。相反，Fdez-Agüera 将 $F_{ms}$ 建模为一个几何级数： \(E = F_{ss} E_{ss} + \sum_{k=1}^{\infty} F_{ss} E_{ss} (1 - E_{avg})^k F_{avg}^k\) 其中 $F_{avg}$ 和 $E_{avg}$ 定义为： \(F_{avg} = F_0 + \frac{1}{21} (1 - F_0), E_{avg} = E_{ss}, F_{avg} = F_0 + \frac{1}{21} (1 - F_0),\)\(L_o(v) = (F_0 f_a + f_b) \text{radiance} + \frac{1 - F_{avg}}{(1 - E_{avg}) (1 - E_{ss})} (F_0 f_a + f_b) F_{avg} \text{irradiance}\)
电解质材质
使用单次散射BRDF在的金属和电介质球体 在较低的粗糙度下，菲涅尔效应太强，而且在较高的粗糙度下仍然有一些变暗。 \(E = 1 = F_{ss} E_{ss} + F_{ms} E_{ms} + E_{diffuse}\) 要比上面的公式多加一个diffuse。
// 根据粗糙度调整菲涅尔效应
vec3 ks = fresnelSchlickRoughness(max(dot(N, V), 0.0), F0, roughness);
vec3 R = (invSkyboxModelMatrix * vec4(reflect(-V, N), 0.0)).xyz;
float mipLevel = floor(roughness * 4.0);
vec3 prefilteredColor = textureLod(prefilterMap, R, 0).rgb;
vec2 envBRDF = texture(brdfLUT, vec2(max(dot(N, V), 0.0), roughness)).rg;
// 单次散射项
vec3 FssEss = ks * envBRDF.x + envBRDF.y;
// 多次散射项
float Ems = (1.0 - (envBRDF.x + envBRDF.y));
vec3 Favg = F0 + (1.0 - F0) / 21.0;
vec3 FmsEms = Ems * FssEss * Favg / (1.0 - Favg * Ems);
// 多加一个diffuse
vec3 kd = albedo * (1.0 - FssEss - FmsEms) * (1.0 - metalness);
vec3 indirectLight = (FssEss * prefilteredColor + (FmsEms + kd) * irradiance) * ao;
结果：
单次散射：
多次散射：
我们可以清楚的从一旁的粗糙金属球看到确实是在高粗糙度下没那么黑了。
参考
-
A Multiple-Scattering Microfacet Model for Real-Time Image-based Lighting: https://www.jcgt.org/published/0008/01/03/paper.pdf ↩
-
Real Shading in Unreal Engine 4:https://cdn2.unrealengine.com/Resources/files/2013SiggraphPresentationsNotes-26915738.pdf ↩ ↩2 ↩3
-
Graphic Rants: Specular BRDF Reference: https://graphicrants.blogspot.com/2013/08/specular-brdf-reference.html ↩
-
Hammersley Points on the Hemisphere http://holger.dammertz.org/stuff/notes_HammersleyOnHemisphere.html ↩
-
Pre-computing the BRDF https://learnopengl.com/PBR/IBL/Specular-IBL ↩