环境光照(Environment Lighting):用于大型、复杂但强度较小的光源,包括天空光(Sky Light)和云散射的光、从场景中的超大型物体反弹的间接光照以及超大型 area light 的直接光照(往往是能包围整个场景的无限远的 area light)。这种光源具有物体位置无关的特性,即物体无论位于场景的哪里,所接受到的环境光照都是一样的(不考虑辐射度传输因素,例如几何相关的 visibility)。
其中,Ambient Light 是最简单的环境光照,相当于给整个场景所有物体提升同样的亮度(同一颜色)。
因为环境光往往是超大型的,而场景中物体相隔不远的情况下,就可以假设场景所有物体所受到的环境光是一样的,例如相隔一条街的物体能受到的天空光照几乎是一样的。当然考虑光的辐射度传输因素的话,不同物体所受到的环境光照将有所不同,但我们可以将其拆解成环境光信息和辐射度传输信息,环境光信息依然可以全场景共享,而辐射度传输信息可以每个物体各有一份。
基于图像的光照(Image Based Lighting,IBL),简单的说就是一类通过 环境贴图(Environment Map) 来保存环境光信息,从而实现基于物理的渲染方法。
为了表示来自四面八方的环境光信息,IBL 使用 Spherical Map 或者 Cube Map 的方式来存储:
IBL 可用于 diffuse/glossy/specular 物体(基本上包含大部分物体了)渲染:这取决于环境贴图,环境贴图分辨率越大,那么所能表示高频信息就越多,从而越适合 specular 物体的渲染;而分辨率越低,所表示的信息更多是低频信息, 则 diffuse 物体的渲染也足以满足,而且还能节省一定存储空间。
IBL 中最常见的算法便是基于 The Split Sum Approximation 的算法。
我们知道,一般的渲染方程如下:
Lr(x,ωr)=∫Ω+Li(x,ωi)fr(x,ωi,ωr)(n⋅ωi)dωi" role="presentation">Lr(x,ωr)=∫Ω+Li(x,ωi)fr(x,ωi,ωr)(n⋅ωi)dωi
在实际渲染的时候,我们当然可以使用蒙特卡洛方法实现该渲染方程,然而这样的开销是巨大的(每个shading point都要做多重采样,而且结果很容易是noisy的)。
为了进一步加快式子,这里有一个很经典的近似公式:
∫Ω+f(x)g(x)dx≈∫ΩGf(x)dx∫ΩG dx⋅∫Ω+g(x)dx" role="presentation">∫Ω+f(x)g(x)dx≈∫ΩGf(x)dx∫ΩG dx⋅∫Ω+g(x)dx
积分域 ΩG" role="presentation">ΩG :在原本 Ω+" role="presentation">Ω+ 的积分域范围内剔除掉 g(x)=0" role="presentation">g(x)=0 的地方而剩余的范围。
想让这个公式近似效果比较精确,那么需要满足以下一种或两种条件:
积分域 ΩG" role="presentation">ΩG 比较小 g(x)" role="presentation">g(x) 比较光滑,即变化不是很大而我们的观察是:
如果 BRDF 是 glossy/specular 的,那么它的 lobe 往往是花瓣状,即只有很小的积分域才能接受环境光。 如果 BRDF 是 diffuse 的,那么它的 lobe 往往是均匀的半球状,即无论哪个方向的环境光打进来, fr" role="presentation">fr 函数的输出几乎没多少变化(甚至是个常数)。于是基于以上理论,Split Sum 方法对渲染方程改造成这样的近似公式来获得渲染的加速:
Lr(x,ωr)≈∫ΩfrLi(x,ωi)dωi∫Ωfr dωi⋅∫Ω+fr(x,ωi,ωr)cos⁡θi dωi" role="presentation">Lr(x,ωr)≈∫ΩfrLi(x,ωi)dωi∫Ωfr dωi⋅∫Ω+fr(x,ωi,ωr)cosθi dωi
对渲染方程拆分成两个部分(环境光积分、BRDF积分)后就可以通过预计算的方式(后面两节会介绍如何预计算)分别减少这些积分的运行时开销,总结这种方法的好处是:
Split Sum 方法和原始蒙特卡洛方法的图像效果几乎一模一样 由于不用对环境贴图进行多重采样,性能开销大大减低了环境光积分:
∫ΩfrLi(x,ωi)dωi∫Ωfr dωi" role="presentation">∫ΩfrLi(x,ωi)dωi∫Ωfr dωi
因为我们已经拥有一张环境贴图(无论是实时的还是预渲染的)来存储环境光信息了,为了计算环境光部分的积分,需要在 Ωfr" role="presentation">Ωfr 范围内做多次光线采样。但是,可以有一个几乎等价但避免运行时多次采样的方式:
预先对纹理进行滤波操作(模糊),我们只需要对滤波后的环境贴图采样一次光线方向就能得到积分。
滤波后的环境贴图实际上称为 Irradiance Environment Map :原来环境贴图中的单位是 radiance L(x,ω)" role="presentation">L(x,ω),贴图积分后的单位则变成了 Irradiance E(x)" role="presentation">E(x)。
当然,BRDF Lobe 的形状越尖锐,即环境光积分范围越小,就需要使用模糊程度更低的环境贴图;反之 BRDF Lobe 的形状越粗壮,即环境光积分范围越大,就需要使用模糊程度更高的环境贴图。
我们可以使用 Mipmap技术 来生成不同Level的环境贴图,通过三线性插值(Trilinear Interportion)的方式来得出任何模糊程度且任何2D位置的环境光滤波结果。
IBL 往往用于支持静态环境光照,这样就可以进行预过滤环境贴图;当然也可以支持动态环境光照,只是需要实时过滤出动态的环境贴图,有一定的开销。
回顾 Microfacet BRDF 的组成:
Fresnel项(举例 Schlick's approximation)F=F0+(1−F0)(1−(h⋅ωr))5=F0+(1−F0)(1−(cos⁡θvh))5" role="presentation">F=F0+(1−F0)(1−(h⋅ωr))5=F0+(1−F0)(1−(cosθvh))5
NDF项(举例 Beckmann NDF)D(h)=1πα2(n⋅h)4⋅exp⁡((n⋅h)2−1a2(n⋅h)2)=1πα2cos4⁡θh⋅exp⁡(cos2⁡θh−1a2cos2⁡θh)" role="presentation">D(h)=1πα2(n⋅h)4⋅exp((n⋅h)2−1a2(n⋅h)2)=1πα2cos4θh⋅exp(cos2θh−1a2cos2θh)
可以预想到,该BRDF 的积分结果依赖于三个参数:
F0" role="presentation">F0(Fresnel项系数) α" role="presentation">α (粗糙度 roughness) θv" role="presentation">θv (反射方向与法线的夹角,实际上决定了 θvh" role="presentation">θvh、θh" role="presentation">θh)我们可以把 F0" role="presentation">F0 项拆出来,让BRDF 积分拆成两个积分,但是这些积分都减少了一个依赖的参数:
(1)∫Ω+fr(p,ωi,ωo)cos⁡θi dωi≈F0∫Ω+frF(1−(1−cos⁡θvh)5)cos⁡θi dωi(2)+∫Ω+frF(1−cos⁡θvh)5cos⁡θi dωi(3)=F0A+B" role="presentation">(1)∫Ω+fr(p,ωi,ωo)cosθi dωi≈F0∫Ω+frF(1−(1−cosθvh)5)cosθi dωi(2)+∫Ω+frF(1−cosθvh)5cosθi dωi(3)=F0A+B
这样,对于相同的材质(相同的BRDF),我们就可以针对剩余两个依赖的参数 α" role="presentation">α 、θv" role="presentation">θv 建立一张 二维查询表,其两个通道分别代表 A 与 B 的值:
SH 是 PRT 的前置知识,单独提出来讲一下。
对于任意函数 f(x)" role="presentation">f(x) ,不管是连续的还是不连续的,我们可以展开成一系列基函数(每项都带某个系数)的线性组合:
f(x)=∑ici⋅Bi(x)" role="presentation">f(x)=∑ici⋅Bi(x)
例如,多项式展开可以看成是一系列基函数(多项式)的线性组合:f(x)=c0+c1⋅x1+c2⋅x2+..." role="presentation">f(x)=c0+c1⋅x1+c2⋅x2+...
再例如,傅立叶变换也可以将 f(x)" role="presentation">f(x) 表示成另一系列基函数(各种频率的正弦谐波)的线性组合:
而 球谐(Spherical Harmonics,SH) 便是定义在球面上的一系列2D基函数,它与2D傅里叶序列有点相似,但非常适合球面函数 f(ω)" role="presentation">f(ω)(即参数为单位球面向量)。
SH 是分阶数的:
在第0阶(l=0)有1个基函数(m=0) 在第1阶(l=1)有3个基函数(m=-1,0,1) ... 在第n阶(l=n)有2n+1个基函数(m=-n,...,0,...,n)。图中蓝色意味着正值,黄色意味着负值。
实际上,阶数越低的基函数所代表的信息就越是低频。如果要完全复原一个任意函数,我们需要无穷阶的 SH,而如果我们只要复原一个任意函数的近似(换句话说只重建出该函数的低频信息),那么我们完全可以只需要前几阶的 SH(包含 l=0,l=1,...,l=n 每一阶的所有基函数)去拟合。
实时渲染中常用二阶SH(其包含系数总共为1+3=4个)和三阶SH(其包含系数总共为1+3+5=9个)来表示。
SH 基函数的形式比较复杂,只建议看一眼就跳过。
SH 的实数形式的基函数如下:
ylm(θ,φ)={2Re⁡(Ylm)m>02Im⁡(Ylm)m<0Yl0m=0={2Klmcos⁡mφPlm(cos⁡θ)m>02Klmsin⁡|m|φPl|m|(cos⁡θ)m<0Kl0Pl0(cos⁡θ)m=0" role="presentation">ylm(θ,φ)={2Re(Ylm)m>02Im(Ylm)m<0Yl0m=0={2KlmcosmφPlm(cosθ)m>02Klmsin|m|φPl|m|(cosθ)m<0Kl0Pl0(cosθ)m=0
其中,
Klm=(2l+1)(l−|m|)!4π(l+|m|)!" role="presentation">Klm=(2l+1)(l−|m|)!4π(l+|m|)!
plm" role="presentation">plm 为勒让德多项式(Legendre polynomials):
P00=1" role="presentation">P00=1 Pmm=(1−2m)Pm−1m−1" role="presentation">Pmm=(1−2m)Pm−1m−1 Pm+1m=(2m+1)zPmm" role="presentation">Pm+1m=(2m+1)zPmm Plm=(2l−1)zPl−1m−(l+m−1)Pl−2ml−m" role="presentation">Plm=(2l−1)zPl−1m−(l+m−1)Pl−2ml−mz为球面向量对应球面坐标的z值,该多项式不依赖x值y值。
以下是前3阶的SH实数形式的基函数:
l=0" role="presentation">l=0:
Y00=s=Y00=121π" role="presentation">Y00=s=Y00=121π
l=1" role="presentation">l=1:
Y1,−1=py=i12(Y1−1+Y11)=34π⋅yrY1,0=pz=Y10=34π⋅zrY1,1=px=12(Y1−1−Y11)=34π⋅xr" role="presentation">Y1,−1=py=i12(Y1−1+Y11)=34π⋅yrY1,0=pz=Y10=34π⋅zrY1,1=px=12(Y1−1−Y11)=34π⋅xr
l=2" role="presentation">l=2:
Y2,−2=dxy=i12(Y2−2−Y22)=1215π⋅xyr2Y2,−1=dyz=i12(Y2−1+Y21)=1215π⋅yzr2Y2,0=dz2=Y20=145π⋅−x2−y2+2z2r2Y2,1=dsz=12(Y2−1−Y21)=1215π⋅zxr2Y2,2=dx2−y2=12(Y2−2+Y22)=1415π⋅x2−y2r2" role="presentation">Y2,−2=dxy=i12(Y2−2−Y22)=1215π⋅xyr2Y2,−1=dyz=i12(Y2−1+Y21)=1215π⋅yzr2Y2,0=dz2=Y20=145π⋅−x2−y2+2z2r2Y2,1=dsz=12(Y2−1−Y21)=1215π⋅zxr2Y2,2=dx2−y2=12(Y2−2+Y22)=1415π⋅x2−y2r2
SH 与一般的基函数相比,具有以下性质:
基函数之间具有正交性(orthonormal)∫ΩBi(ω)⋅Bj(ω)dω=1(i=j)∫ΩBi(ω)⋅Bj(ω)dω=0(i≠j)" role="presentation">∫ΩBi(ω)⋅Bj(ω)dω=1(i=j)∫ΩBi(ω)⋅Bj(ω)dω=0(i≠j)
通过 投影(Projection) 可以很方便得到 SH 系数(SH coefficients)fi=∫Ωf(ω)⋅Bi(ω)dω" role="presentation">fi=∫Ωf(ω)⋅Bi(ω)dω
通过系数向量组与基函数组的点积(前n阶的基函数/系数共有 n*n 个)可以很方便重建球面函数f(ω)≈∑i=1n2fiBi(ω)" role="presentation">f(ω)≈∑i=1n2fiBi(ω)
product projection:c(ω)=a(ω)b(ω)" role="presentation">c(ω)=a(ω)b(ω) ,已知 a" role="presentation">a 的形式而 b" role="presentation">b 未知,有c→=M⋅b→" role="presentation">c→=M⋅b→
其中,c→" role="presentation">c→ 为 c(ω)" role="presentation">c(ω) 的SH系数向量组,b→" role="presentation">b→ 为 b(ω)" role="presentation">b(ω) 的SH系数向量组,M" role="presentation">M 为 n2" role="presentation">n2 X n2" role="presentation">n2 的矩阵,其元素为
Mij=∫Ωa(ω)Bi(ω)Bj(ω)dω" role="presentation">Mij=∫Ωa(ω)Bi(ω)Bj(ω)dω
这样我们就可以预计算好 M" role="presentation">M ,等到知道 b→" role="presentation">b→ 后,乘起来就能得到 c→" role="presentation">c→
推导:令Mi(ω)=a(ω)Bi(ω)" role="presentation">令Mi(ω)=a(ω)Bi(ω)
(4)ci=∫Ωa(ω)b(ω)Bi(ω)dω(5)=∫ΩMi(ω)b(ω)dω(6)=∫Ω∑j=1n2(MijBj(ω))b(ω)dω(7)=∑i=1n2Mij∫ΩBj(ω)b(ω)dω(8)=∑j=1n2Mij⋅bj" role="presentation">(4)ci=∫Ωa(ω)b(ω)Bi(ω)dω(5)=∫ΩMi(ω)b(ω)dω(6)=∫Ω∑j=1n2(MijBj(ω))b(ω)dω(7)=∑i=1n2Mij∫ΩBj(ω)b(ω)dω(8)=∑j=1n2Mij⋅bj
支持插值,对 SH 系数的插值相当于对重建的函数值的插值。
旋转不变性(rotational invariance),对函数 f" role="presentation">f 的旋转 RSH" role="presentation">RSH 等价于对 f(ω)" role="presentation">f(ω) 的自变量的旋转 R3D" role="presentation">R3D:
RSH(f(ω))=f(R3D(ω))" role="presentation">RSH(f(ω))=f(R3D(ω))
预计算辐射度传输(Precomputed Radiance Transfer,PRT) 是一类基于预计算 radiance transfer(辐射度传输)的渲染方法。
注:PRT 核心并不是预计算环境光照,而是预计算辐射度传输。
所谓 radiance transfer,可以理解成是一种关于光路传输分布的几何信息,例如 shadow,ambient occlusion,物体表面互反射等效果都属于此。也就是说 PRT 要求几何信息是静态的,如果应用到物体上,就需要物体本身不会发生形变;如果应用到几何场景上,就需要场景的几何不会发生变化。
PRT 方法在 2002年 SIGGRAPH 会议被 Peter-Pike Sloan 首先提出,并从此掀起了一波 PRT 方法的研究浪潮,当时论文所实现的 PRT 算法便是大量运用了球谐(SH)方式。实际上 PRT 不仅可以通过 SH 去表示光路传输分布,还可以有其它方式去表示(如 Spherical Gaussian、Wavelet)。
PRT(SH lighting & SH transfer)效果(分别为 diffuse 物体的无阴影情况、diffuse 物体的阴影&互反射情况、glossy 物体的无阴影情况、glossy 物体的阴影&互反射情况。):
可以看到 PRT 其实不像 IBL 直接应用环境光照,而是额外考虑了物体几何对环境光照的遮蔽效果。
PRT 首先将渲染方程分为两个球面函数,即环境光部分(lighting function)和传输部分(transfer function):
L(r)=∫Ω+Le(ω)⋅ρ(r,w)⋅V(ω)max(0,n⋅ω)dω=∫Ω+Le(ω)⋅T(r,w)dω" role="presentation">L(r)=∫Ω+Le(ω)⋅ρ(r,w)⋅V(ω)max(0,n⋅ω)dω=∫Ω+Le(ω)⋅T(r,w)dω
由于环境光是场景所有物体共享的,所以往往只需要存储一份全局信息,我们甚至可以使用 IBL 方法,但是在 PRT 采用了一种更加压缩的环境光照表达方式,即 SH Lighting。
球谐光照(SH Lighting):使用低阶的 SH 来表示环境光照,这样只需要预先存储少数个 SH 系数(diffuse 情况下 2 阶SH Lighting 只需 RGB 每个通道各 4 个系数,即共 12 个系数)就可以重建出环境光照。
传统的环境光照信息往往使用环境贴图(Environment Map)表示,这往往需要一个庞大的二维数组存储各个 texel 的值,十分耗费空间,而且采样和纹理I/O也有一定开销。而纹理本质上也是一个函数(信号),输入二维坐标,输出对应纹素的 RGBA 值,因此可以用 SH 表示。
SH Lighting 只预先存储若干个 SH 系数,而不必存储一整张环境贴图,大大减少存储空间和 I/O 开销。 SH 所能表示的 lighting 信息是低频的,因此 SH lighting 只适用于 diffuse 和 glossy 的物体而不适用于 specular 的物体。 环境光照动态性会有所限制(最多只能支持环境光照的 rotation)。传输部分则因为是几何相关的,即每个位置的 transfer 信息往往是不一样的,为此可以每个顶点存储一份 transfer 信息,PRT 同样使用了 SH 来表达 transfer 信息。
PRT 中每个顶点预先存储若干个 SH 系数(diffuse情况下)或者矩阵(glossy情况下),而顶点之间的位置则通过插值的方式来得到对应的 transfer 信息。
预计算 radiance transfer 可以在实时渲染中以极低代价实现单个物体的自阴影、互反射的局部几何传输效果。 物体自身不能有局部的几何变化(但可以整体 scale, rotation 和 transition)。如果要将 PRT 应用到场景中,则就要求场景所有物体不能发生任何几何变化(包括 scale,rotation,transition)。由于场景的几何复杂度往往比单个物体的几何复杂度要高得多,并且渲染效果也明显的多,因此要求的 SH 阶数应当更高(精度更高),导致更沉重的空间需求,所以下文以单个物体的 PRT 方法为例讲解。
Diffuse 物体的渲染方程:
L(r)=∫Ω+Le(ω)⋅ρ⋅V(ω)max(0,n⋅ω)dω" role="presentation">L(r)=∫Ω+Le(ω)⋅ρ⋅V(ω)max(0,n⋅ω)dω
Le" role="presentation">Le 代表环境光; ρ" role="presentation">ρ 为 BRDF 项;V" role="presentation">V 为 Visibility,表示不被遮挡的程度,往往表现为自遮挡产生的阴影现象。
由于物体是 diffuse 的,因此它的 BRDF 将是一个常数 ρ" role="presentation">ρ (无论从哪个方向观察都得到相同的BRDF值),也因此 diffuse 物体的 L" role="presentation">L 将是一个常量值,而不受参数 r" role="presentation">r 影响 对于 lighting 部分 L(ω)" role="presentation">L(ω),原本需要对 Environment Map 进行查询,而现在可以换成使用 SH 函数去表示:Le(ω)≈∑liBi(ω)" role="presentation">Le(ω)≈∑liBi(ω) 对于 transfer function 部分 T(ω)=V(ω)max(0,n⋅ω)" role="presentation">T(ω)=V(ω)max(0,n⋅ω) , 也可以换成使用 SH 函数去表示:T(ω)≈∑TjBj(ω)" role="presentation">T(ω)≈∑TjBj(ω)代入渲染方程整理后得:
(9)L(r)≈ρ∫Ω+∑i=1n2liBi(ω)∑j=1n2TjBj(ω)dω(10)=ρ∑i=1n2∑j=1n2liTj∫Ω+Bi(ω)Bj(ω)dω(11)=ρ∑i=1n2liTi∫Ω+Bi(ω)Bi(ω)dω(12)=ρ∑i=1n2liTi" role="presentation">(9)L(r)≈ρ∫Ω+∑i=1n2liBi(ω)∑j=1n2TjBj(ω)dω(10)=ρ∑i=1n2∑j=1n2liTj∫Ω+Bi(ω)Bj(ω)dω(11)=ρ∑i=1n2liTi∫Ω+Bi(ω)Bi(ω)dω(12)=ρ∑i=1n2liTi
这样,我们只需要预计算出 li" role="presentation">li 和 Ti" role="presentation">Ti,就能让运行时的diffuse物体渲染速度大大提升。
预计算过程:
对整个环境光信息(环境贴图),预计算 lighting 的 SH 系数向量组 l→" role="presentation">l→,其中一个元素为:li=∫Ω+Le(ω)Bi(ω)dω" role="presentation">li=∫Ω+Le(ω)Bi(ω)dω
对每个顶点,预计算 light transfer 的 SH 系数向量组 T→" role="presentation">T→,其中一个元素为:Ti=∫Ω+V(ω)max(0,n⋅ω)Bi(ω)dω" role="presentation">Ti=∫Ω+V(ω)max(0,n⋅ω)Bi(ω)dω
此时,我们可以理解成系数 li" role="presentation">li 代表了环境光照(lighting)的信息,而系数 Ti" role="presentation">Ti 代表了在某个顶点上光路传输(light transfer)的信息;整个模型预计算完成之后将会得到一个 lighting 系数向量和模型顶点数量对应的 transfer 系数向量。
运行时渲染过程:
在 vertex shading 阶段计算顶点的 SH 颜色 L(r)≈ρ∑i=1n2liTi=ρ(l→⋅T→)" role="presentation">L(r)≈ρ∑i=1n2liTi=ρ(l→⋅T→)
在 pixel/fragment shading 阶段得到插值后的 SH 颜色即为该像素的颜色
Glossy 物体的渲染方程:
L(r)=∫Ω+Le(ω)⋅ρ(r,ω)⋅V(ω)⋅max(0,n⋅ω)dω" role="presentation">L(r)=∫Ω+Le(ω)⋅ρ(r,ω)⋅V(ω)⋅max(0,n⋅ω)dω
对于 glossy 物体,不同视角观察物体表面同一点会有不同的光照。因此 glossy 的 BRDF 将是四维的函数(参数不仅包含ω" role="presentation">ω,还包含r" role="presentation">r),这次将 BRDF 算入 transfer,则 transfer function 将额外增加一个二维参数 r" role="presentation">r
对于环境光(lighting)即 L(ω)" role="presentation">L(ω),原本需要对 Environment Map 进行查询,而现在可以换成使用 SH 函数去表示:
Le(ω)≈∑liBi(ω)" role="presentation">Le(ω)≈∑liBi(ω)
对于传输函数(transfer function) T(r,ω)=ρ(r,ω)V(ω)max(0,n⋅ω)" role="presentation">T(r,ω)=ρ(r,ω)V(ω)max(0,n⋅ω) ,换成使用 SH 函数去表示:
T(r,ω)≈∑i=1n2Ti(r)Bi(ω)≈∑i=1n2∑j=1n2TijBj(r)Bi(ω)" role="presentation">T(r,ω)≈∑i=1n2Ti(r)Bi(ω)≈∑i=1n2∑j=1n2TijBj(r)Bi(ω)
由于额外多了一个二维的参数,transfer 系数将是一个矩阵而非之前 diffuse 情况下的系数向量
代入渲染方程整理后得:
(13)L(r)≈∫Ω+∑i=1n2liBi∑j=1n2∑k=1n2TjkBj(ω)Bk(r)dω(14)=∑i=1n2∑j=1n2∑k=1n2liTjkBk(r)∫Ω+Bi(ω)Bj(ω)dω(15)=∑i=1n2∑k=1n2liTikBk(r)" role="presentation">(13)L(r)≈∫Ω+∑i=1n2liBi∑j=1n2∑k=1n2TjkBj(ω)Bk(r)dω(14)=∑i=1n2∑j=1n2∑k=1n2liTjkBk(r)∫Ω+Bi(ω)Bj(ω)dω(15)=∑i=1n2∑k=1n2liTikBk(r)
预计算过程:
对整个环境光信息,预计算 lighting 的 SH 系数向量组 l→" role="presentation">l→,其中一个元素为:li=∫Ω+L(ω)Bi(ω)dω" role="presentation">li=∫Ω+L(ω)Bi(ω)dω 对每个顶点,预计算 light transfer 矩阵 T" role="presentation">T,其中一个元素为:Tij=∫Ω+Bj(ω)Ti(ω)dω" role="presentation">Tij=∫Ω+Bj(ω)Ti(ω)dω 、Ti(r)=∫Ω+T(r,ω)Bi(ω)dω" role="presentation">Ti(r)=∫Ω+T(r,ω)Bi(ω)dω由于 Glossy 物体需要每个顶点存储三个矩阵 T(RGB颜色每个通道对应1个矩阵),这使得空间开销略大,因此不太实用。
运行时渲染过程:
在 vertex shding 阶段计算顶点的 SH 颜色 L(r)≈∑i=1n2∑j=1n2liTijBj(r)=Tl→⋅B(r)→" role="presentation">L(r)≈∑i=1n2∑j=1n2liTijBj(r)=Tl→⋅B(r)→ 在 pixel/fragment shading 阶段得到插值后的 SH 颜色即为该像素的颜色原论文实际上还有个 BRDF 卷积步骤。但是本身用SH表示环境光信息已经足够低频(模糊)了,对这个卷积步骤,个人先存个疑。
而前面的渲染方程中:
Lx(r)=Lx0(r)=∫Ω+Le(ω)⋅ρ(r,ω)⋅V(ω)⋅max(0,n⋅ω)dω" role="presentation">Lx(r)=Lx0(r)=∫Ω+Le(ω)⋅ρ(r,ω)⋅V(ω)⋅max(0,n⋅ω)dω
只考虑几何传播关系的因素只有 V(ω)⋅max(0,n⋅ω)" role="presentation">V(ω)⋅max(0,n⋅ω),即只可实现阴影或者无阴影(令 V=1" role="presentation">V=1)的效果,并不能实现互反射(interreflection)的效果。
所谓 interreflection,实际上就是某个点朝外围环境看时,V=0" role="presentation">V=0 意味着有几何部分遮挡住了直接的环境光,但实际上遮挡还意味着这些几何部分上的出射光会反射给该点,这就形成了该点间接的环境光照。
(迭代)增加 1 次互反射的渲染方程:
Lxn(r):=Lx0(r)+∫Ω+Lx′n−1(ω)⋅ρ(r,ω)⋅(1−V(ω))⋅max(0,n⋅ω)dω" role="presentation">Lxn(r):=Lx0(r)+∫Ω+Lx′n−1(ω)⋅ρ(r,ω)⋅(1−V(ω))⋅max(0,n⋅ω)dω
Lxn" role="presentation">Lxn 意味着在某个点 x" role="presentation">x 的 radiance,而 x′" role="presentation">x′ 则是从 x" role="presentation">x 朝 ω" role="presentation">ω 方向进行 raycast 而得到相交的其它点(x′" role="presentation">x′ 依赖于 ω" role="presentation">ω),而 n" role="presentation">n 则代表了迭代了多少次。
由于 diffuse 与出射方向无关,因此去掉 L" role="presentation">L 的参数 ω" role="presentation">ω。
Lxn:=Lx0+ρ∫Ω+Lx′n−1⋅(1−V(ω))⋅max(0,n⋅ω)dω" role="presentation">Lxn:=Lx0+ρ∫Ω+Lx′n−1⋅(1−V(ω))⋅max(0,n⋅ω)dω
将上述式子左右两边同时分离出 l→" role="presentation">l→ (因为无论点在 x" role="presentation">x 还是别的 x′" role="presentation">x′,其所在的环境光照 lighting SH 系数都是一样的 ),则意味着我们实现 Interreflection 的效果只需要在所有顶点 x" role="presentation">x 的 transfer 向量做以下迭代:
Txn:=Tx0+ρ∫Ω+Tx′n−1⋅(1−V(ω))⋅max(0,n⋅ω)dω" role="presentation">Txn:=Tx0+ρ∫Ω+Tx′n−1⋅(1−V(ω))⋅max(0,n⋅ω)dω
PRT 的一个问题是如果 lighting 部分是预计算的,那就只适用于静态环境光下的静态物体渲染;环境光或者物体只要有变化,PRT 就不得不进行重新预计算;但得益于 SH 的旋转不变性,我们至少可以让 SH Lighting 适用于动态旋转的情形而不必重新预计算。
注:环境光旋转和物体旋转在 PRT 渲染中是等价的(只是说看相对于哪个东西来看待旋转而已)。当然实践中,我们也可以对输入向量来进行反向旋转,来获得等同于环境光正向旋转的效果;只是如果考虑每次对输入向量进行旋转操作比较费的话,可以考虑一次性旋转完 SH lighting,就避免了输入向量旋转的操作了。
现在给定旋转 R" role="presentation">R ,然后有SH投影函数 P" role="presentation">P(输入一个球面向量,输出第 l" role="presentation">l 层band的 SH 系数向量组)。假设我们有 2l+1" role="presentation">2l+1 个任意的球面向量,我们需要想办法求出能等价于旋转 ni" role="presentation">ni 的矩阵 M" role="presentation">M 来旋转 SH 投影函数 P" role="presentation">P :
MP(ni)=P(R(ni)),i∈[−l,l]" role="presentation">MP(ni)=P(R(ni)),i∈[−l,l]
整理得:
M[P(n−l),…,P(nl)]=[P(R(n−l)),…,P(R(nl))]" role="presentation">M[P(n−l),…,P(nl)]=[P(R(n−l)),…,P(R(nl))]
记 A=[P(n−l),...,P(nl)]" role="presentation">A=[P(n−l),...,P(nl)] ,如果矩阵 A" role="presentation">A 是可逆的,则:
M=[P(R(n−l)),…,P(R(nl))]A−1" role="presentation">M=[P(R(n−l)),…,P(R(nl))]A−1
这样,无论 R" role="presentation">R 怎么变化,我们都可以即时根据 R" role="presentation">R 求出 M" role="presentation">M ,然后把最新的 M" role="presentation">M 应用到所有预计算好的 P(ω)" role="presentation">P(ω) 上(即预计算好的SH系数上),而不必重新预计算 P(R(ω))" role="presentation">P(R(ω)) 。
也就是要做 SH 的快速旋转实现,就需要对所有层进行处理,对于第 l" role="presentation">l 层 band 处理如下:
选取共 2l+1" role="presentation">2l+1 个 单位向量 ni" role="presentation">ni,这些向量 ni" role="presentation">ni 投影在该层 band 上的基函数就能得到系数 P(ni)" role="presentation">P(ni)。共 2l+1" role="presentation">2l+1 个 2l+1" role="presentation">2l+1 维向量 P(ni)" role="presentation">P(ni) 构成了矩阵 A" role="presentation">A,并求出 A−1" role="presentation">A−1如何选取单位向量:要保证投影后构成的 A" role="presentation">A 矩阵可逆。
给定旋转 R" role="presentation">R ,对所有 ni" role="presentation">ni 依次做旋转 R(ni)" role="presentation">R(ni) ,这些旋转后的单位向量同样投影在该层 band 上的基函数就能得到系数 P(R(ni))" role="presentation">P(R(ni))。 共 2l+1" role="presentation">2l+1 个 2l+1" role="presentation">2l+1 维向量 P(R(ni))" role="presentation">P(R(ni)) 构成了矩阵 S" role="presentation">S
求出该层 band 上球谐函数的旋转矩阵 M=SA−1" role="presentation">M=SA−1
用 M" role="presentation">M 乘以该层 band 上的 SH 系数向量组就可以得到旋转后的 SH 系数向量组。
l=0" role="presentation">l=0 时只有1个系数,故不需要处理;l=1" role="presentation">l=1 时需要处理3个系数,其中的 A、S、M" role="presentation">A、S、M 将会是 3X3 矩阵;l=2" role="presentation">l=2 时需要处理5个系数,其中的 A、S、M" role="presentation">A、S、M 将会是 5X5矩阵......
最后,将每一层 band 的结果重新拼接起来即可得到完整的旋转后的 SH 系数结果。
// 伪代码:三阶SH旋转 SHCoeffientsAfterRotation(Array coeffients, Rotation rotation){ Array res; // 处理 l = 0 res[0] = coffients[0]; // 处理 l = 1 // step 1 n1 = [1,0,0],n2 = [0,1,0],n3 = [0,0,1]; A = [ // SHProject(n,l,m):向量n投影到SH中l层第m个基函数,输出对应系数 [SHProject(n1,1,-1),SHProject(n2,1,-1),SHProject(n3,1,-1)], [SHProject(n1,1,0),SHProject(n2,1,0),SHProject(n3,1,0)], [SHProject(n1,1,1),SHProject(n2,1,1),SHProject(n3,1,1)] ]; A_inv = InvMatrix(A); // step 2 rn1 = rotation*n1,rn2 = rotation*n2,rn3 = rotation*n3; S = [ [SHProject(rn1,1,-1),SHProject(rn2,1,-1),SHProject(rn3,1,-1)], [SHProject(rn1,1,0),SHProject(rn2,1,0),SHProject(rn3,1,0)], [SHProject(rn1,1,1),SHProject(rn2,1,1),SHProject(rn3,1,1)] ]; // step 3 M = S*A_inv; // step 4 coeff_l1 = M*[coeffients[1],coeffients[2],coeffients[3]]; res[1] = coeff_l1[0]; res[2] = coeff_l1[1]; res[3] = coeff_l1[2]; // 处理 l = 2 // step 1k = 1/sqrt(2); n1 = [1,0,0],n2 = [0,0,1],n3 = [k,k,0],n4 = [k,0,k],n5 = [0,k,k]; A = [[SHProject(n1,2,-2),SHProject(n2,2,-2),SHProject(n3,2,-2),SHProject(n4,2,-2),SHProject(n5,2,-2)],[SHProject(n1,2,-1),SHProject(n2,2,-1),SHProject(n3,2,-1),SHProject(n4,2,-1),SHProject(n5,2,-1)],[SHProject(n1,2,0),SHProject(n2,2,0),SHProject(n3,2,0),SHProject(n4,2,0),SHProject(n5,2,0)],[SHProject(n1,2,1),SHProject(n2,2,1),SHProject(n3,2,1),SHProject(n4,2,1),SHProject(n5,2,1)],[SHProject(n1,2,2),SHProject(n2,2,2),SHProject(n3,2,2),SHProject(n4,2,2),SHProject(n5,2,2)] ]; A_inv = InvMatrix(A); // step 2 rn1 = rotation*n1,rn2 = rotation*n2,rn3 = rotation*n3,rn4 = rotation*n4,rn5 = rotation*n5; S = [[SHProject(rn1,2,-2),SHProject(rn2,2,-2),SHProject(rn3,2,-2),SHProject(rn4,2,-2),SHProject(rn5,2,-2)],[SHProject(rn1,2,-1),SHProject(rn2,2,-1),SHProject(rn3,2,-1),SHProject(rn4,2,-1),SHProject(rn5,2,-1)],[SHProject(rn1,2,0),SHProject(rn2,2,0),SHProject(rn3,2,0),SHProject(rn4,2,0),SHProject(rn5,2,0)],[SHProject(rn1,2,1),SHProject(rn2,2,1),SHProject(rn3,2,1),SHProject(rn4,2,1),SHProject(rn5,2,1)],[SHProject(rn1,2,2),SHProject(rn2,2,2),SHProject(rn3,2,2),SHProject(rn4,2,2),SHProject(rn5,2,2)] ]; // step 3 M = S*A_inv; // step 4 coeff_l2 = M*[coeffients[4],coeffients[5],coeffients[6],coeffients[7],coeffients[8]]; res[4] = coeff_l2[0]; res[5] = coeff_l2[1]; res[6] = coeff_l2[2]; res[7] = coeff_l2[3]; res[8] = coeff_l2[4]; return res; }
这样我们对 lighting 的 SH 系数做 rotation,就可以实现支持环境光旋转或物体旋转的实时 SH Lighting。
但是注意,不应当对 transfer 的 SH 系数做 rotation,因为一个物体的 transfer function 只考虑自身的几何遮蔽(自身与自身的遮蔽关系),而非像物体与环境光那样存在相对关系。
SH rotation 效果图:
前面的 PRT 章节中,我们知道 diffuse 物体的 transfer 迭代公式如下:
Txn:=Tx0+ρ∫Ω+Tx′n−1⋅(1−V(ω))⋅max(0,n⋅ω)dω" role="presentation">Txn:=Tx0+ρ∫Ω+Tx′n−1⋅(1−V(ω))⋅max(0,n⋅ω)dω
要计算出 diffuse 物体的 infinite-bounce 的 transfer,理论上需要进行无限次迭代。然而我们可以利用一些假设来化简掉无限迭代,从而实现一次迭代算出 infinite-bounce 的 transfer 向量。
一个观察是点 x" role="presentation">x 与其发出射线能相交到的点 x′" role="presentation">x′ 的距离往往是比较近的(万一距离远了 x′" role="presentation">x′ 在积分计算中的贡献反正也比较少),因此它们的几何环境也是比较接近的。因此我们假设 x" role="presentation">x 和 x′" role="presentation">x′ 的 transfer 基本一致,即 x=x′" role="presentation">x=x′。
因此原式可以移除无关参数,进一步化简为:
Tn:=T0+Tn−1⋅ρ∫Ω+(1−V(ω))⋅max(0,n⋅ω)dω" role="presentation">Tn:=T0+Tn−1⋅ρ∫Ω+(1−V(ω))⋅max(0,n⋅ω)dω
令 A=ρ⋅∫Ω+(1−V(ω))⋅max(0,n⋅ω)dω" role="presentation">A=ρ⋅∫Ω+(1−V(ω))⋅max(0,n⋅ω)dω
由于
diffuse BRDF为 ρ=albedoπ" role="presentation">ρ=albedoπ,而一般 albedo 都会小于 1 ∫Ω+(1−V(ω))⋅max(0,n⋅ω)<=π" role="presentation">∫Ω+(1−V(ω))⋅max(0,n⋅ω)<=π因此可得 A<1" role="presentation">A<1
原递推式可得:
T1=T0+T0⋅A" role="presentation">T1=T0+T0⋅A T2=T0+T1⋅A" role="presentation">T2=T0+T1⋅A ... Tn=T0+Tn−1⋅A" role="presentation">Tn=T0+Tn−1⋅A递归代入后得:
Tn=T0⋅(1+A+A2+...+An)" role="presentation">Tn=T0⋅(1+A+A2+...+An)
等比数列求和得:
Tn=T0⋅1−An1−A" role="presentation">Tn=T0⋅1−An1−A
当 n 趋近于无穷时,又 A<1" role="presentation">A<1,则最终收敛得到了结论
limn−>∞Tn=limn−>∞(T0⋅1−An1−A)=T0⋅11−A" role="presentation">limn−>∞Tn=limn−>∞(T0⋅1−An1−A)=T0⋅11−A
无论 IBL 还是 PRT 都属于实现环境光照的方案,它们的区别在于:
IBL 是一种从预计算环境光照出发的环境关照渲染方案: 采用环境贴图:存储占用空间较大,同时也占采样 I/O。 能保留高频信息,常用于 diffuse/glossy/specular 物体的渲染。 PRT 是一种从预计算 transfer function 出发的环境光照渲染方案: 采用 SH lighting:存储开销和重建环境光照的开销极低。 只能保留低频信息,常用于 diffuse/glossy 物体的渲染。 效果上,相比 IBL 额外考虑了物体几何对环境光照的遮蔽效果。但只考虑了物体局部 transfer 效果,没有考虑完整场景 transfer 效果。不过在其它 PRT 方案中有支持完整场景的 transfer 效果。 物体不可局部形变,材质不可动态:若发生变化,那么其 transfer 就需要更新。全局光照(Global Illumination,GI),严格意义上指完整的物理的间接光照效果,实际上只要有体现间接光效果(即使是部分效果或者 tricky 的效果)的光照都可称为 GI,因此它总体上分为两派:
有一派是通过一套统一标准的渲染流程把任何物体的 GI 计算出来,往往计算量极大但效果更加物理更加真实,典型的例子便是 path tracing 另一派是把 GI 的效果看成一个个部分来组成,这样我们可以选择其中一些 GI 效果组合使用,适应不同的性能(往往计算量相对低些,尤其是实时渲染)的同时也能带来能接受的 GI 效果(虽然往往不是严谨的物理正确),例如:在单个物体的角度上实现环境光照的局部自遮蔽/自反射效果(例如 PRT),在整个场景的角度来实现间接光照在全局场景 bounce 来 bounce 去的效果(例如 VXGI)。相关知识
使用预计算实时全局光照优化照明
(HDRP)全局光照技术初探(一)
环境光照
光照策略
基于Labview的光照环境控制器智能监控系统
一种基于物联网的农林病虫害监控系统的制作方法
基于改进YOLOv8算法的实时细粒度植物病害检测
基于VR的虚拟花卉真实感光照环境仿真交互应用研究
基于气象信息的膜下滴灌棉花实时灌溉预报
基于机器视觉的植物病虫害实时识别方法
网址: 基于预计算的实时环境光照(Real https://m.huajiangbk.com/newsview733386.html
上一篇: 幼儿园环境设计任务书.doc |
下一篇: 采光好,学生成绩可提高21%!学 |