Learn OpenGL：PBR

发表于 2021-09-25 分类于 CG ， OpenGL Waline：

本文字数： 18k 阅读时长 ≈ 16 分钟

Learn OpenGL: PBR

理论

PBR也就是基于物理的渲染(Physically Based Rendering)，物理指的是物理学。基于物理的渲染满足三个条件：

基于微平面(Microfacet)的表面模型
能量守恒
应用基于物理的BRDF

微表面模型是站在微观角度下看待平面，没有任何一个平面是完全光滑的，所有的表面都是由细小的镜面组成。这些细小的镜面排列越是一致，则表面越是光滑；越粗糙，排列也就越混乱。

这种粗糙程度的量化来说就是：微表面的平均去向与半程向量方向一致的概率。

根据人眼观察到的现象，越是光滑的表面，高光就越是亮而且集中，如下图所示，粗糙度越大则高光的面积就越大，但是高光也越暗

所以是概率越大所以才导致的高光越大？

关于光线的反射与折射

这幅图来自于The Comprehensive PBR Guide。指的是每种材料都是由细小的粒子组成，光线入射之后部分会发生直接反射，部分会进入表面之下。进入表面之下的部分光线会经过几次反射之后再离开表面，而一些光线则会不断反射直到被物体吸收。

直接在表面发生反射的光就成为镜面光（Specular），而经过内部折射几次之后才反射出来的光称为漫反射光（Diffuse）。

能量守恒

反射和折射的能量之和不能超过入射的能量。反射部分形成镜面光照，折射部分进入表面被吸收形成漫反射光照

反射率方程

\[ L_o(p,\omega _o)=\underset{\Omega}{\int}f _r(p,\omega _i,\omega _o)L_i(p,\omega _i)n \cdot \omega _i d\omega _i \]

反射率方程概括了在半球面\(\Omega\)内，打到点\(p\)上的所有光线辐射率之和。

\(n\)是点\(p\)处的法向量；\(f_r\)（BRDF方程）描述了表面的反射规律
\(L_i\)和\(L_o\)分别表示入射光和出射光的辐照度
\(w_i\)和\(w_o\)分别表示入射光和出射光的方向

BRDF

BRDF也就是双向反射分布函数(Bidirectional Reflective Distribution Function)

输入：入射方向\(\omega_i\)，出射方向\(\omega_o\)，平面法线\(n\)，微平面粗糙程度\(a\)
输出：每束光线对一个给定了材质属性的平面上最终反射出来的光线所作出的贡献程度

一般使用的近似模型是Cook-Torrance模型： \[ f_r=k_df_d+k_sf_s \] 左半部分是漫反射，右半部分是镜面反射，\(k_d\)和\(k_s\)分别表示入射光线中被折射和被反射的比率，\(k_d+k_s=1\)

漫反射部分

\(f_d=f_{lambert}=\frac{c}{\pi}\)，\(c\)是表面颜色

镜面反射部分

\[ \begin{aligned} k_s &= F \\ f_s &= f_{cook-torrance} = \frac{DFG}{4(\omega_o\cdot n)(\omega_i \cdot n)} \end{aligned} \]

\(D\)是法线分布函数(Normal Distribution Function)：
Trowbridge-Reitz GGX： \[ NDF_{GGXTR}(n,h,\alpha)=\frac{\alpha^2}{\pi((n\cdot h)^2(\alpha^2-1)+1)^2} \]
\(F\)是菲涅尔系数(Fresnel Equation)：描述反射光所占的比率（视线垂直于表面时，反射较弱）。\(F_0\)是基础反射率
Schlick简化版： \[ F=F_0+(1-F_0)(1-cos\theta)^5 \] 虚幻加速版： \[ F=F_0+(1-F_0)2^{(-5.55473(v\cdot h)-6.98316)(v\cdot h)} \]
\(G\)是几何函数(Geometry Function)：描述微表面之间相互遮蔽的比率
GGX与Schlick-Beckmann近似的结合体，称为Schlick-GGX： \[ G_{SchlickGGX}(n,v,k)=\frac{n\cdot v}{(n\cdot v)(1-k)+k} \] 令\(\alpha=roughness^2\)
对于直接光照：\(k=k_{direct}=\frac{(\alpha+1)^2}{8}\)
对于IBL：\(k=k_{IBL}=\frac{\alpha^2}{2}\)

编写PBR材质-金属工作流

每一个表面参数都可以用纹理保存。

金属工作流引入了金属度，非金属默认0.04，金属则表示表面颜色

反照率Albedo：指定表面颜色
法线Normal：指定表面法线
金属度Metallic：表达金属质地，可以用灰度图
粗糙度Roughness：越粗糙，对应的高光越大而且糊。对应还有光滑度(Smoothness)贴图，糙度=1.0-光滑度
环境光遮蔽AO：可以用来指定表面的阴影，比如砖块的裂缝就应该比较暗

光照

最终的反射方程 \[ L_o(p, w_o)=\underset{\Omega}\int{(k_d\frac{c}{\pi}+k_s\frac{DFG}{4(w_o\cdot n)(w_i\cdot n)})L_i(p,w_i)n\cdot w_i d w_i} \] 因为\(k_s\)就是\(F\)，都是指反射所占比例，应该只乘一次，所以方程变为 \[ L_o(p, w_o)=\underset{\Omega}\int{(k_d\frac{c}{\pi}+\frac{DFG}{4(w_o\cdot n)(w_i\cdot n)})L_i(p,w_i)n\cdot w_i d w_i} \]

直接光照

变量声明

#version 330 core
out vec4 FragColor;
in vec2 TexCoords;
in vec3 WorldPos;
in vec3 Normal;

// material parameters
uniform vec3  albedo;
uniform float metallic;
uniform float roughness;
uniform float ao;

// lights
uniform vec3 lightPositions[4];
uniform vec3 lightColors[4];

uniform vec3 camPos;

const float PI = 3.14159265359;

D、F、G方程声明

// 正态分布函数, N:法线, H:半程向量, roughness:糙度
float DistributionGGX(vec3 N, vec3 H, float roughness)
{
    float a      = roughness*roughness;
    float a2     = a*a;
    float NdotH  = max(dot(N, H), 0.0);
    float NdotH2 = NdotH*NdotH;

    float nom   = a2;
    float denom = (NdotH2 * (a2 - 1.0) + 1.0);
    denom = PI * denom * denom;

    return nom / denom;
}

// 菲涅尔方程的Fresnel-Schlick近似, cosTheta:视线与表面法线的夹角余弦, F0:以0°入射时的反射百分比
// 返回光线被反射的百分比
vec3 fresnelSchlick(float cosTheta, vec3 F0)
{
    return F0 + (1.0 - F0) * pow(1.0 - cosTheta, 5.0);
}  

float GeometrySchlickGGX(float NdotV, float roughness)
{
    float r = (roughness + 1.0);
    float k = (r*r) / 8.0;

    float nom   = NdotV;
    float denom = NdotV * (1.0 - k) + k;

    return nom / denom;
}

// 几何遮蔽函数, N:法线，V:视线, L:光线, roughness:糙度
float GeometrySmith(vec3 N, vec3 V, vec3 L, float roughness)
{
    float NdotV = max(dot(N, V), 0.0);
    float NdotL = max(dot(N, L), 0.0);
    float ggx2  = GeometrySchlickGGX(NdotV, roughness);
    float ggx1  = GeometrySchlickGGX(NdotL, roughness);

    return ggx1 * ggx2;
}

主函数

void main()
{       
    vec3 N = normalize(Normal);
    vec3 V = normalize(camPos - WorldPos);

    vec3 F0 = vec3(0.04); // 以0°入射时的反射百分比
         F0 = mix(F0, albedo, metallic); // 非金属一直是0.04，视觉上正确

    // 计算反射方程的结果
    vec3 Lo = vec3(0.0);
    for(int i = 0; i < 4; ++i) 
    {
        // calculate per-light radiance
        vec3 L = normalize(lightPositions[i] - WorldPos);
        vec3 H = normalize(V + L);
        float distance    = length(lightPositions[i] - WorldPos);
        float attenuation = 1.0 / (distance * distance);
        vec3  radiance    = lightColors[i] * attenuation;        

        // cook-torrance brdf
        float NDF = DistributionGGX(N, H, roughness);        // 正态
        float G   = GeometrySmith(N, V, L, roughness);       // 几何
        vec3  F   = fresnelSchlick(max(dot(H, V), 0.0), F0); // 菲涅尔

        vec3 kS = F; // 反射占比与菲涅尔方程结果相等
        vec3 kD = vec3(1.0) - kS;  // 折射占比
             kD *= 1.0 - metallic; // 金属不折射光线（漫反射为0），设置折射占比为0

        vec3  nominator   = NDF * G * F;
        float denominator = 4.0 * max(dot(N, V), 0.0) * max(dot(N, L), 0.0) + 0.001; // 防止除零错误  
        vec3  specular    = nominator / denominator;

        // 累加出射光线的辐射率，因为已知有限个光源，不需要积分
        float NdotL = max(dot(N, L), 0.0);                
        Lo += (kD * albedo / PI + specular) * radiance * NdotL; 
    }   
    // 加上环境光
    vec3 ambient = vec3(0.03) * albedo * ao;
    vec3 color   = ambient + Lo;
    
    // Reinhard色调映射（HDR）
    color = color / (color + vec3(1.0));
    
    // gamma矫正（PBR要求所有输入都是线性的，所以之前都是在线性空间中计算）
    color = pow(color, vec3(1.0/2.2));  

    FragColor = vec4(color, 1.0);
}

使用纹理的话，之前的参数都可以采样获得

着色器采样

// 采样获得各参数
[...]
uniform sampler2D albedoMap;
uniform sampler2D normalMap;
uniform sampler2D metallicMap;
uniform sampler2D roughnessMap;
uniform sampler2D aoMap;

void main()
{
    vec3  albedo    = pow(texture(albedoMap, TexCoords).rgb, 2.2); // 从sRGB转到线性空间
    vec3  normal    = getNormalFromNormalMap();
    float metallic  = texture(metallicMap, TexCoords).r;
    float roughness = texture(roughnessMap, TexCoords).r;
    float ao        = texture(aoMap, TexCoords).r;
    [...]
}

IBL

基于图像的光照(Image based lighting, IBL)，光源不是一个一个可遍历的，而是来自整体周围环境，也就是立方体贴图的每个像素。

对于单独一个一个的光源可以把反射方程的积分改用加法计算，但是这里只能对半球进行积分。为了提高效率，采用预计算的方法提前计算积分结果，渲染时直接采样即可

对原反射方程进行变换，分成漫反射和镜面反射两部分单独进行预计算 \[ \begin{aligned} L_o(p, w_o)&=\underset{\Omega}\int{(k_d\frac{c}{\pi} + k_s\frac{DFG}{4(w_o\cdot n)(w_i\cdot n)})L_i(p,w_i)n\cdot w_i d w_i} \\ &=\underset{\Omega}\int{(k_d\frac{c}{\pi})L_i(p,w_i)n\cdot w_i d w_i} + \underset{\Omega}\int{(k_s\frac{DFG}{4(w_o\cdot n)(w_i\cdot n)})L_i(p,w_i)n\cdot w_i d w_i} \end{aligned} \]

HDR

为了让环境光照更真实，需要用HDR保存（颜色范围超出[0,1]），而之前使用的天空盒（立方体贴图）只是LDR。在积分之前先要从hdr纹理获取环境数据

载入.hdr纹理

用std_image载入HDR纹理

stbi_set_flip_vertically_on_load(true);
int width, height, nrComponents;
float *data = stbi_loadf("newport_loft.hdr", &width, &height, &nrComponents, 0); // 结果是一个浮点数组
unsigned int hdrTexture;
if (data)
{
    glGenTextures(1, &hdrTexture);
    glBindTexture(GL_TEXTURE_2D, hdrTexture);
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB16F, width, height, 0, GL_RGB, GL_FLOAT, data); 

    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

    stbi_image_free(data);
}
else
{
    std::cout << "Failed to load HDR image." << std::endl;
}

得到hdrTexture，数据就是hdr里的浮点数

采样hdr数据

.hdr文件是将一个360°的全景图保存到了二维平面上，也就是从球坐标系\((\theta, \phi)\)转换到了二维坐标系\((u,v)\)（也称为等距柱状投影图）

给定对立方体贴图采样的方向向量\((x, y, z)\)，对应到球面上的一个点\((\theta, \phi)\)，满足等式： \[ \begin{aligned} x &= cos(\theta)cos(\phi) \\ y &= sin(\theta) \\ z &= cos(\theta)sin(\phi) \end{aligned} \\ \] 其中\(\theta\)是俯仰角(pitch)，满足\(\theta \in [-\frac{\pi}{2},\frac{\pi}{2}]\)，\(\phi\)是方位角(yaw)，满足\(\phi \in [-\pi,\pi]\)

可以算出： \[ \begin{aligned} \phi &= atan(\frac{z}{x}) \\ \theta &= asin(y) \end{aligned} \] 又因为\((u,v)\)的范围都是\([0,1]\)，可以线性映射得到： \[ \begin{aligned} u &= \frac{\phi}{2\pi} + 0.5 \\ v &= \frac{\theta}{\pi} + 0.5 \end{aligned} \] 至此完成从采样方向向量\((x, y, z)\)到纹理坐标\((u,v)\)的转换，也就可以像是对立方体贴图采样一样对HDR贴图采样了

顶点着色器

#version 330 core
layout (location = 0) in vec3 aPos;

out vec3 localPos;

uniform mat4 projection;
uniform mat4 view;

void main()
{
    localPos = aPos;  // 局部坐标作为采样向量
    gl_Position =  projection * view * vec4(localPos, 1.0);
}

片段着色器

#version 330 core
out vec4 FragColor;
in vec3 localPos;

uniform sampler2D equirectangularMap;

const vec2 invAtan = vec2(0.1591, 0.3183); // (1/2PI, 1/PI)
vec2 SampleSphericalMap(vec3 v)
{
    vec2 uv = vec2(atan(v.z, v.x), asin(v.y)); // (phi, theta)
    uv *= invAtan;
    uv += 0.5;
    return uv;
}

void main()
{       
    vec2 uv = SampleSphericalMap(normalize(localPos)); // make sure to normalize localPos
    vec3 color = texture(equirectangularMap, uv).rgb;

    FragColor = vec4(color, 1.0);
}

转换成立方体贴图

可以对hdr采样之后，再利用帧缓冲对象，进行六次渲染：

把HDR当成立方体贴图进行采样，并设置摄像机矩阵保证只渲染完整的一个面（fov为90°，正对当前面），然后把目标立方体贴图的一个面绑定为渲染缓冲的颜色附件，最后渲染一个立方体。这样执行六次，从而把对应的面都渲染到立方体贴图上。

得到立方体贴图称为envCubemap

漫反射积分

得到envCubeMap之后就可以开始进行积分，这里是漫反射部分 \[ \begin{aligned} L_o(p, w_o) & =\underset{\Omega}\int{(k_d\frac{c}{\pi})L_i(p,w_i)n\cdot w_i d w_i} \\ L_o(p, w_o) & =k_d\frac{c}{\pi}\underset{\Omega}\int{L_i(p,w_i)n\cdot w_i d w_i} \end{aligned} \]

对于点\(p\)，漫反射部分的积分只跟立体角\(w_i\)有关，而立体角是球表面积与半径平方的比值，可以用俯仰角\(\theta\)（pitch）和方位角\(\phi\)（yaw）表示 \[ dw = \frac{dA}{r^2}=\frac{(rd\theta)(rsin\theta d\phi)}{r^2}=sin\theta d\theta d\phi \] 因为是半球，其中\(\theta \in [0,\frac{\pi}{2}]\)，\(\phi \in [0,2\pi]\)

反射方程变为： \[ L_o(p,{\phi}_o,{\theta}_o)=k_d\frac{c}{\pi}\int_{\phi=0}^{2\pi}\int_{\theta=0}^{\frac{\pi}{2}}L_i(p,{\phi}_i,{\theta}_i)cos{\theta}sin{\theta}d{\phi}d{\theta} \] 对这个积分用均匀采样的蒙特卡洛方法进行计算： \[ \begin{aligned} & \frac{1}{\pi}\int_0^{2\pi}\int_0^{\frac{\pi}{2}}L_i(p,w_i)ncos\theta sin\theta d\theta d\phi\\ \approx &\frac{1}{\pi}\frac{1}{N_1N_2}\sum_i^{N_1}\sum_j^{N-2}\frac{L_i(p,\phi_j,\theta_i)cos\theta_i sin\theta_i}{p(\theta_i,\phi_j)}\\ =& \frac{1}{\pi}\frac{1}{N_1N_2}\sum_i^{N_1}\sum_j^{N-2}\frac{L_i(p,\phi_j,\theta_i)cos\theta_i sin\theta_i}{\frac{1}{2\pi*0.5\pi}}\\ =& {\pi}\frac{1}{N_1N_2}\sum_i^{N_1}\sum_j^{N-2}{L_i(p,\phi_j,\theta_i)cos\theta_i sin\theta_i}\\ \end{aligned} \] 因为\(\pi\)是常数，所以也直接参与计算保存起来

为了保存积分结果，首先需要创建一个立方体贴图。因为积分中用到了均值，丢失了大部分高频信息，所以分辨率可以低一些。接着对之前转化而来的立方体贴图渲染六次，分别绑定新的立方体贴图六个面。

这样渲染得到新的立方体贴图称为irradianceMap

用到的顶点着色器跟之前的一样，而片段着色器进行数值积分，采样步长是sampleDelta

片段着色器

#version 330 core
out vec4 FragColor;
in vec3 worldPos;

uniform samplerCube environmentMap; // HDR转换来的立方体贴图envCubeMap

const float PI = 3.14159265359;

void main()
{		
    vec3 N = normalize(worldPos);

    vec3 irradiance = vec3(0.0);   
    
    // 计算切线空间
    vec3 up    = vec3(0.0, 1.0, 0.0);
    vec3 right = normalize(cross(up, N));
    up = normalize(cross(N, right));
       
    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)
        {
            // 球面到笛卡尔（在切线空间中）
            vec3 tangentSample = vec3(sin(theta) * cos(phi),  sin(theta) * sin(phi), cos(theta));
            // 切线空间到世界空间
            vec3 sampleVec = tangentSample.x * right + tangentSample.y * up + tangentSample.z * N; 
            // 这里乘了额外的系数cos和sin
            irradiance += texture(environmentMap, sampleVec).rgb * cos(theta) * sin(theta);
            nrSamples++;
        }
    }
    irradiance = PI * irradiance * (1.0 / float(nrSamples)); // 得到平均采样辐照度
    
    FragColor = vec4(irradiance, 1.0);
}

其中在累加之前乘了额外的系数：

\(sin\theta\)是因为俯仰角\(\theta\)比较大的时候，靠近球顶，采样面积会变小，所以用\(sin\theta\)进行平衡
\(cos\theta\)是因为俯仰角较大的光线的权重应该比较小

作为环境光使用

直接光照用的shader里的环境光是常量，现在可以用公式进行计算了：采样irradianceMap得到预积分结果并乘以漫反射率比率。

这个比率可以由菲尼尔公式得到，但是为了防止粗糙非金属表面的菲涅尔过强，需要给原公式引入粗糙度：

fresnelSchlickRoughness

vec3 fresnelSchlickRoughness(float cosTheta, vec3 F0, float roughness)
{
    return F0 + (max(vec3(1.0 - roughness), F0) - F0) * pow(1.0 - cosTheta, 5.0);
}

计算环境光

uniform samplerCube irradianceMap;
[...]
vec3 kS = fresnelSchlickRoughness(max(dot(N, V), 0.0), F0, roughness); 
vec3 kD = 1.0 - kS; // 折射占比，也即漫反射比率
vec3 irradiance = texture(irradianceMap, N).rgb;
vec3 diffuse    = irradiance * albedo;
vec3 ambient    = (kD * diffuse) * ao;

镜面反射积分

现在计算镜面反射部分 \[ \begin{aligned} L_o(p, w_o) =& \underset{\Omega}\int{(k_s\frac{DFG}{4(w_o\cdot n)(w_i\cdot n)})L_i(p,w_i)n\cdot w_i d w_i} \\ {f_r(p,w_i,w_o)=\frac{DFG}{4(w_o\cdot n)(w_i\cdot n)}}\ {\Longrightarrow}&\ \underset{\Omega}\int{f_r(p,w_i,w_o)L_i(p,w_i)n\cdot w_i d w_i} \end{aligned} \] 积分跟入射\(w_i\)和出射\(w_o\)都有关，实时渲染里肯定不能对每种\(w_o\)的取值都进行积分。这里采用虚幻的「分割求和近似法」转换成一个卷积 \[ L_o(p,w_o)\approx \underset{\Omega}{\int}L_i(p,w_i)dw_i * \underset{\Omega}{\int}f_r(p,w_i,w_o)n{\cdot}w_i dw_i \] 这个卷积有两个独立的积分

第一部分积分

\[ \underset{\Omega}{\int}L_i(p,w_i)dw_i \]

只跟\(w_i\)有关，跟漫反射积分一样可以预计算后保存到立方体贴图上

因为越粗糙的反射效果越模糊，采样向量越分散，所以不同的粗糙度积分结果会不同。因为要准备多张贴图而且一个比一个模糊，正好可以用mipmap实现，每个级别的mipmap就对应一个级别的粗糙度。

具体过程是先生成一张新的立方体贴图用来保存积分结果，而且需要开启生成mipmap以及三线性过滤。每个面的分辨率可以是128x128，更光滑的反射需要提高分辨率，生成的新立方体贴图称为prefilterMap

这里不像漫反射一样对半球进行积分，而是围绕反射向量的一个镜面波瓣进行重要性采样。因为越粗糙对应的反射向量范围越大

而且为了快速使积分结果快速向精确值收敛，围绕波瓣范围的同时还需要加入一些随机值，使用低差异序列的过程称为「拟蒙特卡洛积分」

GGX重要性采样

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

    // from spherical coordinates to cartesian coordinates
    vec3 H;
    H.x = cos(phi) * sinTheta;
    H.y = sin(phi) * sinTheta;
    H.z = cosTheta;

    // from tangent-space vector to world-space sample vector
    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);
}

Hammersley低差异序列

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
}
// 样本索引i，总共N个样本
vec2 Hammersley(uint i, uint N)
{
    return vec2(float(i)/float(N), RadicalInverse_VdC(i));
}

prefliter shader

#version 330 core
out vec4 FragColor;
in vec3 localPos;

uniform samplerCube environmentMap; // HDR转换而来的envCubeMap
uniform float roughness;

const float PI = 3.14159265359;

float RadicalInverse_VdC(uint bits);
vec2 Hammersley(uint i, uint N);
vec3 ImportanceSampleGGX(vec2 Xi, vec3 N, float roughness);

void main()
{       
    vec3 N = normalize(localPos);    
    vec3 R = N;
    vec3 V = R;

    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(environmentMap, L).rgb * NdotL;
            totalWeight      += NdotL;
        }
    }
    prefilteredColor = prefilteredColor / totalWeight;

    FragColor = vec4(prefilteredColor, 1.0);
}

要对每种mipmap级别进行积分，算出对应的粗糙度传递给着色器

为防止立方体贴图的面之间出现明显接缝，需要

1	glEnable(GL_TEXTURE_CUBE_MAP_SEAMLESS);

为防止贴图上明亮区域出现点状图案要基于pdf和粗糙度采样mipmap

第二部分积分

\[ \underset{\Omega}{\int}f_r(p,w_i,w_o)n{\cdot}w_i dw_i \]

假设每个方向的入射辐射度都是白色的（\(L(p,x)=1.0\)），对于每种粗糙度和入射角的组合预计算的结果存储为一张2D的LUT，称为BRDF积分图。R和G分别保存的是菲涅尔响应的比例系数和偏差值，横坐标为\(n\cdot w_i\)，纵坐标为粗糙度。

下面对积分做亿点点变换： \[ \begin{aligned} &\ \underset{\Omega}{\int}f_r(p,w_i,w_o)n{\cdot}w_i dw_i \\ =&\ \underset{\Omega}{\int}f_r(p,w_i,w_o)\frac{F(w_o,h)}{F(w_o,h)}n{\cdot}w_i dw_i \\ =&\ \underset{\Omega}{\int}\frac{f_r(p,w_i,w_o)}{F(w_o,h)}{F(w_o,h)}n{\cdot}w_i dw_i \\ {F=F_0+(1-F_0)(1-w_o\cdot h)^5}\ {\Longrightarrow}&\ \underset{\Omega}{\int}\frac{f_r(p,w_i,w_o)}{F(w_o,h)}{(F_0+(1-F_0)(1-w_o\cdot h)^5)}n{\cdot}w_i dw_i \\ {\alpha=(1-w_o\cdot h)^5}\ {\Longrightarrow}&\ \underset{\Omega}{\int}\frac{f_r(p,w_i,w_o)}{F(w_o,h)}{(F_0+(1-F_0)\alpha)}n{\cdot}w_i dw_i \\ =&\ \underset{\Omega}{\int}\frac{f_r(p,w_i,w_o)}{F(w_o,h)}{(F_0*(1-\alpha)+\alpha)}n{\cdot}w_i dw_i \\ =&\ F_0\underset{\Omega}{\int}\frac{f_r(p,w_i,w_o)}{F(w_o,h)}(1-\alpha)n{\cdot}w_i dw_i + \underset{\Omega}{\int}\frac{f_r(p,w_i,w_o)}{F(w_o,h)}(\alpha)n{\cdot}w_i dw_i \\ {\alpha=(1-w_o\cdot h)^5}\ {\Longrightarrow}&\ F_0\underset{\Omega}{\int}\frac{f_r(p,w_i,w_o)}{F(w_o,h)}(1-(1-w_o\cdot h)^5)n{\cdot}w_i dw_i + \underset{\Omega}{\int}\frac{f_r(p,w_i,w_o)}{F(w_o,h)}(1-w_o\cdot h)^5n{\cdot}w_i dw_i \\ {f_r{'}(p,w_i,w_o)=\frac{DG}{4(w_o\cdot n)(w_i\cdot n)}}\ {\Longrightarrow}&\ F_0\underset{\Omega}{\int}{f_r^{'}(p,w_i,w_o)}(1-(1-w_o\cdot h)^5)n{\cdot}w_i dw_i + \underset{\Omega}{\int}{f_r^{'}(p,w_i,w_o)}(1-w_o\cdot h)^5n{\cdot}w_i dw_i \end{aligned} \] 这里又拆成了两个积分，分别是\(F_0\)的比例系数和偏差值，也就是将要保存在LUT中的值，这个LUT又叫做BRDF积分贴图，R和G分量分别就是这两个积分值。采样LUT需要的输入为NdotV和roughness

先生成一张2D纹理，利用帧缓冲把渲染一个屏幕四边形的结果保存到纹理上。

这里渲染过程用的着色器就执行了上面的积分，注意其中的几何函数跟直接光照里的有些不一样，直接光照中\(k_{direct}=(\alpha+1)^2/8\)，而IBL中\(k_{IBL}={\alpha}^2/2\)

着色器

float GeometrySchlickGGX(float NdotV, float roughness)
{
    float a = roughness;
    float k = (a * a) / 2.0; // 这里跟跟直接光照的不一样

    float nom   = NdotV;
    float denom = NdotV * (1.0 - k) + k;

    return nom / denom;
}
// ----------------------------------------------------------------------------
float GeometrySmith(vec3 N, vec3 V, vec3 L, float roughness)
{
    float NdotV = max(dot(N, V), 0.0);
    float NdotL = max(dot(N, L), 0.0);
    float ggx2 = GeometrySchlickGGX(NdotV, roughness);
    float ggx1 = GeometrySchlickGGX(NdotL, roughness);

    return ggx1 * ggx2;
}  

vec2 IntegrateBRDF(float NdotV, float roughness)
{
    vec3 V;
    V.x = sqrt(1.0 - NdotV*NdotV);
    V.y = 0.0;
    V.z = NdotV;

    float A = 0.0;
    float B = 0.0;

    vec3 N = vec3(0.0, 0.0, 1.0);

    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);
}
// ----------------------------------------------------------------------------
void main() 
{
    vec2 integratedBRDF = IntegrateBRDF(TexCoords.x, TexCoords.y);
    FragColor = integratedBRDF;
}

跟漫反射组合

着色器

[...]
uniform samplerCube prefilterMap; // 第一部分积分
uniform sampler2D   brdfLUT; // 第二部分积分
[...]
void main()
{
    [...]
    vec3 R = reflect(-V, N);   
    
    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; // 在与粗糙度对应的mipmap采样  
    vec2 envBRDF  = texture(brdfLUT, vec2(max(dot(N, V), 0.0), roughness)).rg;
    vec3 specular = prefilteredColor * (F * envBRDF.x + envBRDF.y);
    
    // 完整的IBL
    vec3 ambient = (kD * diffuse + specular) * ao; 
}

Vulxin's Blog