基于图像的光照

本文最后更新于：2021年8月23日 凌晨

基于图像的环境光照，是指用一张图片来存储各个方向上环境光照的辐射率，并利用这张图片进行着色计算．本文介绍了基于图像的实时环境光照理论，并提供了利用Vulkan实现这一算法的源码和技术细节．

基于图像的实时环境光照

环境纹理的每一个纹素，都会存储来自某个方向的环境光的辐射率．考虑原始形式的反射方程

这里我们不考虑环境光以外的其他光照，并假设环境光对于任意着色点 $\mathbf{p}$ 都是相同的，因此环境纹理可以记作 $L_i(\boldsymbol{\omega}_i)$．相应地，反射方程可以写作

在本文中，我们所需的点 $\mathbf{p}$ 处的几何信息只有法向量 $\mathbf{n}$，故上式可写作

漫反射

简化 Fresnel 项

如果只考虑漫反射，反射方程为

其中

在实时渲染中，一个常见的思路是对复杂的积分进行预计算．一般地，对任意的 $\mathbf{n}$ 和 ，我们希望得到相应的 $L_o(\mathbf{n},\boldsymbol{\omega}_o)$．一方面，$\boldsymbol{\omega}_o$ 可能的取值有无穷多个．另一方面，如果考虑物体的旋转，$\mathbf{n}$ 的取值也有无穷多个．为所有的 $\mathbf{n}$ 和 $\boldsymbol{\omega}_o$ 预计算出对应的 显然是不可能的．最朴素的想法是对有限个 $\mathbf{n}$ 和 $\boldsymbol{\omega}_o$ 进行预计算，然后对结果进行插值．但是，这要求我们制作一张 4D 的 LUT（lookup table，查找表），因而会耗费大量的空间和时间．更加经济的做法，是通过简化反射方程中需要预计算的被积函数

让包含 $\boldsymbol{\omega}_o$ 的被积项能够提到积分号之外．这样，我们只需为有限个 $\mathbf{n}$ 预计算复杂的积分，并将结果存入一张 2D 的 LUT．关于 $\boldsymbol{\omega}_o$ 的计算，则可以实时去进行．

注意到在被积函数中，只有 Fresnel 项 包含 $\boldsymbol{\omega}_o$．因为 $\mathbf{h}$ 包含 $\boldsymbol{\omega}_i$，该项并不能直接提到积分之外．我们可以用简化版的 Fresnel 项 代替原本的 ，其表达式如下

于是，反射方程简化为

之后，我们只需为有限个 $\mathbf{n}$ 预计算如下积分

注意到该积分表示的是法线为 $\mathbf{n}$ 的着色点接受到的环境光辐照度，因此也被称作辐照度贴图（irradiance map）．

如果环境纹理是一张立方盒贴图，它的每一个纹素都会对应一个方向向量．我们可以让 $\mathbf{n}$ 取遍所有这些方向向量，并将 $I_{\mathrm{diff}}\left(\mathbf{n}\right)$ 的值存入对应的纹素．渲染时，$L_o(\mathbf{n},\boldsymbol{\omega}_o)$ 由下式进行计算

辐照度贴图的预计算

我们可以用 Monte Carlo 积分法计算 ．给定 $\mathbf{n}$，设 $\eta_{\mathbf{n}}$ 是半球 $H^2$ 上的一个概率分布，其概率密度函数为

我们利用 进行重要性采样，由此得到 $I_{\mathrm{diff}}\left(\mathbf{n}\right)$ 的统计量

其中 $N$ 为采样数．下面我们将具体描述如何为 生成样本．我们会提供两种方案：坐标变换法和四元数旋转法．

坐标变换法：

设

切空间的基向量在世界空间下的矩阵 $T_{\mathbf{n}}$可以表示成

其中

设 独立且服从均匀分布 $U(0,1)$，则

服从分布 $\eta_{\mathbf{n}}$．

四元数旋转法：

  • 当 $\mathbf{n}=(0,0,1)$ 时，我们可以使用逆变换采样．设 独立且服从均匀分布 $U(0,1)$，则

    服从分布 $\eta_{(0,0,1)}$．

  • 当 $\mathbf{n}=(0,0,-1)$ 时， 服从分布 $\eta_{(0,0,-1)}$．

  • 当 $\mathbf{n}=(n_x,n_y,n_z)\ne(0,0,\pm 1)$ 时，将向量 $(0,0,1)$ 沿旋转轴

    旋转成向量 $\mathbf{n}$ 的四元数为

    其中 $\alpha$ 是旋转角．因为

    我们有

    由此可以计算出四元数旋转变换 $\mathrm{rot}_q:p\mapsto qpq^{-1}$ 对应的旋转矩阵

    不难看出， 服从分布 $\eta_{\mathbf{n}}$．

镜面反射

下面我们只考虑镜面反射，反射方程为

为化简方程，我们可以将 $L_i(\boldsymbol{\omega}_i)$ 提出积分号外．于是，我们将 $L_o(\mathbf{n},\boldsymbol{\omega}_o)$ 写成如下形式

希望对含 $L_i(\boldsymbol{\omega}_i)$ 和不含 $L_i(\boldsymbol{\omega}_i)$ 的两部分分别进行计算．

分裂和 - 第 1 部分

公式推导

直接计算

是比较复杂的．为了降低预计算的参数空间维度．这里的假设是： lobe 的形状保持不变，始终为光线沿法线入射时的 lobe，即在计算积分时，让 $\boldsymbol{\omega}_o, \mathbf{n}$ 的取值与 $\boldsymbol{\omega}_o$ 的反射向量 $\mathbf{r}(\boldsymbol{\omega}_o)=2(\boldsymbol{\omega}_o\cdot\mathbf{n})\mathbf{n}-\boldsymbol{\omega}_o$ 值相等．或者说，用 $\mathbf{r}(\boldsymbol{\omega}_o)$ 代替被积函数中的 $\boldsymbol{\omega}_o$ 和 $\mathbf{n}$. 记 $\mathbf{n}^*=\mathbf{r}(\boldsymbol{\omega}_o)$，

其中

$C(\boldsymbol{\omega}_i,\mathbf{n}^*)$ 定义为

因为

其中 $\cos\theta_{\mathbf{h}}=\mathbf{n}^* \cdot \mathbf{h}$，可用蒙特卡洛积分法对

进行重要性采样．注意到

我们有

其中

而 $\mathbf{h}^{(k)}$ 是半球 $H^2_{\mathbf{n}^*}$ 上的随机变量，服从概率密度函数为

的概率分布．

值得一提的是，虚幻引擎在实验中发现：即使在得到的估计式去掉 ，效果也不错．同时，这也会让预计算和实时计算更加简单．本文则尝试保留 ，目的在于展示另一种可能的做法．

采样方案

为应用逆变换采样法，先将 $\mathbf{h}^{(k)}$ 在 $\mathbf{n}^*$ 处的切空间上参数化为球坐标 ．设 独立且服从均匀分布 $U(0,1)$．由

可知 是独立的．于是，我们可以分别求出 的边缘分布

使用 cdf 的逆变换，就可得到球坐标下 $\mathbf{h}^{(k)}$ 的采样方案

由此可推出球坐标下 $\boldsymbol{\omega}_i^{(k)}$ 的采样方案

因为上式的球坐标是在切空间中生成的，我们需要先将它们转换到笛卡尔坐标系下，再用坐标变换法转成成世界空间下的坐标．沿用之前定义的 $T_{\mathbf{n}}$，世界空间下 $\boldsymbol{\omega}_i^{(k)}$ 的采样方案可以写作

最终，我们可以导出第一部分分裂积分的估计量

其中

为了简化起见，我们只计算 $F_0$ 为 0 和 1 的情况．确定 $F_0$ 后每给定一个 $\alpha$，我们就可以让 取各个方向，预计算出 的值，得到一张立方盒贴图．对不同的 $\alpha$，我们可以像 Mipmap 一样生成多张立方盒贴图．$\alpha$ 越大，生成的立方盒贴图分辨率就越小．对于其它的 $\alpha$， 的值可由三线性插值得到．而对于 $F_0$ 不为 0 和 1 的情况，则对结果再进行一次线性插值．这些预计算得到的立方盒贴图，被称作 预过滤环境贴图（pre-filtered environment map）．

在渲染时，对于任意的 $\boldsymbol{\omega}_o$ 和 $\alpha$，第一部分分裂积分可估计为

分裂和 - 第 2 部分

公式推导

分裂积分的右边可以化简为

重要性采样

基于Vulkan的算法实现

环境纹理

矩形环境纹理

在 Poly Haven 网站上，我们可以免费下载 CC0 授权的高动态范围（HDR）环境纹理．该网站提供的这些环境纹理都是通过如下 等距矩形投影 得到的

其中 $(x,y,z)$ 和 $(r,\theta,\varphi)$ 分别是方向向量的直角坐标与球坐标，$(u,v)$ 是归一化到 $[0,1]$ 区间的纹理坐标，函数 $\mathrm{arctan2}$ 的定义为

图像为


图1 - $\operatorname{arctan2}$ 的函数图像

取模函数 $\mathrm{mod}$ 的定义为

实际中，纹理的宽度与高度之比通常为 $2:1$．下图就是一张通过等距矩形投影（equirectangular projection）得到的环境纹理．


图2 - 等距矩形投影环境纹理

存储辐射率的所有方向构成了一个球面，我们暂且称之为环境光球面．需要指出的是，一些 3D 软件在解释环境纹理时，会认为环境纹理是贴在环境光球面内部的．以上图为例，如果在 Blender 软件中打开这张环境纹理，我们会发现蓝色的牌子在白色的牌子左边．


图3 - Blender中展示的环境纹理

观察原始的矩形环境纹理，蓝色牌子同样在白色牌子左边．这就说明，Blender 在解释环境纹理时，使用的等距矩形投影其实是下式

立方盒贴图

除了这种矩形环境纹理，立方盒贴图（Cubemap）也常常被用作环境纹理．存储辐射率的所有方向构成了一个球面，立方盒贴图可以看作是把球面上每一点，沿径向投影到与球面中心重合的立方体的表面上．因此，立方盒贴图是由立方体 6 个表面的 6 张正方形贴图共同组成的．因为相对于矩形环境纹理而言，采样立方盒贴图在性能上更有优势，我们在 Vulkan 中也将使用这种立方盒贴图作为环境纹理．

在继续介绍立方盒贴图之前，我们要先讲清楚 Vulkan 的坐标系问题．OpenGL 使用如下世界坐标系


图4 - OpenGL 坐标系

Vulkan 使用的世界坐标系则与 OpenGL 不同．但为了使用基于 OpenGL 坐标系的数学库 GLM，我对 VkViewport 对象进行了如下设置．

1
2
3
4
5
VkViewport viewport{};
viewport.width = renderer.extent.width;
viewport.height = -renderer.extent.height;
viewport.x = 0;   
viewport.y = renderer.extent.height;

这样一来， Vulkan 的坐标系就和与 OpenGL 一样了．这部分内容也可以参考 Sascha Willems 的博客，里面有更加详细的介绍．

接下来，我们介绍 Vulkan 对于立方盒贴图的规定．首先，与矩形环境纹理贴在环境光球面内部不同，6 张立方盒贴图是贴在立方盒外部的． 下图是图 1 对应的立方盒贴图．


图5 - 立方盒贴图环境纹理

放大观察，会发现此时蓝色牌子跑到了白色牌子的右边（见图5左）．但如果我们对立方盒贴图进行采样，然后将它当作环境光进行渲染，得到的显示结果是蓝色牌子在白色牌子的左边（见图5右）．立方盒贴图渲染程序的源码见 003_Vulkan_Skybox，具体实现的细节会在之后进行介绍．


图6 - 立方盒贴图环境纹理（左） 渲染结果（右）

事实上，Vulkan 对立方盒贴图布局的规定可以由下图表明．


图7 - 立方盒贴图布局

以 FACE 0 为例，$+X$ 表示该图贴在 $+X$ 轴所指的立方体的面上，也即立方体的右面．左下角的坐标轴 $(-z,y)$，标明了这张贴图在世界空间的朝向：向右是 $-z$ 方向，向上是 $y$ 方向．右上角的坐标轴 $(u,v)$ 则标明了这张贴图在纹理空间的朝向．这也印证了我们之前所说的：立方盒贴图是贴在立方体外部的．

在 Vulkan 中，立方盒贴图的 6 张图对应 VkImage 对象的 6 层．需要注意的是，6 张图片的顺序必须与图 6 所示的 FACE 0 到 FACE 5 的顺序一致．

从磁盘读写环境纹理

因为环境纹理存储的是辐射率，所以一般情况下图片格式均为浮点类型．最常见的图片格式是 HDR 和 EXR．

HDR 图片的每个纹素占用 32 比特，采用 R8G8B8E8 的编码方式，即每个色彩通道是 8 比特，三个色彩通道共用 8 比特的指数．我们使用只含头文件的轻量级开源库 stb 对它进行读写．读 HDR 可以使用 stbi_loadf 函数，返回的指针所指内存的布局为 R32G32B32A32 ，即 4 通道浮点格式．写 HDR 可以使用 stbi_write_hdr 函数，将布局为 R32G32B32A32 的数据存为 HDR．

EXR 图片的格式比较复杂．对于环境纹理而言，最常用的是 B32G32R32 布局的 3 通道浮点格式，和 B16G16R16 布局的 3 通道半精度浮点格式．我们使用只含头文件的轻量级开源库 tinyexr 对它进行读写．读 EXR 可以使用 LoadEXR 函数，对于 B32G32R32、B16G16R16、A32B32G32R32、A16B16G16R16 格式的 EXR 图片，返回的指针所指内存的布局均为 R32G32B32A32 ，即 4 通道浮点格式．写 EXR 可以用 SaveEXRImageToFile 函数，将布局为 R32G32B32A32 的数据存为指定格式的 EXR 图片．

矩形纹理转立方盒贴图

完整程序源码见 Equirectangular To Cubemap．

将矩形纹理转成立方盒贴图，思路就是把矩形纹理贴到一个立方体上，然后分别对立方体的 6 个面进行拍屏，使立方体的表面正好占据整个屏幕空间．

下面的代码会输出投影矩阵和观察矩阵的乘积，分别表示从右、左、上、下、前、后架设摄像机拍摄立方体．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
glm::mat4 captureProjection =
	glm::perspective((float)M_PI / 2.0f, 1.0f, 0.1f, 10.0f);
glm::mat4 captureViews[] = {
    glm::lookAt(glm::vec3(2.0f, 0.0f, 0.0f), glm::vec3(0.0f,  0.0f,  0.0f), glm::vec3(0.0f, 1.0f,  0.0f)),
    glm::lookAt(glm::vec3(-2.0f, 0.0f, 0.0f), glm::vec3(0.0f,  0.0f,  0.0f), glm::vec3(0.0f, 1.0f,  0.0f)),
    glm::lookAt(glm::vec3(0.0f, 2.0f, 0.0f), glm::vec3(0.0f, 0.0f,  0.0f), glm::vec3(0.0f,  0.0f, -1.0f)),
    glm::lookAt(glm::vec3(0.0f, -2.0f, 0.0f), glm::vec3(0.0f,  0.0f,  0.0f), glm::vec3(0.0f,  0.0f,  1.0f)),
    glm::lookAt(glm::vec3(0.0f, 0.0f, 2.0f), glm::vec3(0.0f,  0.0f,  0.0f), glm::vec3(0.0f, 1.0f,  0.0f)),
    glm::lookAt(glm::vec3(0.0f, 0.0f, -2.0f), glm::vec3(0.0f,  0.0f, 0.0f), glm::vec3(0.0f, 1.0f,  0.0f))
};
for (int i = 0; i < 6; i++)
{
    std::cout << glm::to_string(captureProjection * captureViews[i]) << std::endl;
}

我们可以将输出的矩阵复制到顶点着色器中进行着色．下面是顶点着色器的代码．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
//shader.vert
#version 450

layout(push_constant) uniform PushConstant {
    int index;
} constant;
layout(location = 0) in vec3 aLocalPos;
layout(location = 0) out vec3 localPos;

const mat4 mProjViews[6] = {
    {{0.000000, 0.000000, -1.010101, -1.000000}, {0.000000, 1.000000, 0.000000, 0.000000}, {-1.000000, 0.000000, 0.000000, 0.000000}, {0.000000, 0.000000, 1.919192, 2.000000}},
    {{0.000000, 0.000000, 1.010101, 1.000000}, {0.000000, 1.000000, 0.000000, 0.000000}, {1.000000, 0.000000, 0.000000, 0.000000}, {0.000000, 0.000000, 1.919192, 2.000000}},
    {{1.000000, 0.000000, 0.000000, 0.000000}, {0.000000, 0.000000, -1.010101, -1.000000}, {0.000000, -1.000000, 0.000000, 0.000000}, {0.000000, 0.000000, 1.919192, 2.000000}},
    {{1.000000, 0.000000, 0.000000, 0.000000}, {0.000000, 0.000000, 1.010101, 1.000000}, {0.000000, 1.000000, 0.000000, 0.000000}, {0.000000, 0.000000, 1.919192, 2.000000}},
    {{1.000000, 0.000000, 0.000000, 0.000000}, {0.000000, 1.000000, 0.000000, 0.000000}, {0.000000, 0.000000, -1.010101, -1.000000}, {0.000000, 0.000000, 1.919192, 2.000000}},
    {{-1.000000, 0.000000, 0.000000, 0.000000}, {0.000000, 1.000000, 0.000000, 0.000000}, {0.000000, 0.000000, 1.010101, 1.000000}, {0.000000, 0.000000, 1.919192, 2.000000}}
};

void main() {
    localPos = aLocalPos;
    gl_Position = mProjViews[constant.index] * vec4(aLocalPos, 1.0);
}

这里我们不打算依靠 UV 映射为立方体去贴纹理，而是根据世界坐标去贴纹理，所以需要将世界坐标传入片元着色器．片元着色器的代码如下．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
//shader.frag
#version 450

layout(binding = 0) uniform sampler2D equirectangularMap;
layout(location = 0) in vec3 localPos;
layout(location = 0) out vec4 outColor;

const vec2 invAtan = vec2(0.1591, 0.3183);

vec2 SampleSphericalMap(vec3 v) {
    vec2 uv = vec2(atan(v.x, v.z), acos(v.y));
    uv *= invAtan;
    uv.x = 0.5 - uv.x;
    return uv;
}
void main() {
    vec2 uv = SampleSphericalMap(normalize(localPos)); 
    vec3 color = texture(equirectangularMap, uv).rgb;
    outColor = vec4(color, 1.0);
}

片元着色器中，SampleSphericalMap 函数基本上就是之前提到的投影

但是，具体形式上有所差别．一是因为如果不在意贴图是否发生了旋转，那么取模函数 $x\mapsto\mathrm{mod}(x,2\pi)$ 可以用平移函数 $x\mapsto x+\pi$ 进行替换．二是因为我们在 Vulkan 中使用的坐标系与本文在数学推导时采用的坐标系不同，在将本文公式转化为 GLSL 代码时，应该先把公式中所有的向量 $\mathbf{v}=(x,y,z)$ 用 $\mathbf{v}’=(z,x,y)$ 进行替换．

Vulkan中的立方盒贴图

在 Vulkan 中，在为立方盒贴图创建 VkImage 对象时，需要如下创建信息．

1
2
3
4
VkImageCreateInfo imageInfo{};
//...
imageInfo.arrayLayers = 6;
imageInfo.flags = VK_IMAGE_CREATE_CUBE_COMPATIBLE_BIT;

其他涉及到VkImage的操作，也同样要将图片的层数layerCount设置为 6．比如说，如果使用了 staging buffer，在进行vkCmdCopyBufferToImage时，要设置

1
2
3
VkBufferImageCopy region{};
//...
region.imageSubresource.layerCount = 6;

在使用vkCmdPipelineBarrier时，要设置

1
2
3
VkImageMemoryBarrier barrier{};
//...
barrier.subresourceRange.layerCount = 6;

在使用vkCmdBlitImage创建 Mipmap 的时候，也要设置

1
2
3
4
VkImageBlit blit{};
//...
blit.srcSubresource.layerCount = 6;
blit.dstSubresource.layerCount = 6;

在为立方盒贴图创建 VkImageView 对象时，需要如下创建信息．

1
2
3
4
VkImageViewCreateInfo viewInfo{};
//...
viewInfo.viewType = VK_IMAGE_VIEW_TYPE_CUBE;
viewInfo.subresourceRange.layerCount = 6;

立方盒贴图的渲染

完整程序源码见 Vulkan Skybox．

在渲染立方盒贴图时，我们要消除摄像机观察矩阵中的平移，因此会对观察矩阵做如下处理．

1
ubo.view = glm::mat4(glm::mat3(camera.viewMatrix()));

着色器的代码就是简单地渲染一个立方体．只需要注意，在采样立方盒贴图时，需要使用一个一个vec3类型的方向向量作为纹理坐标．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
//shader.vert
#version 450

layout(binding = 0) uniform UniformBufferObject {
    mat4 view;
    mat4 proj;
} ubo;

layout(location = 0) in vec3 inPosition;
layout(location = 0) out vec3 fragTexCoord;

void main() {
    gl_Position = ubo.proj * ubo.view * vec4(inPosition, 1.0);
    fragTexCoord = inPosition;
}
1
2
3
4
5
6
7
8
9
10
//shader.frag
#version 450

layout(binding = 1) uniform samplerCube texSampler;
layout(location = 0) in vec3 fragTexCoord;
layout(location = 0) out vec4 outColor;

void main() {
    outColor = texture(texSampler, fragTexCoord);
}

预计算

辐照度贴图

完整程序源码见 Irradiance Map．

对于重要性采样所需的随机数，我们使用 Sobol 低差异序列来加速收敛．

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
//shader.frag
#version 450

layout(binding = 1) uniform samplerCube cubeMap;

layout(location = 0) in vec3 localPos;
layout(location = 0) out vec4 outColor;

#define PI 3.14159265359
#define SAMPLE_NUM 262144

uvec2 Sobol(uint n) {
    uvec2 p = uvec2(0u);
    uvec2 d = uvec2(0x80000000u);
    for(; n != 0u; n >>= 1u) {
        if((n & 1u) != 0u)
            p ^= d;        
        d.x >>= 1u; 2 Van der Corput
        d.y ^= d.y >> 1u;
    }
    return p;
}

uint OwenHash(uint x, uint seed) {
    x ^= x * 0x3d20adeau;
    x += seed;
    x *= (seed >> 16) | 1u;
    x ^= x * 0x05526c56u;
    x ^= x * 0x53a22864u;
    return x;
}

uint OwenScramble(uint p, uint seed) {
    p = bitfieldReverse(p);
    p = OwenHash(p, seed);
    return bitfieldReverse(p);
}

vec2 sobol2d(uint index, uint seed) {
    uvec2 p = Sobol(index);
    p.x = OwenScramble(p.x, 0xe7843fbfu);
    p.y = OwenScramble(p.y, 0x8d8fb1e0u);
    return vec2(p) / float(0xffffffffu);
}

vec3 genNu(uint index, uint seed) {
    vec2 psi = sobol2d(index, seed);
    float a = sqrt(psi.x);
    float b = 2.0 * PI * psi.y;
    return vec3(a * sin(b), sqrt(1.0 -psi.x), a * cos(b));
}

vec3 genSampleDirection(vec3 nu, vec3 N) {
    vec3 T = cross(vec3(0.0, 1.0, 0.0), N);
    vec3 B = cross(N, T);
    return B * nu.x + N *nu.y + T * nu.z;
}

void main() {
    vec3 irradiance = vec3(0.0);
    for(uint i = 0; i < SAMPLE_NUM; i++){
        uint seed = uint(float(0x00004000u) * gl_FragCoord.x) +
            uint(float(0x40000000u) * gl_FragCoord.y);
        vec3 nu = genNu(i, OwenHash(i, seed));
        vec3 dir = genSampleDirection(nu, normalize(localPos));
        irradiance += texture(cubeMap, dir).rgb;
    }
    irradiance /= float(SAMPLE_NUM);
    outColor = vec4(irradiance, 1.0);
}

$2^{18}=262144$ 次采样下得到的辐照度贴图如下．


图8 - 辐照度贴图

预过滤环境贴图

完整程序源码见 Pre-filtered Map．

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
//shader.frag
#version 450

layout(binding = 0) uniform samplerCube cubeMap;

layout(location = 0) in vec3 localPos;
layout(location = 0) out vec4 outColor;

#define PI 3.14159265359
#define SAMPLE_NUM 524288
#define ALPHA 0.999999
#define F0 1

uvec2 Sobol(uint n) {
    uvec2 p = uvec2(0u);
    uvec2 d = uvec2(0x80000000u);
    for(; n != 0u; n >>= 1u) {
        if((n & 1u) != 0u)
            p ^= d;
        d.x >>= 1u;
        d.y ^= d.y >> 1u;
    }
    return p;
}

uint OwenHash(uint x, uint seed) {
    x ^= x * 0x3d20adeau;
    x += seed;
    x *= (seed >> 16) | 1u;
    x ^= x * 0x05526c56u;
    x ^= x * 0x53a22864u;
    return x;
}

uint OwenScramble(uint p, uint seed) {
    p = bitfieldReverse(p);
    p = OwenHash(p, seed);
    return bitfieldReverse(p);
}

vec2 sobol2d(uint index, uint seed) {
    uvec2 p = Sobol(index);
    p.x = OwenScramble(p.x, 0xe7843fbfu);
    p.y = OwenScramble(p.y, 0x8d8fb1e0u);
    return vec2(p) / float(0xffffffffu);
}

vec3 genNu(uint index, uint seed) {
    vec2 psi = sobol2d(index, seed);
    float a_sqr = ALPHA * ALPHA;
    float denominator = 1.0 - (1.0 - a_sqr) * psi.x;
    float numerator_part = 2.0 * ALPHA * sqrt(psi.x * (1.0 - psi.x));
    float Nu_z = numerator_part * cos(2.0 * PI * psi.y);
    float Nu_x = numerator_part * sin(2.0 * PI * psi.y);
    float Nu_y = 1.0 - (1.0 + a_sqr) * psi.x;
    return vec3(Nu_x, Nu_y, Nu_z) / denominator;
}

vec3 genSampleDirection(vec3 nu, vec3 N) {
    vec3 T = cross(vec3(0.0, 1.0, 0.0), N);
    vec3 B = cross(N, T);
    return B * nu.x + N * nu.y + T * nu.z;
}

void main() {
    vec3 N = normalize(localPos);
    vec3 radiance = vec3(0.0);
    float weight_sum = 0.0;
    for(uint i = 0; i < SAMPLE_NUM; i++){
        uint seed = uint(float(0x00004000u) * gl_FragCoord.x) +
            uint(float(0x40000000u) * gl_FragCoord.y);
        vec3 nu = genNu(i, OwenHash(i, seed));
        vec3 wi = genSampleDirection(nu, N);
        float NdotWi = max(dot(N, wi), 0.0);
#if (F0 == 1)
        float c = 0.5 / (ALPHA + (2.0 - ALPHA) * NdotWi);
#else
        float c_base = 1 - dot(N, normalize(N + wi));
        float c = 0.5 * c_base * c_base * c_base * c_base * c_base /
            (ALPHA + (2.0 - ALPHA) * NdotWi);
#endif
        float weight = NdotWi * c;
        weight_sum += weight;
        radiance += weight * texture(cubeMap, wi).rgb;
    }
    radiance /= weight_sum;
    outColor = vec4(radiance, 1.0);
}

我们选取 $\alpha=0.0^2, 0.2^2, 0.4^2, 0.6^2, 0.8^2, 1.0^2$ 这 6 个值进行预计算．当 $F_0=1$ 时，$2^{19}=524288$ 次采样下得到的预过滤贴图如下．


图9 - 预过滤贴图

BRDF 2D LUT

计算机图形
计算机图形 PBR

本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议 ，转载请注明出处！