引言
对真实感的追求,让渲染技术在 “更详细地模拟光与物体的交互” 的方向上发展。PBR(基于物理的渲染) 是所有这些尝试通过对光的物理模拟来实现照片级效果的技术的统称。
目前,模拟光的最佳模型是一个被称为 the rendering equation(渲染方程) 的方程。渲染方程尝试描述的是:在给出场景中与特定点的所有交互的入射光的情况下,如何得到这一点发出的光。我们在随后会了解其细节并介绍正确的术语。需要注意的是——我们并不会尝试解决完整的渲染方程,而是使用以下简化的版本:
L
o
(
p
,
ω
o
)
=
∫
Ω
f
r
(
p
,
ω
i
,
ω
o
)
L
i
(
p
,
ω
i
)
(
n
⋅
ω
i
)
d
ω
i
L_o(p,\omega_o) = \int_\Omega f_r(p,\omega_i,\omega_o)L_i(p,\omega_i)(\bf{n} \cdot \omega_i) \it{d} \omega_i
Lo(p,ωo)=∫Ωfr(p,ωi,ωo)Li(p,ωi)(n⋅ωi)dωi
要理解这个方程,我们首先需要了解“光”的原理,然后,我们需要了解一些常用术语。为了让你大致了解该公式的含义,简单来说,我们可以说该公式描述了:在给定了一个像素所受到的所有光以及如何混合它们的函数的情况下,这个像素的颜色是什么。
物理术语
为了能正确理解渲染方程,我们需要了解一些物理量的含义。而这些量中最重要的是 辐射率(Radiance)(公式中的
L
L
L)。
它是一个不好理解的东西,因为它是一些其他物理量的结合。因此,在正式定义它之前,我们来看一些其他的物理量。
辐射通量(Radiant Flux):它是光源发出的总能量的量度,以瓦特表示(瓦特=焦耳/秒)。我们将用
Φ
\Phi
Φ表示通量。
任何光源都会发射能量,发射能量的数额是波长的函数。
立体角(Solid angle):它衡量了当观察者从一个点观察一个物体时,物体看起来有多大。为此,我们将物体的轮廓投影到以我们观察的点为中心的单位球体的表面上,我们得到的形状的面积就是“立体角”。在下图中,你可以看到立体角ω是淡蓝色多边形在单位球面上的投影:
辐射强度(Radiant Intensity):它指每个立体角上的通量。假设你有一个向所有方向发射光的光源,那么有多少光(通量)实际上是朝特定方向发射的?辐射强度就是用来回答这个问题,它是在一个方向上通过一个给定立体角的通量。描述它的公式是
I
=
d
Φ
d
ω
I=\frac{d\Phi}{d\omega}
I=dωdΦ,其
Φ
\Phi
Φ表示辐射通量而ω表示立体角。
辐射率(Radiance):终于,我们可以给出辐射率的公式了:
L
=
d
Φ
2
d
ω
d
A
c
o
s
θ
L=\frac{d\Phi^2}{d\omega dAcos\theta}
L=dωdAcosθdΦ2
其中:
Φ
\Phi
Φ是辐射通量,
A
A
A是光所影响的这个“小平面”的面积,ω是光所穿过的立体角,
θ
\theta
θ是这个“小平面”的法线与光线方向的夹角(因此,
A
c
o
s
θ
Acos\theta
Acosθ就是这个“小平面”在光线方向上的投影面积)。
我们喜欢这个公式是因为我们感兴趣的所有物理学部分它都包含了。我们可以用辐射率来描述“一条光线”:它是通过一个无限小的立体角,到达一个无限小的区域的通量——这样就描述了光线的行为。因此,当我们讨论辐射率时,我们讨论的是“从某个方向到某个区域的一定量的光”。
当我们对一个点进行渲染时,进入该点的所有入射光我们都感兴趣,即所有照射到以该点本身为中心的半球的辐射率的总和,这个总和的名字叫做 辐照度(Irradiance)。辐照度 和 辐射率 是我们主要关注的物理量,我们将研究它们以实现PBR。
渲染方程
现在我们可以回到渲染方程并尝试完全理解它了:
L
o
(
p
,
ω
o
)
=
∫
Ω
f
r
(
p
,
ω
i
,
ω
o
)
L
i
(
p
,
ω
i
)
(
n
⋅
ω
i
)
d
ω
i
L_o(p,\omega_o) = \int_\Omega f_r(p,\omega_i,\omega_o)L_i(p,\omega_i)(\bf{n} \cdot \omega_i) \it{d} \omega_i
Lo(p,ωo)=∫Ωfr(p,ωi,ωo)Li(p,ωi)(n⋅ωi)dωi
我们已经明白了 L L L就是辐射率,而它是关于世界上的一个点和某一方向加上立体角的函数(从现在开始我们将始终使用无限小的立体角,因此可以简单地将其视为方向向量)。这个方程描述的是:如何计算从一个点 p p p朝着 ω o \omega_o ωo发出的辐射率 L o ( p , ω o ) L_o(p,\omega_o) Lo(p,ωo),而这个结果就是我们为屏幕上的像素着色所需的全部内容。
为了计算它,我们需要像素所在位置的表面的法线 n \bf{n} n,和来自于场景的辐照度。为了获得辐照度,我们对所有入射的辐射率 L i ( p , ω i ) L_i(p,\omega_i) Li(p,ωi)相加,因此方程前面有积分符号 ∫ Ω \int_\Omega ∫Ω。请注意积分的域 Ω \Omega Ω是一个以我们正在计算的点为中心的半球,它朝向表面,也就是说你可以沿着表面法线的方向从该点移动到半球的顶部。
点积 n ⋅ ω i \bf{n} \cdot \omega_i n⋅ωi 考虑了光线的入射角。如果光线垂直于表面,它将更多地集中在光照区域;而如果角度较小,它将传播到更大的区域,最终传播得太多以至于实际上不可见。
现在我们可以明白,这个方程其实就是为了计算射出的辐射率,而你需要提供所有的入射辐射率,它们还会加权每条入射光线与表面法线的夹角的余弦。
接下来我们还需要介绍的就是 f r ( p , ω i , ω o ) f_r(p,\omega_i,\omega_o) fr(p,ωi,ωo) ,它就是BRDF。这个函数将“位置”、“入射光线”、“出射光线”作为函数的输入,而函数会输出 “入射光有多少比例会作用于最终出射的辐射率”。对于一个绝对的镜面反射比如镜子,它的BRDF会是这样:函数在那条精确与入射光线有相同反射角的出射光线上返回1,其余所有的出射光线函数都会返回0。有重要的一点要注意:一个基于物理的BRDF必须遵守能量守恒定律,即 ∀ ∫ Ω f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω o ≤ 1 \forall \int_\Omega f_r(p,\omega_i,\omega_o)(\bf{n} \cdot \omega_{\it{i}}) \it{d} \omega_o \le 1 ∀∫Ωfr(p,ωi,ωo)(n⋅ωi)dωo≤1 ,这意味着所有出射光的总量不能超过入射光的量。
翻译为代码
现在我们已经掌握了所有需要的知识,那么我们该如何将它们应用上,来实际写出渲染代码呢?我们有两个主要问题:
- 首先,我们该如何表示出场景中所有的辐射率?
- 其次,我们该如何快速求解积分以便能够在实时引擎中使用它?
第一个问题很简单——使用环境贴图(Environment Maps)。这个贴图(CubeMap,尽管使用球形的贴图更准确)上的每个像素都编码了一个点周围一个特定方向上的入射辐射率。如果我们想象CubeMap上的每个像素都是一个发光源而RGB颜色是通量,那么我们就可以将
L
(
p
,
ω
)
L(p,\omega)
L(p,ω) 近似为从CubeMap上采样的值(其中
p
p
p是CubeMap的中心)。也就是说
L
(
p
,
ω
)
≈
t
e
x
C
U
B
E
(
c
u
b
e
m
a
p
,
ω
)
L(p,\omega)\approx texCUBE(cubemap,\omega)
L(p,ω)≈texCUBE(cubemap,ω)。显然,让场景中的每个点都有一个立方体贴图会消耗太多内存!因此我们牺牲了一些质量:在场景中创建一定数量的立方体贴图,并让每个点都选择最接近的那个。为了减少误差,我们可以使用立方体贴图的世界位置来校正用于采样的向量,使其更准确。这样,我们最终得到了一种估算辐射率的方法,即:
L
(
p
,
ω
)
≈
t
e
x
C
U
B
E
(
c
u
b
e
m
a
p
,
ω
p
)
L(p,\omega)\approx texCUBE(cubemap,\omega_p)
L(p,ω)≈texCUBE(cubemap,ωp)
其中 ω p \omega_p ωp 是经由点和立方体贴图在世界中的位置校正后的采样向量。
对于第二个问题——如何求解积分,就有点棘手了。因为在某些情况下,我们是无法足够快地求解的。但如果在某些情况下,比如BRDF返回的结果碰巧只取决于入射的辐射率而和出射的方向无关,或者更简单些——不取决于任何东西而是个常数,那么此时我们就可以做一些漂亮的优化了。所以让我们看看在 Lambert’s BRDF 中的情况,在这里它就是一个常数,即所有入射光作用于出射光的比例都是个固定值。
Lambert
在 Lambert’s BRDF 中,
f
r
(
p
,
ω
i
,
ω
o
)
=
c
π
f_r(p,\omega_i,\omega_o) = \frac{c}{\pi}
fr(p,ωi,ωo)=πc,其中 c 是表面的颜色。如果我们将其代入渲染方程,就将得到:
L
o
(
p
,
ω
o
)
=
∫
Ω
c
π
L
i
(
p
,
ω
i
)
(
n
⋅
ω
i
)
d
ω
i
=
c
π
∫
Ω
L
i
(
p
,
ω
i
)
(
n
⋅
ω
i
)
d
ω
i
\begin{aligned} L_o(p,\omega_o) & = \int_\Omega\frac{c}{\pi}L_i(p,\omega_i)(\bf{n} \cdot \omega_i) \it{d} \omega_i \\ & = \frac{c}{\pi}\int_\Omega L_i(p,\omega_i)(\bf{n} \cdot \omega_i) \it{d} \omega_i \\ \end{aligned}
Lo(p,ωo)=∫ΩπcLi(p,ωi)(n⋅ωi)dωi=πc∫ΩLi(p,ωi)(n⋅ωi)dωi
现在,积分只取决于
ω
i
\omega_i
ωi 没有别的,这意味着我们可以预先计算它(例如用蒙特卡罗积分)并将结果存储到另一个CubeMap中。这些值将按照
ω
o
\omega_o
ωo 的方向来存储,意思是给定一个方向,我们可以对CubeMap进行采样并获得该方向上的光。这将整个渲染方程简化成了一次对一个预计算的CubeMap的采样:
L
(
p
,
ω
o
)
=
t
e
x
C
U
B
E
(
l
a
m
b
e
r
t
C
u
b
e
m
a
p
,
ω
o
p
)
L(p,\omega_o) = texCUBE(lambertCubemap,\omega_{op})
L(p,ωo)=texCUBE(lambertCubemap,ωop)
其中 ω o p \omega_{op} ωop 是经由点和立方体贴图在世界中的位置校正后的出射方向。
现在,我们已经有所有元素了,终于可以写一个shader了!先看看结果吧:
对于一个只是采样单个贴图的shader来说,效果还不错吧?请注意随着环境的变化整个光源是如何变化的。渲染所使用的CubeMap并非由原CubeMap经过 convolved(译注:可能是类似下采样的操作) 得到的,而是看起来模糊得多,如下图所示:
(左图是环境的入射辐射率图,右图是Lambert版渲染方程下的出射辐射率图)
下面让我们展示shader代码。请注意,为了简单起见,我没有使用蒙特卡罗积分,而是简单地将积分离散化了。理论上当样本数目无限时二者没有区别,但现实中我这样做会比蒙特卡罗积分方法产生更多的带状瑕疵。在我的测试里这已经足够好了,因为我已经将CubeMap的分辨率降低到每面 32x32了。但是如果你想尝试使用它,请记住这一问题。
我们所需的第一个shader是用来生成模糊的环境贴图的,它常被称为convolved envmap(卷积化的环境贴图),因为它是环境贴图与卷积核 ( n ⋅ ω i ) (\bf{n} \cdot \omega_i) (n⋅ωi)的卷积结果。
由于在shader中我们将在球面坐标上积分,所以我们将以球面坐标的形式改写公式:
L
o
(
p
,
θ
o
,
ϕ
o
)
=
c
π
∫
Φ
∫
Θ
L
i
(
p
,
θ
i
,
ϕ
i
)
c
o
s
(
θ
i
)
s
i
n
(
θ
i
)
d
θ
i
d
ϕ
i
L_o(p,\theta_o,\phi_o)=\frac{c}{\pi} \int_\Phi \int_\Theta L_i(p,\theta_i,\phi_i)cos(\theta_i)sin(\theta_i)d\theta_id\phi_i
Lo(p,θo,ϕo)=πc∫Φ∫ΘLi(p,θi,ϕi)cos(θi)sin(θi)dθidϕi
译者补充:
首先,这个球面坐标可以看作是放在模型表面一点 p p p 上,而球面坐标的天顶方向就是表面的法线方向。
然后是关于球面坐标的符号:
θ \theta θ 表示极角, ϕ \phi ϕ 表示方位角。
θ o \theta_o θo 就是出射的极角, ϕ o \phi_o ϕo 就是出射的方位角,由于 Lambert’s BRDF 是无关出射光线的方向的,所以可以看到此处的计算和这二者无关。
这里主要关注的是入射光线的极角 θ i \theta_i θi 和方位角 ϕ i \phi_i ϕi,公式主要是对其求了积分 ∫ Φ ∫ Θ L i ( p , θ i , ϕ i ) d θ i d ϕ i \int_\Phi \int_\Theta L_i(p,\theta_i,\phi_i)d\theta_id\phi_i ∫Φ∫ΘLi(p,θi,ϕi)dθidϕi,即对所有方向上的辐射率求和。
而在这个积分公式中还有两项: c o s ( θ i ) cos(\theta_i) cos(θi) 与 s i n ( θ i ) sin(\theta_i) sin(θi)。
其中 c o s ( θ i ) cos(\theta_i) cos(θi) 自不必说,入射光线方向越接近于表面切线方向时,表面得到的光照自然越小。问题是为什么还会乘 s i n ( θ i ) sin(\theta_i) sin(θi) ?这个问题就是下一段作者详细讨论的。
你或许已经注意到了,在公式中有个额外的 s i n ( θ i ) sin(\theta_i) sin(θi) 。这是由于:“积分”实际上可以看作是在微小的均匀的步长上计算的,当我们使用立体角来积分时,一切都是OK的,因为立体角本来就是在积分域上均匀分布的。然而当我们换成球面坐标时,我们会在 θ \theta θ 为0时得到相对更多的采样,而在 θ \theta θ 为 π 2 \frac{\pi}{2} 2π 时得到相对更少的采样。你可以在建模软件中创建一个球然后在线框模式观察它,然后你就明白我的意思了。
译者的理解:
θ \theta θ 和 ϕ \phi ϕ 的变化产生了立体角,即 d ω d\omega dω 由 d θ d\theta dθ 和 d ϕ d\phi dϕ 共同作用产生。
然而问题是,当 θ \theta θ 在不同的值上时,立体角的变化并不一致。
如图所示,在 θ \theta θ 不同的位置,同样的 d θ d\theta dθ 与 d θ d\theta dθ 会产生不同的 d ω d\omega dω 。
因此乘算 s i n ( θ i ) sin(\theta_i) sin(θi) 是对立体角分布上的补偿,即 d ω i = s i n ( θ ) d θ d ϕ d\omega_i=sin(\theta)d\theta d\phi dωi=sin(θ)dθdϕ
这个公式中的双重积分是通过为每个积分使用蒙特卡洛估算法来求解的。
而我们将使用下面离散的方程:
L
o
(
p
,
θ
o
,
ϕ
o
)
=
c
π
2
π
N
1
π
2
N
2
∑
N
1
∑
N
2
L
i
(
p
,
θ
i
,
ϕ
i
)
c
o
s
(
θ
i
)
s
i
n
(
θ
i
)
=
π
c
N
1
N
2
∑
N
1
∑
N
2
L
i
(
p
,
θ
i
,
ϕ
i
)
c
o
s
(
θ
i
)
s
i
n
(
θ
i
)
\begin{aligned} L_o(p,\theta_o,\phi_o) & = \frac{c}{\pi} \frac{2\pi}{N_1} \frac{\pi}{2N_2} \sum^{N_1} \sum^{N_2} L_i(p,\theta_i,\phi_i)cos(\theta_i)sin(\theta_i)\\ & = \frac{\pi c}{N_1N_2} \sum^{N_1} \sum^{N_2} L_i(p,\theta_i,\phi_i)cos(\theta_i)sin(\theta_i)\\ \end{aligned}
Lo(p,θo,ϕo)=πcN12π2N2π∑N1∑N2Li(p,θi,ϕi)cos(θi)sin(θi)=N1N2πc∑N1∑N2Li(p,θi,ϕi)cos(θi)sin(θi)
(译者注:这个公式将原先的积分变成了离散的求和。其中 N 1 N_1 N1是方位角划分的步数, N 2 N_2 N2是极角划分的步数)
下面终于可以将其变为shader代码了:
...
float4 PixelShaderFunction(VertexShaderOutput input) : COLOR
{
//这个函数用来计算“特定表面法线方向”上的出射辐射率,结果将存在CubeMap的对应位置上
//cubeFace代表的是CubeMap中哪一个面
//InterpolatedPosition代表“特定表面法线方向”对应在这一个CubeMap的面上的位置,范围(-1~1)
//下面根据cubeFace来换算normal:“特定表面法线方向”在世界空间下的方向
float3 normal = normalize( float3(input.InterpolatedPosition.xy, 1) );
if(cubeFace==2)
normal = normalize( float3(input.InterpolatedPosition.x, 1, -input.InterpolatedPosition.y) );
else if(cubeFace==3)
normal = normalize( float3(input.InterpolatedPosition.x, -1, input.InterpolatedPosition.y) );
else if(cubeFace==0)
normal = normalize( float3( 1, input.InterpolatedPosition.y,-input.InterpolatedPosition.x) );
else if(cubeFace==1)
normal = normalize( float3( -1, input.InterpolatedPosition.y, input.InterpolatedPosition.x) );
else if(cubeFace==5)
normal = normalize( float3(-input.InterpolatedPosition.x, input.InterpolatedPosition.y, -1) );
//下面计算球面坐标的轴在世界空间下的方向,
//其中up是天顶方向,right是方位角基准方向
float3 up = float3(0,1,0);
float3 right = normalize(cross(up,normal));
up = cross(normal,right);
float3 sampledColour = float3(0,0,0); //逐渐累加的采样得到的颜色值,最后会根据采样次数平均化
float index = 0; //采样次数,用于最后的平均化
for(float phi = 0; phi < 6.283; phi += 0.025) //采样所有的方位角方向
{
for(float theta = 0; theta < 1.57; theta += 0.1)//采样顶半球所有的极角方向
{
//下面计算sampleVector:本次采样的方向在世界空间下的方向
float3 temp = cos(phi) * right + sin(phi) * up;
float3 sampleVector = cos(theta) * normal + sin(theta) * temp;
//从环境贴图中采样,并套入本篇上面讨论的公式
sampledColour += texCUBE( diffuseCubemap_Sampler, sampleVector ).rgb *
cos(theta) * sin(theta);
index ++;
}
}
//最终结果平均化:
return float4( PI * sampledColour / index), 1 );
}
...
我省略了顶点着色器、变量定义。
使用正常的环境贴图CubeMap作为输入,为 convolved cubemap 的每一个面都运行此shader,就得到了“辐照度图”。
然后,我们就可以使用“辐照度图”作为下一个着色器(模型着色器)的输入了:
...
float4 PixelShaderFunction(VertexShaderOutput input) : COLOR
{
float3 irradiance= texCUBE(irradianceCubemap_Sampler, input.SampleDir).rgb;
float3 diffuse = materialColour * irradiance;
return float4( diffuse , 1);
}
...
非常简单且计算迅速。
关于PBR的文章的第一部分到此结束