LearnOpenGL学习笔记（十四） - PBR与IBL

本文最后更新于 2025年1月7日 下午

终局之战。

PBR

理论

PBR，即基于物理的渲染（Physically Based Rendering），为一系列物理正确的渲染技术集合，以科学的方式描述材质与光照。

PBR遵循三大条件：

  • 基于微平面（Microfacet）的表面模型
  • 能量守恒
  • 应用基于物理的BRDF

微平面模型

根据微平面理论，任何表面在达到微观尺度后都可以用被称为微平面的微小镜面表示。表面越粗糙，组成该表面的微平面的指向分布变化越大。

image-20241224123832698

指向分布变化越大，入射光线在表面就越趋向向不同的方向发散，导致分布更宽泛的镜面反射。

我们使用粗糙度（Roughness）属性以统计学的方式描述表面的微平面分布混乱程度。我们可以基于一个平面的粗糙度来计算出众多微平面中，朝向方向沿着某个向量hh方向的比例。该向量hh为此表面的半程向量。

h=l+v∣∣l+v∣∣h=\frac{l+v}{||l+v||}，其中ll为入射光方向，vv为视线方向。微平面的朝向方向与半程向量的方向越是一致，镜面反射的效果就越是强烈越是锐利。

能量守恒

出射光线的能量永远不能超过入射光线的能量（发光面除外）

因此，粗糙度增大时，尽管镜面反射的区域增大了，但相应地，镜面反射的亮度也变小了。

能量守恒的一个关键方面在于对折射光和反射光做出了明确的区分。反射光会直接被表面反射，不进入表面内部；折射光则回进入表面内部，与表面内部的物质粒子进行碰撞，随后被全部吸收或部分吸收。

如果考虑为部分吸收，未被吸收的部分在经过内部碰撞、发散后会在略远的地方重新离开表面。一般情况下不考虑此部分，而次表面散射会考虑，使得诸如皮肤、蜡质等材质更为真实。

在通常情况下，根据能量守恒，我们只是单纯地用1减去镜面反射分量，就得到了漫反射分量。

然而，对于金属(Metallic)表面，其所有折射光均会被吸收。因此，在PBR管线中，金属表面不存在漫反射光，只存在镜面反射光。

反射光与折射光它们二者之间是互斥的关系。无论何种光线，其被材质表面所反射的能量将无法再被材质吸收。

1
2
float kS = calculateSpecularComponent(...); // 反射/镜面 部分
float kD = 1.0 - ks;                        // 折射/漫反射 部分

反射率方程

Lo(p,ωo)=∫Ωfr(p,ωi,ωo)Li(p,ωi)n ⋅ωidωiL_o(p,\omega_o) = \int_\Omega f_r(p,\omega_i,\omega_o)L_i(p,\omega_i)n\ ·\omega_id\omega_i

已知的：

  • pp：入射点坐标
  • ωo\omega_o：出射方向
  • ωi\omega_i：入射方向
  • nn：法向量

未知的：

  • 辐射率（Radiance）以LL表示，代表来自单一方向上的光线强度。亦可理解为一个拥有辐射通量Φ\Phi的光源在单位面积AA，单位立体角ω\omega上辐射的总能量。

L=d2ΦdAdωcos⁡θL = \frac{d^2\Phi}{dAd\omega\cos\theta}

在这里，我们不把光源视为点，而是平面，所以光源有单位面积的概念。

image-20241224134419511

辐射率受到**入射光线与平面法线夹角θ\theta的余弦值cos⁡θ\cos\theta（dot(lightDir, N）**的影响。光线垂直于平面时，强度最高。因此，反射率方程中的n ⋅ωin\ · \omega_i就代表着影响辐射率的参数。

而辐射率本身，在黎曼和代码中可以从光源处获得。也可以通过环境贴图测算所有入射方向上的辐射率。在代码中，辐射率用一个函数L表示。

当我们把立体角和面积视为无限小，此时辐射率表示单束光线穿过空间中一个点的通量。此时，我们可以计算出作用于单个点（即片段）上的单束光线的辐射率。

这种情况下，立体角变为方向向量，面A变为点P。由此，我们可以在着色器中使用辐射率计算单束光线对每个片段的作用了。

  • 辐射通量（Radiant Flux）以Φ\Phi表示，代表一个光源输出的能量，单位瓦特。

  • 辐射强度（Radiant Intensity）以II表示，代表单位球面上，一个光源向每单位立体角所投送的辐射通量。

I=dΦdωI = \frac{d\Phi}{d\omega}

立体角为投射到单位球体上的一个截面的面积。可以将其理解为带有体积的方向。

image-20241224133553743

简单理解：一束光，其通过ωi\omega_i方向入射到点p的能量大小可以以n⋅ωin ·\omega_i表示。入射后，反射光向四面八方反射出去。此时，点p的观察者向ωo\omega_o方向看去，发现ωo\omega_o方向的辐照度就是反射率方程右侧式子未经积分的结果。

**辐照度（Irradiance）**指所有投射到点p上的光线总和。

然而，实际情况下肯定不止一束光线投射到点p，入射光应当也是来自四面八方的。因此，我们需要统计来自于以点p为球心的半球领域（Hemisphere）内所有方向上的入射光。这个半球领域以Ω\Omega表示。

image-20241224141410017

我们将单束光线计算得到的辐照度，对半球领域内的所有入射方向进行积分，得到的就是真正的辐照度。

反射率方程没有解析解，所以要通过离散的方式，按一定的步长将其分散求解。然后根据步长大小将结果平均化。该方法被称为黎曼和（Reimann Sum）。

1
2
3
4
5
6
7
8
9
10
11
int steps = 100;
float sum = 0.0f;
vec3 P    = ...;
vec3 Wo   = ...;
vec3 N    = ...;
float dW  = 1.0f / steps;
for(int i = 0; i < steps; ++i)
{
    vec3 Wi = getNextIncomingLightDir(i);
    sum += Fr(p, Wi, Wo) * L(p, Wi) * dot(N, Wi) * dW;
}

BRDF

双向反射分布函数（Bidirectional Reflective Distribution Function, BRDF）用于基于表面材质属性对入射辐照率进行加权。

可以注意到，入射辐照率有两个加权参数，BRDF和入射方向与法向量夹角余弦值。理解了这个以后反射率方程就更好理解了。

BRDF接受入射方向ωi\omega_i，出射方向ωo\omega_o，平面法线nn和粗糙度α\alpha作为输入。它可以近似求出每束光线对一个给定了材质属性的平面上最终反射出的光线的贡献程度。

例如，如果一个平面绝对光滑，那么除入射光线输出1.0外，其余光线全部输出0.0。

在现代实时渲染管线中，我们一般使用Cook-Torrance BRDF模型。它兼有漫反射和镜面反射两部分：

fr=kdflambert+fcooktorrancef_r = k_df_{lambert}+f_{cook_torrance}

此处：

  • kdk_d表示入射光线中被折射部分能量的比率
  • ksk_s表示被反射部分的比率
  • flambertf_{lambert}表示漫反射部分，flambert=cπf_{lambert} = \frac{c}{\pi}，其中cc为表面颜色
  • fcook−torrancef_{cook-torrance}表示镜面反射部分。如下：

fcook−torrance=DFG4(ωo⋅n)(ωi⋅n)f_{cook-torrance=\frac{DFG}{4(\omega_o·n)(\omega_i·n)}}

其中，DFG为三个函数，分别为：

  • 法线分布函数（Normal Distribution Function）：估算在给定粗糙度的条件下，朝向方向与半程向量一致的微平面数量。是主要受粗糙度影响的值。
  • 几何函数（Geometry Function）：当一个平面比较粗糙时，其上的部分微表面可能挡住其余微表面，从而减少表面反射的光线，形成阴影。这类阴影被称为**“自成阴影”**。
  • 菲涅尔方程（Fresnel Equation）：描述不同表面角下表面反射光线的比率。

例如，一个不算光滑的表面，当视线与其逐渐平行时，镜面反射的强度逐渐增大。

法线分布函数

该函数（NDF）从统计学上近似表示与半程向量取向一致的微平面比率。其返回的值越小，则说明在此范围内与半程向量取向一致的微平面越少，该范围就越暗。

常用的NDF为Trowbridge-Reitz GGX。其GLSL实现如下：

1
2
3
4
5
6
7
8
9
10
float D_GGX_TR(vec3 N, vec3 H, float a)
{
    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;
}

NDFGGXTR(n,h,α)=α2π((α2−1)(n⋅h)2+1)2NDF_{GGXTR}(n,h,\alpha) = \frac{\alpha^2}{\pi((\alpha^2-1)(n · h)^2+1)^2}

其中，H为半程向量，α为粗糙度。

几何函数

该函数从统计学近似求得微平面自成阴影的比率。其接受粗糙度作为输入参数，粗糙度越高，自成阴影概率更高。

GSchlickGGX(n,v,k)=n⋅v(n⋅v)(1−k)+kG_{SchlickGGX}(n,v,k) = \frac{n · v}{(n·v)(1-k)+k}

其中，k为α的重映射。根据情况不同，二者关系也不同。

针对直接光照：kdirect=(α+1)28k_{direct} = \frac{(\alpha+1)^2}{8}

针对IBL光照：kIBL=α22k_{IBL}=\frac{\alpha^2}{2}

对于几何部分，有一部分的遮蔽是由视线本身无法到达被遮蔽的微表面导致（几何遮蔽，Geometry Obstruction），也有一部分由微表面被自成阴影纳入范围导致（几何阴影，Geometry Shadowing）。

为了将上面两个部分都纳入公式考虑范围，我们可以采用Smith’s method来计算几何函数：

G(n,v,l,k)=Gsub(n,v,k)Gsub(n,l,k)G(n,v,l,k) = G_{sub}(n,v,k)G_{sub}(n,l,k)

其中，GsubG_{sub}就是指SchlickGGX法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
float GeometrySchlickGGX(float NdotV, float k)
{
    float nom = NdotV;
    float denom = NdotV * (1.0 - k) + k;
    return nom / denom;
}

float GeometrySmith(vec3 N, vec3 V, vec3 L, float k)
{
    float NdotV = max(dot(N,V),0.0);
    float NdotL = max(dot(N,L),0.0);
    float ggx1 = GeometrySchlickGGX(NdotV,k);
    float ggx2 = GeometrySchlickGGX(NdotL,k);
}

菲涅尔方程

此方程描述被反射光线对比被折射光线的比率，其返回值随观察角度不同而不同。

垂直观察时，任何材质表面都有基础反射率（Base Reflectivity）。但如果以近乎90度的角度观察表面，反光就明显的多。

我们使用Fresnel-Schlick近似法求取菲涅尔方程近似解：

FSchlick(h,v,F0)=F0+(1−F0)(1−(h⋅v))5F_{Schlick}(h,v,F_0)=F_0+(1-F_0)(1-(h·v))^5

其中，F0F_0表示平面基础反射率，由**折射指数（Indices of Refraction, IOR）**计算得出。

该近似法仅对非金属表面有定义。对于金属，需要使用不同结构的菲涅尔，这十分不方便。因此，我们与计算出F0F_0（即垂直观察时的基础反射率）结果，然后根据观察角进行插值，就可以对任何材质使用同一公式了。

F0F_0以RGB三原色表示，因为金属材质表面的反射光有时带有色彩。所以我们可以看到，非金属表面的F0三个分量相等，而金属则不相等。

**金属度（Metalness）**参数用于描述一个材质表面是金属还是非金属的。通过控制金属度，可以细微调整材质表面的视觉效果。

Fresnel-Schlick近似的GLSL实现如下：

1
2
3
4
5
6
7
8
vec3 F0 = vec(0.04); // 0.04为常见电介质的基础反射率平均值
F0 = mix(F0, surfaceColor,rgb, metalness); // 金属度越大，基础反射率越接近表面颜色。

vec3 fresnelSchlick(float cosTheta, vec3 F0)
{
    // 其中cosTheta为h dot v
    return F0 + (1.0 - F0) * pow(1.0 - cosTheta, 5.0);
}

Cook-Torrance反射率方程

Lo(p,ωo)=∫Ω(Kdcπ+ksDFG4(ωo⋅n)(ωi⋅n)Li(p,ωi)dωi)L_o(p,\omega_o) = \int_\Omega(K_d\frac{c}{\pi}+k_s\frac{DFG}{4(\omega_o·n)(\omega_i·n)}L_i(p,\omega_i)d\omega_i)

PBR材质

一个典型的PBR材质应当包含下列贴图：

  • 反照率（Albedo）：指定金属材质表面颜色（即基础反照率F0），类似于漫反射纹理（仅包含表面颜色）。
  • 法线：法线贴图，用于制造凹凸不平的视觉效果。
  • 金属度：逐纹素指定金属质地。
  • 粗糙度：逐纹素指定材质粗糙程度。有时基于One Minus变为**光滑度（Smoothness）**贴图。
  • AO：环境光遮蔽贴图。

光照

点光源在每个方向的辐射通量均相同。因此，在不考虑衰减的情况下，其辐射通量可用一个三维向量表示。

然而，真正的辐射率计算必定需要考虑衰减，而辐射率L(p,ω)L(p,\omega)也确实需要接受一个位置pp作为输入。

由此，我们得到以下过程：

对于直接点光源，辐射率函数L进行以下操作：

  1. 获取光源颜色值
  2. 将颜色值按光源和某点p的距离衰减
  3. 按照视线与法线夹角余弦值，即n⋅ωin·\omega_i缩放。其中ωi\omega_i就是光源和pp点之间的方向向量。

代码如下：

1
2
3
4
5
vec3 lightColor = vec3(23.47, 21.31, 20.79);
vec3 wi = normalize(lightPos - fragPos);
float cosTheta = max(dot(N, Wi), 0.0);
float attenuation = calculateAttenuation(fragPos, lightPos);
float radiance = lightColor * attenuation * cosTheta;

对于其他光源，平行光不需要计算衰减；聚光灯的辐射强度需要根据聚光灯方向向量进行缩放。

反射率方程中，计算总辐射率需要对点p所在的半球领域进行积分。但实际上，程序中的光源数量是有限的。我们只需要迭代光源数量的次数，便可以计算出所有光源对某点的辐射率贡献，直接相加就能得到总辐射率。

考虑间接光源就需要用到IBL和积分计算了。

PBR着色器示例

片段着色器

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
#version 330 core
// I/O相关
out vec4 FragColor;
in vec2 TexCoords;
in vec3 WorldPos;
in vec3 Normal;

// 相机相关
uniform vec3 camPos;

// 材质属性相关
uniform sampler2D albedoMap;
uniform sampler2D metallicMap;
uniform sampler2D roughnessMap;
uniform sampler2D aoMap;

// 光照相关
const int LIGHT_COUNT = 4;
uniform vec3 lightPositions[LIGHT_COUNT];
uniform vec3 lightColors[LIGHT_COUNT];

// 辅助常量
const float PI = 3.14159265359;

// 函数签名
vec3 fresnelSchlick(float cosTheta, vec3 F0);
float DistributionGGX(vec3 N, vec3 H, float roughness);
float GeometrySchlickGGX(float NdotV, float roughtness);
float GeometrySmithGGX(vec3 N, vec3 V, vec3 L, float roughness);

void main(){
    vec3 albedo     = pow(texture(albedoMap, TexCoords).rgb, 2.2); // 漫反射贴图通常为Gamma空间下制作，需要先转换到线性
    float metallic  = texture(metallicMap, TexCoords).r;
    float roughness = texture(roughnessMap, TexCoords).r; // 金属度和粗糙度贴图一般位于线性空间
    float ao        = texture(aoMap, TexCoords).r; // 一般而言AO贴图也需要转换到线性，这里图方便没转

    vec3 N = normalize(Normal); // 法线
    vec3 V = normalize(camPos - WorldPos); // 视线方向，即ViewDir
    // -----直接光照部分-----
    vec3 Lo = vec3(0.0); // 最终计算结果，即总辐射率
    for(int i=0; i<LIGHT_COUNT; ++i){
        // 计算Li项
        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; // 单个光源在入射方向的辐射率
        // 计算BRDF项
        // 1. 菲涅尔
        vec3 F0 = vec3(0.04);
        F0 = mix(F0, albedo, metallic);
        vec3 F = fresnelSchlick(max(dot(H, V),0.0), F0);
        // 2. 法线分布
        float NDF = DistributionGGX(N, H, roughness);
        // 3. 几何
        float G = GeometrySmith(N, V, L, roughness);
        // 4. Cook-Torrance高光项
        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;
        // 5. Cook-Torrance漫反射项
        vec3 diffuse = albedo / PI;
        // 6. 总结
        vec3 kS = F; // 菲涅尔方程表示反射光线在所有光线内的占比
        vec3 KD = 1.0 - kS; // 所谓漫反射，就是折射光线在表面内部多次弹射后再次（在入射点附近）弹出表面的光，根据能量守恒可以求得
        kD *= 1.0 - metallic; // 金属不存在折射，只有反射
        float NdotL = max(dot(N, L), 0.0);
        Lo += (kD * diffuse + specular) * radiance * NdotL;
    }
    // 计算环境光照
    vec3 ambient = vec3(0.03) * albedo * ao; // AO（环境光遮蔽）作用于环境光
    vec3 color = ambient + Lo;
    // Reinhard色调映射，从LDR映射到HDR
    color = color / (color + vec3(1.0));
    // Gamma矫正
    color = pow(color, vec3(1.0 / 2.2));
}

// -----BRDF部分-----
// ---菲涅尔方程---
// 计算菲涅尔方程的Schlick近似
// cosTheta一般为H dot V，表示法线和视线的夹角。夹角越小，此值越大，菲涅尔项越弱
// F0为垂直视角反射率，材质独有属性，一般绝缘体为0.04，金属需要根据情况插值。则有F0 = mix(0.04, albedo, metallic)
vec3 fresnelSchlick(float cosTheta, vec3 F0){
    // clamp用于将第一参数的值超出第二-第三参数区间的部分设置为边缘值，此处用于防止亮点或暗点
    return F0 + (1.0 - F0) * pow(clamp(1.0 - cosTheta, 0.0, 1.0), 5.0);
}
// ---法线分布函数---
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;
}
// ---几何函数---
float GeometrySchlickGGX(float NdotV, float roughtness){
    float r = (roughness + 1.0);
    float k = (r*r) / 8.0; // 此处未用IBL，否则k计算方式将有所不同
    float nom = NdotV;
    float denom = NdotV * (1.0 - k) + k;
    return nom / denom;
}
float GeometrySmithGGX(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, roughtness);
    return ggx1 * ggx2;
}

注意：

  • 此着色器中，我们在线性空间计算光照。
  • 线性空间钟，逆平方衰减量更加物理准确。
  • Cook-Torrance BRDF的镜面反射分量系数就是菲涅尔项，它代表反射光线占总光线的比例。因此，最终计算时无需重复乘以kS。
  • 由于PBR基于物理，所以所有光照输入都与真实的物理值相仿，很容易使得最终计算的Lo超过1.0而被截断，所以需要将颜色值进行色调映射到HDR范围
  • PBR要求所有数值均处于线性空间，所以要在颜色计算完毕后进行Gamma矫正计算。

IBL

漫反射辐照度

基于图像的光照（Image based lighting， IBL）将环境立方体贴图上的每个像素视为光源，并在渲染方程中直接使用这些光源。相比于直接光源，IBL能捕捉近乎全部的环境光照，是全局光照的粗略近似。因此，在PBR中使用IBL能大大提高光照真实度。

前文的反射率方程：

Lo(p,ωo)=∫Ω(kdcπ+ksDFG4(ωo⋅n)(ωi⋅n))Li(p,ωi)n⋅ωidωiL_o(p,\omega_o) = \int_\Omega(k_d\frac{c}{\pi}+k_s\frac{DFG}{4(\omega_o·n)(\omega_i·n)})L_i(p,\omega_i)n·\omega_id\omega_i

我们的目标是计算半球Ω\Omega上所有入射光方向ωi\omega_i的积分。在使用直接光照时，因为仅存在N个直接光源，所以对于着色点pp，只存在N个入射方向，可以把积分化为N次循环。

然而，使用IBL时，来自周围环境的每个方向的入射光都有可能具有辐射度。因此，我们必须实现下面两个需求：

  • 给定任何入射向量ωi\omega_i，都可以获取来自该方向上场景的辐射度。
  • 解决积分的方法必须快速、实时。

对于要求1，我们可以将环境贴图作为立方体贴图进行采样，即vec3 radiance = texture(_cubemapEnvironment, w_i).rgb

对于要求2，我们可以预先计算好所有需要的数值，将它们存储在纹理中，需要时进行采样即可。

我们知道，反射率方程的积分内多项式由漫反射项和镜面反射项构成。将二者拆开：

Lo(p,ωo)=∫Ω(kdcπ)Li(p,ωi)n⋅ωidωi+∫Ω(ksDFG4(ωo⋅n)(ωi⋅n))Li(p,ωi)n⋅ωidωiL_o(p,\omega_o) = \int_\Omega(k_d\frac{c}{\pi})L_i(p,\omega_i)n·\omega_id\omega_i + \int_\Omega(k_s\frac{DFG}{4(\omega_o·n)(\omega_i·n)})L_i(p,\omega_i)n·\omega_id\omega_i

对于漫反射项，我们知道kd、c、πk_d、c、\pi都是常数项。所以可以将其移出积分，得到新BRDF漫反射项kdcπ∫ΩLi(p,ωi)n⋅ωidωik_d\frac{c}{\pi}\int_\Omega L_i(p,\omega_i)n·\omega_id\omega_i

接下来我们着重来看漫反射部分。

从上面的漫反射项可以看出，该积分只依赖于入射方向。因此，我们进行的预处理就是：给定一个立方体贴图，计算其在每个采样方向，即纹素上的漫反射积分结果。该预处理通过卷积实现。

卷积的特性是，对数据集中的一个条目做一些计算时，要考虑到数据集中的所有其他条目。这里的数据集就是场景的辐射度或环境贴图。因此，要对立方体贴图中的每个采样方向做计算，我们都会考虑半球Ω上的所有其他采样方向。

然而，半球内存在的方向实在太多，如果一个个考虑必将消耗非常多的性能。因此，我们对半球上的方向进行离散采样，并对其辐射度取平均值，来计算每个输出采样方向的积分。

卷积后，我们得到的输出采样方向积分究竟是什么？

预计算的立方体贴图，在每个采样方向ωo\omega_o上存储其积分结果，可以理解为场景中所有能够击中以ωo\omega_o为法线的表面的间接漫反射光的预计算总和。是对全局光照的近似。辐照度指投射到某点的所有光线的总和。所以完成积分计算后的贴图才叫辐照度贴图。

完成积分后，我们会得到一张新的立方体贴图。它有点类似于对原本的立方体环境贴图进行模糊卷积处理。这张新贴图被称为辐照度贴图。通过该贴图，给定一个入射方向，我们可以直接采样到该方向的预计算辐照度。

然而，这种方式默认材质表面始终处于立方体贴图的中心，也就是说我们默认立方体贴图是作为天空盒的。如果一个场景很复杂时，就会显得不太真实（例如，室内的镜面反射和室外的应当不同，但在我们的方法里却相同）。

为了解决这一问题，可以使用反射探针技术。

环境立方体贴图和其辐照度贴图

左侧是环境立方体贴图，可以被认为是辐射率贴图；右侧为辐照度贴图。

.hdr

先前我们提到，PBR工作流中的一切数值都不应当被局限在LDR范围内。因此，直接点光源的颜色也必定会大于一。在IBL中，类似地，我们通过辐照度贴图采样到的值也必定在HDR范围而非LDR范围内。也就是说，我们的环境立方体贴图必须为HDR贴图。

辐射率贴图又被称为辐射度文件，扩展名为**.hdr**。该文件存储了一张完整的环境辐射度贴图（以等距柱状投影贴图的形式，而非六面立方体贴图），六个面的颜色数据均为浮点数，允许0-1范围外的值。它使用8 bit存储每个通道，并用alpha通道存放指数。

使用stb_image库可以方便地将.hdr文件加载为浮点数组。如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#include "stb_image.h"
// ...
stbi_set_flip_vertically_on_load(true); // 由于OpenGL的默认纹理倒转，需要翻转预处理
int width, height, nrComponents; // 长、宽、通道
float* data = stbi_load("name.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_TEXTRUE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTRUE_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;
}

上述代码执行后，将自动把.hdr文件存储到浮点数组中，保持HDR范围。

随后，我们尝试将等距圆柱投影贴图转换为立方体贴图，以便更方便地采样。

1
2
3
4
5
6
7
8
9
10
// 顶点着色器，原样渲染立方体
#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);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// 片段着色器
#version 330 core
out vec4 FragColor;
in vec3 localPos;
uniform sampler2D equirectangularMap; // 等距柱状纹理贴图
const vec3 invAtan = vec2(0.1591, 0.3183);
// 通过变换UV，使得采样方向转换到等距柱状空间
vec2 SampleSphericalMap(vec3 v){
    vec2 uv = vec2(atan(v.z, v.x), asin(v.y));
    uv *= invAtan;
    uv += 0.5;
    return uv;
}
void main(){
    vec2 uv = SampleSphericalMap(normalize(localPos));
    vec3 color = texture(equirectangularMap, uv).rgb;
    FragColor = vec4(color, 1.0);
}

我们知道，在片段着色器运行时，不处于[0,1]范围内的uv，以及我未通过深度测试的片段将被剔除。也就是说，只有我们在视口内看到的东西才能被渲染。

我们通过SampleSphericalMap变换UV，也只是对视口内存在的UV进行变换。然而，要生成立方体贴图的六个面，我们就必定需要把六个面都渲染到纹理，然后将六个纹理合成为立方体贴图。

当然，我们也可以通过预计算将.hdr转换为cubemap，例如使用Blender的烘焙等。

随后，我们对同一立方体渲染六次，每次面对立方体的一个面，用FBO记录结果：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// 创建FBO、RBO
unsigned int captureFBO, captureRBO;
glGenFrameBuffers(1, &captureFBO);
glGenRenderBuffers(1, &captureRBO); // 帧缓冲需要绑定附件才能运作，RBO用于深度缓冲，纹理附件用于颜色缓冲
glBindFrameBuffer(GL_FRAMEBUFFER, captureFBO);
glBindRenderBuffer(GL_RENDERBUFFER, captureRBO);
glRenderbufferStorage(GL_RENDERBUFFER, GL_DEPTH_COMPONENT24, 512, 512); // 尺寸为512*512，代表cubemap一个面的分辨率
glFramebufferRenderBuffer(GL_FRAMEBUFFER, GL_DEPTH_ATTACHMENT, GL_RENDERBUFFER, captureRBO);  // 将RBO附加到FBO
// 生成立方体贴图
unsigned int envCubemap;
glGenTextures(1, &envCubemap);
glBindTexture(GL_TEXTURE_CUBE_MAP, envCubemap);
for(unsigned int i=0;i<6;i++){
    glTexImage2D(GL_TEXTURE_CUBE_MAP_POSITIVE_X+i, 0, GL_RGB16F, 512, 512, 0, GL_RGB, GL_FLOAT, nullptr); //为cubemap的每个面分配内存，故最后参数为nullptr
}
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

随后，在渲染循环中，将等距柱状2D纹理捕捉到立方体贴图面上。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
// FOV90度，分别LookAt六个面
glm::mat4 captureProjection = glm:;perspective(glm::radians(90.0f),1.0f,0.1f,10.0f);
glm::mat4 captureViews[] = {
   glm::lookAt(glm::vec3(0.0f, 0.0f, 0.0f), glm::vec3( 1.0f,  0.0f,  0.0f), glm::vec3(0.0f, -1.0f,  0.0f)),
   glm::lookAt(glm::vec3(0.0f, 0.0f, 0.0f), glm::vec3(-1.0f,  0.0f,  0.0f), glm::vec3(0.0f, -1.0f,  0.0f)),
   glm::lookAt(glm::vec3(0.0f, 0.0f, 0.0f), glm::vec3( 0.0f,  1.0f,  0.0f), glm::vec3(0.0f,  0.0f,  1.0f)),
   glm::lookAt(glm::vec3(0.0f, 0.0f, 0.0f), glm::vec3( 0.0f, -1.0f,  0.0f), glm::vec3(0.0f,  0.0f, -1.0f)),
   glm::lookAt(glm::vec3(0.0f, 0.0f, 0.0f), glm::vec3( 0.0f,  0.0f,  1.0f), glm::vec3(0.0f, -1.0f,  0.0f)),
   glm::lookAt(glm::vec3(0.0f, 0.0f, 0.0f), glm::vec3( 0.0f,  0.0f, -1.0f), glm::vec3(0.0f, -1.0f,  0.0f))
};
equirectangularToCubemapShader.use();
equirectangularToCubemapShader.setInt("equirectangularMap",0);
equirectangularToCubemapShader.setMat4("projection", captureProjection);
glActiveTexture(GL_TEXTURE0);
glBindTexture(GL_TEXTURE_2D, hdrTexture);
glViewport(0,0,512,512); //与纹理附件、RBO分辨率对应
glBindFrameBuffer(GL_FRAMEBUFFER, captureFBO);
// 对于每个面：
for(unsigned int i=0;i<6;i++){
    equirectangularToCubemapShader.setMat4("view",captureViews[i]);
    // 将已绑定CUBEMAP（此时六个面分配了内存但没有数据）的每个面作为纹理附件，从而渲染到纹理
    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_CUBE_MAP_POSITIVE_X + i, envCubemap, 0);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTJ_BUFFER_BIT);
    // 此时渲染的立方体面将被绘制到CUBEMAP的某个面上
    renderCube();
}
glBindFramebuffer(GL_FRAMEBUFFER, 0);

.随后，使用天空盒着色器渲染完成的Cubemap：

1
2
3
4
5
6
7
8
9
10
11
12
// 天空盒顶点着色器
#version 330 core
layout(location = 0) in vec3 aPos;
uniform mat4 projection;
uniform mat4 view;
out vec3 localPos;
void main(){
    localPos = aPos;
    mat4 rotView = mat4(mat3(view)); //将平移变换移除
    vec4 clipPos = projection * rotView * vec4(localPos, 1.0);
    gl_Position = clipPos.xyww; //确保深度值始终为1
}
1
2
3
4
5
6
7
8
9
10
11
12
#version 330 core
out vec4 FragColor;
in vec3 localPos;
uniform samplerCube environmentMap;
void main(){
    vec3 envColor = texture(environmentMap, localPos).rgb;
    // HDR色调映射，转换到LDR空间
    envColor = envColor / (envColor + vec3(1.0));
    // .hdr为线性空间，需要进行伽马矫正
    envColor = pow(envColor, vec3(1.0/2.2));
    FragColor = vec4(envColor, 1.0);
}

Cubemap卷积

我们先前提到，对辐射度贴图（即环境立方体贴图或.hdr文件）进行积分，会得到辐照度贴图。对辐照度贴图使用一个方向向量采样，得到的值就是以该方向向量为法线的材质表面接收到的总光线强度的近似。

那么，这个卷积操作具体该如何实现呢？

首先要明确，均匀采样辐射度并离散计算辐照度是可行的，但它的性能消耗依然比较大。所以我们才要对辐射度贴图进行与计算。

与SSAO类似，我们随机均匀采样的范围为一个半球。为了生成辐照度贴图，我们需要将环境光照求卷积，转换为立方体贴图。假设对于每个片段，表面的半球朝向法向量 N ，对立方体贴图进行卷积等于计算朝向 N的半球 Ω 中每个方向 wi 的总平均辐射率。

半球到底是什么？

我们知道，.hdr文件存储的辐射度贴图以等距柱状投影形式呈现。等距柱状投影和立方体贴图都是对球面贴图投影，用于把球面展开为单张贴图或六张贴图。

因此，真正的辐射度贴图应当是球面的形式。

其次，我们知道每个半球都有一个朝向向量。把辐射度贴图球面一切为二，使得其中一个半球的朝向向量为wi。我们便得到了wi所代表的半球。

具体卷积步骤为：

  1. 对于立方体贴图的每个纹素，在纹素所代表的方向（可以采样到该纹素的方向向量的反方向）的半球Ω\Omega内生成固定数量的采样向量，然后让这些向量对原立方体贴图进行采样。
  2. 对采样的结果取平均值，赋值给新立方体贴图的同一纹素位置。

观察BRDF漫反射分量：

kdcπ∫ΩLi(p,ωi)n⋅ωidωik_d\frac{c}{\pi}\int_\Omega L_i(p,\omega_i)n·\omega_id\omega_i

该积分对立体角进行积分。而立体角极难处理，所以我们用球坐标代替。

image-20241231121744834

对于围绕大圆的Head角ϕ\phi，采样范围为[0,2π][0,2\pi]。对于从半球顶点出发的Pitch角θ\theta，采样范围为[0,12π][0,\frac{1}{2}\pi]。

反射积分方程的漫反射项更新为下：

Lo(p,ϕo,θo)=kdcπ∫ϕ=02π∫θ=012πLi(p,ϕi,θi)cos(θ)sin(θ)dϕdθL_o(p,\phi_o,\theta_o)=k_d\frac{c}{\pi}\int_{\phi=0}^{2\pi}\int_{\theta=0}^{\frac{1}{2}\pi}L_i(p,\phi_i,\theta_i)cos(\theta)sin(\theta)d\phi d\theta

为了求解积分，我们需要在半球内采集固定数量的离散样本并对结果取平均值。分别给每个球坐标轴指定样本数量n1、n2n1、n2以求积分黎曼和。上述积分式转化为：

Lo(p,ϕo,θo)=kdcπ1n1n2∑ϕ=0n1∑θ=0n2Li(p,ϕi,θi)cos(θ)sin(θ)dϕdθL_o(p,\phi_o,\theta_o) = k_d\frac{c}{\pi}\frac{1}{n1n2}\sum_{\phi=0}^{n1}\sum_{\theta=0}^{n2}L_i(p,\phi_i,\theta_i)cos(\theta)sin(\theta)d\phi d\theta

当我们离散地对两个球坐标轴采样时，每个采样近似代表了半球上的一小块区域。

下面的片段着色器对已转换为立方体贴图的.hdr文件进行卷积，从而得到辐照度立方体贴图。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#version 330 core
out vec4 FragColor;
in vec3 localPos;
uniform samplerCube environmentMap; // 转换为为立方体贴图的.hdr文件
const float PI = 3.14159265359;
void main(){
    vec3 normal = normalize(localPos);
    vec3 irradiance = vec3(0.0);
    //---卷积操作---
    // 构建以法线为中心的局部坐标系，类似于切线空间，便于采样
    vec3 up = vec3(0.0,1.0,0.0);
    vec3 right = cross(up, normal);
    up = cross(normal, 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;
            // 对辐射度贴图采样并累加
            irradiance += texture(environmentMap, sampleVec).rgb * cos(theta) * sin(theta);
            nrSamples++;
        }
    }
    irradiance = PI * irradiance * (1.0 / float(nrSamples));
    FragColor = vec4(irradiance, 1.0);
}

在辐射度立方体贴图生成后，便可以调用此着色器生成辐照度立方体贴图了。配套的C++代码如下：

1
2
3
4
5
6
7
8
9
10
11
unsigned int irradianceMap;
glGenTextures(1, &irradianceMap);
glBindTexture(GL_TEXTURE_CUBE_MAP, irradianceMap);
for(unsigned int i=0;i<6;i++){
    glTexImage2D(GL_TEXTURE_CUBE_MAP_POSITIVE_X+i,0,GL_RGB16F,32,32,0,GL_RGB,GL_FLOAT,nullptr); // 为每个面分配内存
    glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
	glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
	glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
	glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
}

辐照度图不需要太高的分辨率，因为它的高频细节较少。这里用32*32。

需要注意的是，这里我们不使用先前预计算的辐射度贴图，而是重新将立方体渲染六次并渲染到CubeMap纹理，再用卷积着色器对其操作。

1
2
3
glBindFramebuffer(GL_FRAMEBUFFER, captureFBO);
glBindRenderbuffer(GL_RENDERBUFFER, captureRBO);
glRenderbufferStorage(GL_RENDERBUFFER, GL_DEPTH_COMPONENT24, 32,32);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
irradianceShader.use();
irradianceShader.setInt("environmentMap",0);
irradianceShader.setMat4("projection",captureProjection);
glActiveTexture(GL_TEXTURE0);
glBindTexture(GL_TEXTURE_CUBE_MAP, envCubemap);
glViewport(0,0,32,32);
glBindFramebuffer(GL_FRAMEBUFFER, captureFBO);
for(unsigned int i=0;i<6;i++){
    irradianceShader.setMat4("view",captureViews[i]);
    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_CUBE_MAP_POSITIVE_X+i, irradianceMap, 0);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    renderCube();
}
glBindFramebuffer(GL_FRAMEBUFFER, 0);

应用到PBR

辐照度贴图表示所有周围间接光累计的反射率的漫反射部分的积分。先前的PBR代码中，我们设置vec3 ambient = vec3(0.03)，现在我们要用辐照度贴图的采样值代替：

1
2
3
4
5
6
vec3 kS = fresnelSchlick(max(dot(N, V), 0.0), F0);
vec3 kD = 1.0 - kS;
// 此处N就是法线
vec3 irradiance = texture(irradianceMap, N).rgb;
vec3 diffuse    = irradiance * albedo;
vec3 ambient    = (kD * diffuse) * ao;

我们希望粗糙度高的表面菲涅尔反射弱。然而现在的方法并没有把菲涅尔项与粗糙度结合。通过改进Schlick菲涅尔近似，可以解决这一问题：

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

最终的环境光代码为：

1
2
3
4
5
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;

镜面反射IBL

理论

反射方程的镜面反射项如下：

∫Ω(ksDFG4(ωo⋅n)(ωi⋅n))Li(p,ωi)n⋅ωidωi=∫Ωfr(p,ωi,ωo)Li(p,ωi)n⋅ωidωi\int_\Omega(k_s\frac{DFG}{4(\omega_o·n)(\omega_i·n)})L_i(p,\omega_i)n·\omega_id\omega_i = \int_\Omega f_r(p,\omega_i,\omega_o)L_i(p,\omega_i)n·\omega_id\omega_i

对于镜面反射项，其中DFGDFG三项均与视线方向vv有关。因此，如果要进行穷举遍历，存在的可能性将极其庞大，无法实时计算。为了解决这一问题，我们可以预先计算镜面部分的卷积，在后续着色时采样即可。该方案被称为分割求和近似法。

在分割求和近似法中，我们将方程的镜面项分割为两个独立部分，然后对其分别求卷积，并在着色器中求和，得到真正的镜面辐照度量。

我们将镜面反射项如此拆分：

Lo(p,ωo)=∫ΩLi(p,ωi)dωi ∗ ∫Ωfr(p,ωi,ωo)n⋅ωidωiL_o(p,\omega_o)=\int_\Omega L_i(p,\omega_i)d\omega_i\ *\ \int_\Omega f_r(p,\omega_i,\omega_o)n·\omega_id\omega_i

因为有两个部分，所以我们要对这两个部分分别进行卷积，生成贴图。

第一部分为预滤波环境贴图（PrefilterMap），它是积分中与漫反射项共用的一个系数。它类似于辐照度图，是考虑了粗糙度的环境卷积贴图。之所以考虑粗糙度，是因为随着粗糙度的增加，镜面反射会更模糊。

而粗糙度越大，贴图越模糊，需要的精度就越小。因此，高粗糙度的预滤波环境贴图的分辨率可以变得较小。因此，我们可以把不同级别的粗糙度指的预卷积结果存储在Mipmap中。如图：

image-20250102145713279

第二部分为BRDF积分贴图（BRDFLUT），它是一张查找贴图（Look Up Texture, LUT），存储BRDF对每个粗糙度和入射角的响应结果。贴图的横坐标是BRDF输入，粗糙度为纵坐标；R通道存储菲涅尔响应系数，G通道存储菲涅尔响应偏差值。

image-20250103141657325

通过结合采样预滤波环境贴图和BRDFLUT我们就可以获得镜面积分的结果：

1
2
3
4
float lod = getMipLevelFromRoughness(roughness); // 根据粗糙度，获取指定的预滤波环境贴图Mipmap层级
vec3 prefilteredColor = textureCubeLod(PrefilteredEnvMap, refVec, lod); // 从指定的Mipmap层级中取出对应粗糙度的预滤波环境贴图
vec2 envBRDF = texture2D(BRDFIntegrationMap, vec2(NdotV, roughness)).xy; // 使用(NdotV, roughness)作为UV，采样BRDFLUT，获取BRDF项的数值
vec3 indirectSpecular = prefilteredColor * (F * envBRDF.x + envBRDF.y); // 得到间接光照的镜面反射项数值。其中，envBRDF的R通道，即.x，为菲涅尔响应系数，应当与材质基础反射率F相乘；G通道，即.y，为菲涅尔响应偏差值

操作

预滤波HDR环境贴图

先前我们已经明确，预滤波环境贴图类似于漫反射IBL中的辐照度图，但引入了粗糙度作为变量（粗糙度越大，“辐照度图”就越“模糊”），需要根据不同的粗糙度生成不同精度/模糊程度的辐照度图，并存储在Mipmap中。

为了生成Mipmap，我们首先需要在为预滤波环境立方体贴图分配内存后，通过glGenerateMipmap为其Mipmap分配内存：

1
2
3
4
5
6
7
8
9
10
11
12
unsigned int prefilterMap;
glGenTexture(1, &prefilterMap);
glBindTexture(GL_TEXTURE_CUBE_MAP, prefilterMap);
for(unsigned int i=0;i<6;i++){
    glTexImage2D(GL_TEXTURE_CUBE_MAP_POSITIVE_X+i, 0, GL_RGB16F, 128, 128, 0, GL_RGB, GL_FLOAT, nullptr);
}
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR);  // 注意这里
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
glGenerateMipmap(GL_TEXTURE_CUBE_MAP); // 为Mipmap分配内存

漫反射IBL中，我们使用标准半球均匀采样入射向量，因为对于漫反射，其入射向量的方向均匀且随机。然而，对于镜面反射来说，根据粗糙度不同，入射向量的分布也有所不同。如下：

image-20250106123555674

可能出射的反射光构成的形状被称为镜面波瓣（Specular Lobe）。随着粗糙度增加，镜面波瓣的大小增加；随着入射方向改变，镜面波瓣的形状也会发生变化。我们在卷积进行采样时，应当选取镜面波瓣内的向量，此过程被称为重要性采样。

镜面波瓣始终指向微表面的半程向量。

蒙特卡洛积分

蒙特卡洛积分建立在大数定律的基础上，不为近乎无限的样本值x求积分，而是简单地从总体中随机挑选样本N生成采样值并求平均。样本数N越大，结果越接近积分的精确结果。公式表示：

O=∫abf(x)dx=1N∑i=0N−1f(x)pdf(x)O = \int^b_af(x)dx=\frac{1}{N}\sum^{N-1}_{i=0}\frac{f(x)}{pdf(x)}

我们从[a,b]范围内采样N个样本，将它们相加并除样本数以取平均。公式中，pdfpdf为概率密度函数，表示特定样本在整个样本集上发生的概率。

当每次取到的样本均服从pdf(x)，则此蒙特卡洛估算无偏（Unbias），即，随着样本数的增加，蒙特卡洛积分结果必定逐渐收敛到精确解。

当生成的样本并不完全服从pdf(x)，而是有特定的倾向，则此蒙特卡洛估算有偏。有偏蒙特卡洛积分有更快的收敛速度，对性能敏感的应用程序来说十分合适。只要偏差的倾向较为合理，最终收敛的结果也不会有太多偏差。

使用低差异序列（Low-Discrepancy Sequence）进行蒙特卡洛积分可以进一步提升收敛速度。低差异序列生成的随机样本相较于完全随机样本，具有更加均匀的分布。此过程被称为拟蒙特卡洛积分（Quasi-Monte Carlo）。

image-20250106130826690

总结：蒙特卡洛积分时一种以高效的离散方式求积分的直观方法。为了进一步提升蒙特卡洛积分的计算速度，我们采用下面两种方法：1. 使用有偏蒙特卡洛积分，将采样范围按照实际情况进行限制。2. 使用低差异序列生成采样样本。

在镜面IBL引入蒙特卡洛的关键在于，借助有偏蒙特卡洛的范围限制，将材质表面的反射特性纳入采样考虑范围内，即重要性采样。

我们接下来将使用重要性采样来预计算间接镜面反射项。

低差异序列

我们将使用Hammersley序列作为拟蒙特卡洛过程的低差异序列。该序列通过把十进制数字的二进制表示镜像翻转到小数点右边得到。代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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
}
// N为总样本数，i为样本索引
vec2 Hammersley(uint i, uint N)
{
    return vec2(float(i)/float(N), RadicalInverse_VdC(i));
}
GGX重要性采样

首先，生成随机低差异序列：

1
2
3
4
const uint SAMPLE_COUNT = 4096u;
for(uint i = 0u; i<SAMPLE_COUNTL; i++){
    vec2 Xi = Hammersley(i, SAMPLE_COUNT);
}

随后，进行重要性采样。先对采样向量进行偏移，使其朝向特定粗糙度的镜面波瓣方向。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
vec3 ImportanceSampleGGX(vec2 Xi, vec3 N, float roughness) {
    // 经验公式，平方粗糙度视觉效果更佳
    float a = roughness * roughness;
    // 生成半程向量 H 的方位角 phi。
    // Xi.x 用于随机生成 phi 的均匀分布。
    float phi = 2.0 * PI * Xi.x;
    // 根据 GGX 分布公式生成极角的余弦值 cosTheta。
    // 使用 GGX 分布的逆变换采样公式：
    // cosTheta = sqrt((1 - Xi.y) / (1 + (a² - 1) * Xi.y))
    // Xi.y 为均匀分布的随机数，用于控制分布方向。
    float cosTheta = sqrt((1.0 - Xi.y) / (1.0 + (a * a - 1.0) * Xi.y));
    // 计算 sinTheta（极角的正弦值）。
    // 用于后续将球面坐标（phi, theta）转换为笛卡尔坐标。
    float sinTheta = sqrt(1.0 - cosTheta * cosTheta);
    // 创建半程向量 H（微表面法线）的局部空间坐标。
    // 使用球面坐标公式将 (phi, cosTheta, sinTheta) 转换为 3D 笛卡尔坐标。
    vec3 H;
    H.x = cos(phi) * sinTheta; // x 分量
    H.y = sin(phi) * sinTheta; // y 分量
    H.z = cosTheta;            // z 分量（指向半球上的点）
    // 确定一个基向量 up，与法线 N 足够正交。
    // 如果 N 接近 z 轴（垂直向上），选择 x 轴作为 up；否则选择 z 轴。
    vec3 up = abs(N.z) < 0.999 ? vec3(0.0, 0.0, 1.0) : vec3(1.0, 0.0, 0.0);
    // 计算切线空间的两个正交基向量：tangent 和 bitangent。
    // tangent 是 N 和 up 的叉积，确保与 N 和 up 都正交。
    vec3 tangent = normalize(cross(up, N));      // 切线向量
    vec3 bitangent = cross(N, tangent);          // 副切线向量
    // 将半程向量 H 从局部空间转换到世界空间。
    // 将 H 的坐标（局部坐标系下）组合成一个世界坐标向量 sampleVec。
    vec3 sampleVec = tangent * H.x + bitangent * H.y + N * H.z;
    // 返回归一化后的采样向量 sampleVec。
    // 采样向量是在 GGX 分布下，基于法线 N 和粗糙度生成的半程向量。
    return normalize(sampleVec);
}

在镜面IBL中，微表面法线与半程向量H的朝向一直。因为镜面反射定律要求入射光方向和反射光方向对称于微表面法线，而半程向量的方向正好满足这种对称性。

通过调用该函数，我们便可以获得一个采样向量，该向量大体围绕着预估的微表面半程向量（法线）。最终的计算预滤波环境贴图的着色器如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#version 330 core
out vec4 FragColor;
in vec3 localPos;

uniform samplerCube environmnetMap;
uniform float roughness;

const float PI = 3.141591265359;

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

void main(){
    // 天空盒永远位于世界原点，故其世界坐标等效于局部坐标
    // localPos归一化后得到从原点指向该片段的单位向量
    // 天空盒本质上可以看作球体，指向片段的向量就是该片段的法线
    vec3 N = normalize(localPos);
    // 由于我们在卷积环境贴图时事先不知道视角方向，因此假设视角方向——也就是镜面反射方向——总是等于输出采样方向，尽管该情况掠射镜面反射效果不完美，但也足够
    vec3 R = N;
    vec3 V = N;
    
    const uint SAMPLE_COUNT = 1024u; // 采样样本数N
    float totalWeight = 0.0; // 累加权重，表示光照的贡献值累加，用于最后取有权平均值
    vec3 prefilteredColor = vec3(0.0);
    // 开始循环采样
    for(uint i=0;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); // 使用反射公式推导,V相当于镜面反射方向，H相当于法线，此公式可以求出入射光方向L
        float NdotL = max(dot(N, L),0.0); //N dot L可以表示表面法线和光照方向的夹角大小。NdotL值越大，光照贡献越大。
        if(NdotL>0.0){
            prefilteredColor += texture(environmentMap, L).rgb * NdotL;
            totalWeight += NdotL;
        }
    }
    prefilteredColor = prefilteredColor / totalWeight;
    FragColor = vec4(prefilteredColor, 1.0);
}
捕获预过滤Mipmap级别

预滤波卷积着色器需要在不同的Mipmap级别上运行：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
prefilterShader.use();
prefilterShader.setInt("environmentMap", 0);
prefilterShader.setMat4("projection", captureProjection);
glActiveTexture(GL_TEXTURE0);
glBindFrameBuffer(GL_FRAMEBUFFER, envCubemap);
unsigned int maxMipLevels = 5;
for(unsigned int mip = 0; mip<maxMipLevels; mip++){
    // 计算Mipmap尺寸
    unsigned int mipWidth = 128 * std::pow(0.5,mip);
    unsigned int mipHeight = 128 * std::pow(0.5,mip);
    // 深度缓冲用RBO
    glBindRenderBuffer(GL_RENDERBUFFER, captureRBO);
    glRenderbufferStorage(GL_RENDERBUFFER, GL_DEPTH_COMPONENT24, mipWidth, mipHeight);
    // 记得调整视口尺寸
    glViewport(0, 0, mipWidth, mipHeight);
    float roughness = (float)mip/(float)(maxMipLevels - 1);
    prefilterShader.setFloat("roughness", roughness);
    // 渲染到立方体贴图的六个面
    for(unsigned int i=0;i<6;i++){
        prefilterShader.setMat4("view", captureViews[i]);
        // glFramebufferTexture2D的最后一个参数可以指定要渲染的目标mip级别
        glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_CUBE_MAP_POSITIVE_X + i, prefilterMap, mip);
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
        renderCube();
    }
}
glBindFramebuffer(GL_FRAMEBUFFER, 0);

通过glsl的textureLod函数，可以从指定Mip级别的纹理中采样。

预过滤卷积伪像

在使用上面的预过滤环境贴图时，会遇到下面两个问题：

  • 立方体贴图接缝

默认情况下，OpenGL不会在立方体贴图的面之间进行插值。在使用高粗糙度的预滤波环境贴图时，由于其分辨率较低，面之间的接缝会格外明显。

我们可以通过glEnable(GL_TEXTURE_CUBE_MAP_SEAMLESS)解决这一问题。

  • 亮点

镜面反射的光强度变化大，高频细节多，所以需要进行大量采样。尽管我们已经进行了数千次采样，但对某些级别的预滤波环境贴图Mipmap来说依然不够，导致明亮区域周边出现亮点：

image-20250106150230700

一种方案是提高样本数量，但会降低性能表现，治标不治本。

另一种方案是基于积分的PDF和粗糙度采样环境贴图的Mipmap。低分辨率的环境贴图，高频细节也少，采样突变的情况更少：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
float D = DistributionGGX(NdotH, roughness); // 计算NDF
float pdf = (D * NdotH / (4.0 * HdotV)) +0.0001; // 计算PDF（当前采样点的概率密度），用于权重归一化
float resolution = 512.0; // 环境贴图分辨率
// 计算每个纹素在立方体贴图中的立体角
// 	一个立方体有6个面，每个面包含resolution^2个纹素
// 	立方体贴图（即球面）覆盖了4PI个立体角
float saTexel = 4.0*PI/(6.0*resolution*resolution); 
// 立体角的大小与点的分布密度，即PDF成反比。如果一个采样点的概率密度很高，那么它出现的概率更大，所占的立体角就更小
// 概率密度越高，点更集中，单个采样点在球面上占的立体角更小。反之亦然。
float saSample = 1.0/(float(SAMPLE_COUNT)*pdf+0.0001); // 每个采样点对应的立体角
// 根据每个采样点的立体角 saSample和每个像素的立体角 saTexel，计算两者的比例
// 如果二者一致，说明采样点与像素刚好一对一对应，那么就不需要对环境贴图进行缩放
// 否则，意味着一个采样点对应着多个纹素，需要进行缩放（模糊处理）
float mipLevel = roughness == 0.0 ? 0.0 : 0.5 * log2(saSample / saTexel); // 当前粗糙度对应的Mipmap层级

进行该操作的前提是开启环境贴图的三线性过滤，以及生成环境贴图的Mipmap：

1
2
glBindTexture(GL_TEXTURE_CUBE_MAP, envCubemap);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR); 
1
glGenerateMipmap(GL_TEXTURE_CUBE_MAP);
预计算BRDF

回顾镜面部分的反射率方程：

Lo(p,ωo)=∫ΩLi(p,ωi)dωi∗∫Ωfr(p,ωi,ωO)n⋅ωidωiL_o(p,\omega_o)=\int_\Omega L_i(p,\omega_i)d\omega_i*\int_\Omega f_r(p,\omega_i,\omega_O)n·\omega_id\omega_i

目前，我们以及完成了左半部分的计算，得到了入射辐射率的卷积——辐照度。因此，我们可以将左半部分视为1（即纯白的环境光，或辐射度恒定为1.0），进行右侧计算。

我们作如下化简：

∫Ωfr(p,ωi,ωo)n⋅ωidωi=∫Ωfr(p,ωi,ωo)F(ωo,h)F(ωo,h)n⋅ωidωi\int_\Omega f_r(p,\omega_i,\omega_o)n·\omega_id\omega_i = \int_\Omega \frac{f_r(p,\omega_i,\omega_o)}{F(\omega_o,h)}F(\omega_o,h)n·\omega_id\omega_i

使用Fresnel-Schlick公式替换右侧的FF得到

∫Ωfr(p,ωi,ωo)F(ωo,h)(F0+(1−F0)(1−ωo⋅h)5)n⋅ωidωi\int_\Omega\frac{f_r(p,\omega_i,\omega_o)}{F(\omega_o,h)}(F_0+(1-F_0)(1-\omega_o·h)^5)n·\omega_id\omega_i

使用α\alpha替换(1−ωi⋅h)5(1-\omega_i·h)^5，得到：

∫Ωfr(p,ωi,ωo)F(ωo,h)(F0+(1−F0)α)n⋅ωidωi=∫Ωfr(p,ωi,ωo)F(ωo,h)(F0∗(1−α)+α)n⋅ωidωi\int_\Omega\frac{f_r(p,\omega_i,\omega_o)}{F(\omega_o,h)}(F_0+(1-F_0)\alpha)n·\omega_id\omega_i = \int_\Omega\frac{f_r(p,\omega_i,\omega_o)}{F(\omega_o,h)}(F_0*(1-\alpha)+\alpha)n·\omega_id\omega_i

拆分得：

∫Ωfr(p,ωi,ωo)F(ωo,h)(F0∗(1−α))n⋅ωidωi+∫Ωfr(p,ωi,ωo)F(ωo,h)(α)n⋅ωidωi\int\limits_{\Omega} \frac{f_r(p, \omega_i, \omega_o)}{F(\omega_o, h)} (F_0 * (1 - \alpha)) n \cdot \omega_i d\omega_i + \int\limits_{\Omega} \frac{f_r(p, \omega_i, \omega_o)}{F(\omega_o, h)} (\alpha) n \cdot \omega_i d\omega_i

还原：

F0∫Ωfr(p,ωi,ωo)(1−(1−ωo⋅h)5)n⋅ωidωi+∫Ωfr(p,ωi,ωo)(1−ωo⋅h)5n⋅ωidωiF_0 \int\limits_{\Omega} f_r(p, \omega_i, \omega_o)(1 - {(1 - \omega_o \cdot h)}^5) n \cdot \omega_i d\omega_i + \int\limits_{\Omega} f_r(p, \omega_i, \omega_o) {(1 - \omega_o \cdot h)}^5 n \cdot \omega_i d\omega_i

注意，BRDF项中得F项与原本分母的F项进行了约分，后续的BRDF项，即f，不计算F项。

由此，我们可以对BRDF求卷积，以nn和ωo\omega_o的夹角以及粗糙度作为输入，并将卷积结果存储在LUT中，得到BRDF积分贴图。

BRDF卷积着色器在2D平面上执行计算，使用其2D纹理坐标作为卷积输入。它根据BRDF几何函数和Fresnel-Schlick近似来处理采样向量。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
float GeometrySchlickGGX(float NdotV, float roughness)
{
    float a = roughness;
    float k = (a * a) / 2.0; // 在IBL中，k值不同

    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, flaot 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;
}
1
2
3
4
5
6
7
8
9
10
unsigned int brdfLUTTexture;
glGenTextures(1, &brdfLUTTexture);

// pre-allocate enough memory for the LUT texture.
glBindTexture(GL_TEXTURE_2D, brdfLUTTexture);
glTexImage2D(GL_TEXTURE_2D, 0, GL_RG16F, 512, 512, 0, GL_RG, GL_FLOAT, 0);
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);
1
2
3
4
5
6
7
8
9
10
11
glBindFramebuffer(GL_FRAMEBUFFER, captureFBO);
glBindRenderbuffer(GL_RENDERBUFFER, captureRBO);
glRenderbufferStorage(GL_RENDERBUFFER, GL_DEPTH_COMPONENT24, 512, 512);
glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, brdfLUTTexture, 0);

glViewport(0, 0, 512, 512);
brdfShader.use();
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
RenderQuad();

glBindFramebuffer(GL_FRAMEBUFFER, 0);
完成IBL反射

PBR着色器关键部分如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
#version 330 core
out vec4 FragColor;
in vec2 TexCoords;
in vec3 WorldPos;
in vec3 Normal;

uniform vec3 albedo;
uniform float metallic;
uniform float roughness;
uniform float ao;

uniform samplerCube irradianceMap;
uniform samplerCube prefilterMap;
uniform sampler2D brdfLUT;

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

uniform vec3 camPos;

const float PI = 3.14159265359;
// ----------------------------------------------------------------------------
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;
}
// ----------------------------------------------------------------------------
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;
}
// ----------------------------------------------------------------------------
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;
}
// ----------------------------------------------------------------------------
vec3 fresnelSchlick(float cosTheta, vec3 F0)
{
    return F0 + (1.0 - F0) * pow(clamp(1.0 - cosTheta, 0.0, 1.0), 5.0);
}
// ----------------------------------------------------------------------------
vec3 fresnelSchlickRoughness(float cosTheta, vec3 F0, float roughness)
{
    return F0 + (max(vec3(1.0 - roughness), F0) - F0) * pow(clamp(1.0 - cosTheta, 0.0, 1.0), 5.0);
}   
// ----------------------------------------------------------------------------
void main()
{		
    vec3 N = Normal;
    vec3 V = normalize(camPos - WorldPos);
    vec3 R = reflect(-V, N); 

    // 结合金属度，计算基本反射率F0
    vec3 F0 = vec3(0.04); 
    F0 = mix(F0, albedo, metallic);

    // 计算直接光源的辐射率
    vec3 Lo = vec3(0.0);
    for(int i = 0; i < 4; ++i) 
    {
        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; // 灯光在入射方向的辐射率

        // 计算BRDF
        float NDF = DistributionGGX(N, H, roughness); // 法线分布函数
        float G   = GeometrySmith(N, V, L, roughness); // 结合了几何遮蔽和几何阴影的几何函数。其中K值因为使用了IBL所以有所不同
        vec3 F    = fresnelSchlick(max(dot(H, V), 0.0), F0); // 菲涅尔方程
        
        vec3 numerator    = NDF * G * F;
        float denominator = 4.0 * max(dot(N, V), 0.0) * max(dot(N, L), 0.0) + 0.0001;
        vec3 specular = numerator / denominator;
        
        // 镜面反射系数与菲涅尔项相同，代表镜面反射光线的占比
        vec3 kS = F;
        // 基于能量守恒，计算漫反射光线占比
        vec3 kD = vec3(1.0) - kS;
        // 金属度越高，漫反射就越少
        kD *= 1.0 - metallic;	                
            
        // NdotL表示光线垂直照射到表面的程度。越大，入射光越垂直，对总辐射率的贡献越大
        float NdotL = max(dot(N, L), 0.0);        

        // 计算此片段在此光源的作用下的出射辐射率Lo
        Lo += (kD * albedo / PI + specular) * radiance * NdotL; // specular项中已包含了菲涅尔
    }   
    
    // ---使用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;    
    // 采样BRDFLUT，以获取BRDF项
    vec2 brdf  = texture(brdfLUT, vec2(max(dot(N, V), 0.0), roughness)).rg;
    // 计算镜面反射
    vec3 specular = prefilteredColor * (F * brdf.x + brdf.y);
    // 得到环境光
    vec3 ambient = (kD * diffuse + specular) * ao;
    
    vec3 color = ambient + Lo;

    // HDR色调映射
    color = color / (color + vec3(1.0));
    // 伽马矫正
    color = pow(color, vec3(1.0/2.2)); 

    FragColor = vec4(color , 1.0);

#图形学 #OpenGL #学习笔记
LearnOpenGL学习笔记（十四） - PBR与IBL
http://example.com/2025/01/07/LearnOpenGL学习笔记（十四）-PBR与IBL/
作者
Yoi
发布于
2025年1月7日
许可协议