好久不见！最近在忙别的工作，这几天难得有时间，赶紧把 202 系列作业做个收尾ε٩(๑> ₃ <)۶з！

0x00 To begin with

这篇文章将会包含以下内容：

  • 部分课程内容回顾
  • GAMES202 作业 4

For reference👇：

master 分支上是 202 全部作业汇总，working 分支用于提交代码。需要原始作业可下载 master 分支的包

0x01 课程回顾 Recap

能量守恒的 BSDF Energy-Preserving BSDFs

Cook-Torrance 提出的微表面 BRDF 只考虑光线在表面经过一次弹射后出射到观察方向的能量，然而忽略了多次弹射最终出射到观察方向的能量，这就造成了能量损失。尤其是当材质的 roughness 较高时，G 项减小，多次弹射的比例增加，能量损失越严重。

Missing Energy of Microfacet BRDF, Revisiting Physically Based Shading at Imageworks, SIGGRAPH 2017

2017 年 Kulla 和 Conty 受到 Disney Principle BRDF 的启发，提出了一种改进的方法，即通过引入一个能量补偿项，将多次弹射出射到观察方向的能量补充回来，就能近似地保持能量守恒。

有了渲染方程，我们可以很容易地求得微表面 BRDF 在观察方向 $\mu_{o}$ 上出射的总能量：

$$ E(\mu_{o})=\int_{0}^{2\pi}\int_{0}^{1}f(\mu_{o},\mu_{i},\phi)\mu_{i}d\mu_{i}d\phi $$

这里记 $\mu_{i}=\sin \theta,\theta=n\cdot \mu_{i}$，之后就可以将 $f_{r}\cos \theta \sin \theta d\theta d\phi$ 写为 $f_{r}\mu_{i} d\mu_{i}d\phi$。然而原 PPT 里提到的是记为 $\cos \theta$，不过就积分结果来说都一样

$E(\mu_{o})$ 的范围为 0-1，对于不同的 $\mu_{o}$，$E(\mu_{o})$ 的值也会发生变化。微表面 BRDF 的 $E(\mu_{o})$ 会随着 roughness 的增加而衰减。为了弥补这部分能量，Kulla 和 Conty 设计了一个新的 BRDF lobe：

$$ f_{ms}=\frac{(1-E(\mu_{o}))(1-E(\mu_{i}))}{\pi(1-E_{avg})},E_{avg}=2\int_{0}^1E(\mu)\mu d\mu $$

其中 $f_{ms}$ 为弥补损失的 BRDF，$E_{avg}$ 为在所有观察方向下出射能量的平均值。同样计算这个 BRDF 在观察方向 $\mu_{o}$ 上的积分：

$$ \begin{aligned} E_{ms}(\mu_{o})&=\int_{0}^{2\pi}\int_{0}^{1}f_{ms}(\mu_{o},\mu_{i},\phi)\mu_{i}d\mu_{i}d\phi \\ &=2 \int_{0}^1 \frac{(1-E(\mu_{o}))(1-E(\mu_{i}))}{\pi(1-E_{avg})}\mu_{i}d\mu_{i}\\ &=2 \frac{1-E(\mu_{o})}{1-E_{avg}}\int_{0}^1(1-E(\mu_{i}))\mu_{i}d\mu_{i}\\ &= \frac{1-E(\mu_{o})}{1-E_{avg}}(1-E_{avg})\\ &= 1-E(\mu_{o}) \end{aligned} $$

将其加到原来的微表面 BRDF 上即可弥补损失：

$$ E(\mu_{o})+E_{ms}(\mu_{o})\equiv 1 $$

Micro Energy Compensation, Revisiting Physically Based Shading at Imageworks, SIGGRAPH 2017

计算 $E(\mu)$ 和 $E_{avg}$ 都需要积分，我们可以在离线端预先计算好这两项，将其按需要的变量打表，存储在纹理上，然后渲染时对应查询即可。PPT 中指出 32×32 的大小的纹理就足够了，作业中则使用了更大的 128×128。

颜色表面 Colored Surfaces

在上一节中，表面将所有入射的光全部反射，即光线不存在被吸收或者透射的情况。然而，绝大多数表面是带有颜色的，即其会吸收某种成分的光，这种情况下的能量损失是符合物理的，即 $E(\mu)<1$。

Kulla 和 Conty 受到 Jakob 等人的启发，多次弹射下能量是 diffuse 的，所以可以使用余弦加权的 Fresnel 项代表平均被反射的能量：

$$ F_{avg}=2\int_{0}^1F(\mu)\mu d\mu $$

我们已经知道所有反射出射的能量比例为 $E_{avg}$，那么因为微表面的遮挡而没有出射的能量比例则为 $1-E_{avg}$，这部分能量（即发生了二次 bounce）中又有 $F_{avg}$ 的比例能够被反射，所以经过 $n$ 次反射出射到观察方向的能量可以归纳为：

$$ \begin{aligned} E_{0}&=F_{avg}E_{avg}\\ E_{1}&=F_{avg}(1-E_{avg})F_{avg}E_{avg}\\ \dots\\ E_{n}&=F_{avg}^n (1-E_{avg})^n F_{avg}E_{avg} \end{aligned} $$

很显然这是一个等比数列，将其求和后即可得到所有经过多次弹射最终出射到观察方向的总能量：

$$ E=F_{avg}E_{avg}\sum_{k=0}^{\infty}F_{avg}^k(1-E_{avg})^k=\frac{F_{avg}E_{avg}}{1-F_{avg}(1-E_{avg})} $$

$E_{avg}$ 可以通过预计算得到，但是 $F_{avg}$ 无法这么做。所以 Kulla 和 Conty 分别对介电质和导体列出了两种前人用到的数值近似方法：

Dielectric Fresnel, Revisiting Physically Based Shading at Imageworks, SIGGRAPH 2017

Conductor Fresnel, Revisiting Physically Based Shading at Imageworks, SIGGRAPH 2017

如此，在渲染时即可采样 $E_{avg}$ 和根据表面性质近似 $F_{avg}$ 计算出颜色表面的出射能量比例，然后再将其与上节中没有能量吸收的 $f_{ms}$ 相乘即可。

0x02 问题修复 TroubleShooting

做作业时遇到一些出现问题的地方，这里一并提一下：

  1. 找不到 vec.h
  2. 未定义标识符“M_PI”
  3. uniform 变量 uEavgFLut 命名错误

第一个问题是配置问题，用 vscode 在 task.json 的 args 字段加上：

1 "-I",
2"${workspaceFolder}/assignment4/lut-gen/ext",

第二个问题把所有 M_PI 改成 PI。

第三个问题在 KullaContyMaterial.js 中，调用父类构造函数时传入了命名和 Shader 中不对应的 uniform 变量，把 uEavgFLut 改成 uEavgLut 即可。

0x03 预积分 $E(\mu)$

蒙特卡洛求解 Monte Carlo

在 Emu_MC.cpp 的 IntegrateBRDF 中完成对 $E(\mu)$ 的积分，使用蒙特卡洛方法在半球上均匀采样入射光方向，求解这个积分。积分内部为 BRDF 以及 cos 项：

$$ f_{r}(i,o)=\frac{F(i,h)G(i,o,h)D(h)}{4(n\cdot i)(n\cdot o)},\cos \theta=n\cdot i $$

F，G，D 的实现框架中均已给出，按需调用即可。

 1Vec3f IntegrateBRDF(Vec3f V, float roughness, float NdotV) {
 2    float A = 0.0;
 3    float B = 0.0;
 4    float C = 0.0;
 5    const int sample_count = 1024;
 6    Vec3f N = Vec3f(0.0, 0.0, 1.0);
 7    
 8    samplePoints sampleList = squareToCosineHemisphere(sample_count);
 9    for (int i = 0; i < sample_count; i++) {
10      // calculate (fr * ni) / p_o here
11      Vec3f L = sampleList.directions[i];
12      float pdf = sampleList.PDFs[i];
13
14      float NdotL = dot(N,L);
15      Vec3f H = normalize(V+L);
16
17      // here we assume all lights are reflected, no absorption!
18      // deal with color later
19      float F = 1.0;
20      float NDF = DistributionGGX(N, H, roughness);
21      float G = GeometrySmith(roughness, NdotV, NdotL);
22
23      float nom = F * NDF * G;
24      float denom = std::max(4.0 * NdotL * NdotV,0.001);
25      float fr = nom / denom;
26
27      // accumulate each component
28      A = B = C +=fr*NdotL/pdf;
29
30    }
31
32    return {A / sample_count, B / sample_count, C / sample_count};
33}

运行程序即可得到这样一张 LUT：

可以看到随着粗糙度降低（自上而下），噪声越来越明显。这是因为低粗糙度的微表面法线 $m$ 集中在宏观法线 $n$ 附近，但通过随机采样入射光方向 $i$ 计算得到的半程向量 $h$ 则不然，因此导致低粗糙度下积分值的方差显著增加。但我们可以通过对微表面法线 $m$ 进行重要性采样来改善这个问题，即使得采样的 $m$ 更多地位于 $n$ 附近。

重要性采样求解 Importance Sampling

最终的任务仍然是采样 $i$ 和计算 $\frac{f_{r}(i,o,h)(i,n)}{pdf_{i}(i)}$，但现在 $i$ 并不由直接采样得到，而是根据某个微表面法线分布模型重要性采样得到的微表面法向 $m$ 间接求解。

框架中使用的 NDF 为 GGX 分布，我们需要根据该分布生成采样的微表面法向 $m$，再由 $m$ 和 $o$ 反求入射方向 $i$。

由微表面法线分布模型的性质，可以计算出采样 $m$ 的概率密度：

$$ pdf_{m}(m)=D(m)(m\cdot n) $$

通过该概率密度，可以计算出该分布下对应的采样点：

$$ \begin{aligned} \theta_{m}&=\arctan\left( \frac{\alpha \sqrt{ \xi_{1} }}{\sqrt{ 1-\xi_{1} }} \right)\\ \phi_{h}&=2\pi \xi_{2} \end{aligned} $$

其中，$\alpha=roughness^{2},\xi_{1},\xi_{2}\in[0,1)$。

在 Emu_IS.cpp 的 ImportanceSampleGGX 函数中，完成重要性采样的逻辑，并返回切线空间的微表面法线 $h$。

 1Vec3f ImportanceSampleGGX(Vec2f Xi, Vec3f N, float roughness) {
 2    float a = roughness * roughness;
 3
 4    // in spherical space - Bonus 1
 5    float Xi1 = Xi.x;
 6    float Xi2 = Xi.y;
 7    float theta = atan(a*sqrt(Xi1)/sqrt(1-Xi1));
 8    float phi = 2.0 * PI * Xi2;
 9
10    // from spherical space to cartesian space - Bonus 1
11    float x = cos(phi)*sin(theta);
12    float y = sin(phi)*sin(theta);
13    float z = cos(theta);
14
15    // tangent coordinates - Bonus 1
16    Vec3f up = abs(N.z) < 0.999 ? Vec3f(0.0, 0.0, 1.0) : Vec3f(1.0, 0.0, 0.0); // N is fixed 001 so not select 001
17    Vec3f tangent = normalize(cross(up,N));
18    Vec3f bitangent = cross(N,tangent);
19
20    // transform H to tangent space - Bonus 1
21    Vec3f H = normalize(tangent*x + bitangent*y + N*z);
22    return H;
23}

接下来计算 $\frac{f_{r}(i,o,h)(i,n)}{pdf_{i}(i)}$ 的值，需要先将采样 $m$ 的概率密度转化为采样 $i$ 的概率密度：

$$ pdf_{i}(i)=pdf_{m}(m)\left\lvert \left\lvert \frac{\partial\omega_{m}}{\partial \omega_{i}} \right\rvert \right\rvert =\frac{D(m)(m\cdot n)}{4(i\cdot m)} $$

所以最终需要累加的权重为：

$$ \begin{aligned} weight(i)&=\frac{f_{r}(i,o,h)(i,n)}{pdf_{i}(i)} \\ &=\frac{F(i,m)G(i,o,m)D(m)}{4(n\cdot i)(n\cdot o)}\cdot(n\cdot i)\cdot \frac{4(i\cdot m)}{D(m)(n\cdot m)} \\ &=\frac{(o\cdot m)G(i,o,h)}{(o\cdot n)(m\cdot n)} \end{aligned} $$

这里暂且考虑所有能量都被反射，即 F 为 1；最后留在分子上的实际上是 $i\cdot m$，但是 $i\cdot m=o\cdot m$，所以可以替换

在 Emu_IS.cpp 的 IntegrateBRDF 函数中每次采样累加该权重即可：

 1Vec3f IntegrateBRDF(Vec3f V, float roughness) {
 2
 3    const int sample_count = 1024;
 4    Vec3f N = Vec3f(0.0, 0.0, 1.0);
 5    Vec3f Emu = Vec3f(0.0f);
 6    for (int i = 0; i < sample_count; i++) {
 7        Vec2f Xi = Hammersley(i, sample_count);
 8        Vec3f H = ImportanceSampleGGX(Xi, N, roughness);
 9        Vec3f L = normalize(H * 2.0f * dot(V, H) - V);
10
11        float NoL = std::max(L.z, 0.0f);
12        float NoH = std::max(H.z, 0.0f);
13        float VoH = std::max(dot(V, H), 0.0f);
14        float NoV = std::max(dot(N, V), 0.0f);
15        
16        // calculate (fr * ni) / p_o here - Bonus 1
17        float w = VoH * GeometrySmith(roughness,NoV,NoL)/(NoV*NoH);
18        Emu += Vec3f(w);
19        // Split Sum - Bonus 2
20    }
21
22    return Emu / sample_count;
23}

完成重要性采样之后得到的 $E(\mu)$ 就没有噪点了：

0x04 预积分 $E_{avg}$

预积分 $E_{avg}$ 的过程比较简单，甚至因为 $E(\mu),\mu$ 都是定值直接返回 $2E(\mu)\mu$ 即可。

 1Vec3f IntegrateEmu(Vec3f V, float roughness, float NdotV, Vec3f Ei) {
 2    return Ei * NdotV * 2.0f;
 3    Vec3f Eavg = Vec3f(0.0f);
 4    const int sample_count = 1024;
 5    Vec3f N = Vec3f(0.0, 0.0, 1.0);
 6
 7    for (int i = 0; i < sample_count; i++) 
 8    {
 9        Vec2f Xi = Hammersley(i, sample_count);
10        Vec3f H = ImportanceSampleGGX(Xi, N, roughness);
11        Vec3f L = normalize(H * 2.0f * dot(V, H) - V);
12
13        float NoL = std::max(L.z, 0.0f);
14        float NoH = std::max(H.z, 0.0f);
15        float VoH = std::max(dot(V, H), 0.0f);
16        float NoV = std::max(dot(N, V), 0.0f);
17
18        // calculate Eavg here - Bonus 1
19        Eavg += Ei * NoV * 2.0f;
20    }
21
22    return Eavg / sample_count;
23}

求解 $E_{avg}$ 的过程仅依赖于采样到的 $E(\mu)$ ，所以 Eavg_IS.cpp 和 E_avg_MC.cpp 的实现是一样的。

运行程序可以得到这样一张 LUT：

如果正确实现了 $E(\mu)$ 的重要性采样，上图就会变得 Smooth 一点：

0x05 实时渲染 Realtime Rendering

PBR 材质的实现就不贴了，公式已经给出了，可以直接搬离线端的实现

在实时端需要采样预积分的 $E(\mu)$ 和 $E_{avg}$ 的 LUT，然后仍然是套公式计算能量守恒的 BRDF：

$$ f_{r}=f_{micro}+f_{add}*f_{ms} $$

在 MultiScatterBRDF 函数中，需要完成对 $f_{add},f_{ms}$ 的计算：

$$ \begin{aligned} f_{add}&=\frac{F_{avg}E_{avg}}{1-F_{avg}(1-E_{avg})}\\ f_{ms}(\mu_{o},\mu_{i})&=\frac{(1-E(\mu_{o}))(1-E(\mu_{i}))}{\pi(1-E_{avg})} \end{aligned} $$

 1vec3 MultiScatterBRDF(float NdotL, float NdotV)
 2{
 3  vec3 albedo = pow(texture2D(uAlbedoMap, vTextureCoord).rgb, vec3(2.2));
 4
 5  vec3 E_o = texture2D(uBRDFLut, vec2(NdotL, uRoughness)).xxx;
 6  vec3 E_i = texture2D(uBRDFLut, vec2(NdotV, uRoughness)).xxx;
 7
 8  vec3 E_avg = texture2D(uEavgLut, vec2(0, uRoughness)).xxx;
 9  // copper
10  vec3 edgetint = vec3(0.827, 0.792, 0.678);
11  vec3 F_avg = AverageFresnel(albedo, edgetint);
12  
13  // calculate fms and missing energy here
14  vec3 f_add = F_avg * E_avg / (1.0 - F_avg * (1.0 - E_avg));
15  vec3 f_ms = (1.0 - E_o) * (1.0 - E_i) / (PI * (1.0 - E_avg));
16  return f_add*f_ms;
17  
18}

完成后即可达到如图所示的效果：

上：Kulla-Conty BSDF 下：Microfacet BRDF

0x06 Split Sum

仍然是计算 $E(\mu)$，结果保存到一张与 roughness 和曲率相关的 LUT 中，我们可以把 F 项从 BRDF 中分离出来（这个思路在 Lecture5 Real-time Environment Mapping 1h 15min 左右有提及）：

$$ \begin{aligned} E&=\int_{\Omega^+}f_{r}(n\cdot i) d\omega_{i} \\ &=\int_{\Omega^+} \frac{f_{r}}{F}(n\cdot i) Fd\omega_{i}\\ &=\int_{\Omega^+} \frac{f_{r}}{F}(n\cdot i) (F_{0}+(1-F_{0})(1-\cos \theta)^5)d\omega_{i}\\ &=\int_{\Omega^+} \frac{f_{r}}{F}(n\cdot i) (F_{0}+(1-F_{0})\alpha)d\omega_{i}\\ &=\int_{\Omega^+} \frac{f_{r}}{F}(n\cdot i) (F_{0}(1-\alpha)+\alpha)d\omega_{i}\\ &=F_{0}\int_{\Omega^+} \frac{f_{r}}{F}(n\cdot i) (1-\alpha)d\omega_{i}+\int_{\Omega^+} \frac{f_{r}}{F}(n\cdot i) \alpha d\omega_{i} \end{aligned} $$

记 $\alpha=(1-\cos \theta)^5$

更改一下 IntegrateBRDF，将这两个 Split 的结果分别存入 LUT 的 GB 通道，R 通道仍然存原来的 $E(\mu)$：

然后在实时端重载一个 MultiScatterBRDF，增加一个 F0 作为参数，用两个 Split 的积分还原 $E(\mu),E_{avg}$：

调用这个新的函数并传入 F0 即可看到使用 Split Sum 计算出来的结果：

上：Split Sum 下：Microfacet BRDF