PBR中引入IBL——镜面反射篇

PBR中引入IBL——镜面反射篇

回顾上一篇：漫反射

• Cook-Torrance反射方程分解
  $$L_o(p,{\omega }_o)=\underset{\Omega }{\int }(k_d\frac{c}{\pi })L_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i+\underset{\Omega }{\int }(k_s\frac{DFG}{4({\omega }_o\cdot n)({\omega }_i\cdot n)})L_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$$
  $$=L_o^{\prime }(p,{\omega }_o)+L_o^{\prime \prime }(p,{\omega }_o)$$

• 漫反射部分为：
  $$L_o^{\prime }(p,{\omega }_o)=k_d\frac{c}{\pi }\underset{\Omega }{\int }{L}_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$$

• 对于已给出的或计算出的环境立方体贴图envCubeTexture，在方向${\omega }_i$上采样，得到的是此方向上的radiance $L_i(p,{\omega }_i)$

• 我们的目的是得到$L_o^{\prime }(p,{\omega }_o)$，所以关键在计算$\underset{\Omega }{\int }{L}_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$

• 计算$\underset{\Omega }{\int }{L}_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$时，不可能从半球 Ω 的每个可能的方向采样envCubeTexture。不过我们可以对有限数量的方向采样以近似求解，在半球内均匀间隔取方向可以获得一个相当精确的积分结果，从而离散地计算积分

• 用有限数量的样本近似求解积分的方法叫黎曼和：${\int }_a^{b}f(x)dx={\sum }_{k=1}^{n}f({\epsilon }_k)\Delta x_k$，用黎曼和的方法求出积分结果后再进行卷积，在这里指的就是对纹理envCubeTexture的采样结果的进行卷积

• 对每个片段实时计算积分显然不现实，我们通过预计算，将积分结果提前存入纹理irradianceMap中

• 对纹理irradianceMap，在方向${\omega }_o$上采样，可以理解为场景中所有能够击中面向${\omega }_o$的表面的间接漫反射光的总和，即$L_o^{\prime }(p,{\omega }_o)$

镜面反射的计算思路

• 这一篇的目标是计算$L_o^{\prime \prime }(p,{\omega }_o)$，计算的关键同样是积分的计算
• 根据上一篇的经验，这一部分的处理思路也是先分解渲染方程，然后预计算基础环境纹理envCubeTexture的卷积，创建纹理保存卷积结果

镜面反射部分方程的分解

$$L_o^{\prime \prime }(p,{\omega }_o)=\underset{\Omega }{\int }(k_s\frac{DFG}{4({\omega }_o\cdot n)({\omega }_i\cdot n)})L_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$$
$$=\underset{\Omega }{\int }{f}_r(p,{\omega }_i,{\omega }_o)L_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$$
$$=\underset{\Omega }{\int }{L}_i(p,{\omega }_i)d{\omega }_i\ast \underset{\Omega }{\int }{f}_r(p,{\omega }_i,{\omega }_o)n\cdot {\omega }_{i}d{\omega }_i$$

• $L_o^{\prime \prime }(p,{\omega }_o)$被分成了前后两部分，所以要进行两次卷积，生成两个新的纹理

生成预计算纹理

• 对基础纹理envCubeTexture按照$\underset{\Omega }{\int }{L}_i(p,{\omega }_i)d{\omega }_i$进行卷积，将卷积的结果保存至新建的纹理中，这张纹理被称为预滤波环境贴图prefilterMap
• 在生成prefilterMap时，要考虑粗糙度的影响。因为随着粗糙度的增加，参与基础纹理envCubeTexture卷积的采样向量会更分散，导致反射更模糊，所以对于卷积的每个粗糙度级别，我们将按顺序把模糊后的结果存储在mipmap中
• 对基础纹理envCubeTexture按照$\underset{\Omega }{\int }{f}_r(p,{\omega }_i,{\omega }_o)n\cdot {\omega }_{i}d{\omega }_i$进行卷积，将卷积后的结果保存至2d查找纹理（Look Up Texture, LUT）中，这张纹理被称为BRDF积分贴图brdfLUT
• 我们假设每个方向的入射辐射度都是白色，在给定粗糙度、光线 ωi 法线 n 夹角 n⋅ωi 的情况下，预计算 BRDF 的响应结果。Epic Games 将预计算好的 BRDF 对每个粗糙度和入射角的组合的响应结果存储在一张 2d查找纹理上，即brdfLUT。这个纹理存储是菲涅耳响应的系数（R 通道）和偏差值（G 通道），它提供了分割版镜面反射积分的第二个部分
• 然后分别对这两张纹理采样，将采样结果相乘就得到了$L_o^{\prime \prime }(p,{\omega }_o)$，对吗？
• 也对也不对，因为其中还有极其复杂的细节。。。

必要的预备知识

镜面波瓣

• 镜面反射依赖于表面的粗糙度，反射光线可能比较松散，也可能比较紧密，但是一定会围绕着反射向量，除非表面极度粗糙
• 所有可能出射的反射光构成的形状称为镜面波瓣
• 在微表面模型里给定入射光方向，则镜面波瓣指向微平面的半程向量的方向

蒙特卡洛积分（Monte Carlo Integration）

$${\int }_a^{b}f(x)dx\approx \frac{1}{N}\sum _{i=0}^{N-1}\frac{f(x)}{pdf(x)}$$

• pdf (probability density function，概率密度函数)它的含义是特定样本在整个样本集上发生的概率

• 蒙特卡洛积分告诉我们：要计算等式左边的定积分，只用简单地从总体中随机挑选N个样本，并按照等式右边的式子进行计算，就能得到定积分的结果；并且随着样本数量的增加，结果会越来越接近定积分的精确结果

• 蒙特卡洛积分在计算机图形学中应用非常普遍，因为它是一种以高效的离散方式对连续的积分求近似而且非常直观的方法

重要性采样

• 重要性采样是蒙特卡洛积分的一种采样策略
• 蒙特卡洛积分中pdf(x)的意义在于：采样的样本不局限于均匀采样，我们可以根据被积分函数的形状去选择更高效的pdf(x)来指导采样点的生成，这也是重要性采样的核心思想
• 一个区间被采样到的概率越高，则这次采样在求和中的权重就越小，这也符合我们的直觉，假如不降低高频样本的权重，采样次数越多，估计出来的结果就越大，而不是收敛于期望

低差异序列（Low Discrepancy Sequence）

• 目的是生成分布均匀的随机数

• 对于一切需要采样的算法来说，分布均匀的随机数就意味着更加优秀的样本分布；在计算蒙特卡洛积分的过程中采样无处不在，所以好的样本分布直接影响积分过程的收敛速度

• 伪随机数组成的二维点集
  

• 低差异序列点集
  

预滤波环境贴图prefilterMap

生成prefilterMap的mipmap贴图

• 由于镜面反射受粗糙度的影响，所以对于每个粗糙度级别，我们将按顺序把模糊后的结果存储在预滤波贴图prefilterMap的mipmap中
• 为了确保为其 mip 级别分配足够的内存，一个简单方法是调用 glGenerateMipmap
• 注意将prefilterMap的缩小过滤器设置为 GL_LINEAR_MIPMAP_LINEAR 以启用三线性过滤

计算积分

• 根据蒙特卡洛积分算法，求解积分需要对样本进行采样，这里的样本指的就是入射光线${\omega }_i$
• 根据重要性采样的思想，把采样的范围锁定在镜面波瓣内才会有意义
• 由低差异序列可知，不同种类的随机数也会影响蒙特卡洛积分的结果
• 因此使用低差异序列生成的随机数，并采用重要性采样的思想，在镜面波瓣内进行采样，最终可以得到不错的结果

实现低差异序列

• 低差异序列也有很多种，这里使用 Hammersley 序列， Hammersley序列又是基于 Van Der Corput 序列
• Van Der Corput 序列，该序列是把十进制数字的二进制表示镜像翻转到小数点右边而得

float RadicalInverse_VdC(uint bits)
{
    bits = (bits << 16u) | (bits >> 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);
    return float(bits) * 2.3283064365386963e-10; // / 0x100000000 
}

vec2 Hammersley(uint i, uint N) 
{
    return vec2(float(i)/float(N), RadicalInverse_VdC(i)); 
}

实现重要性采样

• 这里的重要性采样的目标是：使采样基础纹理的样本的反射光线在镜面波瓣内；为了实现这个目标，我们需要将重要性采样的结果作为半程向量使用，然后用反射光线和半程向量计算出采样基础纹理的样本，这也是生成prefilterMap的关键步骤
• 重要性采样的过程：开始一个大循环，生成一个低差异序列值，用该序列值在切线空间中生成样本向量，将样本向量变换到世界空间

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);
       vec3 H;
       H.x = cos(phi) * sinTheta;
       H.y = sin(phi) * sinTheta;
       H.z = cosTheta;
       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);
       vec3 sampleVec = tangent * H.x + bitangent * H.y + N * H.z;
       return normalize(sampleVec);
}

const uint SAMPLE_COUNT = 4096u;
for(uint i = 0u; i < SAMPLE_COUNT; ++i)
{ 
    vec2 Xi = Hammersley(i, SAMPLE_COUNT);
    vec3 H = ImportanceSampleGGX(Xi, N, roughness);
}

prefilterMap的着色器

• 利用重要性采样生成的半程向量H和观察向量V，计算采样基础纹理的样本L，即入射光线${\omega }_i$
• Epic将观察方向V近似为光线的出射方向
• 在立方体贴图中，顶点坐标又可以认为是法线N和观察方向V
• 计算L的方向跟施密特正交化的方法很像，但不是生成正交向量，而是V关于H的反射向量：vec3 L = normalize(2.0 * dot(V, H) * H - V);

void main() 
{ 
    vec3 N = normalize(localPos); 
    vec3 V = N; 
    const uint SAMPLE_COUNT = 1024u; 
    float totalWeight = 0.0; 
    vec3 prefilteredColor = vec3(0.0); 
    for(uint i = 0u; i < SAMPLE_COUNT; ++i) 
    { 
        vec2 Xi = Hammersley(i, SAMPLE_COUNT); 
        vec3 H = ImportanceSampleGGX(Xi, N, roughness); 
        vec3 L = normalize(2.0 * dot(V, H) * H - V); 
        float NdotL = max(dot(N, L), 0.0); 
        if(NdotL > 0.0) 
        { 
            prefilteredColor += texture(envCubeTexture, L).rgb * NdotL; 
            totalWeight += NdotL; 
        } 
    } 
    prefilteredColor = prefilteredColor / totalWeight; 
    FragColor = vec4(prefilteredColor, 1.0); 
}

BRDF积分贴图brdfLUT

整理积分方程

• 首先通过变形把$\underset{\Omega }{\int }{f}_r(p,{\omega }_i,{\omega }_o)n\cdot {\omega }_{i}d{\omega }_i$中菲涅尔系数$F_0$移到积分外部
• 菲涅尔函数$F_{Schlick}(h,v,F_0)=F_0+(1-F_0)(1-h\cdot v{)}^5$，其中的v就是${\omega }_o$，$F_0$是个常数，因此又可以写成$F({\omega }_o,h)=F_0+(1-F_0)(1-{\omega }_o\cdot h{)}^5$

$$\underset{\Omega }{\int }{f}_r(p,{\omega }_i,{\omega }_o)n\cdot {\omega }_{i}d{\omega }_i=\underset{\Omega }{\int }{f}_r(p,{\omega }_i,{\omega }_o)\frac{F({\omega }_o,h)}{F({\omega }_o,h)}n\cdot {\omega }_{i}d{\omega }_i$$
$$=\underset{\Omega }{\int }\frac{f_r(p,{\omega }_i,{\omega }_o)}{F({\omega }_o,h)}F({\omega }_o,h)n\cdot {\omega }_{i}d{\omega }_i$$
$$=\underset{\Omega }{\int }\frac{f_r(p,{\omega }_i,{\omega }_o)}{F({\omega }_o,h)}(F_0+(1-F_0)(1-{\omega }_o\cdot h{)}^5)n\cdot {\omega }_{i}d{\omega }_i$$
用$\alpha$代替$(1-{\omega }_o\cdot h{)}^5$，得：
$$=\underset{\Omega }{\int }\frac{f_r(p,{\omega }_i,{\omega }_o)}{F({\omega }_o,h)}(F_0+(1-F_0)\alpha )n\cdot {\omega }_{i}d{\omega }_i$$
$$=\underset{\Omega }{\int }\frac{f_r(p,{\omega }_i,{\omega }_o)}{F({\omega }_o,h)}(F_0+\alpha -F_0\alpha )n\cdot {\omega }_{i}d{\omega }_i$$
$$=\underset{\Omega }{\int }\frac{f_r(p,{\omega }_i,{\omega }_o)}{F({\omega }_o,h)}(F_0(1-\alpha )+\alpha )n\cdot {\omega }_{i}d{\omega }_i$$
$$=\underset{\Omega }{\int }\frac{f_r(p,{\omega }_i,{\omega }_o)}{F({\omega }_o,h)}(F_0(1-\alpha ))n\cdot {\omega }_{i}d{\omega }_i+\underset{\Omega }{\int }\frac{f_r(p,{\omega }_i,{\omega }_o)}{F({\omega }_o,h)}\alpha \cdot n\cdot {\omega }_{i}d{\omega }_i$$
用$f_r$中得菲涅尔函数与分子进行约分，把$F_0$提到积分外，再替换$\alpha$，得：
$$=F_0\underset{\Omega }{\int }{f}_r^{\prime }(p,{\omega }_i,{\omega }_o)(1-(1-{\omega }_o\cdot h{)}^5)n\cdot {\omega }_{i}d{\omega }_i+\underset{\Omega }{\int }{f}_r^{\prime }(p,{\omega }_i,{\omega }_o)(1-{\omega }_o\cdot h{)}^5\cdot n\cdot {\omega }_{i}d{\omega }_i$$
最终这个式子的前半部分表示$F_0$的比例，后半部分表示偏差，我们把这两个部分分别计算和存储，所以就成了开始时候提到的2d纹理

brdfLUT着色器

//...
const uint SAMPLE_COUNT = 1024u; 
for(uint i = 0u; i < SAMPLE_COUNT; ++i) 
{ 
    vec2 Xi = Hammersley(i, SAMPLE_COUNT); 
    vec3 H = ImportanceSampleGGX(Xi, N, roughness); 
    vec3 L = normalize(2.0 * dot(V, H) * H - V); 
    float NdotL = max(L.z, 0.0); //??
    float NdotH = max(H.z, 0.0); //??
    float VdotH = max(dot(V, H), 0.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; 
    } 
} 
A /= float(SAMPLE_COUNT); 
B /= float(SAMPLE_COUNT); 
return vec2(A, B); 

• 注意：几何函数$G_{SchlickGGX}(n,v,k)=\frac{n\cdot v}{(n\cdot v)(1-k)+k}$中，k的值不再是$k_{direct}=\frac{(\alpha +1{)}^2}{8}$，而是$k_{I}BL=\frac{{\alpha }^2}{2}$，$\alpha$表示粗糙度

完成IBL的镜面反射

• 首先，使用反射向量采样预过滤的环境贴图，获取表面的间接镜面反射。请注意，我们会根据表面粗糙度在合适的 mip 级别采样，以使更粗糙的表面产生更模糊的镜面反射。
• 然后我们用已知的材质粗糙度和视线-法线夹角作为输入，采样 brdfLUT

uniform samplerCube prefilterMap; 
uniform sampler2D brdfLUT;
void main()
{
//1.
vec3 R = reflect(-V, N); 
const float MAX_REFLECTION_LOD = 4.0; 
vec3 prefilteredColor = textureLod(prefilterMap, R, roughness * MAX_REFLECTION_LOD).rgb;
//2.
vec3 F = FresnelSchlickRoughness(max(dot(N, V), 0.0), F0, roughness); 
vec2 envBRDF = texture(brdfLUT, vec2(max(dot(N, V), 0.0), roughness)).rg; 
vec3 specular = prefilteredColor * (F * envBRDF.x + envBRDF.y);//直接用间接光菲涅尔项F代替F0
}

• 至此IBL的镜面反射计算完成
  

完成IBL

• 上一篇完成的漫反射 + 这一篇完成的镜面反射

vec3 F = FresnelSchlickRoughness(max(dot(N, V), 0.0), F0, roughness);

vec3 kS = F; 
vec3 kD = 1.0 - kS; 
kD *= 1.0 - metallic; 

vec3 irradiance = texture(irradianceMap, N).rgb; 
vec3 diffuse = irradiance * albedo; 

const float MAX_REFLECTION_LOD = 4.0; 
vec3 prefilteredColor = textureLod(prefilterMap, R, roughness * MAX_REFLECTION_LOD).rgb; 
vec2 envBRDF = texture(brdfLUT, vec2(max(dot(N, V), 0.0), roughness)).rg; 
vec3 specular = prefilteredColor * (F * envBRDF.x + envBRDF.y); 

vec3 ambient = (kD * diffuse + specular) * ao;

• 这里的specular 没有乘以 kS，因为已经乘过菲涅耳系数了
  

终极合并——全局光照

• 直接光照 + 间接光照（IBL）=全局光照

vec3 N = Normal;
vec3 V = normalize(camPos - WorldPos);
vec3 R = reflect(-V, N);
vec3 F0 = vec3(0.04); 
vec3 Lo = vec3(0.0);
F0 = mix(F0, albedo, metallic);
vec3 L = normalize(lightPositions[i] - WorldPos);
vec3 H = normalize(V + L);
//...
vec3 kS = F;
Lo += (kD * albedo / PI + specular) * radiance * NdotL;//直接光照

vec3 F = fresnelSchlickRoughness(max(dot(N, V), 0.0), F0, roughness);
//...
vec3 ambient = (kD * diffuse + specular) * ao;//间接光照（IBL）

vec3 color = ambient + Lo;//全局光照
color = color / (color + vec3(1.0));// HDR tonemapping
color = pow(color, vec3(1.0/2.2)); // gamma correct
FragColor = vec4(color , 1.0);