PBR理论基础3.1:基于图像的光照(下)

 

接上文:PBR理论基础3:基于图像的光照(上)

四、蒙特卡洛积分

蒙特卡洛积分在前篇计算漫反射项时就已经应用过了

对于定积分  A=\int_{a}^{b} f(x) d x  ,如果无法简单求得原函数,一个尝试去求得其近似解的方法就是蒙特卡洛积分:

  • 在积分区间 [a, b] 上进行均匀采样得到样本 X_{1}, \cdots, X_{N},就可以求和得到 F(N) \approx \frac{b-a}{N} \sum_{i=1}^{N} f\left(X_{i}\right),随着样本数量 N 的增大,其结果也会逐渐收敛到积分的精确解

这个方法和黎曼和的求法非常类似:\frac{b-a}{N} 就是图中矩形的宽,而 f\left(X_{i}\right) 就是对应矩形的高

4.1 概率密度函数

不过上面的计算没有考虑一个东西,那就是对于样本 X_{1}, \cdots, X_{N},它不一定是均匀分布的。对于前面章节漫反射的预计算,我们可以认为任意方向的光源对出射向量 \omega_o 的贡献都是一致的,这就意味着以下两点

  • 积分仅依赖于 \omega_i,而不依赖于 \omega_o
  • 生成的样本均匀分布,就意味着随着样本数量的不断增加,我们最终将收敛到积分的精确解(无偏估计)
对于完美的漫反射模型,已知确定的出射光(入射光)时,其入射光(出射光)在半球中均匀分布

但是镜面反射就不一定了,任意方向的光源对出射向量 \omega_o 的贡献并不一致:极端情况对于绝对光滑的平面,有且只有一条光线对 \omega_o 有 1.0 的贡献,其它光线对其贡献都为 0.0。这时如果再对入射光方向 \omega_i 进行随机采样,必然会得出错误的结果,因为你无论怎么采样,都几乎不可能正好命中到那条有 100% 贡献的入射光,这时得出的辐射率结果就会为 0,毕竟这些向量都与积分无关。就算你运气好,其中一次采样刚好命中了那条有 100% 贡献的入射光,你会发现随着采样次数变多,得出来的结果也会越加趋近于 0,也就是收敛到一个错误值

怎么算都不对,这就是因为没有考虑概率密度函数,由于这个那个例子的有偏性,随机采样可能永远不会收敛到精确解,或者说是会收敛到一个错误解

那么这个问题怎么破呢?这就要引用概率密度函数:

众所周知:

  1. 假设某一连续型随机变量 X 的样本空间为 D,其概率密度分布函数为 p(x),数学期望就为 E(X)=\int_{D} x p(x) d x
  2. 若另一连续随机变量 Y 满足 Y = f(x),则 Y 的数学期望为 E(Y)=\int_{D} f(x) p(x) d x

这样对于上面绝对光滑平面镜面反射的例子,就可以认为对于那条有 100% 贡献的入射光 \omega_i,其概率密度函数 p(x) 满足 p(x) = 1,其它无贡献的入射光 p(x) = 0。换成用概率密度函数的蒙特卡洛积分描述就是:

  • 在积分区间 [a, b] 上进行采样得到样本 X_{1}, \cdots, X_{N}(这个采样要根据概率密度函数指定规则,并非均匀采样),最终积分结果就可近似为:F_{N}=\frac{1}{N} \sum_{i=1}^{N} \frac{f\left(X_{i}\right)}{p\left(X_{i}\right)}

这里也可以发现前面推导的 F(N) \approx \frac{b-a}{N} \sum_{i=1}^{N} f\left(X_{i}\right) 其实就是上面的一种特殊情况:均匀分布,即 p(x) 满足  p(x)=\left\{\begin{array}{ll} \frac{1}{b-a} & a \leq x \leq b \\ 0 & \text {other } \end{array}\right.

4.2 重要性采样

还有一个问题需要解决,那就是对于非均匀分布的样本,如何对其采样才可得到准确的收敛结果呢?

可以再举一个例子:现在要统计全国大学生的平均身高,有以下2个方案:

  1. 直接挨家挨户调查全国所有大学生每一个人的身高,然后求其平均
  2. 随便进入一所大学,统计这所大学里面所有学生的身高,然后求其平均作为近似解

可见,选择方案②确实是一个很不错的选择,它就是随机采样的一个案例,但是现在有这样一个问题:这所大学是文科大学,也因此男女比例是 1:2,要知道女性的平均身高是要低于男性的,这样若只是单纯的随机抽样必然会得出错误的结果。对于这个问题,该如何设计采样方案,才能使得最终的结果尽可能的客观真实呢?

没错这就是一道初中题,但是他和我们现在要解决的问题却是一致的

我们已知 f(x) 在概率分布 p(x) 的期望为  E=\int_{D} f(x) p(x) d x,这样就是要找一个采样分布函数 q(x) 满足 E=\int_{D} f(x) p(x) d x=\int_{D} q(x) \frac{p(x)}{q(x)} f(x) d x ,也就是函数 f(x) 在概率分布 p(x) 的期望可以被看作函数 \frac{p(x)}{q(x)} f(x) 在概率分布 q(x) 上的分布,并且这个 q(x) 是支持均匀采样的,这样就可以得到:

  • 在分布区间 q(x) 上进行采样得到样本 x_{1}, \cdots, x_{N},最终积分结果就可近似为:F_{N}=\frac{1}{N} \sum_{i=1}^{N} \frac{p\left(x_{i}\right)}{q\left(x_{i}\right)} f\left(x_{i}\right),其中样本 x_{1}, \cdots, x_{N} 为实际均匀样本 X_{1}, \cdots, X_{N} 关于 q(x) 的映射 q(X_{1}), \cdots, q(X_{N})\frac{p\left(x_{i}\right)}{q\left(x_{i}\right)} 就是重要性权重

回到上面那个统计身高的例子,假设男生 =1, 女生 =0,那么就可以得到对应的概率密度函数 p(x)=\left\{\begin{array}{ll} 1/3 & x=1 \\ 2/3 & x=0 \end{array}\right.,一个计算 q(x) 的方向是,尽可能保证 \frac{p(x)}{q(x)} f(x) 和需要积分的函数 f(x) 尽可能接近,此时 q(x)=\left\{\begin{array}{ll} 1 & x=1 \\ 2 & x=0 \end{array}\right.,完美,用简单的话来讲就是:每次抽样到男生的时候,都对他的身高加权,权值为女生的两倍。又或者说将男生分成一模一样的两份,和女生混合后作为样本从中抽取求平均

 

 

五、镜面IBL预滤波项

对于反射方程的镜面反射项,可以进行卷积近似:

\small L_o(p,\omega_o) = \int\limits_{\Omega} f_r(p, \omega_i, \omega_o) L_i(p,\omega_i) n \cdot \omega_i d\omega_i\\ = \int\limits_{\Omega}L_i(p,\omega_i) d\omega_i * \int\limits_{\Omega} f_r(p, \omega_i, \omega_o) n \cdot \omega_i d\omega_i

其中  \small \int\limits_{\Omega}L_i(p,\omega_i) d\omega_i  这一部分被称为预滤波项,和上一章的辐照度图一样,这里也可以用一样的方式预处理,不过这次在确定入射光线 \small \omega_i 时,反射光线 \small \omega_o 的分布并不均匀(反过来也一样),并且同时还受表面的粗糙度影响

就如上图,确定入射光线 \small \omega_i 时,所有的反射光线 \small \omega_o 构成的形状称为镜面波瓣,随着粗糙度的增加,镜面波瓣的大小增加;随着入射光方向不同,形状也会发生变化,因此对于影响光线分布的两个因素:粗糙度(Roughness) 和光方向,就需要做额外的考量

5.1 低差异序列(Low Discrepancy Sequence)

确定粗糙度时,由于不同方向上光线密度不同,需要引用上一节中的重要性采样来确保积分结果的正确收敛以及得到更快的收敛速度,当然还可以使用一个黑科技,那就是低差异序列:

低差异序列很简单,可以理解为完全随机采样和均分采样的一个折中,一个经典的低差异序列是 Hammersley 序列,关于低差异序列和随机,体量不小可以后面开个专篇讲,这里暂时了解下就好了

  • 对于一个在 [0,1]^{n} 区域内的点集,任意选取一个空间中的区域 B,此区域内点的数量 A 和点集个数的总量 N 的比值和此区域的体积 \lambda_s 的差的绝对值的最大值,就是这个点集的差异(Discrepancy) ,即有  D_{N}(P)=\sup _{B \in J}\left|\frac{A(B)}{N}-\lambda_{s}(B)\right|,分布越均匀的点集,任意区域内的点集数量占点总数量的比例也会越接近于这个区域的体积
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
}
// ----------------------------------------------------------------------------
vec2 Hammersley(uint i, uint N)
{
    return vec2(float(i)/float(N), RadicalInverse_VdC(i));
}  

上面就是 Hammersley 序列的具体实现,除此之外在 openGL.cn 的网站上,也有 Hammersley 序列的非位运算版本

5.2 基于GGX的重要性采样

这一章中,了解了几个 NDF,其中一个重要的 GGX - NDF 如下:

G G X(n, h, \alpha)=\frac{\alpha^{2}}{\pi\left(\cos \theta^{2}\left(\alpha^{2}-1\right)+1\right)^{2}},使用倾斜角 \theta 和航向角 \phi 来表示就是 p(\theta, \phi)=\frac{\alpha^{2} \cos \theta \sin \theta}{\pi\left(\cos \theta^{2}\left(\alpha^{2}-1\right)+1\right)^{2}}

为了求得累积分布函数(Cumulative Distribution Function),也就是是概率密度函数的积分,需要先求得 \theta 和 \phi 的边缘概率密度函数:

  • p_{h}(\theta)=\int_{0}^{2 \pi} p_{h}(\theta, \phi) d \phi=\frac{2 \alpha^{2} \cos \theta \sin \theta}{\left(\left(\alpha^{2}-1\right) \cos ^{2} \theta+1\right)^{2}}
  • p_{h}(\phi)=\int_{0}^{\frac{\pi}{2}} p_{h}(\theta, \phi) d \theta=\int_{0}^{\frac{\pi}{2}} \frac{\alpha^{2} \cos \theta \sin \theta}{\pi\left(\left(\alpha^{2}-1\right) \cos ^{2} \theta+1\right)^{2}} d \theta = \frac{1}{2\pi}

然后就可以得到 CDF:

  • P_{h}(\phi)=\frac{\phi}{2 \pi}
  • P_{h}(\theta)=\int_{0}^{\theta} \frac{2 \alpha^{2} \text { costsint }}{\left(\left(\alpha^{2}-1\right) \cos ^{2} t+1\right)^{2}} d t =\frac{\alpha^{2}}{\left(\alpha^{2}-1\right)^{2} \cos ^{2} \theta+\left(\alpha^{2}-1\right)}-\frac{1}{\alpha^{2}-1}

定义前面依据低差异序列生成的随机向量 (x, y) 对应航向角 \phi 和倾斜角 \theta,那么就有

  • \phi=2 \pi x
  • \theta=\arccos \sqrt{\frac{1-y}{y\left(\alpha^{2}-1\right)+1}}

这基于逆变换算法:对于概率分布函数 p(x),以及它的累积分布函数 P(x),根据性质可以知道它的值域是一定是分布在 [0,1] 区间内的,这样通过它的反函数就可以实现输入 [0,1] 之间的随机数,得到符合概率密度分布的采样点。这个反函数也正对应着前一节的 q(x),openGL.cn 中的参考代码如下:

vec3 ImportanceSampleGGX(vec2 Xi, vec3 N, float a)
{
    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);
}

5.3 粗糙度(Roughness) 和光方向考量

好了,接下来就是如何处理不同粗糙度和光方向

粗糙度:关于粗糙度的处理非常的暴力,那就是对于每个特定粗糙度级别都算一遍卷积的结果,并把它们存在预滤波环境贴图的 mipmap 中,对于中间的粗糙度可以对两张 mipmap 进行插值

预滤波环境贴图,来源 openGL.cn

光方向:规定镜面反射方向 n、以及视角方向 \small \omega_o 相等,此时也满足 \small \omega_o = \small \omega_i,这样掠射角方向可能会得出不太正确的结果,但是也算是一种折中

由于低级的 mipmap 精度不够,在这个过程中需要注意的是:

  • 在立方体贴图边缘插值时,需要考虑相邻的立方体面
  • 在采样环境贴图的时,最好对同级的环境贴图 mipmap 进行采样,否则在采样精度不足的情况下可能会出现噪点

这里更多的是实践

 

 

六、预计算 BRDF 项

接下来就是计算 BRDF 项:

\small \int\limits_{\Omega} f_r(p, \omega_i, \omega_o) n \cdot \omega_i d\omega_i =\int\limits_{\Omega} \frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)} n \cdot \omega_i d\omega_i

它受到3个参数影响:\small n \cdot \omega_o、粗糙度 \alpha 以及菲涅尔系数 F_0,但是一个好消息是,由于不需要考虑光照辐射度 L(已经在前一节被抽离出来了),因此在这3个值确定时,得出来的结果也是一个定值

离线打表的手微微颤抖,这不就又可以打表嘛,直接开个3维常量数组,把结果算好打进去就可以了。但如果只需要存2维数据的话,是不是就可以将这些数据存到纹理中而不是作为常量呢?很显然是可以的!

6.1 BRDF项查找纹理

一张非常著名的 2D 查找纹理(Look Up Texture, LUT) 横空出世:其中

  • 横向对应 \small n \cdot \omega_o,也就是 cos\theta,范围刚好是 [0, 1],
  • 纵向对应粗糙度常数 \alpha,范围也是 [0, 1]
  • 纹理 R 通道为积分比例项,G 通道为积分偏差项(至于比例项和偏差项是什么可以看后面的推导)

那现在要解决的问题就是:对于3个输入参数:\small n \cdot \omega_o、粗糙度 \alpha 以及菲涅尔系数 F_0,需要先干掉一个参数 F_0,步骤如下:

  1. 尝试把 NDF 中的菲涅尔项 F 移出去:\int\limits_{\Omega} \frac{f_r(p, \omega_i, \omega_o)}{F(\omega_o, h)} F(\omega_o, h) n \cdot \omega_i d\omega_i
  2. 对于右边的 F(\omega_o, h) 将其展开:\int\limits_{\Omega} \frac{f_r(p, \omega_i, \omega_o)}{F(\omega_o, h)} (F_0 + (1 - F_0){(1 - \omega_o \cdot h)}^5) n \cdot \omega_i d\omega_i
  3. 用 \alpha 替换 {(1 - \omega_o \cdot h)}^5,可以得到:\int\limits_{\Omega} \frac{f_r(p, \omega_i, \omega_o)}{F(\omega_o, h)} (F_0 * (1 - \alpha) + \alpha) n \cdot \omega_i d\omega_i

最后拆掉,并且将常量移出积分,就可以得到:

\int\limits_{\Omega} f_r(p, \omega_i, \omega_o) n \cdot \omega_i d\omega_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 \\ = F_0 \int\limits_{\Omega} \frac{DG}{4(\omega_o \cdot n)(\omega_i \cdot n)}(1 - {(1 - \omega_o \cdot h)}^5) n \cdot \omega_i d\omega_i \\ + \int\limits_{\Omega} \frac{DG}{4(\omega_o \cdot n)(\omega_i \cdot n)} {(1 - \omega_o \cdot h)}^5 n \cdot \omega_i d\omega_i

其中   \int\limits_{\Omega} \frac{DG}{4(\omega_o \cdot n)(\omega_i \cdot n)}(1 - {(1 - \omega_o \cdot h)}^5) n \cdot \omega_i d\omega_i 为 F_0 的比例项\int\limits_{\Omega} \frac{DG}{4(\omega_o \cdot n)(\omega_i \cdot n)} {(1 - \omega_o \cdot h)}^5 n \cdot \omega_i d\omega_i 为偏差项,搞定!

6.2 完成IBL反射

按理说到这里就要结尾了,不过还是有一个东西要提一下

在计算阴影遮蔽函数(几何函数G)时,对于 G_{SchlickGGX}(n, v, k)=\frac{n \cdot v}{(n \cdot v)(1-k)+k},其一般在针对光照和针对 IBL 时其 k 项不同:这里应用第二个

  • 针对光照时:k_{direct} = \frac{(\alpha + 1)^2}{8}
  • 针对 IBL 的重映射(Remapping):k_{IBL} = \frac{\alpha^2}{2}

除此之外,生成的 LUT 应该是浮点数存储,一般都为16位精度浮点格式

……

OK,搞定所有的离线计算后,计算实时光照时就非常简单了:

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;   
vec2 envBRDF = texture(brdfLUT, vec2(max(dot(N, V), 0.0), roughness)).rg;
vec3 specular = prefilteredColor * (F * envBRDF.x + envBRDF.y);

//因为NDF已经乘过了菲涅耳系数,所以这里不需要KS
vec3 ambient = (kD * diffuse + specular) * ao; 

如果顺利的话,就可以得到一个非常不错的效果:

当然这些都是 learnopengl.cn 上的效果,后面会有一篇文章,尝试在 Unity 中实现 PBR(造轮子

 

参考文章:

 

  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值