(翻译)Physically Based Rendering i

2023-06-08  本文已影响0人  Dragon_boy

光照

光照环境的准确性和相干性在营造可信的视觉效果时至关重要,在调查已有的渲染引擎后和相关文献后,不难发现,相干性基本达不到。

例如虚幻中,使用者可以调节以流明为单位的点光源的亮度,但平行光却使用了另外一种不知名的单位。为了匹配点光源5000流明的亮度,平行光需要设为10的亮度,这样子不匹配单位,使用者容易混乱。如果完全使用单一的自创单位,的确能避免问题,不过复用光照设备就比较困难了。例如,外景中,10亮度的平行光作为太阳,其它所有的光源都要与之匹配,再迁到室内,相同的数值可能就太亮,过曝了。

因此,我们需要一开始就让所有的光照在默认情况下就符合常理,同时给予高度的自由以供修改,我们会提供多种光源,两类,直接光和间接光:

直接光包括精确光源,光域光,区域光。间接光包括基于图像的光源(IBL),同时涵盖局部和平行光照探针。

局部光照探针对移动设备来说花销太大,我们会主要讨论平行光照探针。

单位

接下来会讨论如何去实现多种不同的光源,这是会涉及到的各种单位:

光度术语 标记 单位
发光功率 \Phi 流明(lm)
发光强度 I 坎德拉(cd)或\frac{lm}{sr}sr是球面度)
照度 E 勒克司(lx)或\frac{lm}{m^2}
亮度 L 尼特(nt)或\frac{cd}{m^2}
辐射功率 \Phi_e 瓦特(W)
光视效能 \eta 流明每瓦特(\frac{lm}{W})
发光效率 V 百分比(\%)

为了得到恰当的相干照明,我们必须使用那些可以代表不同灯光强度比率的单位,强度的跨度很广,比如家用灯泡800流明,白天太阳光120000流明。

最简单的方式就是使用物理灯光单位,这样可以复用灯光设备,使用物理灯光单位也意味着可以用基于物理的摄像机。

下面展示了我们为各种灯光类型选用的单位:

灯光类型 单位
平行光 照度(lx\frac{lm}{m^2})
点光 发光功率(lm)
聚光 发光功率(lm)
光域光 发光强度(cd)
遮罩光域光 发光功率(lm)
区域光 发光功率(lm)
基于图像的光 亮度(\frac{cd}{m^2})

关于辐射功率的单位

虽说超市里卖的灯泡包装上标注的基本都是流明亮度,但我们日常生活中还是常常用瓦特来表明灯光的亮度。瓦特的量只是说明灯泡会使用多少能量,而不是表明亮度,这对我们也有帮助,毕竟节能灯越来越多了,功率不等于亮度。

然而啊,美术估计还是适应用功率来衡量灯光的亮度,我们得提供这样的选择。流明(发光功率)和瓦特(辐射功率)的转换如下:

\Phi=\Phi_e \eta \tag{52}

上述等式中,\eta是灯光的光视效能,代表每单位瓦特的流明量,其最大的可能值是683\frac{lm}{W},我们同样可以用发光效率V来表示:

\Phi=\Phi_e683 \times V\tag{53}

下图是个转化索引:

灯光类型 效能\eta 效率V
白炽灯 14-35 2-5%
LED 28-100 4-15%
荧光灯 60-100 9-15%

灯光单位验证

使用物理灯光单位的一个重大优势在于我们可以用物理手段检验我们的公式,我们可以用相关设备去测量。

照度

照度可以用入射光计量表来测量,比如下图的世光L-478D:

计量表使用一个白色的圆顶来接收光到达表面的照度,我们需要让圆顶恰当朝向某个要测量的方向,比如,垂直面对太阳和水平面对太阳,结果肯定不一样。

亮度

表面的光亮度可以用亮度计量计测量,也经常被称为点亮度计,它使用一块直板测量单一方向的入射光。把我们的入射光计量计用世光5°反光镜改装一下就行。

发光强度

灯光的发光强度无法直接测量,但如果我们知道光源和测量设备之间的距离,可以通过测量到的照度计算得到,公式如下:

I=E \cdot d^2 \tag{54}

直接光照

我们已经为各种灯光类型定义了灯光单位,不过对于光照等式的结果的单位还没什么着落,使用物理灯光单位意味着我们需要在shader里计算亮度值,因此每时每刻我们都需要所有的光照函数中计算亮度L_{out}。亮度取决于照度E和BSDF f(v, l)

L_{out}=f(v, l)E \tag{55}

平行光

平行光最重要的目的就是模拟室外环境,充当太阳月亮的角色。现实生活其实并不存在真正意义上的平行光,不过距离遥远的光源可以看作是平行光。

这一近似对漫反射计算很友好,但高光可能就不太对了。寒霜引擎的做法是将太阳平行光看作是一个圆盘区域光,不过我们认为这样做得不偿失。

之前我们说过为平行光选择的单位是照度(lx),某种程度上是因为太阳和天空的照度很容易在晚上查到,这也可以简化亮度公式:

L_{out}=f(v,l)E_{\perp}\big<n\cdot l\big>\tag{56}

其中E_{\perp}是垂直于光源的表面所接受的光源照度,如果用平行光代表太阳,E_{\perp}就是垂直于太阳光方向的表面所接受的太阳照度。

下表提供了部分太阳和天空照度的参考值,在加利福尼亚三月的某个晴空万里的日子测量的:

灯光 10am 12pm 5:30pm
Sky_{\perp} + Sun_{\perp} 120000 130000 90000
Sky_ {\perp} 20000 25000 9000
Sun_ {\perp} 100000 105000 81000

动态平行光的实时开销很小的,实现如下:

vec3 l = normalize(-lightDirection);
float NoL = clamp(dot(n, l), 0.0, 1.0);

// lightIntensity is the illuminance
// at perpendicular incidence in lux
float illuminance = lightIntensity * NoL;
vec3 luminance = BSDF(v, l) * illuminance;

下图的场景展示了模拟正午阳光的效果:

精确光

我们的引擎提供两种精确光源,点光和聚光,这两种光源实际上在物理学里不太准确,因为:

  1. 它们本身没有尺寸概念

  2. 它们并不遵循平方反比原则

第一个问题可以用区域光来解决,不过实际上还是会尽量用开销小的点光。

第二个问题也好解决,对于某一精确光源,感知到的强度会随与观察者之间的距离增加,根据距离的平方成比例衰减。对于平方衰减的精确光,公式55里的E的等式如下,其中d是点到光源的距离:

E=L_{in}\big<n \cdot l\big>=\frac{I}{d^2}\big<n\cdot l\big>\tag{57}

点光和聚光的区别在于E的计算,特别是如何根据发光功率\Phi计算发光强度l

点光

点光代表空间中的某一位置:

点光的发光功率是通过其发光强度在光源的立体角范围内积分所计算得到的,公式如下,由此可通过发光功率得到发光强度:

\Phi=\int_{\Omega}Idl=\int_{0}^{2\pi}\int_{0}^{\pi}Id\theta d\phi=4\pi I\tag{58}
I=\frac{\Phi}{4\pi}

有之前推导出的公式,我们可以得到点光关于发光功率的亮度等式:

L_{out}=f(v,l)\frac{\Phi}{4\pi d^2}\big<n \cdot l\big>\tag{59}

下图是点光加上平方衰减照亮某个简单场景的效果,衰减做了夸大处理:

聚光

聚光由位置,方向,还有两个锥形角\theta_{inner}\theta_{outer}确定,这两个角的作用是用来定义聚光角衰减的,因此需要同时考虑两种衰减。

等式60描述了聚光的发光功率的计算过程,和点光类似,用了\theta_{outer},范围是[0..\pi]

\Phi=\int_{\Omega}Idl=\int_0^{2\pi}\int_0^{\theta_{outer}}Id\theta d\phi=2\pi (1-cos\frac{\theta_{outer}}{2})I\tag{60} $$$$ I=\frac{\Phi}{2\pi(1-cos\frac{\theta_{outer}}{2})}

因为这个等式在物理上正确,这会让聚光有点难以使用:改变外角会改变照度。下图展示了不同外角下的效果,左侧为55°,右侧为15°,可以看到外角变小,照度等级提升:

因为外角和照度有所关联了,美术就不好去单独调节了,因此需要提供一个参数来关闭这种联系,下面的等式搞好达成目的:

\Phi=\pi I \tag{61}
I=\frac{\Phi}{\pi}

这样搞就不会有外角变化,照度也变化的问题了:

那如何让它在物理上正确呢?如果聚光的反射体可以被完美吸收光的漫反射遮罩替代就行。

聚光的亮度公式可以用两种方式表示:

L_{out}=f(v,l)\frac{\Phi}{\pi d^2}\big<n\cdot l\big>\lambda(l)\tag{62}

L_{out}=f(v,l)\frac{\Phi}{2\pi (1-cos\frac{\theta_{outer}}{2})d^2}\big<n \cdot l\big>\lambda(l)\tag{63}

\lambda( l )是聚光角衰减因数:

\lambda(l)=\frac{l \cdot spotDirection - cos\theta_{outer}}{cos\theta_{inner} - cos\theta_{outer}}\tag{64}

衰减函数

要让光源基于物理,遵循平方反比原则的衰减是必须的。不过这一简单的数学公式在具体实现时并不实用:

  1. 因为距离的问题,可能会除以0。

  2. 光源的范围理论上是无限的,这也就意味着每个像素需要计算场景中所有光源对它的影响。

第一个问题很好解决,可以假设精确光源其实是有一定大小的,我们可以假定它至少是个半径1cm的球,这样可以优化公式:

E=\frac { I } { max( d^2, 0.01^2 ) } \tag{65}

第二个问题的话,可以通过设定每个光源的范围来解决,这样子做的话,美术可以快速确定光照区域,引擎也可以据此信息快速地剔除光源,我们也可以有意识地调整光源范围来帮助引擎进行优化。

数学上,光的照度需要在接近边界时平缓过渡到0,下面是相应的优化后的的公式,r是光源的范围:

E=\frac{I}{max(d^2,0.01^2)}\Big< 1-\frac{d^4}{r^4}\Big >^2\tag{66}

下图是精确光源的实现,注意灯光强度是以cd为单位的发光强度I,计算自CPU端传来的发光功率。这段代码并没有优化,有些计算可以挪到CPU进行,例如灯光范围的平方反比,聚光的缩放和锥形角。

float getSquareFalloffAttenuation(vec3 posToLight, float lightInvRadius) {
    float distanceSquare = dot(posToLight, posToLight);
    float factor = distanceSquare * lightInvRadius * lightInvRadius;
    float smoothFactor = max(1.0 - factor * factor, 0.0);
    return (smoothFactor * smoothFactor) / max(distanceSquare, 1e-4);
}

float getSpotAngleAttenuation(vec3 l, vec3 lightDir,
        float innerAngle, float outerAngle) {
    // the scale and offset computations can be done CPU-side
    float cosOuter = cos(outerAngle);
    float spotScale = 1.0 / max(cos(innerAngle) - cosOuter, 1e-4)
    float spotOffset = -cosOuter * spotScale

    float cd = dot(normalize(-lightDir), l);
    float attenuation = clamp(cd * spotScale + spotOffset, 0.0, 1.0);
    return attenuation * attenuation;
}

vec3 evaluatePunctualLight() {
    vec3 l = normalize(posToLight);
    float NoL = clamp(dot(n, l), 0.0, 1.0);
    vec3 posToLight = lightPosition - worldPosition;

    float attenuation;
    attenuation  = getSquareFalloffAttenuation(posToLight, lightInvRadius);
    attenuation *= getSpotAngleAttenuation(l, lightDir, innerAngle, outerAngle);

    vec3 luminance = (BSDF(v, l) * lightIntensity * attenuation * NoL) * lightColor;
    return luminance;
}

光域光

精确光挺好用,就是控制选项太少了。工业光照领域就设计了一套光照系统来满足人们的各种需求,考虑了以下几点:

我们目前为止介绍的光照系统已经满足前两点了,我们还需要找到定义空间中灯光分布的方法。灯光分布在室内场景,部分室外场景甚至是道路照明中都很重要,下图展示了美术可以控制的灯光分布在场景中的效果:

光域光使用光域配置文件来描述强度分布,有两种流行格式,IES(美国照明工程协会)和EULUMDAT(欧洲流明数据格式),我们只考虑IES。有很多工具和引擎支持IES配置文件,如虚幻,寒霜,Renderman,Maya,Killzone。此外,IES灯光配置文件通常供应给电灯和灯具制造商。光域配置文件在规划灯具时非常有用,光源会被遮挡,这样会让光源向特定方向发光,造造形。

IES配置文件存储了球面各个角度的发光强度,这一球面坐标系统会展开在光域网上,方便查看,可以用类似IESviewer的工具可视化。下图展示了皮克斯在Renderman里使用的一个X箭头样式的IES文件,以光域网形式展开:

简要来说,IES配置文件存储了环绕光源的各个角度下的发光强度,单位为坎德拉。对于每个测量好的水平角度,会提供对应的一系列不同垂直角度下的发光强度,因此往往是水平对称的测算光源。上述的X箭头文件就是个很好的例子,垂直方向强度不一,但水平对称。垂直角的范围在0-180°,水平角的范围在0-360°。

下图展示了Pixar在Renderman里所使用的提系列IES配置文件,用我们的工具渲染了一下:

IES配置文件可以直接用到任何精确光源上,为了实现这一点,我们需要处理IES文件,然后生成一个纹理形式的光域文件。出于对性能的考虑,我们会生成一张1D文件,用来存储在某一特定的垂直角度下,所有水平角度的发光强度的平均值,每个像素对应一个垂直角度。实际上要是想最为正确地表示光域灯,我们是得使用2D纹理的,但是大多数的光远基本是在水平面上对称的,用1D也不是不行。纹理中存储的值通过定义在IES文件内的最大光强的倒数进行了标准化处理,这就让我们可以用任意浮点数类型来存储,或者考虑下性能,牺牲下精度,存储在8位图中,这样做我们也可以把光域文件当作是个遮罩:

光域配置文件作为遮罩

对于任意精确光源,美术通过调节灯光的发光功率调节发光强度,然后除由IES文件计算得到的光强度。IES文件本身包含一个发光强度,但这只是针对纯灯泡的,现实情况还得考虑灯具的遮挡。为了测量灯具的强度,我们使用来自配置文件里的强度在单位球里进行蒙特卡洛积分计算。

光域配置文件

发光强度来自于文件本身,所有采样自1D纹理的值会乘上最大强度,我们还提供了一个乘数来方便调节。

光域配置可以当作简单的衰减来实时使用,公式如下:

L_{out}=f(v,l)\frac{I}{d^2}\big<n\cdot l\big>\Psi(l)\tag{67}

\Psi(l)就是光域衰减函数,它依赖于光线方向和光源朝向,光源朝向和聚光类似。

光域衰减函数实现很简单,在精确光的实现里加一个衰减因数就行:

float getPhotometricAttenuation(vec3 posToLight, vec3 lightDir) {
    float cosTheta = dot(-posToLight, lightDir);
    float angle = acos(cosTheta) * (1.0 / PI);
    return texture2DLodEXT(lightProfileMap, vec2(angle, 0.0), 0.0).r;
}

vec3 evaluatePunctualLight() {
    vec3 l = normalize(posToLight);
    float NoL = clamp(dot(n, l), 0.0, 1.0);
    vec3 posToLight = lightPosition - worldPosition;

    float attenuation;
    attenuation  = getSquareFalloffAttenuation(posToLight, lightInvRadius);
    attenuation *= getSpotAngleAttenuation(l, lightDirection, innerAngle, outerAngle);
    attenuation *= getPhotometricAttenuation(l, lightDirection);

    float luminance = (BSDF(v, l) * lightIntensity * attenuation * NoL) * lightColor;
    return luminance;
}

光强在CPU端计算,并且取决于是否要将光域作为遮罩使用。

float multiplier;
// Photometric profile used as a mask
if (photometricLight.isMasked()) {
    // The desired intensity is set by the artist
    // The integrated intensity comes from a Monte-Carlo
    // integration over the unit sphere around the luminaire
    multiplier = photometricLight.getDesiredIntensity() /
            photometricLight.getIntegratedIntensity();
} else {
    // Multiplier provided for convenience, set to 1.0 by default
    multiplier = photometricLight.getMultiplier();
}

// The max intensity in cd comes from the IES profile
float lightIntensity = photometricLight.getMaxIntensity() * multiplier;

灯光参数化

我们的目的是让灯光参数更加直观易用,本着这样的原则,我们决定将灯光颜色从灯光强度中分离出来,定义为一个线性RGB颜色参数。

完整的灯光参数列表如下:

参数 定义
Type 平行,点,聚,区域
Direction 用于平行光,聚光,光域点光,线性和管状光
Color 灯光发出的颜色,线性RGB。在工具里可以以sRGB或者色温的形式存在
Intensity 灯光亮度,单位取决于灯光类型
Falloff radius 最大影响范围
Inner angle 聚光锥形内角
Outer angle 聚光锥形外角
Length 区域光的长度,用于创建线性或管状光
Radius 区域光的半径,用于创建球形或管状光
Photometric profile 存储光域灯光配置的纹理,只对精确光有效
Masked profile 布尔值,用于判断IES配置是否作为遮罩使用
Photometric multiplier 光域灯光的亮度乘数(如果IES作为遮罩使用则关闭)

注意,为了简化实现,所有的发光效率值会在送往shader前转化为发光强度(cd)。灯光类型可由参数的值判断,比如点光的长度,半径,内角和外角都是0。

色温

现实生活中的人造光源基本由色温去定义,单位是开尔文(K),光源的色温定义是:理想黑体辐射源向光源辐射可见光后的温度。方便起见,工具里需要提供调节光源颜色的选项,作为色温(通常范围是1000K-12500K)。

为了从温度中计算RGB值,我们可以用普朗克轨迹,如下图,这一轨迹是发光黑体在色度空间(CIE 1931)里随温度变化的颜色变化路径。

从这一轨迹中计算RGB值最简单的方法就是使用如下的等式(色度空间CIE 1960(UCS)),其中T是温度,uv是UCS里的坐标:

u(T)=\frac{0.860117757+1.54118254\times 10^{-4}T+1.28641212\times 10^{-7}T^2}{1+8.42420235\times 10^{-4}T+7.08145163\times 10^{-7}T^2}\tag{68}
v(T)=\frac{0.317398726+4.22806245\times 10^{-5}T+4.20481691\times 10^{-8}T^2}{1-2.89741816\times 10^{-5}T+1.61456053\times 10^{-7}T^2}

在1000K到15000K范围内,这一近似结果大概是9\times10^{-5}。根据CIE 1960空间结果,可以计算得到CIES 1931空间(xyY)的值:

x=\frac{3u}{2u-8v+4}\tag{69}
y=\frac{2v}{2u-8v+4}

上述公式针对黑体色温,与标准光源的色温相关。如果我们向计算D系列中标准CIE光源的精确色度坐标,可以使用下列等式:

x=\left\{ \begin{aligned} &0.244063+0.09911\frac{10^3}{T}+2.9678\frac{10^6}{T^2} −4.6070\frac{10^9}{T^3} &\ 4000K\le T\le7000K\\ &0.237040+0.24748\frac{10^3}{T}+1.9018\frac{10^6}{T^2}−2.0064\frac{10^9}{T^3} &\ 7000K\le T\le25000K \end{aligned} \right. \tag{70}
y=−3x^2+2.87x−0.275

根据xyY空间结果,可以转换得到CIE XYZ空间结果:

X=\frac{xY}{y}\tag{71}
Z=\frac{(1-x-y)Y}{y}

根据需要,我们令Y=1,这样可以用一个三维矩阵将XYZ空间的结果转换到线性RGB空间:

\left[ \begin{matrix} R\\ G\\ B \end{matrix} \right] = M^{-1} \left[ \begin{matrix} X\\ Y\\ Z \end{matrix} \right] \tag{72}

M由目标RGB颜色空间的基色计算得到,下面的等式展示了使用针对sRGB空间的反转矩阵进行转化的过程:

\left[ \begin{matrix} R\\ G\\ B \end{matrix} \right] = \left[ \begin{matrix} 3.2404542 \ −1.5371385 \ 0.4985314\\ −0.9692660\ 1.8760108 \ 0.0415560\\ 0.0556434 \ −0.2040259\ 1.0572252 \end{matrix} \right] \left[ \begin{matrix} X\\ Y\\ Z \end{matrix} \right] \tag{73}

这些等式的结果是一个在sRGB空间中的线性RGB三元组。由于我们比较关注结果的色度,我们必须加入一个标准化步骤来避免截断大于1的值,搞花颜色:

\hat C_{linear}=\frac{C_{linear}}{max(C_{linear})}\tag{74}

最后需要执行一步光电转换来得到可展示的值:

C_{sRGB}= \left\{ \begin{aligned} &12.92 \times \hat C_{linear} \ &\hat C_{linear}\le 0.0031308\\ &1.055\times \hat C_{linear}^{\frac{1}{2.4}}-0.055 \ &\hat C_{linear}>0.0031308 \end{aligned} \right. \tag{75}

下图展示了1000K-12500K的相关色温,所有颜色都将CIE \ D_{65}作为白平衡点:

同理,1000K-12500K的CIE D系列标准发光体:

没有标准化的相关色温:

下表将不同光源的相关色温以色板形式表现了,这些颜色与D_{65}白点相关,因此我们看到的颜色与你的显示器有关。

温度(K) 光源 颜色
1,700-1,800 火柴焰 \color{ #ff7d00}{█████}
1,850-1,930 蜡烛焰 \color{ #ff8701}{█████}
2,000-3,000 日出日落时的太阳 \color{ #ffa64c}{█████}
2,500-2,900 家用钨丝灯泡 \color{ #ffb05e}{█████}
3,000 千瓦钨丝灯 \color{ #ffb86f}{█████}
3,200-3,500 石英灯 \color {#ffbf7b}{█████}
3,200-3,700 荧光灯 \color{ #ffbf7b}{█████}
3,275 2千瓦钨丝灯 \color{ #ffc180}{█████}
3,380 5千瓦,10千瓦钨丝灯 \color{#ffc486}{█████}
5,000-5,400 正午太阳 \color{ #ffe9d7}{█████}
5,500-6,500 日光(太阳+天空) \color{#fff3f1}{█████}
5,500-6,500 太阳穿过云层 \color{#fff3f1}{█████}
6,000-7,500 阴天 \color{#faf6ff}{█████}
6,500 显示器的白点 RGB \color{#fff8fe}{█████}
7,000-8,000 户外阴影区 \color{ #ebecff}{█████}
8,000-10,000 局部多云 \color{#d6e0ff}{█████}

预曝光灯光

使用基于物理的渲染和物理灯光单位会面临一个问题,该如何在代码里处理范围跨度巨大的值?假设在shader中全精度计算后得到结果,我们还是得将结果存储在某一大小的缓冲中,比如RGB16F,最简单的方法就是存储前线应用摄像机曝光,如下:

fragColor = luminance * camera.exposure;

这一操作需要在光照计算时使用单精度浮点值,不过出于性能的考虑,我们还是希望使用半精度浮点,但这对照度和亮度值的适配不太好,可能会溢出。解决方案是让光源就先自己曝光下,就不在光照计算后进行了,这一操作刚好可以在CPU里进行,节省开销,当然GPU也可以。

// The inputs must be highp/single precision,
// both for range (intensity) and precision (exposure)
// The output is mediump/half precision
float computePreExposedIntensity(highp float intensity, highp float exposure) {
    return intensity * exposure;
}

Light getPointLight(uint index) {
    Light light;
    uint lightIndex = // fetch light index;

    // the intensity must be highp/single precision
    highp vec4 colorIntensity  = lightsUniforms.lights[lightIndex][1];

    // pre-expose the light
    light.colorIntensity.w = computePreExposedIntensity(
            colorIntensity.w, frameUniforms.exposure);

    return light;
}

实际中预曝光操作如下:

基于图像的灯光

现实生活中,有直接来自光源的光,也有经环境反射的间接光,在某种意义上一个物体周围的环境也看做是光源,这些环境光可以用图像编码,称为基于图像的光照或者间接光照。

不过这种基于图像的光照有一些先置,比如得提前准备好需要的环境图像,而且得预处理一下才能够用于光照计算,图像可以从现实取得,或者由引擎离线或实时生成。

这里我们只关注平行环境探针,也就是假设光源来自无限远的地方。

环境入射的光照称之为日照强度(E),反射的光称之为辐射率(L_{out})。

辐射率L_{out}关于材质模型BRDF的计算公式如下(\Theta代表各种参数,如金属度,粗糙度等):

L_{out}(n,v,\Theta)=\int_{\Omega}f(l,v,\Theta)L_{\perp}(l)\big<n\cdot l\big>dl\tag{76}

注意这里只考虑了宏观情况,本质上我们需要对来自所有方向的点光应用BRDF,然后把结果编码进IBL。

IBL类型

现代渲染引擎包含四种常用的IBL类型:

另外我们得区分静态和动态IBL,比如实现动态昼夜循环需要动态计算平行灯光探针。平面反射和屏幕空间反射本身都是动态的。

IBL单位

灯光必须使用物理单位,IBL会使用亮度单位\frac{cd}{m^2}

高动态范围图像并不好处理,摄像机并不能存储环境亮度,而是个和环境亮度相关的值,与设备型号有关,因此我们会提供一个乘数参数来进行调节,尽量去达到环境本身的亮度值。

为了给IBL恰当地重构HDRI的亮度,拍照时得记录些额外信息:

处理灯光探针

之前讲过,IBL的辐射率是由沿表面半球积分得到的,这显然很耗,那我们就得把灯光谈真转化到适合实时计算的形式。

下面这俩就是加速探针计算会用到的技术:

平行灯光探针

漫反射BRDF积分

根据兰伯特BRDF(并不取决于\vec{l},\vec{v},\theta,因此L_d(n,v,\theta)=L_d(n,\sigma)),得到辐射率:

f_{d}(\sigma)=\frac{\sigma}{\pi}\\ \begin{aligned} L_d(n,\sigma) &=\int_{\Omega}f_d(\sigma)L_{\perp}(l)\big<n\cdot l\big>dl \\ &= \frac{\sigma}{\pi}\int_{\Omega}L_{\perp}(l)\big<n\cdot l\big>dl \\ &= \frac{\sigma}{\pi}E_d(n) \ 其中辐照度E_d(n)=\int_{\Omega}L_{\perp}(l)\big<n\cdot l\big>dl \end{aligned}

或者离散域:

E_d(n)=\sum_{\forall i\in image}L_{\perp}(s_i)\big<n \cdot s_i\big>\Omega_s

\Omega_s是采样i的立体角(针对立方体贴图,利用\frac{2\pi}{6 \cdot width \cdot height}得到)。

E_d可以预先计算好存储在立体贴图中。image可以是立方贴图或者等距投影贴图。\frac{\sigma}{\pi}独立于IBL,单独计算,实时添加。

基于图像的环境:

使用兰伯特BRDF后的基于图像的光照贴图:

辐照度可以由球谐分解近似得到结果,这样可以就可以实时计算。

球谐分解在概念上和傅里叶变换类似,即在频域标准坐标系内遍历信号。有一些特质:

实际操作中只需要4或9个系数就足够了。

9个:

4个:

实际操作中,我们会使用\big<cos\theta\big>预卷积L_{\perp},并通过基础缩放系数K_l^m来预缩放这些系数,实现如下:

vec3 irradianceSH(vec3 n) {
    // uniform vec3 sphericalHarmonics[9]
    // We can use only the first 2 bands for better performance
    return
          sphericalHarmonics[0]
        + sphericalHarmonics[1] * (n.y)
        + sphericalHarmonics[2] * (n.z)
        + sphericalHarmonics[3] * (n.x)
        + sphericalHarmonics[4] * (n.y * n.x)
        + sphericalHarmonics[5] * (n.y * n.z)
        + sphericalHarmonics[6] * (3.0 * n.z * n.z - 1.0)
        + sphericalHarmonics[7] * (n.z * n.x)
        + sphericalHarmonics[8] * (n.x * n.x - n.y * n.y);
}

高光BRDF积分

IBL辐照度和BRDF计算得到的辐射率L_{out}如下:

L_{out}(n,v,\Theta)=\int_{\Omega}f(l,v,\Theta)L_{\perp}(l)\big<n \cdot l\big>\partial l\tag{77}

我们通过f(l,v, \Theta)\big<n \cdot l\big>得到L_{\perp}的卷积结果,即使用BRDF作为算子对环境滤波。粗糙度越高,高光反射越模糊。

f扩开,组合起来,得到:

L_{out}(n,v,\Theta)=\int_{\Omega}D(l,v,\alpha)F(l,v,f_0,f_{90})V(l,v,\alpha)L_{\perp}(l)\big<n \cdot l\big>\partial l\tag{78}

这一表达式依赖于积分内的v,\alpha,f_0,f_{90},因此计算相当耗能,不适合在移动设备上实时计算。

简化BRDF积分

因为实在是找不到计算L_{out}的近似和简易方法,我们会用另一个简化的等式替代:\hat I,并假设v=n,即观察方向等于法线方向,当然,这样去近似的化会破坏掉卷积中和观察方向相关的效果,如视距模糊。

这样子简化会影响结果中常量值大小,我们可以用一个缩放因数K来修正。

I是个可以拆开的积分乘法,即I(g()f())=I(g())I(f())

如此:

I(f(\Theta)L_{\perp})\approx \tilde I(f(\Theta)L_{\perp}) \tag{79}
\tilde I(f(\Theta)L_{\perp}) = K\times \hat I(f(\Theta)L_{\perp})
K=\frac{I(f(\Theta))}{\hat I(f(\Theta))}

由上可知,当L_{\perp}是常量时,\tilde I等于I,就可以得到正确的结果:

\begin{aligned} \tilde I(f(\Theta)L_{\perp}^{constant}) &=L_{\perp}^{constant}\hat I(f(\Theta)) \frac{I(f(\Theta))}{\hat I(f(\Theta))} \\ &=L_{\perp}^{constant} I(f(\Theta))\\ &=I(f(\Theta)L_{\perp}^{constant}) \end{aligned}

同样也可以证明在v=n时结果也正确,因为I=\hat I:

\begin{aligned} \tilde I(f(\Theta)L_{\perp}) &=I(f(\Theta)L_{\perp}) \frac{I(f(\Theta))}{I(f(\Theta))}\\ &=I(f(\Theta)L_{\perp}) \end{aligned}

最后,代入L_{\perp}=\bar {L_{\perp}}+(L_{\perp}-\bar {L_{\perp}})=\bar {L_{\perp}}+\Delta L_{\perp}

\begin{aligned} \tilde I(f(\Theta)L_{\perp}) &=\tilde I[f(\Theta)(\bar {L_{\perp}}+\Delta L_{\perp})]\\ &=K\times \hat I[f(\Theta)(\bar {L_{\perp}}+\Delta L_{\perp})]\\ &=K\times [\hat I(f(\Theta)\bar {L_{\perp}})+\hat I(f(\Theta)\Delta {L_{\perp}})]\\ &=K\times \hat I(f(\Theta)\bar {L_{\perp}})+K\times \hat I(f(\Theta)\Delta {L_{\perp}})\\ &=\tilde I(f(\Theta)\bar {L_{\perp}})+\tilde I(f(\Theta)\Delta {L_{\perp}})\\ &=I(f(\Theta)\bar {L_{\perp}})+\tilde I(f(\Theta)\Delta {L_{\perp}}) \end{aligned}

上面的结果显示平均辐照度结果正确,即I(f(\Theta)\bar {L_{\perp}})

我们可以把上述的近似当作是把辐射率L_{\perp}拆分了两部分,平均值\bar {L_{\perp}}和差值\Delta L_{\perp},先计算平均部分,再加上差值部分:

approximation(L_{\perp})=correct(\bar{L_{\perp}})+simplified(L_{\perp}-\bar{L_{\perp}}\tag{80})

下面是几个积分式子的公式:

\hat I(f(n,\alpha)L_{\perp})=\int_{\Omega}f(l,n,\alpha)L_{\perp}(l)\big<n\cdot l\big>\partial l\tag{81}
\hat I(f(n,\alpha))=\int_{\Omega}f(l,n,\alpha)\big<n\cdot l\big>\partial l
I(f(n,v,\alpha))=\int_{\Omega}f(l,n,v,\alpha)\big<n\cdot l\big>\partial l

上述三个等式可以提前计算好,并且存储在查找表中。

离散域

在离散域中,等式81变为:

\hat I(f(n,\alpha)L_{\perp})\equiv \frac{1}{N}\sum_{\forall i\in image}f(l_i,n,\alpha)L_{\perp}(l_i)\big<n\cdot l\big>\tag{82}
\hat I(f(n,\alpha))\equiv\frac{1}{N}\sum_{\forall i\in image}f(l_i,n,\alpha)\big<n\cdot l\big>
I(f(n,v,\alpha))\equiv\frac{1}{N}\sum_{\forall i\in image}f(l_i,n,v,\alpha)\big<n\cdot l\big>

然而,实际操作中,我们需要考虑概率分布函数使用重要性采样,添加一项,即\frac{4\big<v \cdot h\big>}{D\big<h_i,\alpha\big>\big<n\cdot h\big>}

\hat I(f(n,\alpha)L_{\perp}) \equiv\frac{4}{N} \sum_{i}^N f(l_i,n,\alpha) \frac{\big<v\cdot h\big>} {D(h_i,\alpha)\big<n\cdot h\big>}L_{\perp}(l_i)\big<n\cdot l\big>\tag{83}
\hat I(f(n,\alpha))\equiv\frac{4}{N}\sum_{i}^Nf(l_i,n,\alpha)\frac{\big<v\cdot h\big>}{D(h_i,\alpha)\big<n\cdot h\big>}\big<n\cdot l\big>
I(f(n,v,\alpha))\equiv\frac{4}{N}\sum_{i}^Nf(l_i,n,v,\alpha)\frac{\big<v\cdot h\big>}{D(h_i,\alpha)\big<n\cdot h\big>}\big<n\cdot l\big>

对于\hat I,我们假设过v=n,上述式子简化为:

\hat I(f(n,\alpha)L_{\perp})\equiv \frac{4}{N}\sum_{i}^N \frac{f(l_i,n,\alpha)}{D(h_i,\alpha)}L_{\perp}(l_i)\big<n\cdot l\big>\tag{84}
\hat I(f(n,\alpha)) \equiv\frac{4}{N}\sum_{i}^N \frac{f(l_i,n,\alpha)}{D(h_i,\alpha)}\big<n\cdot l\big>
I(f(n,v,\alpha)) \equiv\frac{4}{N}\sum_{i}^N \frac{f(l_i,n,v,\alpha)}{D(h_i,\alpha)} \frac{\big<v\cdot h\big>}{\big<n\cdot h\big>}\big<n\cdot l\big>

由前两个式子可得到LD(n,\alpha)=\frac{\hat I(f(n,\alpha)L_{\perp})}{\hat I(f(n,\alpha))}

LD(n,\alpha)\equiv \frac {\sum_{i}^N \frac{f(l_i,n,\alpha)}{D(h_i,\alpha)}L_{\perp}(l_i)\big<n\cdot l\big>} {\sum_{i}^N \frac{f(l_i,n,\alpha)}{D(h_i,\alpha)}\big<n\cdot l\big>} \tag{85}
I(f(n,v,\alpha)) \equiv\frac{4}{N}\sum_{i}^N \frac{f(l_i,n,v,\alpha)}{D(h_i,\alpha)} \frac{\big<v\cdot h\big>}{\big<n\cdot h\big>}\big<n\cdot l\big> \tag{86}

到此已经可以离线计算这些等事了的结果了,唯一的难题是我们不知道f_0f_{90}的值,我们可以将这些值实时嵌入,不过对等式85无效,我们需要假设f_0=f_{90}=1

我们还需要处理BRDF里的可见性项,设为1即可,即V=1

f(l_i,n,\alpha)=D(h_i,\alpha)F(f_0,f_{90},\big<v\cdot h\big>)V(l_i,v,\alpha) \tag{87}

代入上述式子,得:

LD(n,\alpha)\equiv \frac {\sum_{i}^N V(l_i,v,\alpha)L_{\perp}(l_i)\big<n\cdot l\big>} {\sum_{i}^N \big<n\cdot l\big>} \tag{88}
I(f(n,v,\alpha)) \equiv\frac{4}{N}\sum_{i}^N F(f_0,f_{90},\big<v\cdot h\big>)V(l_i,v,\alpha) \frac{\big<v\cdot h\big>}{\big<n\cdot h\big>}\big<n\cdot l\big> \tag{89}

将菲涅尔项代入等式89:

F(f_0,f_{90},\big<v\cdot h\big>) = f_0(1-F_c(\big<v\cdot h\big>)) + f_{90}F_c(\big<v\cdot h\big>) \tag{90}
F_c(\big<v\cdot h\big>) = (1-\big<v\cdot h\big>)^5
I(f(n,v,\alpha)) \equiv\frac{4}{N}\sum_{i}^N [f_0(1-F_c(\big<v\cdot h\big>)) + f_{90}F_c(\big<v\cdot h\big>)] V(l_i,v,\alpha) \frac{\big<v\cdot h\big>}{\big<n\cdot h\big>}\big<n\cdot l\big> \tag{91}
\begin{aligned} I(f(n,v,\alpha)) &\equiv f_0\frac{4}{N}\sum_{i}^N (1-F_c(\big<v\cdot h\big>))V(l_i,v,\alpha) \frac{\big<v\cdot h\big>}{\big<n\cdot h\big>}\big<n\cdot l\big> \\ &+f_{90}\frac{4}{N}\sum_{i}^N F_c(\big<v\cdot h\big>)V(l_i,v,\alpha) \frac{\big<v\cdot h\big>}{\big<n\cdot h\big>}\big<n\cdot l\big> \end{aligned}

最后,我们将可以离线计算得等式提取出来:

DFG_1(\alpha,\big<n\cdot v\big>)= \frac{4}{N}\sum_{i}^N (1-F_c(\big<v\cdot h\big>))V(l_i,v,\alpha) \frac{\big<v\cdot h\big>}{\big<n\cdot h\big>}\big<n\cdot l\big> \tag{92}
DFG_2(\alpha,\big<n\cdot v\big>)=\frac{4}{N}\sum_{i}^N F_c(\big<v\cdot h\big>)V(l_i,v,\alpha) \frac{\big<v\cdot h\big>}{\big<n\cdot h\big>}\big<n\cdot l\big>
I(f(n,v,\alpha)) = f_0DFG_1(\alpha,\big<n\cdot v\big>)+f_{90}DFG_2(\alpha,\big<n\cdot v\big>)

注意DFG_1DFG_2只依赖于n\cdot v,即积分关于法线对称,当积分是,可以选用任意观察方向。

最终总结下公式:

L_{out}(n,v,\alpha,f_0,f_{90}) \approx[f_0DFG_1(n \cdot v,\alpha)+f_{90}DFG_2(n \cdot v,\alpha)]\times LD(n,\alpha)
DFG_1(\alpha,\big<n \cdot v\big>)= \frac{4}{N}\sum_{i}^N (1-F_c(\big<v \cdot h\big>))V(l_i,v,\alpha) \frac{\big<v \cdot h\big>}{\big<n \cdot h\big>}\big<n \cdot l\big>
DFG_2(\alpha,\big<n \cdot v\big>)=\frac{4}{N}\sum_{i}^N F_c(\big<v \cdot h\big>)V(l_i,v,\alpha) \frac{\big<v \cdot h\big>}{\big<n \cdot h\big>}\big<n \cdot l\big>
LD(n,\alpha)\equiv \frac {\sum_{i}^N V(h_i,\alpha)L_{\perp}(l_i)\big<n \cdot l\big>} {\sum_{i}^N \big<n \cdot l\big>}

DFG_1DFG_2可视化

DFG_1DFG_2可以预计算存储在(n \cdot v,\alpha)坐标系的2D纹理中,双线性采样。纹理样例如下(Y轴\alpha,X轴cos\theta):

DFG_1DFG_2的范围在[0,1]内,然而8位图的精度不够,可惜并不是所有移动设备支持16位或浮点纹理,采样器数量也有限制,如果不考虑纹理形式,用近似分析去模拟结果可能更好。纹理格式可以选用OpenGL ES3.0的RG16F。

下面是应用近似分析的方法得到的可视化结果:

LD可视化

LD是个只取决于\alpha参数(与粗糙度有关)的卷积,可以存储在mip-map立方体贴图中。lod层级与\alpha的关系:

\alpha=perceptualRoughness^2\\ lod_{\alpha}=\alpha^{\frac{1}{2}}=perceptualRoughness

例子:

间接高光和间接漫反射可视化

例子如下,移除了直接光:

IBL计算实现

vec3 ibl(vec3 n, vec3 v, vec3 diffuseColor, vec3 f0, vec3 f90,
        float perceptualRoughness) {
    vec3 r = reflect(n);
    vec3 Ld = textureCube(irradianceEnvMap, r) * diffuseColor;
    float lod = computeLODFromRoughness(perceptualRoughness);
    vec3 Lld = textureCube(prefilteredEnvMap, r, lod);
    vec2 Ldfg = textureLod(dfgLut, vec2(dot(n, v), perceptualRoughness), 0.0).xy;
    vec3 Lr =  (f0 * Ldfg.x + f90 * Ldfg.y) * Lld;
    return Ld + Lr;
}

我们可以用球谐替代光照贴图,用DFG查找表的近似分析。

vec3 irradianceSH(vec3 n) {
    // uniform vec3 sphericalHarmonics[9]
    // We can use only the first 2 bands for better performance
    return
          sphericalHarmonics[0]
        + sphericalHarmonics[1] * (n.y)
        + sphericalHarmonics[2] * (n.z)
        + sphericalHarmonics[3] * (n.x)
        + sphericalHarmonics[4] * (n.y * n.x)
        + sphericalHarmonics[5] * (n.y * n.z)
        + sphericalHarmonics[6] * (3.0 * n.z * n.z - 1.0)
        + sphericalHarmonics[7] * (n.z * n.x)
        + sphericalHarmonics[8] * (n.x * n.x - n.y * n.y);
}

// NOTE: this is the DFG LUT implementation of the function above
vec2 prefilteredDFG_LUT(float coord, float NoV) {
    // coord = sqrt(roughness), which is the mapping used by the
    // IBL prefiltering code when computing the mipmaps
    return textureLod(dfgLut, vec2(NoV, coord), 0.0).rg;
}

vec3 evaluateSpecularIBL(vec3 r, float perceptualRoughness) {
    // This assumes a 256x256 cubemap, with 9 mip levels
    float lod = 8.0 * perceptualRoughness;
    // decodeEnvironmentMap() either decodes RGBM or is a no-op if the
    // cubemap is stored in a float texture
    return decodeEnvironmentMap(textureCubeLodEXT(environmentMap, r, lod));
}

vec3 evaluateIBL(vec3 n, vec3 v, vec3 diffuseColor, vec3 f0, vec3 f90, float perceptualRoughness) {
    float NoV = max(dot(n, v), 0.0);
    vec3 r = reflect(-v, n);

    // Specular indirect
    vec3 indirectSpecular = evaluateSpecularIBL(r, perceptualRoughness);
    vec2 env = prefilteredDFG_LUT(perceptualRoughness, NoV);
    vec3 specularColor = f0 * env.x + f90 * env.y;

    // Diffuse indirect
    // We multiply by the Lambertian BRDF to compute radiance from irradiance
    // With the Disney BRDF we would have to remove the Fresnel term that
    // depends on NoL (it would be rolled into the SH). The Lambertian BRDF
    // can be baked directly in the SH to save a multiplication here
    vec3 indirectDiffuse = max(irradianceSH(n), 0.0) * Fd_Lambert();

    // Indirect contribution
    return diffuseColor * indirectDiffuse + indirectSpecular * specularColor;
}

多次散射预积分

之前讲过应对能量损失的做法。能量补偿项的缩放值取决于r,定义如下:

r=\int_{\Omega}D(l,v)V(l,v)\big<n\cdot l\big>\partial l \tag{93}

也可以利用重要性采样:

r=\frac{4}{N}\sum_i^NV(l_i,v,\alpha)\frac{\big<v\cdot h\big>}{\big<n\cdot h\big>\big<n\cdot l\big>}\tag{94}

这和DFG一样,只是没有菲涅尔项。如果进一步假设f_{90}=1,我们可以重构式子:

L_{out}(n,v,\alpha,f_0) \approx[(1-f_0)DFG_1^{multiscatter}(n\cdot v,\alpha) +f_{0}DFG_2^{}multiscatter(n\cdot v,\alpha)]\times LD(n,\alpha)
DFG_1^{multiscatter}(\alpha,\big<n\cdot v\big>)= \frac{4}{N}\sum_{i}^N F_c(\big<v\cdot h\big>)V(l_i,v,\alpha) \frac{\big<v\cdot h\big>}{\big<n\cdot h\big>}\big<n\cdot l\big>
DFG_2^{multiscatter}(\alpha,\big<n\cdot v\big>)=\frac{4}{N}\sum_{i}^N V(l_i,v,\alpha) \frac{\big<v\cdot h\big>}{\big<n\cdot h\big>}\big<n\cdot l\big>
LD(n,\alpha)\equiv \frac {\sum_{i}^N V(l_i,n,\alpha)L_{\perp}(l_i)\big<n\cdot l\big>} {\sum_{i}^NV(l_i,n,\alpha) \big<n\cdot l\big>}

实现里只需简要修改下:

vec2 dfg = textureLod(dfgLut, vec2(dot(n, v), perceptualRoughness), 0.0).xy;
// (1 - f0) * dfg.x + f0 * dfg.y
vec3 specularColor = mix(dfg.xxx, dfg.yyy, f0);

总结

重要性采样和预滤波IBL的区别:

假定v=n所产生的反射错误(下),缺少拉伸的反射效果:

在0.0625粗糙度下的立方体贴图LOD,刚好在两个级别间,会有交叉渐变:

在0.125粗糙度下的立方体贴图LOD,刚好在级别1,由于v=n,在极端角度下会很明显:

\alpha=0时,由于纹理缩放操作出现的波纹条样错误:

清漆

在采样IBL时,清漆层会作为一个额外的高光项参与计算,因为我们不能合理的沿半球积分,这一高光项要朝向观察方向。下面是相关实现:

// clearCoat_NoV == shading_NoV if the clear coat layer doesn't have its own normal map
float Fc = F_Schlick(0.04, 1.0, clearCoat_NoV) * clearCoat;
// base layer attenuation for energy compensation
iblDiffuse  *= 1.0 - Fc;
iblSpecular *= sq(1.0 - Fc);
iblSpecular += specularIBL(r, clearCoatPerceptualRoughness) * Fc;

各向异性

使用了“弯曲反射矢量”技术来对各向异性光照进行模拟。如下图:

实现如下:

vec3 anisotropicTangent = cross(bitangent, v);
vec3 anisotropicNormal = cross(anisotropicTangent, bitangent);
vec3 bentNormal = normalize(mix(n, anisotropicNormal, anisotropy));
vec3 r = reflect(-v, bentNormal);

这一技术可以使用负数anisotropy,当值为负,高光方向会出现在副切线方向而不是切线方向。

vec3 anisotropicDirection = anisotropy >= 0.0 ? bitangent : tangent;
vec3 anisotropicTangent = cross(anisotropicDirection, v);
vec3 anisotropicNormal = cross(anisotropicTangent, anisotropicDirection);
vec3 bentNormal = normalize(mix(n, anisotropicNormal, anisotropy));
vec3 r = reflect(-v, bentNormal);

正负值示例:

布料

布料材质模型的IBL更加复杂,难点在于它使用了不同的NDF。之前提到过的DFG项不能用于布料BRDF,由于我们设计的布料BRDF里并不包含菲涅尔项,我们可以生成一个简单的DG项,如下图:

实现如下:

float diffuse = Fd_Lambert() * ambientOcclusion;
#if defined(SHADING_MODEL_CLOTH)
#if defined(MATERIAL_HAS_SUBSURFACE_COLOR)
diffuse *= saturate((NoV + 0.5) / 2.25);
#endif
#endif

vec3 indirectDiffuse = irradianceIBL(n) * diffuse;
#if defined(SHADING_MODEL_CLOTH) && defined(MATERIAL_HAS_SUBSURFACE_COLOR)
indirectDiffuse *= saturate(subsurfaceColor + NoV);
#endif

vec3 ibl = diffuseColor * indirectDiffuse + indirectSpecular * specularColor;

透明和半透明光照

透明

我们得先理解材质不透明度是如何应用的,观察一扇窗,可以看见漫反射是透明的,而且高光反射越量,窗户看起来更透明。如下图:

为了恰当的实现不透明度,我们会使用预乘alpha。给定不透明度\alpha_{opacity}和漫反射颜色\sigma(线性,没有预乘),我们可以计算得到一个片段的有效不透明度:

\begin{aligned} color&=\sigma * \alpha_{opacity}\\ opacity&=\alpha_{opacity} \end{aligned}

更确切地描述下就是,源色的RGB定义了像素所发出光的量,alpha定义了该像素会遮挡多少光的量。使用以下混合函数:

\begin{aligned} Blend_{src}&=1\\ Blend_{dst}&=1-src_{\alpha} \end{aligned}

GLSL实现:

// baseColor has already been premultiplied
vec4 shadeSurface(vec4 baseColor) {
    float alpha = baseColor.a;

    vec3 diffuseColor = evaluateDiffuseLighting();
    vec3 specularColor = evaluateSpecularLighting();    

    return vec4(diffuseColor + specularColor, alpha);
}

半透明

半透明材质可以分为两类:

体半透明常用于粒子系统光照,例如云或烟。表面半透明用于模拟穿透散射的材质,例如蜡,大理石,皮肤等。

遮蔽

遮蔽是一个重要的变灰因数,可以在任意缩放级别创建阴影效果:

寒霜引擎中,漫反射微型遮蔽会预烘培进漫反射贴图中,高光微型遮蔽会预烘培进反射贴图中。在我们的系统中,微型遮蔽会烘培进基础颜色贴图中。

中缩放环境遮蔽会预烘培进环境遮蔽贴图中,作为一个材质参数存在。

大缩放环境遮蔽会使用基于屏幕空间的技术进行计算,如SSAO,HBAO等。

注意,为了避免过黑,推荐使用min(AO_{medium},AO_{large})

漫反射遮蔽

介绍一个构建环境遮蔽的方法。L_a是环境照度函数,利用球谐编码,接着定义一个可见性函数V,如果方向l上不存在遮挡,则V(l)=1,否则为0。

渲染等式的环境项如下:

L(l,v)=\int_{\Omega}f(l,v)L_a(l)V(l)\big<n\cdot l\big>dl \tag{95}

我们可以将可见性项分离出来近似:

L(l,v)=\Big(\pi\int_{\Omega}f(l,v)L_a(l)dl\Big) \Big(\frac{1}{\pi}\int_{\Omega}V(l)\big<n\cdot l\big>dl\Big) \tag{96}

只有在平行光L_a是常量,f是兰伯特项时,近似才会很准确。如果两个函数在大多数球面上相对平滑,那么近似结果还是可信的。

上述近似的左边一项是IBL的预计算漫反射部分,右边一项是范围在0-1之间的标量,由它可得漫反射环境遮蔽项:

AO=1-\frac{1}{\pi}\int_{\Omega}V(l)\big<n\cdot l\big>dl \tag{97}

由于我们使用的是预计算漫反射项,我们无法实时精确计算,为了解决这一问题,我们会通过应用环境遮蔽因数来部分地重构入射光。

实际操作中,烘培环境遮蔽数据会存储在灰度图中,通常是低精度格式。需要注意的是我们的材质模型里的环境遮蔽属性是用来创建宏观级别的漫反射环境遮蔽效果的。这样的近似结果在物理上并不准确,但也是折中考虑了质量和性能。

下图是有无漫反射环境遮蔽效果的对比图,左无右有:

GLSL实现:

// diffuse indirect
vec3 indirectDiffuse = max(irradianceSH(n), 0.0) * Fd_Lambert();
// ambient occlusion
indirectDiffuse *= texture2D(aoMap, outUV).r;

注意环境遮蔽只应用于非直接光照。

高光遮蔽

高光微型遮蔽可以由f_0推导得到,这一点基于现实物体反射率一定高于2%的准则,因此0-2%的值可以看作是预烘焙高光遮蔽,用来平滑地消减菲涅尔项。

float f90 = clamp(dot(f0, 50.0 * 0.33), 0.0, 1.0);
// cheap luminance approximation
float f90 = clamp(50.0 * f0.g, 0.0, 1.0);

之前提到的对于环境遮蔽地推导假定在兰伯特表面,并且只对非直接漫反射光照合理,表面可及性的信息缺失在重构非直接高光光照时会有重大缺陷,通常表现为漏光。

下面提出一个由漫反射遮蔽推到高光遮蔽的经验方法,结果并不基于物理,但效果不错。这一方法的目的是返回粗糙表面未经修改的漫反射遮蔽项,对于光滑表面,消除法向入射光的遮蔽影响,增加极端入射余角处的遮蔽影响。

float computeSpecularAO(float NoV, float ao, float roughness) {
    return clamp(pow(NoV + ao, exp2(-16.0 * roughness - 1.0)) - 1.0 + ao, 0.0, 1.0);
}

// specular indirect
vec3 indirectSpecular = evaluateSpecularIBL(r, perceptualRoughness);
// ambient occlusion
float ao = texture2D(aoMap, outUV).r;
indirectSpecular *= computeSpecularAO(NoV, ao, roughness);

注意高光遮蔽因数只应用于间接光照。

水平高光遮蔽

当为使用法线贴图的表面计算高光IBL部分时,可以用一个指向表面的反射矢量来决定其存在。如果反射矢量直接用于着色计算,那些不该被照亮的地方也会被当作照亮,这也是产生漏光的可能性,可以用以下方法消除影响。

核心概念是屏蔽来自表面后方的光,即反射矢量和法线的点积为负,那么反射矢量会指向表面。实现如下:

// specular indirect
vec3 indirectSpecular = evaluateSpecularIBL(r, perceptualRoughness);

// horizon occlusion with falloff, should be computed for direct specular too
float horizon = min(1.0 + dot(r, n), 1.0);
indirectSpecular *= horizon * horizon;

法线贴图

法线贴图作用有二:用低面数模型表示高面数模型,增加表面细节。

有无法线贴图对比,左无右有:

法线贴图:

如果我们想叠加法线贴图的话,如下图,增加裂痕的细节贴图:

对于存储在切线空间的法线贴图,很显然普通的线性或叠加混合并不可用。

重定向法线贴图

给定基础法线贴图和细节法线贴图,将细节图的法线基矢量旋转到基础贴图的法线基矢量。

实现如下:

vec3 t = texture(baseMap,   uv).xyz * vec3( 2.0,  2.0, 2.0) + vec3(-1.0, -1.0,  0.0);
vec3 u = texture(detailMap, uv).xyz * vec3(-2.0, -2.0, 2.0) + vec3( 1.0,  1.0, -1.0);
vec3 r = normalize(t * dot(t, u) - u * t.z);
return r;

这里假定了法线信息未被压缩,范围在[0..1]。

这一操作开销较大,通常会离线使用。下面是一个对比,左混法线,右应用:

UDN混合

UDN混合是偏导混合技术的一个变体,性能不错,可实时进行,不过会损失一些细节。

vec3 t = texture(baseMap,   uv).xyz * 2.0 - 1.0;
vec3 u = texture(detailMap, uv).xyz * 2.0 - 1.0;
vec3 r = normalize(t.xy + u.xy, t.z);
return r;

体积效果

指数高度雾

图像处理管线

OETF步骤中,会应用目标颜色空间的光电转换函数。

基于物理的摄像机

图像变换的第一步是使用基于物理的摄像机恰当地曝光场景的亮度。

曝光设定

因为我们在光照管线里所使用的是光度单位,到达摄像机的光能量用亮度L表示,单位是cd/m^2。摄像机传感器接收的光亮度范围在10^{-5}-10^9,范围太大,我们得先映射到合适范围,再去处理。

这一重映射由摄像机在一定时间内曝光传感器完成。为了最大程度使用传感器的范围限制,场景中灯光范围会位于中间灰度附近,介于黑白之间。曝光有三个相关设置:

光圈:标记为N,单位为光圈系数(f-stops)f,这一设置控制了摄像机系统的光圈开闭程度。光圈系数即镜头焦距入瞳直径之间的比率,值越高,表明光圈越小,值越小,表明光圈越大。光圈还可以控制景深。

快门速度:标记为t,单位为秒s,控制了光圈维持打开状态的时间(也控制了传感器快门的时机)。另外,还可以控制动态模糊。

感光度:标记为S,单位为ISO,用于量化进入传感器的光。此外,还可以控制噪声程度。

曝光度

我们用一个曝光度的变量来涵盖上述三种设置,标记为EV。

EV表示为一个以2为底的对数缩放,1EV的差分称之为1级。一正级(+1EV)表示两倍亮度,一负级(-1EV)表示一半的亮度。

EV=log_2(\frac{N^2}{t})\tag{98}

这只是关于光圈和快门速度的函数,不包含感光度。按照惯例,曝光度是针对100感光度定义的,即EV_{100}

对任意给定感光度,曝光度等式:

EV_S=EV_{100}+log_2(\frac{S}{100})\tag{99}

由此可得EV_{100}:

EV_{100}=EV_S-log_2(\frac{S}{100})=log_2(\frac{N^2}{t})-log_2(\frac{S}{100})\tag{100}

曝光度和亮度

摄像机可以测量场景的平均亮度,并且转换为曝光度以实现自动曝光。EV和亮度的关系:

EV=log_2(\frac{L\times S}{K})\tag{101}

K是一个反光计量计常量,不同设备不一。常用的有两个,12.5和14。这里选用K=12.5

针对EV_{100}得:

EV=log_2(L\frac{100}{12.5})\tag{102}

在引擎中,我们可以先测定某一帧的平均亮度,从而实现自动曝光。测量亮度最简单的方法是降采样亮度缓冲到1像素,接着读取剩余的值。这一方法并不稳定,会被极端值所影响,因此大多数游戏里会使用亮度直方图来移除极端值。

由给定曝光度EV可以计算得到亮度:

L=2^{EV_{100}} \times \frac{12.5}{100}=2^{EV_{100}-3}\tag{103}

曝光度和照度

EV和照度关系:

EV=log_2(\frac{E\times S}{C})\tag{104}

C是入射光测量计常量,不同设备不一,针对平面传感器,常用值为250,针对半球传感器,常用的值有320和340。

针对EV_{100},得:

EV=log_2(E\frac{100}{C})\tag{105}

对于给定的曝光度,对于平面传感器,可得照度:

E=2^{EV_{100}}\times 2.5\tag{106}

对于半球传感器:

E=2^{EV_{100}}\times 3.4\tag{107}

曝光补偿

对摄像机来说,曝光度是那三个设定的组合,但对用户来说,曝光度通常被认为是光强度,因此摄像机会提供一个曝光补偿项来辅助调节。曝光补偿EC作为偏移量存在:

EV_{100}^{'}=EV_{100}-EC \tag{108}

这里的补偿作为负数,是因为我们把光圈系数作为EC的单位来调整最终的曝光。EV提升即调小镜头光圈,高曝光度会让图片变暗。

曝光

为了将场景亮度转换为标准亮度,我们必须使用光度曝光,用于调整到达摄像机传感器的光通量。光度曝光H,单位为勒克斯秒:

H=\frac{q\cdot t}{N^2}L \tag{109}

L是场景亮度,t是快门速度,N是光圈,q是镜头光晕衰减(一般情况下,q=0.65)。上述定义未考虑传感器感光度,我们可以使用三种方法与感光度关联:基于饱和度的速度,基于噪声的速度,标准输出感光度。

这里介绍基于饱和度的速度,得到最大可能曝光H_{sat}

H_{sat}=\frac{78}{S_{sat}}\tag{110}

组合109和110可以计算得到最大亮度L_{max}:

L_{max}=\frac{N^2}{q \cdot t}\frac{78}{S}\tag{111}

由此可以计算得到标准化入射光亮度:

L'=L\frac{1}{L_{max}}\tag{112}

L_{max}可以用公式98简化,S=100,q=0.65

\begin{aligned} L_{max}&=\frac{N^2}{t}\frac{78}{q \cdot S}\\ L_{max}&=2^{EV_{100}}\frac{78}{q \cdot S}\\ L_{max}&=2^{EV_{100}}\times 1.2 \end{aligned}

曝光应用:

// Computes the camera's EV100 from exposure settings
// aperture in f-stops
// shutterSpeed in seconds
// sensitivity in ISO
float exposureSettings(float aperture, float shutterSpeed, float sensitivity) {
    return log2((aperture * aperture) / shutterSpeed * 100.0 / sensitivity);
}

// Computes the exposure normalization factor from
// the camera's EV100
float exposure(float ev100) {
    return 1.0 / (pow(2.0, ev100) * 1.2);
}

float ev100 = exposureSettings(aperture, shutterSpeed, sensitivity);
float exposure = exposure(ev100);

vec4 color = evaluateLighting();
color.rgb *= exposure;

自动曝光

上述的曝光需要用户手动调节,但在拍摄动态场景时,亮度变化剧烈,手动调整不现实。我们已经介绍过如何由亮度计算曝光度,应用在测光表上也是可以的。

测量场景亮度有两种方法:

注意上述两种方法会在与颜色相乘后寻找平均亮度,这样子做会不太对,我们可以保留一个亮度缓冲,存储相乘操作前的每个像素的亮度。开销大。

这两个方法也限制了平均测光过程,因为每个像素的权重是始终一致的。摄像机提供三种测光模式:

点测光:只有图像中央的小圈范围内的像素会参与最终的曝光,小圈尺寸通常是全尺寸的1%-5%。

中心权重测光:屏幕中心的亮度值权重更大。

多重域或矩阵测光:每种设备不一。这一模式会优先曝光场景中最重要的部分,操作是网格切分图像并分类每个单元,还可以进一步与一直的数据进行比对来实现更恰当的曝光操作。

点测光

每个亮度值的权重w的计算公式如下:

w(x,y)= \left\{ \begin{aligned} &1 \ \ |p_{x,y}-s_{x,y}|\le s_r&\\ &0 \ \ |p_{x,y}-s_{x,y}|> s_r& \end{aligned} \right. \tag{113}

其中p是像素的位置,s是点的中心,s_r是点的半径。

中心权重测光

w(x,y)=smooth(|p_{x,y}-c|\times \frac{2}{width})\tag{114}

其中c是时间中心,smooth是平滑函数,如smoothstep()

适应

为了平滑测量结果,我们可以用以下等式:

L_{avg}=L_{avg}+(L-L_{avg})\times(1-e^{-\Delta t\cdot \tau})\tag{115}

其中\Delta t是链子前一帧的间隔时间,\tau是控制适应率的常量。

发光

因为曝光度缩放在视觉上几乎是连续的,所以常常用来当作一个灯光单位。这也意味着我们可以让用户使用曝光补偿作为单位来调节光强度,光强度因此与曝光设置相关,我们应尽量避免这种使用方法,不过可以应用在发光表面创建发光效果。

根据发光颜色c和当前曝光度EV_{100},我们可以计算发光值得亮度:

EV_{bloom}=EV_{100}+EC \tag{116}
L_{bloom}=c\times 2^{EV_{bloom}-3}

实现如下:

vec4 surfaceShading() {
    vec4 color = evaluateLights();
    // rgb = color, w = exposure compensation
    vec4 emissive = getEmissive();
    color.rgb += emissive.rgb * pow(2.0, ev100 + emissive.w - 3.0);
    color.rgb *= exposure;
    return color;
}

灯光路径

引擎所使用的灯光路径,或者说渲染方法,可能会存在严重的性能问题,在光源数量上也有很大的限制。传统的渲染方法有两种,前向和延迟。

我们使用的渲染方法要遵循以下限制:

另外,我们想要支持:

许多现代渲染引擎使用的延迟渲染可以很轻易地支持成百上千的光源,但很耗带宽。如果使用我们默认地PBR材质模型,我们地G-buffer会每像素使用160-192位空间,也就是所需带宽量很高。

前向渲染向来不适合处理大量光源。通用的步骤是对于每个可见光渲染一次场景,然后叠加结果,还有一种方法是为每个物体绑定一些光源,不过对于大型物体来说不太现实(比如建筑,道路等)。

分块着色(tiled)可以同时用于前向和延迟渲染,核心在于将屏幕按网格进行切分,对于每个单元块,找到影响其中像素的一系列光源,这项技术可以减少重绘次数,也可以为大型物体进行着色计算。这项技术会面临深度不连续的问题,额外增加许多不必要的工作。

集群前向渲染(clustered):

场景分块:

集群前向渲染

我们进一步介绍另一个方法,集群着色,它是分块渲染的眼神,在三维层面分块。集群切分在视图空间进行,将视锥切分为3D网格。

视锥首先沿深度轴切分:

接着进行屏幕切分,体素化视锥。我们把每一簇称之为锥素(froxel,估计是frustum+voxel,锥体+体素)。

在渲染一帧前,每个光源会绑定到范围交叉的任何锥素上,这一绑定过程的结果就是对于每个锥素会有一个光源列表。在渲染过程中,我们会计算一个片段所属的锥素ID,由此可得影响该片段的灯光列表。

深度并不是线性切分的,而是指数。在场景中,靠近近平面的像素多于远平面,这样可以在更为重要的区域提升光源绑定的精度,

指数切分中每次深度切分使用的世界空间距离量:

简单的指数体素化并不够,上面的图表很清晰地表明了沿切分方向的空间分布,但对于接近近平面的分布出错了。我么可以测试一下更小范围内的分布:

由上可知,前半切分都很接近近平面。16次切分,前5米就分了8次,对于点光源和聚光来说(球体和锥体),精度太高,过于浪费了,得不偿失。

解决方案是根据场景和远近平面来调整第一次切分的锥素大小,这样的话,我们可以更好地分布余下的锥素。下图展示了在0.1m-0.5m之间使用一个特定锥素后的结果:

实现笔记

灯光绑定可以在GPU或CPU上进行。

GPU灯光绑定

这一实现需要计算着色器的支持。光源存储在SSBO中,并且传输到进行绑定操作的计算着色器中。

锥体体素化由第一个计算着色器执行,可以只进行一次,灯光绑定由其它计算着色器逐帧执行。

计算着色器的线程模型非常适配这项任务,我们可以很简单地调用和锥素数量一致的工作组(XYZ工作组数量对应锥素网格分辨率),每个工作区会作为一个线程,遍历所有灯光来进行绑定。

CPU灯光绑定

在不支持计算着色器的设备上进行,算法和GPU实现不同,引擎会将每个光源进行类似光栅化的操作,作为锥素存在。例如,给定点光源的位置和半径,很容易计算它所交叉的锥素列表。

着色

逐锥素灯光列表可以以SSBO或纹理形式传递到片元着色器中。

从深度到锥素

给定近平面n,远平面f,最大深度切分数量m,[0..1]范围的线性深度值z,可以为给定位置计算簇的索引:

zToCluster(z,n,f,m)=floor(max(log2(z \frac{m}{-log2( \frac{n}{f})}+m,0))\tag{117}

这一公式有分辨率相关问题,我们可以引入sn,一个特殊的近平面值,定义了初次切分的锥素的扩展(第一次切分的锥素范围[n..sn],其余[sn..f])。

zToCluster(z,sn,f,m)=floor(max(log2(z \frac{m-1}{-log2(\frac{sn}{f})}+m,0))\tag{118}

下面的等式可以根据gl.FragCoord.z计算线性深度值:

linearZ(z)=\frac{n}{f+z(n-f)}\tag{119}

可以预先计算c0c1来简化公式:

c1=\frac{f}{n}\tag{120}
c0=1-c1
linearZ(z)=\frac{1}{z\cdot c0 + c1}

这一简化很重要,因为我们会将线性Z值传入log2中,恰好可以避免除0的问题。

计算给定片段的锥素索引的实现如下:

define MAX_LIGHT_COUNT 16 // max number of lights per froxel

uniform uvec4 froxels; // res x, res y, count y, count y
uniform vec4 zParams;  // c0, c1, index scale, index bias

uint getDepthSlice() {
    return uint(max(0.0, log2(zParams.x * gl_FragCoord.z + zParams.y) *
            zParams.z + zParams.w));
}

uint getFroxelOffset(uint depthSlice) {
    uvec2 froxelCoord = uvec2(gl_FragCoord.xy) / froxels.xy;
    froxelCoord.y = (froxels.w - 1u) - froxelCoord.y;

    uint index = froxelCoord.x + froxelCoord.y * froxels.z +
            depthSlice * froxels.z * froxels.w;
    return index * MAX_FROXEL_LIGHT_COUNT;
}

uint slice = getDepthSlice();
uint offset = getFroxelOffset(slice);

// Compute lighting...

下面是预先计算的一些变量:

froxels[0] = TILE_RESOLUTION_IN_PX;
froxels[1] = TILE_RESOLUTION_IN_PX;
froxels[2] = numberOfTilesInX;
froxels[3] = numberOfTilesInY;

zParams[0] = 1.0f - Z_FAR / Z_NEAR;
zParams[1] = Z_FAR / Z_NEAR;
zParams[2] = (MAX_DEPTH_SLICES - 1) / log2(Z_SPECIAL_NEAR / Z_FAR);
zParams[3] = MAX_DEPTH_SLICES;

从锥素到深度

给定锥素索引i,特殊近平面sn,远平面f和最大切分数量m,可以计算得到给定锥素的最小深度:

clusterToZ(i\ge 1,sn,f,m)2^{(i-m)\frac{-log2(\frac{sn}{f})}{m-1}}\tag{121}

i=0,z为 0。结果范围在[0..1]间,乘以f得到世界空间距离。

在计算着色器中应该用exps代替pow,商可以预先计算好,作为常量传入。

坐标系统

世界空间

Filamen使用Y朝上的右手坐标系(红+X,绿+Y,蓝+Z)

摄像机空间

Filament的摄像机看向它的局部-Z轴。

立方体贴图坐标系统

Filament的立方体贴图遵循OpenGL的惯例:

镜像

为了简化反射渲染,IBL立方体贴图X轴镜像存储。

等距投影环境贴图

为了将等距投影贴图转换为水平或垂直排布的立方体贴图,我们需要将+Z面的位置放在源方格环境贴图的中心。

环境贴图的世界空间定向和天空盒

挡在Filament里定义一个天空盒或IBL时,立方体填土的-Z面应面对世界空间的+Z轴,对于镜像过的,则是-Z对-Z。

上一篇 下一篇

猜你喜欢

热点阅读