Shading From Environment Lighting

Environment Lighting

环境光照是用一张图来记录的是在场景中向任何方向看去的光照情况,其认为光照来自无限远处,因此对场景中的每个点,环境光照是相同的
SphericalMapCubeMap

Image-Based Lighting(IBL)

工业界进行Environment Lighting的方法通常是Image-Based Lighting(IBL)

假设现在已经得到环境光照的映射结果图,计算Shading就是解Rendering Equation,因为Shading不考虑阴影,所以此处的Rendering Equation中不需要考虑Visibility项
IBLRednderingEquation
求解上述问题通常使用Monte Carlo积分。蒙特卡罗积分可参考GAMES101光线追踪部分

但是,蒙特卡罗积分需要大量的采样,通常来说,在Shader中需要采样的方法都因其较缓慢而不会优先采用,但是近年来随着Temporal方面的研究,基于采样的方法在RTR中的应用越来越多(当然,通常能避免采样肯定是最好的)

The Split Sum

那么如何避免采样呢?
首先,考察一下BRDF的情况
BRDFobservation
可见,Glossy BRDF的积分域比较小,而Diffuse BRDF的值在积分域上没有变化
由此可见,BRDF至少会满足Small Support或Smooth两个条件的其中一个,于是就可以用到上一节中提到的经典的积分近似方法

Ωf(x)g(x)dxΩGf(x)dx/ΩGdxΩg(x)dx\int_{\Omega}f(x)g(x)\mathrm{d}x\approx\int_{\Omega_{G}}f(x)\mathrm{d}x/\int_{\Omega_{G}}\mathrm{d}x\cdot\int_{\Omega}g(x)\mathrm{d}x

于是可以将渲染方程拆分如下:

Lo(p,ωo)ΩfrLi(p,ωi)dωi/ΩfrdωiΩ+fr(p,ωi,ωo)cosθidωiL_o(p,\omega_o)\approx\int_{\Omega_{f_r}}L_i(p,\omega_i)\mathrm{d}\omega_i\bigg/\int_{\Omega_{f_r}}\mathrm{d}\omega_i\cdot\int_{\Omega_{+}}f_r(p,\omega_i,\omega_o)\cos\theta_i\mathrm{d}\omega_i

注:先前在实时阴影的部分有如下的拆分形式

Lo(p,ωi)Ω+V(p,ωi)dωi/Ω+dωiΩ+Li(p,ωi)fr(p,ωi,ωo)cosθidωiL_o(p,\omega_i)\approx\int_{\Omega^+}V(p,\omega_i)\mathrm{d}\omega_i\bigg/\int_{\Omega^+}\mathrm{d}\omega_i\cdot\int_{\Omega^+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)\cos\theta_i\mathrm{d}\omega_i

在RTR中,拆分渲染方程的形式可以根据需要选用

现在回过头考察积分拆分后的前半边
SplitSumFirstHalf
实际上,该项就是对Environment Light贴图定义的环境光照情况进行Filtering,其Filter时的Kernel大小由BRDF的覆盖范围决定

EnvironmentLightingFiltering
由于已提前获得环境光照贴图,因此该项可以在渲染前就预计算完毕。类似MipMap的思想,此处也可以使用不同大小的Kernel生成多张图像,在渲染中不同的BRDF查询对应大小的Kernel生成的图,也可以进行类似三线性插值的操作

上述对环境光照贴图进行Filter的原理可以由下图表示
SplitSum1

通常来说,计算Shading Point的环境光照结果需要在该点BRDF进行采样,在其Lobe上采样多条路径对其结果加权平均。而上述的Filtering方法相当于对于Shading Point镜面反射方向,将其周围一小片区域的环境光照加权平均后写回该点,那么在渲染时只需要查询其镜面反射方向的环境光贴图而无需进行采样。如此便可获得较好的近似结果

上面所述是The Split Sum中的第一部分,仅避免拆分后的积分前半边积分的采样是不够的,还需要避免后半边的采样
SplitSumSecondHalf

首先复习一下GAMES101微表面模型中对BRDF的定义
MicrofacetBRDF
BRDF可以拆分成Fresnel项、自遮挡项和微表面法线分布项(NDF),其中Fresnel项和NDF最为关键
对于Fresnel项,GAMES101介绍了Schlick’s Approximation
Schlick
在Schlick’s Approximation中,菲涅尔项取决于基础反射率R0R_0和入射角θ\theta

而对于NDF,可以由Beckmann Distribution近似表示
Beckmann
其中,α\alpha可用于定义Roughness,θh\theta_h为入射角

使用上述方式拆分BRDF来进行预计算,就可以将原先五维的预计算降至三维,但是三维依然太高
使用下式改写后半边积分
SplitSum2
该式实际上是先除去Fresnel项再乘上Fresnel项,将乘上的Fresnel项用Schlick’s Approximation近似,然后将拆成两项的和,拆出基础反射率R0R_0
于是左边积分对于基础反射率R0R_0的依赖被消除,对该积分的预计算被降成二维(Roughness和入射角)
预计算的二维结果可以直接写成一张Texture
SplitSum2Tex

将积分写成求和形式如下,通常在RTR领域较少用积分表示,更多习惯于使用求和的写法,这也是为什么这个方法会叫做"The Split Sum",而不是"The Split Integeral"
The Split Sum Approximation
TheSplitSumApproximation

Shadow Form Environment Lighting

目前,没有很好的方法从环境光贴图中生成阴影。首先,如果将从环境光贴图生成阴影看成是一个多光源问题,则需要生成很多Shadow Map;如果将其看成是一个采样问题,则很难获得渲染方程中的Visibility项,并且其很可能不符合渲染方程近似方法的两个适用条件。

工业界目前的解决方案通常是选取场景中一个(或多个)最重要的光源来生成阴影。

Real-Time Environment Lighting(& Global Illumination)

背景知识

傅里叶级数

Fourier

滤波

图像及其频谱
ImageAndFrequency

低通滤波
LPF

卷积
Convolution

一个更通用的理解是:任何先相乘再积分的操作(Product Integral)都可以被认为是滤波

Basis Function

将一个函数描述成一系列函数的线性组合,如下:

f(x)=iciBi(x)f(x)=\sum_ic_i\cdot B_i(x)

其中,Bi(x)B_i(x)就被称为基函数
常见的基函数系列有:傅里叶级数、n次多项式系列(1,x,x2,x3...1,x,x^2,x^3...)等

Spherical Harmonics

球谐函数是一系列2D基函数,其中每个基函数都是定义在球上的(三维空间中的方向可以用θ,ϕ\theta,\phi描述,因此是二维)
球谐函数的理解可以类比1D傅里叶级数

球谐函数可视化如下:
SphericalHarmonicsVisualization
其中每一行基函数描述的频率相同,每一阶基函数的数量为2l+12l+1,前nn阶基函数总数为n2n^2

球谐函数的每个基函数都是利用Legendre Polynomial表示的

如果一个定义在球面上的函数ff用球谐函数表示,其每一项基函数的系数都可以用下式计算:

ci=Ωf(ω)Bi(ω)dωc_i=\int_\Omega f(\omega)B_i(\omega)\mathrm{d}\omega

也就是上面提到的Product Integral
求其系数的过程也被称为投影,即原函数在该基函数上的投影

注:如果用前4阶球谐函数的16个基函数表示一个函数,则最高只能保留l=3l=3,即第4层基函数所对应的频率,而丢失更高频的信息

下面利用球谐函数进行Diffuse物体的Shading
首先,如前所述,对于环境光预先进行滤波在对Shading Point的镜面反射方向进行一次Query,就约等于对Shading Point镜面反射方向邻近角度进行多次Query后求加权平均
PrefilteringSingleQuery

球谐函数的良好性质:
SHproperties
注:SH基函数的旋转都可以用其同阶的基函数的线性组合来描述

SH for Diffuse Case

Diffuse物体的BRDF包含高频信息较少,利用球谐函数的前3阶就可以恢复出效果较好的结果(其BRDF在3阶之后的基函数上的投影接近于0)
SH_BRDF
又如前所述,两个函数进行Product Integral得到的函数的频率取决于两者中频率较低的函数。于是Diffuse物体的BRDF相当于一个低通滤波器。也就是说即便使用一个带有高频信息的光源照亮Diffuse的物体,Diffuse物体表面也不会保留多少高频信息
综上所述,Lighting项也完全不需要描述太高的频率。对于Diffuse物体,其光照也只需要用前3阶球谐函数就可以恢复出很好的结果
Lighting_SH

Percomputed Radiance Transfer(PRT)

在实时渲染中,渲染方程通常被拆成三个部分来理解,如下:
RenderingEquation3Parts
其中,Lighting项、Visibility项、BRDF项都可以描述成定义在球面上的函数
注:Environment Lighting中Visibility项与光源的情况无关,因为默认光从无限远处打过来;BRDF项原本是一个4维函数,但是此处表示的是给定出射方向的BRDF,因此只剩下2维可变

利用上述方式理解渲染方程,就是对于某个Shading Point,将对其Lighting项、Visibility项、BRDF项求Product Integral就是该Shading Point的颜色
但是这种方式如果暴力求解,其计算量很大
引入预计算方法,提前计算好与Lighting无关的部分,则渲染时的计算速度就能大大提高

首先,假设场景中只有光源可以运动,其他物体全都静止。则对于渲染方程,还有另一种理解方式,即分为Lighting和Light Transport两部分:
RenderingEquationPRT
对于Lighting项,可以拆成一系列球谐函数基函数表示:

L(i)liBi(i)L(i)\approx\sum l_i B_i(i)

而Light Transport项对于当前的Shading Point是固定不变的(场景中只有光源可以运动,其他物体全都静止),则Light Transport项在渲染前就可以进行预计算,并且Light Transport项仍是一个球面函数,也可以表示成一系列球谐基函数

PRT(Diffuse Case)

对于Diffuse物体,其BRDF是一个常数,将其提出渲染方程积分之外。将其Lighting项表示成一系列球谐基函数的和,然后将求和操作提出积分之外(在RTR中,通常认为积分与求和的顺序总是可以互换)。此时,剩下部分的积分可以看作是Light Transport项(除了被提出的BRDF)在基函数上的投影,因此该积分也可预计算,于是有了下面的形式,此时对于一个Shading Point的渲染方程求解被转化为求Lighting和Light Transport的向量积,并且由于包含了Visibility项,其渲染的结果是包含Shadow的:
PRT_DiffuseCase

上述预计算的代价就是,其场景(除了光源)不可变化。(光源的旋转变换可以对应到基函数的旋转变换,因此只需再乘上一个旋转变换的系数即可)

对于上面的预计算,有另一种理解方式:先计算出用各个基函数描述环境光照来照亮场景时的情况作为基本情况,在渲染过程中将当前光照照亮场景的情况表示成前面计算完成的基本情况的组合
PRTanotherobservation

注:PRT通常有两种方式,可以先对三角形的每个顶点预计算出渲染方程的结果,对其内部的Fragment进行插值;也可以先根据顶点的Lighting向量和Light Transport向量,插值出内部点的两种向量,在渲染过程中对其内部的Fragment求向量积

换一种思路进行PRT,可将Lighting和Light Transport分别各自进行预计算,如下:
PRTseparately
此时变成了一个双重求和,看似变成了O(n2)O(n^2)的复杂度与前述的推导不符,但是要记得球谐函数具有正交性,其同阶的不同基函数相互投影为0。因此,上式的积分只有在BpB_pBqB_q是相同基函数时,其积分结果为1,否则为0,就相当于在向量积后面乘上一个对角矩阵

PRT(Glossy Case)

对于Glossy材质,其BRDF与观察方向有关,因此其投影到基函数之后的Light Transport项是一个出射方向oo的函数,而对于每一个出射方向oo都可以将其投影至基函数也写成一系列基函数的和,因此每一个oo对应一个基函数系数vector,将所有oo合并到一起则Light Transport项预计算结果变成了一个矩阵(而不是向量),此时渲染方程的计算就变成了Lighting向量乘上Light Transport矩阵
(其实我不太明白,出射方向oo不是二维的吗,那如此形成的Light Transport项不应该是三维的吗,为什么是一个矩阵)
PRTglossy

由于Glossy材质的BRDF比Diffuse材质的BRDF略复杂些,因此有时只用前3阶SH是不够的,通常会用到4阶、5阶,在科研领域可能会用到10阶甚至更高

关于计算量:(以用前4阶SH为例)

  1. Diffuse Case:
    AT each point: dot-product of size 16
  2. Glossy Case:
    At each point: vector(16) * matrix(16*16)

PRT Multi-Bounce

PRT甚至可以预计算光线多次弹射的情况来实现全局光照效果
对于光线弹射的情况,可以借鉴正则表达式的写法,表示如下:
PRT-MultiBounce
其中,LELE就是光线直接打到观察点,中间的D,S,G分别表示光线打到Diffuse,Specular,Glossy物体
由上述表达式可以见,无论中间的弹射过程多么复杂,Lighting和Light Transport始终是分开的,因此可以对Light Transport中不同的光线弹射情况进行预计算,这虽然增加了预计算的开销,但是实时渲染的速度与其预计算的复杂程度无关,因此可以用几乎同样快速的方式渲染出全局光照

那么如何进行考虑多次弹射的Light Transport预计算呢?回想前文所述,PRT可以理解成用SH基函数描述的光源照亮场景,将各基函数光源渲染情况作为基本情况,然后再用SH基函数的组合描述实际光源,将实际光源照亮场景的渲染结果转换为用上述基本情况的组合来描述。
可见,预计算Lighting和Light Transport实际上就是进行渲染,只不过是用SH基函数描述的光源。那么对于考虑多次弹射的情况,其实就是在上述预渲染的过程中利用Ray Tracing等方法计算出SH基函数光源照亮场景且考虑多次弹射的基本情况,再用这种基本情况来组合描述实际光源多次弹射照亮场景的情况

Follow Up Works

FollowUpWorks

More Basis Functions

MoreBasisFunctions

Wavelet

HaarWavelet
小波相比SH,可以全频率地描述一个函数

用小波描述环境光照通常使用Cube Map
每一面不断进行小波变换,留下高频信息,将低频信息放在左上角,最终记录各频率对应的小波系数
WaveletTransform
于是我们会发现,高频信息实际上很少,大部分系数为0,于是我们得到了一个环境光照的高质量压缩
注:JPEG格式利用了类似这样思想的压缩方式,称为离散余弦变换(Discrete Cosine Transform, DCT)

用相同的存储量,相比SH,Wavelet可以恢复出较高频率的信息,但是,小波像SH那样不支持快速旋转
SHvsWavelet