前言
本篇主要是对《Coding Labs :: Physically Based Rendering》中内容的总结与实践。
本篇将对一些原理进行概况，但可能理解并不准确。
随后，将实践原文中的代码并观察效果。
渲染方程
概括上讲，渲染方程所描述的是：在知道一位置的所有入射光的情况下，这个位置某一方向的出射光是什么？
L
o
(
p
,
ω
o
)
=
∫
Ω
f
r
(
p
,
ω
i
,
ω
o
)
L
i
(
p
,
ω
i
)
(
n
⋅
ω
i
)
d
ω
i
L_o(p,\omega_o) = \int_\Omega f_r(p,\omega_i,\omega_o)L_i(p,\omega_i)(\bf{n} \cdot \omega_i) \it{d} \omega_i
Lo(p,ωo)=∫Ωfr(p,ωi,ωo)Li(p,ωi)(n⋅ωi)dωi
方程中基础符号的意思如下：
| 符号 | 意思 |
|---|---|
| p p p | 位置 |
| n \bf{n} n | 表面的法线方向（宏观） |
| i i i 与 o o o | 入射 与 出射。 |
| ω i \omega_i ωi 与 ω o \omega_o ωo | 入射方向 与 出射方向。 |
| L i L_i Li 与 L o L_o Lo | 入射光 与 出射光。 |
| f r f_r fr | 双向反射分布函数（BRDF），他表示了一条入射光在各个方向上的出射情况。 |
| ∫ Ω \int_\Omega ∫Ω | 在半球面上积分 |
那么，组合起来的意思就是：
| 符号 | 意思 |
|---|---|
| L o ( p , ω o ) L_o(p,\omega_o) Lo(p,ωo) | 在 p p p 位置上， ω o \omega_o ωo 方向上的出射光。可以看作是需要求出的结果。 |
| L i ( p , ω i ) L_i(p,\omega_i) Li(p,ωi) | 在 p p p 位置上， ω i \omega_i ωi 方向上的入射光。可以看作是需要提供的数据。 |
| f r ( p , ω i , ω o ) f_r(p,\omega_i,\omega_o) fr(p,ωi,ωo) | 在 p p p 位置上， ω i \omega_i ωi 方向上的入射光有多少比例作用到 ω o \omega_o ωo 方向上。例如对于一个绝对的镜面反射比如镜子，它的情况会是这样：当 ω o \omega_o ωo 精确为 ω i \omega_i ωi 的反射方向时，函数结果为 1.0，否则函数结果就是 0.0 。 |
| ( n ⋅ ω i ) (\bf{n} \cdot \omega_i) (n⋅ωi) | 表面法线方向点乘入射光方向，也就是乘算上入射角的余弦。 |
| ∫ Ω f r ( p , ω i , ω o ) L i ( p , ω i ) ( n ⋅ ω i ) d ω i \int_\Omega f_r(p,\omega_i,\omega_o)L_i(p,\omega_i)(\bf{n} \cdot \omega_i) \it{d} \omega_i ∫Ωfr(p,ωi,ωo)Li(p,ωi)(n⋅ωi)dωi | 在整个半球上对所有入射光对 ω o \omega_o ωo方向上的贡献进行积分。（为什么不是一个完整的球？因为不需要算背面半球的光） |
如果这一函数能够解答，那么就可以求出场景中某一点到眼睛的方向的出射光，也就是应渲染出的颜色。
CubeMap
要计算渲染方程所遇到的第一个问题是：如何知道场景中一点的所有方向上的入射光？
CubeMap 是一种方案，其将所有方向的入射光存到六个面的立方体中。
在之前的博客《图形API学习工程（21）：使用CubeMap纹理》中已经实践了相关操作。
Lambert
接下来遇到的问题是关于BRDF。
首先，BRDF肯定是取决与材质的。比如对于一个绝对的镜面反射，他的BRDF就很简单：只有在精确是反射的方向上BRDF结果为1，其余都是0。那么在这种情况下就只需要知道那一条完全反射的光线就可以了。
但是对于更粗糙的表面呢？四面八方的入射光都需要考虑，那么就需要进行实时的积分了，而这样的计算在实时中计算速度是比较慢的。
不过对于Lambert的BRDF，问题变得简化：它假设了BRDF的结果只取决于入射光，而与出射光无关（即与观察的视角无关），此时相当于表面是一个非常粗糙的完全没有高光的材质。公式：
L
o
(
p
,
ω
o
)
=
∫
Ω
c
π
L
i
(
p
,
ω
i
)
(
n
⋅
ω
i
)
d
ω
i
=
c
π
∫
Ω
L
i
(
p
,
ω
i
)
(
n
⋅
ω
i
)
d
ω
i
\begin{aligned} L_o(p,\omega_o) & = \int_\Omega\frac{c}{\pi}L_i(p,\omega_i)(\bf{n} \cdot \omega_i) \it{d} \omega_i \\ & = \frac{c}{\pi}\int_\Omega L_i(p,\omega_i)(\bf{n} \cdot \omega_i) \it{d} \omega_i \\ \end{aligned}
Lo(p,ωo)=∫ΩπcLi(p,ωi)(n⋅ωi)dωi=πc∫ΩLi(p,ωi)(n⋅ωi)dωi
（公式里为什么要除以 π \pi π还不理解，待后续学习）
这样，就可以在所有 n \bf{n} n方向上进行积分，将结果存到一张CubeMap上了，这也是本篇要做的事情。
积分
由于原文作者在shader中将在球面坐标上积分，所以下面以球面坐标的形式改写公式：
L
o
(
p
,
θ
o
,
ϕ
o
)
=
c
π
∫
Φ
∫
Θ
L
i
(
p
,
θ
i
,
ϕ
i
)
c
o
s
(
θ
i
)
s
i
n
(
θ
i
)
d
θ
i
d
ϕ
i
L_o(p,\theta_o,\phi_o)=\frac{c}{\pi} \int_\Phi \int_\Theta L_i(p,\theta_i,\phi_i)cos(\theta_i)sin(\theta_i)d\theta_id\phi_i
Lo(p,θo,ϕo)=πc∫Φ∫ΘLi(p,θi,ϕi)cos(θi)sin(θi)dθidϕi
求解它的一种方式是使用蒙特卡洛，但原文作者选择使用离散化的方式：
L
o
(
p
,
θ
o
,
ϕ
o
)
=
c
π
2
π
N
1
π
2
N
2
∑
N
1
∑
N
2
L
i
(
p
,
θ
i
,
ϕ
i
)
c
o
s
(
θ
i
)
s
i
n
(
θ
i
)
=
π
c
N
1
N
2
∑
N
1
∑
N
2
L
i
(
p
,
θ
i
,
ϕ
i
)
c
o
s
(
θ
i
)
s
i
n
(
θ
i
)
\begin{aligned} L_o(p,\theta_o,\phi_o) & = \frac{c}{\pi} \frac{2\pi}{N_1} \frac{\pi}{2N_2} \sum^{N_1} \sum^{N_2} L_i(p,\theta_i,\phi_i)cos(\theta_i)sin(\theta_i)\\ & = \frac{\pi c}{N_1N_2} \sum^{N_1} \sum^{N_2} L_i(p,\theta_i,\phi_i)cos(\theta_i)sin(\theta_i)\\ \end{aligned}
Lo(p,θo,ϕo)=πcN12π2N2π∑N1∑N2Li(p,θi,ϕi)cos(θi)sin(θi)=N1N2πc∑N1∑N2Li(p,θi,ϕi)cos(θi)sin(θi)
这个公式将原先的积分变成了离散的求和。其中 N 1 N_1 N1是方位角划分的步数， N 2 N_2 N2是极角划分的步数
代码
对应的着色器代码如下（原文是hlsl，而下面是glsl）：
#version 450
#extension GL_ARB_separate_shader_objects : enable
layout(location = 0) in vec2 InterpolatedPosition;
//UniformBuffer
layout(binding = 0) uniform UniformBufferObject {
int cubeFace;
} ubo;
//CubeMap纹理
layout(binding = 9) uniform samplerCube ourCubeMap;
out vec4 fColor;
void main()
{
float PI = 3.1415926;
//这个函数用来计算“特定表面法线方向”上的出射辐射率，结果将存在CubeMap的对应位置上
//cubeFace代表的是CubeMap中哪一个面
int cubeFace = ubo.cubeFace;
//InterpolatedPosition代表“特定表面法线方向”对应在这一个CubeMap的面上的位置，范围(-1~1)
//下面根据cubeFace来换算normal：“特定表面法线方向”在世界空间下的方向
vec3 normal = normalize( vec3(InterpolatedPosition.xy, 1) );
if(cubeFace==2)
normal = normalize( vec3(InterpolatedPosition.x, 1, -InterpolatedPosition.y) );
else if(cubeFace==3)
normal = normalize( vec3(InterpolatedPosition.x, -1, InterpolatedPosition.y) );
else if(cubeFace==0)
normal = normalize( vec3( 1, InterpolatedPosition.y,-InterpolatedPosition.x) );
else if(cubeFace==1)
normal = normalize( vec3( -1, InterpolatedPosition.y, InterpolatedPosition.x) );
else if(cubeFace==5)
normal = normalize( vec3(-InterpolatedPosition.x, InterpolatedPosition.y, -1) );
//下面计算球面坐标的轴在世界空间下的方向，
//其中up是天顶方向，right是方位角基准方向
vec3 up = vec3(0,1,0);
vec3 right = normalize(cross(up,normal));
up = cross(normal,right);
vec3 sampledColour = vec3(0,0,0); //逐渐累加的采样得到的颜色值，最后会根据采样次数平均化
float index = 0; //采样次数，用于最后的平均化
for(float phi = 0; phi < 6.283; phi += 0.025) //采样所有的方位角方向
{
for(float theta = 0; theta < 1.57; theta += 0.1)//采样顶半球所有的极角方向
{
//下面计算sampleVector：本次采样的方向在世界空间下的方向
vec3 temp = cos(phi) * right + sin(phi) * up;
vec3 sampleVector = cos(theta) * normal + sin(theta) * temp;
//从环境贴图中采样，并套入本篇上面讨论的公式
sampledColour += texture( ourCubeMap, sampleVector ).rgb *
cos(theta) * sin(theta);
index = index+1;
}
}
//最终结果平均化：
fColor = vec4(( PI * sampledColour / index), 1 );
}
其中，cubeFace是UniformBuffer，在C++程序中控制从0递增到5。且每次渲染后都保存为图片。（关于UniformBuffer已经在《图形API学习工程（6）：创建并使用UniformBuffer》中做了实践，保存图片的功能在《将OpenGL渲染的结果保存为图片》中做了实践）
//六个面保存的图片：
std::vector<std::string> fileNames = {
"right.bmp",
"left.bmp",
"top.bmp",
"bottom.bmp",
"front.bmp",
"back.bmp"
};
//渲染六个面：
for (int i = 0; i < 6; i++)
{
//更新cubeFace的UniformBuffer：
UBData.cubeFace = i;
UpdateUniformBuffer();
//显示（调用 draw call）
display();
//将当前渲染结果保存为图片
ScreenShot(fileNames[i]);
//交换前后缓冲
glfwSwapBuffers(window);
//轮询并处理事件
glfwPollEvents();
}
完整工程代码见 https://codechina.csdn.net/u013412391/cubemaplambertdiffuse
结果
运行程序后，将会依次输出6张贴图。与原始的CubeMap对比如下：
我用原文中的图片做测试，并仿照他的样子拼合，结果如下：
但是原文的结果如下：
目前我并不清楚为什么和原文有差异，比较怀疑的一个原因是原文的CubeMap是HDR，但我这里只是保存了一个低分辨率的图片。