深入理解 PBR/基于图像照明 (IBL)

深入理解 PBR/基于图像照明 (IBL)

话接上文 离线渲染 | 金属工作流,当时我只实现了关于点光源的 PBR,效果不是很好。之前看到了学长的文章

膜力鸭苏蛙可:引擎搭建记录(1) - 基于物理的渲染(PBR)zhuanlan.zhihu.com图标

我才知道,环境光强无敌呀!

所以我自己也来实现一下。

这篇文章,我会把公式的由来,推导和理解都尽量讲清楚,做到知其然和知其所以然。故文章取名“深入理解”。

// 不是非常深入,应该是浅入(保命

同时我会给出离线渲染和实时渲染中的基于图像照明(Image-Based Lighting, IBL)的实现。

目录

  1. 基础
    1. 反射方程
      1. 菲涅尔系数 F
      2. 法线分布函数 D
      3. 几何函数 G
    2. 金属工作流
    3. 基于图像的照明(Image Based Lighting, IBL)
  2. 离线 IBL
    1. 别名方法 Alias Method
    2. 采样
  3. 实时 IBL
    1. irradiance map
    2. specular
      1. pre-filtered environment map
      2. BRDF LUT
    3. 综合

1. 基础

1.1 反射方程

这个与上篇文章高度重复。
因为补充了一些细节,主要是各公式的发展历程,所以这里再讨论一遍

反射方程[1],描述了给定入射光和反射规律的情况下,从 p 点沿着观察方向 \omega_o 的出射光

其中 f_r 是双向反射分布函数(Bidirectional Reflectance Distribution Function, BRDF),用于描述反射规律。示例如下

完美镜面反射 ideal specular
完美漫反射 ideal diffuse
光滑镜面 glossy specular
回射 retro-refective

基于物理的渲染(Physically Based Rendering,PBR)是指使用基于物理原理和微平面理论(microfacet theory)建模的着色/光照模型,以及使用从现实中测量的表面参数来准确表示真实世界材质的渲染理念。

图片来源:https://zhuanlan.zhihu.com/p/53086060

体现在反射方程上,就是要求 f_r 要基于物理原理和微平面理论。

常用的微平面模型为 Cook-Torrance[2] 模型

\begin{aligned} f_r(N,V,L) &= k_d f_d + f_s\\ f_d &= \frac{\rho}{\pi}\\ k_s &= F\\ f_s &= k_s \frac{DG}{4(N\cdot V)(N\cdot L)}=\frac{FDG}{4(N\cdot V)(N\cdot L)}\\ k_d + k_s &= 1\\  \end{aligned}

其中 f_d 是漫反射项, f_s 是镜面反射项,系数 k_dk_s 保证了能量守恒。

以下分别讲讲 FDG 的具体形式。

1.1.1 菲涅尔系数 F

F 是菲涅尔系数[3],用于描述反射率(值的范围是 0 - 1)。

原公式比较复杂,Schlick 给出了简化版本[4]

\begin{aligned} F ( \theta ) & = F _ { 0 } + \left( 1 - F _ { 0 } \right) ( 1 - \cos \theta ) ^ { 5 } \\ F _ { 0 } & = \left( \frac { \eta _ { 1 } - \eta _ { 2 } } { \eta _ { 1 } + \eta _ { 2 } } \right) ^ { 2 } \end{aligned}

虚幻给出了更加速的版本[5]

F ( \mathbf { v } , \mathbf { h } ) = F _ { 0 } + \left( 1 - F _ { 0 } \right) 2 ^ { ( - 5.55473 ( \mathbf { v } \cdot \mathbf { h } ) - 6.98316 ) ( \mathbf { v } \cdot \mathbf { h } ) }

关键是 exp2 这个内建函数,比 pow 快

1.1.2 法线分布函数 D

D 是法线分布函数(Normal Distribution Function),用于描述微平面法线分布。Cook-Torrance[2] 使用的是 Beckmann 分布[6]

k _ { \mathrm { spec } } = \frac { \exp \left( - \tan ^ { 2 } ( \alpha ) / m ^ { 2 } \right) } { \pi m ^ { 2 } \cos ^ { 4 } ( \alpha ) } , \alpha = \arccos ( N \cdot H )

Walter et al.[7] 给出了“更好”的法线分布函数 GGX/Trowbridge-Reitz,这也被迪士尼和虚幻[5]所采用

D ( \mathbf { m } ) = \frac { \alpha _ { g } ^ { 2 } \chi ^ { + } ( \mathbf { m } \cdot \mathbf { n } ) } { \pi \cos ^ { 4 } \theta _ { m } \left( \alpha _ { g } ^ { 2 } + \tan ^ { 2 } \theta _ { m } \right) ^ { 2 } }

1.1.3 几何函数 G

G 是几何函数(Geometry Function),用于描述几何遮蔽(Geometry Obstruction)和几何阴影(Geometry Shadowing)(值的范围是 0-1)。

图片来源:https://learnopengl.com/PBR/Theory

一般还用 Smith shadowing-masking approximation 来简化[8]

G ( \mathbf { i } , \mathbf { o } , \mathbf { m } ) \approx G _ { 1 } ( \mathbf { i } , \mathbf { m } ) G _ { 1 } ( \mathbf { o } , \mathbf { m } )

G_1 可由 D 导出,故一般两者是相对应的。

GGX 对应的 G_1 比较复杂,考虑简化一下。

之前提到的 Schlick[4] 的简化公式是基于 Beckmann 的

G_1 ( v, m ) = \frac { v\cdot m } { (v \cdot m)(1 - k) + k } \quad \text { with } \quad k = \sqrt { \frac { 2 \alpha } { \pi } }

其中 \alpha=\text{roughness}^2

将 Beckmann 改为 GGX ,需要做如下改动[5]

k = \frac{\alpha}{2}

另外,对于解析光源(点光源、方向光、聚光灯),我们将粗糙度重映射[5]

\frac{(\text{roughness}+1)}{2}

注意,只对于解析光源进行该操作,本文只涉及 IBL,所以不需要考虑这个。完整的 PBR 计算里才会用到。

1.2 金属工作流 Metallic Workflow

之前提到的 Schlick 版本的菲涅尔公式,只能用于描述电介质(如玻璃,水)。对于导体(如金属),要使用另一个菲涅尔公式。这就很麻烦。金属工作流则使用参数金属度(metalness)来确定金属性。为了能让菲涅尔简化公式能同时用于金属和非金属,我们可用金属性来对 F_0 进行插值,非金属的 F_0 默认近似为 0.04,金属的 F_0 则为表面颜色。

vec3 F0 = vec3(0.04);
F0      = mix(F0, surfaceColor.rgb, metalness);

1.3 基于图像的照明(Image Based Lighting, IBL)

IBL 相当于一个无限大的球面光源在照射场景。辐射度由 environment map 决定。

图片来源:https://learnopengl.com/PBR/IBL/Diffuse-irradiance

即 environment map 给出了函数

L_i(\omega_i)=\text{texture}(\omega_i)

如果使用的是 2D 的 environment map,那要将方向向量 \omega_i 转为纹理坐标,或者先将该 environment map 转化成 cubemap。

2. 离线 IBL

这里不对离线渲染算法 path tracing[1] 做介绍

在 path tracing 中实现 IBL,关键就是对其进行重要性采样。

pbrt[9] 给出的实现是通过 2D 的二分查找来采样像素,时间复杂度为 O(\log w \log h)

还可以更快,我下边给出时间复杂度为 O(1) 的采样算法。

2.1 别名法 Alias Method

别名法[10]可以实现对有限的离散随机变量快速采样,时间复杂度为 O(1) 。而 environment map 从像素的视角上看,就是一个离散的随机变量。所以我们可以用别名法来加速这个采样过程。

一个像素的采样概率由亮度和像素在球面上的面积(正比于 \sin\theta )决定,如下

p_\text{pixel}(i,j) \varpropto \text{texture}(i,j).y()*\sin\theta

其中 y() 是亮度, \theta 用像素中点对应的 \theta 来近似,经过归一化后即可得到它们的具体概率值。

算法请参考[10],具体实现如下

AliasMethod.cppgithub.com

2.2 采样

首先用别名法采样离散的像素。然后我们继续在像素内均匀采样。这样纹理坐标的概率为

p_\text{texcoord}(u,v)=\frac{p_{\text{pixel}}(i_u,j_v)}{\frac{1}{w*h}}=w*h*p_{\text{pixel}}(i_u,j_v)

在 path tracing 中,我们需要的是 p(\omega_i) ,其与 p_\text{texcoord}(u,v) 的关系[9]

p ( \omega ) = \frac { p ( \theta , \phi ) } { \sin \theta } = \frac { p_\text{texcoord} ( u , v ) } { 2 \pi ^ { 2 } \sin \theta }

具体实现如下

InfiniteAreaLight.cppgithub.com

3. 实时 IBL

对于 Cook-Torrance 模型,我们要解的方程为

\begin{aligned} L _ { o } \left( p , \omega _ { o } \right)  &= \int _ { \Omega } f _ { r } \left( p , \omega _ { i } , \omega _ { o } \right) L _ { i } \left( p , \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d} \omega _ { i }\\  &=\int _ { \Omega }(k_d f_d + f_s) L _ { i } \left( p , \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d} \omega _ { i }\\  &=\int_\Omega k_d f_d L_i \left( p , \omega _ { i } \right) n \cdot \omega _ { i }\mathrm{d}\omega_i +\int_\Omega f_s L_i \left( p , \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d}\omega_i\\  &= L_d(p.\omega_o) + L_s(p,\omega_o) \end{aligned}

方程拆成了两部分,左边是漫反射项,右边是镜面项。对于 IBL,位置假设不变,所以有

L(p,\omega_o)=L(\omega_o)

则方程简化为

L _ { o } \left(p,\omega _ { o } \right) =  L_d(\omega_o) + L_s(\omega_o)

离线渲染中我们可以用蒙特卡洛积分法去求解这个积分值,但计算量很大。

所以实时渲染中面临的问题就是如何快速计算。

主要思路就是预计算,把复杂的积分都先算好。我们会分别预计算漫反射项和镜面项,最终在实时渲染中只需通过简单的纹理采样即可得到结果。

3.1 irradiance map

漫反射项为

\begin{aligned} L_d(\omega_o) &=\int_\Omega k_d f_d L_i \left( p , \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d}\omega_i\\ &=\frac{\rho}{\pi}\int_\Omega k_dL_i \left( p , \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d} \omega _ { i }\\ \end{aligned}

因为 k_d=(1-F(n\cdot \omega_i))(1-\text{metalness}),因此我们不能将其移出积分式。

为了简化,我们对其做了近似

k_d\approx (1-F_\text{roughness}(n\cdot\omega_o))(1-\text{metalness})\overset{\Delta}{=}k_d^*

其中 F_\text{roughness}(n\cdot\omega_o)

\begin{aligned} F_\text{roughness}(n\cdot\omega_o) &=F_0+(\Lambda-F_0)(1-n\cdot\omega_o)^5\\ \Lambda&=\max\{1-\text{roughness},F_0\}\\ \end{aligned}

相比于 Schlick 的菲涅尔公式,就是 \Lambda 替代了原本的 1

这样,漫反射项可简化为

\begin{aligned} L_d(\omega_o) &=\frac{\rho}{\pi}\int_\Omega k_d L_i \left( p , \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d} \omega _ { i }\\ &\approx k_d^*\frac{\rho}{\pi}\int_\Omega L_i \left( p , \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d} \omega _ { i }\\ \end{aligned}

现在的积分项只取决于法向 n (半球 \Omega 是法向所在的半球),注意,并不是取决于视线方向 \omega_o

因此我们可以预计算这个积分值,得到一个 cubemap,称为 irradiance map

图片来源:https://learnopengl.com/PBR/IBL/Diffuse-irradiance

积分方法就是蒙特卡洛积分,我们可以简单的在半球面上均匀采样。

对于均匀采样的蒙特卡洛积分,公式为

\begin{aligned} \frac{1}{\pi}\int_\Omega L_i \left( p , \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d} \omega _ { i }  &=\frac{1}{\pi}\int_0^{2\pi}\int_0^{\frac{\pi}{2}} L_i \left( p , \omega _ { i } \right) n \cos\theta \sin\theta \mathrm{d}\theta\mathrm{d}\phi\\  &\approx \frac { 1 } { \pi } \frac { 1 } { N_1N_2 } \sum _ { i } ^ { N_1 } \sum _ { j } ^ { N_2 } \frac{L _ { i } \left( p , \phi _ { j } , \theta _ { i } \right) \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 } \left( p , \phi _ { j } , \theta _ { i } \right) \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 } \left( p , \phi _ { j } , \theta _ { i } \right) \cos ( \theta_i ) \sin ( \theta_i )\\  \end{aligned}

注意我们把 \frac{1}{\pi} 也放到了 irradiance map 中。这样在使用它时我们就不必除以π。

我们还可以用余弦为权重的半球采样

\begin{aligned}  \frac{1}{\pi}\int_\Omega L_i \left( p , \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d} \omega _ { i }  &\approx\frac{1}{\pi}\frac{1}{N}\sum_i^N\frac{L(p,\omega_i)(n\cdot\omega_i)}{\frac{n\cdot\omega_i}{\pi}}\\  &= \frac{1}{N}\sum_i^N L(p,\omega_i)\\  \end{aligned}

使用 irradiance map 时根据法向 n 采样,并且乘上 k_d \rho 就可以得到漫反射项。

这个预计算可以 CPU 算,当然也可以用 GPU 来算,这里附上 glsl 的实现

Ubpa/RenderLabgithub.com图标

3.2 speular

镜面项为

L_s(\omega_o)=\int_\Omega f_s(\omega_i,\omega_o) L_i \left( \omega _ { i } \right) n \cdot \omega _ { i } d \omega _ { i }

这个值取决于 n\omega_o ,此外还有 BRDF 中的粗糙度和 F_0 (金属度和 albedo 决定),总共有 9 个维度(方向向量用球面坐标表示的话就是 2 个维度),非常多,所以我们没法像漫反射项那样直接用一个 cubemap 解决问题。

首先考虑最简单暴力的方法。我们直接在实时渲染中用蒙特卡洛积分来求解

\int_\Omega f_s(\omega_i,\omega_o) L_i \left( \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d}\omega_i \approx \frac{1}{N}\sum_j^N\frac{f_s(\omega_i^{(j)},\omega_o)L_i(\omega_i)}{p(\omega_i^{(j)})}

但是为了快速计算,我们没法过多采样。但采样较少的话误差又很大。所以这个方法在实时渲染中比较糟糕。

虚幻给出了近似解决的方法 spilt sum approximation[4]

\frac{1}{N}\sum_j^N\frac{f_s(\omega_i^{(j)},\omega_o)L_i(\omega_i^{(j)})}{p(\omega_i^{(j)})}  \approx  \left(\frac {\sum_k^N L_i(\omega_i^{(k)})W(\omega_i^{(k)})} {\sum_k^N W(\omega_i^{(k)})} \right)  \left(\frac{1}{N}\sum_j^N\frac{f_s(\omega_i^{(j)},\omega_o)(n \cdot \omega _ { i }^{(j)})}{p(\omega_i^{(j)})}\right)

其中权重函数 W(\omega_i)=n \cdot \omega_i

虚幻没有给出具体解释,解释可以参考[11],但写的也不是很好。

这里我给出自己的详细解释。

首先我们想把 L 从积分中去除

\begin{aligned}  \int_\Omega f_s(\omega_i,\omega_o) L_i \left( \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d}\omega_i &=L_c(\omega_o) \int_\Omega f_s(\omega_i,\omega_o) n \cdot \omega _ { i } d \omega _ { i }\\  L_c(\omega_o) &= \frac {\int_\Omega f_s(\omega_i,\omega_o) L_i \left( \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d}\omega_i} {\int_\Omega f_s(\omega_i,\omega_o) n \cdot \omega _ { i } d \omega _ { i }}\\  \end{aligned}

接着用法线分布函数来进行重要性采样,有 \begin{aligned} L_c(\omega_o) &\approx \frac {\frac1N \sum_k^N \frac {D(\omega_h^{(k)})F(\omega_o\cdot\omega_h^{(k)})G(\omega_o,\omega_i^{(k)})L_i(\omega_i^{(k)})(n\cdot\omega_i^{(k)})} {4(\omega_o\cdot n)(\omega_i^{(k)} \cdot n)p(\omega_i^{(k)})}}  {\frac1N \sum_k^N \frac {D(\omega_h^{(k)})F(\omega_o\cdot\omega_h^{(k)})G(\omega_o,\omega_i^{(k)})(n\cdot\omega_i^{(k)})} {4(\omega_o\cdot n)(\omega_i^{(k)} \cdot n)p(\omega_i^{(k)})}}\\  &= \frac {\sum_k^N \frac {D(\omega_h^{(k)})F(\omega_o\cdot\omega_h^{(k)})G(\omega_o,\omega_i^{(k)})L_i(\omega_i^{(k)})} {4(\omega_o\cdot n)\frac{D(\omega_h^{(k)})(\omega_h^{(k)}\cdot n)}{4(\omega_o\cdot\omega_h^{(k)})}}}  {\sum_k^N \frac {D(\omega_h^{(k)})F(\omega_o\cdot\omega_h^{(k)})G(\omega_o,\omega_i^{(k)})} {4(\omega_o\cdot n)\frac{D(\omega_h^{(k)})(\omega_h^{(k)}\cdot n)}{4(\omega_o\cdot\omega_h^{(k)})}}}\\   &= \frac {\sum_k^N \frac {F(\omega_o\cdot\omega_h^{(k)})G(\omega_o,\omega_i^{(k)})L_i(\omega_i^{(k)})(\omega_o\cdot\omega_h^{(k)})} {\omega_h^{(k)}\cdot n}}  {\sum_k^N \frac {F(\omega_o\cdot\omega_h^{(k)})G(\omega_o,\omega_i^{(k)})(\omega_o\cdot\omega_h^{(k)})} {\omega_h^{(k)}\cdot n}}\\  &\approx \frac {\sum_k^N \frac {F(\omega_o\cdot\omega_h^{(k)})G_1(\omega_i^{(k)})L_i(\omega_i^{(k)})(\omega_o\cdot\omega_h^{(k)})} {\omega_h^{(k)}\cdot n}}  {\sum_k^N \frac {F(\omega_o\cdot\omega_h^{(k)})G_1(\omega_i^{(k)})(\omega_o\cdot\omega_h^{(k)})} {\omega_h^{(k)}\cdot n}}\\  \end{aligned}

Unity 发现 F 对结果影响不大[12],就去掉了 F

为什么 F 没多大影响,我想有两点原因

  • 对于光滑情形, \boldsymbol{\omega}_h 接近 \mathbf{n} ,所以 F 基本是定值,分子分母可约去
  • 对于非光滑情形,L 已经变得很粗略了,所以做这种近似也影响不大。

去掉 F 后,公式为

L_c(\omega_o) \approx \frac {\sum_k^N \frac {G_1(\omega_i^{(k)})L_i(\omega_i^{(k)})(\omega_o\cdot\omega_h^{(k)})} {\omega_h^{(k)}\cdot n}}  {\sum_k^N \frac {G_1(\omega_i^{(k)})(\omega_o\cdot\omega_h^{(k)})} {\omega_h^{(k)}\cdot n}}

镜面的 BRDF f(\omega_o,\omega_i) 大致会在反射反向 R=\text{reflect}(-\omega_o,\mathbf{n}) 附近有值

图片来源:pbrt, p 509

假设不同方向入射,这个波瓣形状变化不大,则有

f(\omega_o,\omega_i(\mathbf{n}),\mathbf{n}) \approx f(\mathbf{R},\omega_i(\mathbf{R}),\mathbf{R})

其中 \omega_i(\mathbf{n})\mathbf{n} 的夹角等于 \omega_i(\mathbf{R})\mathbf{R} 的夹角

示意图大概为

我忘画 wi 了,自行脑补下...

因此有

L_c(\omega_o)  = \frac {\int_{\Omega(\mathbf{n})} f_s(\omega_i,\omega_o) L_i \left( \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d}\omega_i} {\int_{\Omega(\mathbf{n})} f_s(\omega_i,\omega_o) n \cdot \omega _ { i } \mathrm{d} \omega _ { i }}  \approx \frac {\int_{\Omega(\mathbf{R})} f_s(\omega_i,\mathbf{R}) L_i \left( \omega _ { i } \right) \mathbf{R} \cdot \omega _ { i } \mathrm{d}\omega_i} {\int_{\Omega(\mathbf{R})} f_s(\omega_i,\mathbf{R}) \mathbf{R} \cdot \omega _ { i } \mathrm{d} \omega _ { i }}\overset{\Delta}{=}L_c^*(\mathbf{R})

其中 \boldsymbol{\omega}_o=\mathbf{n}=\mathbf{R},因此在之前化简的基础上做如上近似,可把点积约去

L_c^*(\mathbf{R}) \approx \frac {\sum_k^N L(\omega_i^{(k)})G_1(\omega_i^{(k)})} {\sum_k^N G_1(\omega_i^{(k)})}

可见,虚幻所用的权重函数 W(\omega_i)=n \cdot \omega_i 相当于 G_1(\omega_i)

虚幻之所以觉得 W(\omega_i)=n \cdot \omega_i 更好,其比较对象是 W(\omega_i)=1 。根据推导,我们应该使用 G_1(\omega_i) 。但刚好两者性质差不多(如值域相同,变化趋势相似),所以虚幻这么做确实有其合理性。

现在,整个积分变为

\int_\Omega f_s(\omega_i,\omega_o) L_i \left( \omega _ { i } \right) n \cdot \omega _ { i } \mathrm{d}\omega_i =L_c^*(\mathbf{R}) \int_\Omega f_s(\omega_i,\omega_o) n \cdot \omega _ { i } d \omega _ { i }

总之,以上过程做了很多简化,很不物理。

3.2.1 pre-filter environment map

现在我们要计算

L_c^*(\mathbf{R}) \approx\frac {\sum_k^N L_i(\omega_i^{(k)})(n \cdot \omega_i^{(k)})} {\sum_k^N (n \cdot \omega_i^{(k)})}

之前做了镜面波瓣的相近假设,这样忽略了 \omega_o\mathbf{n} 的具体方向,只考虑了 R 。这样会引入误差,最明显的现象就是在掠射角(grazing angles)看表面时没法得到拖长的反射[5]

图片来源:https://learnopengl.com/PBR/IBL/Specular-IBL

L_c^*(\mathbf{R}) 有两个变量,为粗糙度和 \mathbf{R} 。对于一个特定的粗糙度,变量就只剩 R ,因此我们可以像 irrandiance map 那样,预计算得到一个 cubemap。我们均匀的取粗糙度为 0, 0.25, 0.5, 0.75, 1.0,这样得到 5 个 cubemap。实时渲染时,就用 R 和粗糙度进行三线性插值。

这个预计算有时需要实时更新,那么快速计算就十分重要。这个预计算关键的就是计算积分。为了加速收敛,我们可以用重要性采样,但还是需要不少的样本。Krivanek[13] 的 Pre-filtered importance sampling 可以减少样本数量,收敛提升明显,只引入了一小些偏差。

图片来源:Moving Frostbite to Physically Based Rendering

这样我们得到的 5 个 cubemap,就是 pre-filtered environment map

图片来源:https://learnopengl.com/PBR/IBL/Specular-IBL

这个预计算可以 CPU 算,当然也可以用 GPU 来算,这里附上 glsl 的实现

prefiltergithub.com

从公式上看,这个 cubemap 是 L_c^*(\mathbf{R}),要用 R 采样,虽然在预计算时 R=\omega_o=n,但其实是 \omega_o,\mathbf{n}\leftarrow R ,即 \mathbf{R} 才是决定因素。

3.2.2 BRDF LUT

上边说到,积分改为了

\int_\Omega f_s(\omega_i,\omega_o) L_i \left( \omega _ { i } \right) n \cdot \omega _ { i } d \omega _ { i }\mathrm{d}\omega_i =L_c^*(\mathbf{R}) \int_\Omega f_s(\omega_i,\omega_o) n \cdot \omega _ { i } d \omega _ { i }

其中的 L_c^*(\mathbf{R}) 我们已经估计好了,现在就剩 \int_\Omega f_s(\omega_i,\omega_o) n \cdot \omega _ { i } d \omega _ { i }

这个积分含有变量 \omega_onF_0 (金属度和 albedo 决定)和粗糙度。

由于我们使用的 GGX 是各向同性的,所以我们只需要知道 \omega_on 的夹角 \theta 即可。

这样我们还剩三个变量 \thetaF_0 (金属度和 albedo 决定)和粗糙度。因此我们可以考虑将 F_0 移出积分式

\begin{aligned} \int _ { \Omega } f _ { s } \left(\omega _ { i } , \omega _ { o } \right) n \cdot \omega _ { i } d \omega _ { i }   &= \int _ { \Omega } f _ { s } \left(\omega _ { i } , \omega _ { o } \right) \frac { F \left( \omega _ { o } , h \right) } { F \left( \omega _ { o } , h \right) } n \cdot \omega _ { i } d \omega _ { i }\\  &= \int _ { \Omega } \frac { f _ { s } \left(\omega _ { i } , \omega _ { o } \right) } { F \left( \omega _ { o } , h \right) } \left( F _ { 0 } + \left( 1 - F _ { 0 } \right) \left( 1 - \omega _ { o } \cdot h \right) ^ { 5 } \right) n \cdot \omega _ { i } d \omega _ { i }\\  &=F_0 \int _ { \Omega } \frac { f _ { s } \left(\omega _ { i } , \omega _ { o } \right) } { F \left( \omega _ { o } , h \right) }(1-(1-\omega_o\cdot h)^5) n \cdot \omega_i d\omega_i\\ &\quad + \int _ { \Omega } \frac { f _ { s } \left(\omega _ { i } , \omega _ { o } \right) } { F \left( \omega _ { o } , h \right) }(1-\omega_o\cdot h)^5 n \cdot \omega_i d\omega_i\\  &=F_0*\text{scale}+\text{bias} \end{aligned}

我们将积分拆成了两部分 scale 和 bias,它们都消去了 F ,因此变量只剩 \theta 和粗糙度,为了方便,使用 \cos\theta 做变量(这样,这些变量的范围都是 [0,1] )。

scale 和 bias 可以用根据法线分布函数进行重要性采样的蒙特卡洛积分来解决。

\begin{aligned} \text{scale}  &\approx \frac{1}{N}\sum_k^N\frac {\frac{D(\omega_h^{(k)})G(\omega_o,\omega_i^{(k)})}{4(\omega_o\cdot n)(\omega_i^{(k)}\cdot n)}(1-(1-\omega_o\cdot\omega_h^{(k)})^5)(\omega_i^{(k)} \cdot n)} {p(\omega_i^{(k)})}\\  &= \frac{1}{N}\sum_k^N\frac {\frac{D(\omega_h^{(k)})G(\omega_o,\omega_i^{(k)})}{4(\omega_o\cdot n)(\omega_i^{(k)}\cdot n)}(1-(1-\omega_o\cdot\omega_h^{(k)})^5)(\omega_i^{(k)} \cdot n)} {\frac{D(\omega_h^{(k)})(\omega_i^{(k)}\cdot n)}{4(\omega_o\cdot\omega_h^{(k)})}}\\  &= \frac{1}{N}\sum_k^N\frac {G(\omega_o,\omega_i^{(k)})(\omega_o\cdot \omega_h^{(k)})(1-(1-\omega_o\cdot\omega_h^{(k)})^5)} {(\omega_o\cdot n)(\omega_i^{(k)}\cdot n)}\\ \end{aligned}

bias 也类似。

这两部分都可以用一个 2D 的纹理来表示预计算的结果,它们都是 1 维变量,可以放在同一张纹理里边,scale 放在红色通道, bias 放在绿色通道。这个纹理称为 LUT。

这个 LUT 是由 BRDF 决定的,所以确定的 BRDF 就有确定的 LUT。我们用的 BRDF 对应的 LUT 为

LUT

这个预计算可以直接用这张图。自己也可以用 CPU 算,当然也可以用 GPU 来算,这里附上 glsl 的实现

brdf LUTgithub.com图标

3.3 综合

最终的近似公式为

L(\omega_o)=k_d^*\rho\ \text{IrradianceMap}(\mathbf{n}) +L_c^*(\mathbf{R})(F_0*\text{scale}+\text{bias})

我们拿到了三个纹理 irrandiance map、pre-filtered environment map 和 LUT。我们可以用它们实时地计算 L(\omega_o)

核心代码如下

vec3 F0 = mix(vec3(0.04), albedo, metalness);
vec3 F = F_roughness(wo, normal, F0, roughness);
vec3 kS = F;
vec3 kD = (1 - metalness) * (vec3(1) - kS);

vec3 irradiance = texture(irradianceMap, norm).rgb;

vec3 diffuse = irradiance * albedo;

vec3 R = reflect(-wo, norm);
const float MAX_REFLECTION_LOD = 4.0;
// 用 R 来采样
vec3 prefilteredColor = textureLod(prefilterMap, R, roughness * MAX_REFLECTION_LOD).rgb;
vec2 scale_bias = texture(brdfLUT, vec2(max(dot(norm, wo), 0.0), roughness)).rg;
vec3 specular = prefilteredColor * (F0 * scale_bias.x + scale_bias.y);

vec3 ambient = kD * diffuse + specular;

结果

实时 IBL

离线 IBL

对比图

(之后找模型再渲染几个图来

总结

IBL 的效果真的很好,用 metalness 和 roughness 就可以控制材质,非常方便。

但这些只是近似结果,最大的误差来自于 spilt sum approximation,也就是说,反射效果会有较大的差异。

还有很多很多的变化形式,但总的思想是不变的,都是为了实时渲染时能快速求解渲染方程。


按照国际惯例,最后贴一下我最近写的渲染引擎 RenderLab

RenderLabgithub.com
RenderLab

包含了离线渲染,实时渲染,编辑器的功能,未来还会作为游戏引擎。


本文是我的毕业纪念

因此本大四老咸鱼想在这里简单回忆一下过去

本着对 MMD 的兴趣,大三下选修了图形学

woc,课程作业量超大,一周一个编程实验,总共有 9 个作业,肝的难受

// 附上作业链接 github/Ubpa/CG_LGLiu,感兴趣的同学可以了解下

而我之前很少写代码,学了这门课我才知道,编程真有意思,图形学真有趣

// 没错,我之前不怎么码代码,太蔡了,QAQ,github 一片灰

最后,课程项目是随便做一个东西,我就尝试做了低仿崩崩崩(抄袭,苟命)

崩坏3 仿,从零开始的 7 天 Unitywww.bilibili.com
游戏截图

当时的我还经验浅薄(当然现在也是),什么都不懂还必须得强行上,勉强交上了一份答卷

当时我本打算读研选 AI 方向,但一番抉择过后还是选择图形学吧。

未来两至三年,继续学习,磨练技术(秃头√)


感谢大家的阅读

参考

  1. ^abKajiya, James T. (1986),"The rendering equation"(PDF),Siggraph 1986: 143–150. http://www.cse.chalmers.se/edu/year/2011/course/TDA361/2007/rend_eq.pdf
  2. ^abR. Cook and K. Torrance. "A reflectance model for computer graphics". Computer Graphics (SIGGRAPH '81 Proceedings), Vol. 15, No. 3, July 1981, pp. 301–316. http://inst.eecs.berkeley.edu/~cs283/sp13/lectures/cookpaper.pdf
  3. ^https://en.wikipedia.org/wiki/Fresnel_equations
  4. ^abSchlick, C. (1994)."An Inexpensive BRDF Model for Physically-based Rendering"(PDF).Computer Graphics Forum.13(3): 233–246. http://www.cs.virginia.edu/~jdl/bib/appearance/analytic%20models/schlick94b.pdf
  5. ^abcdeKaris B, Games E. Real shading in unreal engine 4[J]. Proc. Physically Based Shading Theory Practice, 2013, 4. https://cdn2.unrealengine.com/Resources/files/2013SiggraphPresentationsNotes-26915738.pdf
  6. ^Petr Beckmann, André Spizzichino, The scattering of electromagnetic waves from rough surfaces, Pergamon Press, 1963, 503 pp
  7. ^Walter B, Marschner S R, Li H, et al. Microfacet models for refraction through rough surfaces[C]//Proceedings of the 18th Eurographics conference on Rendering Techniques. Eurographics Association, 2007: 195-206.
  8. ^SMITH B. G.: Geometrical shadowing of a random rough surface. IEEE Trans. on Antennas and Propagation (1967), 668–671
  9. ^abFrom Theory to Implementation, pp. 845-851 https://www.pbrt.org/
  10. ^abhttps://en.wikipedia.org/wiki/Alias_method
  11. ^Alternative Take on the Split Sum Approximation for Cubemap Pre-filtering https://zero-radiance.github.io/post/split-sum/
  12. ^https://github.com/Unity-Technologies/ScriptableRenderPipeline/blob/680100c6b7f638fcd4d2c05e8e4ac32cf0d7338c/com.unity.render-pipelines.core/ShaderLibrary/ImageBasedLighting.hlsl#L531
  13. ^Křivánek J, Colbert M. Real‐time shading with filtered importance sampling[C]//Computer Graphics Forum. Wiley/Blackwell (10.1111), 2008, 27(4): 1147-1154. http://dcgi.felk.cvut.cz/publications/2008/krivanek-cgf-rts
编辑于 2019-06-25

文章被以下专栏收录

    图形学 Repeater

    本专栏收集高品质游戏开发方面的原创文章,并不局限于技术文章,尤其欢迎美术类策划类管理类文章。投资市场营销类暂不接受。欢迎投稿但是否采纳全看运气。如无特别声明,本专栏作品采用知识共享署名 4.0 国际许可协议进行许可。