GraphicsLab Project之基于物理的着色系统（Physical based shading）-直接光照

引言

近些年来，基于物理的光照着色系统（Physical based shading）越来越流行。主流的渲染器，游戏引擎都支持了这种着色方式。相比于以前的Phong和Blinn-Phong的光照着色模型，这种着色方式更加真实，也更加容易适应不同场景。接下来的几篇文章，我将从实现的角度分析如何设计实现一个这样的光照系统。主要仿照Unreal4的实现方式来进行实现。

基于物理的着色

那么，在我们进行实际的编码之前，先大致的了解下什么是基于物理的着色。所谓基于物理的着色，就是根据物理世界中光与材质的交互方式，近似的提出一种方案来模拟这种效果。读者需要注意，这种方式是“基于物理的”，并不是“物理精确的”。它不是真实的再现现实世界中光与物质的交互，而是把握它们交互的主要特性，通过数学近似的方法来模拟。一个基于物理的着色，需要有如下三个基本的条件：

1.基于微表面理论建模的反射模型 
2.能量守恒 
3.基于物理的BRDF 
（除此之外，还有什么相互性等等，不过本篇文章主要关注以上三个条件）

渲染方程

我们接下来的实现，将主要围绕一个方程来展开。这个方程总结了想要渲染出一个场景的一般性步骤，这个方程称为渲染方程（针对于实时渲染来说，又可以称之为反射方程）： 
Lo=∫Ωf(pi,wi,wo)Li(pi,wi)n⋅widwi
Lo=∫Ωf(pi,wi,wo)Li(pi,wi)n⋅widwi

下面简要的介绍下公式中每一个部分的意思： 

LoLo：表示的是经过着色之后，从wowo方向观察点pipi时的颜色
f(pi,wi,wo)f(pi,wi,wo)：表示的是点pipi上，从wowo方向反射出去的光照与从wiwi方向的入射的光照的比值，该函数称为BRDF（Bidirectional Reflection Distribution Function） 

Li(pi,wi)Li(pi,wi)：从wiwi方向入射到pipi点上的光照 

n⋅win⋅wi：表示的是光照的反射强度与光照方向与接受光照表面法线方向之间角度的关系，即Lambert Law 

∫Ω…dwi∫Ω…dwi：表示的是在点pipi法线方向上的半球，所有入射光线方向的积分 

 

光与材质的交互

当一束光线照射到材质表面的时候，会分成两个主要的部分，反射的光线和折射的光线。 

 

对于折射的部分，由于物质本身的特性，一部分会被吸收，一部分会在物质本身的粒子之间来回碰撞，直至又从物体表面反射回去。 

 

从物体内部再次反射出去的光线，是不规则的。想要精确的模拟这种情况，需要用到次表面散射的技术(Subsurface scattering)。引入这种技术，需要非常复杂的计算，消耗过多，所以业界一般使用漫反射(Diffuse)来进行模拟。 

所以，当处理一束光照射到物质表面的时候，需要分成了两个不同的反射部分进行处理，分别是漫反射和镜面反射。所以渲染方程就变成了如下的形式： 
Lo=∫Ω(fd+fs)Li(n⋅wi)dwi
Lo=∫Ω(fd+fs)Li(n⋅wi)dwi

也就是将原先的BRDF函数分为了两个不同的部分：漫反射部分和镜面反射部分。 

直接光照处理理论

由于本篇文章处理直接光照，也就是只处理解析光源，比如平行光，点光源这样的光源，所以渲染方程能够再一次的简化成为如下的形式： 
Lo=(fd+fs)Li(n⋅wi)
Lo=(fd+fs)Li(n⋅wi)

这种简化是因为对于单一解析光源来说，物体表面上同一个点只会从一个方向接受到来自光源的光照，所以不需要进行积分计算。对于多个解析光源，我们也只要简单的分别对光源进行上述公式的计算，然后叠加结果即可。

漫反射部分

前面说过，反射方程主要分为了两个部分，其中一个为漫反射，它模拟的是从物体内部反射出去的光线效果。业界有很多模拟漫反射的方法，如经典的Lambertian Reflection模型和Oren-Nayar模型等。在Unreal4中采用的是Lambertian Reflection模型，这种模型十分简单，很容易计算出来，消耗非常小，并且也能够得到很好的结果。我的实现也将采用这个模型。 

Lambertian Reflection模型实际上假设物体内部反射出来的光线，是均匀的分布在半球上的，所以这种模型又被称之为Perfect Diffuse Reflection，如图所示： 

 

而该模型的BRDF如下： 
fd=cπ
fd=cπ

其中c为该物体材质的Albedo属性，和Phong模型中的Diffuse Color/Texture类似。而除以ππ的原因是为归一化BRDF，从而使得BRDF符合能量守恒的条件，不会导致反射出去的光线比入射的光线还要多的情况。 

采用Lambertian Reflection漫反射模型是不是非常的简单，它和以前的Phong十分的类似，唯一不同之处就在于多了一个归一化的操作。而这就是以前经验式的光照模型，与基于物理的光照模型之间的区别。 

至此，我们的渲染方程变成了如下的模样： 
Lo=(kd⋅cπ+fs)Li(n⋅wi)
Lo=(kd⋅cπ+fs)Li(n⋅wi)

（kd为漫反射部分所占的入射光照的比例，后文会详细解释）

镜面反射部分

上面讲解完毕了漫反射部分之后，接下来了解下镜面反射部分。从前面的描述我们了解到，基于物理着色系统的漫反射光照和普通的Phong式的光照差别不大，那么镜面反射了？实际上，基于物理着色系统主要就是提高了镜面反射部分的细节，能够通过镜面反射模拟出更加真实的物体表面属性。 

在详细了解具体的镜面反射光照模型之前，我们先来了解下该模型提出的一个基本假设–微表面模型。 

微表面模型

以前的光照模型，都是假设被照射到的一个平面都是光滑的表面，但是实际上表面并不总是那么光滑的。现实世界中，物体的表面从粗糙到光滑都存在，而粗糙的表面具有更加多的细节。微表面模型就是描述了一个表面粗糙程度的模型。 

我们可以假设，被光照到的一块表面区域，实际上是由很多细小的光滑表面组成，如下图所示： 

 

表面越粗糙，那么微表面分布就越加的崎岖不平；表面越光滑，微表面分布就越加的平整，如下图所示： 

 

所以，我们这里要讲解的镜面反射光照模型就是根据此微表面模型来构想的。根据不同的输入，模拟不同程度的微表面分布情况。 

Cook-Torrance镜面反射光照模型

目前业界经常使用的镜面光照模型就是Cook-Torrance镜面反射光照模型。实际上，Cook-Torrance镜面反射光照模型不是单一的模型，而是一种光照模型框架，它描述了一些基本的构成部分，模型可以自己根据需要来选择合适的公式。 

下面就是Cook-Torrance的整体公式： 
fs=D⋅F⋅G4⋅(n⋅l)(n⋅v)
fs=D⋅F⋅G4⋅(n⋅l)(n⋅v)

D项

D项是一个Normal Distribution Function(NDF)函数，这个函数描述了微表面法线的分布情况。当越多的微表面法线与宏观表面法线相近的时候，高光区域就越亮，亮斑就越小，反之就越暗，亮斑逐渐扩散。这个函数返回的是一个标量。 

F项

F项模拟的是Fresnel效果。从前面的光与材质交互一节我们知道，光照射到表面的时候，分为了两个不同的部分，反射（镜面反射部分）和折射（漫反射部分）。那么这两个部分是如何分布的了，每个部分应该占据多少比例了？这个就是由F项来确定。大体上，当我们观察到的视角向量vv与表面的法相向量nn之间夹角越小的时候，折射部分就越多，反射就越小；反之就折射越少，反射越多。这个效果你可以想象，当你垂直的观察水面的时候，你能够很清晰的看到水底的东西，而很难看到水面反射的东西；当你接近与水面平行的时候，你能够很清晰的看到水面反射的物体，而很难看到水底的物体。所以，F项描述了镜面反射部分所占的比例。同时需要知道，不同波长的光线，照射到同一个物体表面的时候，反射和折射的比例并不总是相同的，所以F项实际上是一个RGB三元分量。 

G项

G项，即Geometry Function，又称之为Shadow and Mask项。当一束光线照射到表面的时候，由于表面的微表面特征，有一些光线会被其他表面遮挡住（Shadow）；同时，当照射到微表面上之后，一些被反射出去的光线也会被其他光线给遮挡住（Mask）。所以就使用这个项来描述那些遮挡住光线的微表面占所有微表面的比值情况，所以这个G项也是一个标量。 

 

函数选取

接下来，我们就依次的讲解Cook-Torrance三个函数D,F,G函数的选择情况。这里还是以Unreal4来举例说明。 

D项

对于D项，我们选取的函数是GGX，该函数如下所示： 
D(n,h,α)=α2π((n⋅h)2(α2−1)+1)2
D(n,h,α)=α2π((n⋅h)2(α2−1)+1)2

其中： 
nn：表面的法向向量 
hh：half向量，该向量可以通过l+v∣∣l+v∣∣l+v∣∣l+v∣∣计算得到 
αα：表面的粗糙值，通过表面的Roughness属性计算得来 

F项

前面我们知道，F项描述的是镜面反射部分的比例。精确的模拟Fresnel效果的方程十分复杂，业界一般使用如下的近似函数来模拟： 
F(n,v,F0)=F0+(1−F0)(1−(n⋅v))5
F(n,v,F0)=F0+(1−F0)(1−(n⋅v))5

其中： 
nn：表面的法向向量 
vv：观察向量 
F0F0：表面的基本反射属性 

这里有一个需要注意的地方，既然F项描述的是镜面反射部分的比例，那么漫反射部分的比例自然就是1 - F了。大部分情况下，的确是这样。但是对于金属材质的物体来说，光照射到物体表面，一部分被反射，一部分被折射，而被折射的部分很快的被材质给吸收了，所以并没有从内部出去。 

所以，对于金属材质来说，只有镜面反射，没有漫反射。也就是说，这里的F项，需要区别对待金属材质和非金属材质。 

为了能够让我们的算法能够更容易的处理这种情况，我们引入了一个属性Metallic，用来描述表面的金属度。同时，最好能够让算法统一的处理金属材质和非金属材质。 

下图，是不同材质的F0F0值: 

 

观察上图，我们得到金属材质的F0F0值普遍高于0.5，而非金属材质的F0F0值都小于0.17。也就说，这个非金属材质的F0F0值集中在一个非常小的范围里面，所以我们就使用一个平均值来近似的表示非金属材质的F0F0属性，这个值为0.04。也就是说，对于非金属材质来说，我们就直接认定F0=0.04F0=0.04。由于金属表面没有漫反射部分，所以它的外观主要由镜面反射部分来表现，因此我们保留它的F0F0值。 

所以，当我们在使用这里的F项函数的时候，需要计算一下它需要使用的F0F0值。 

G项

对于G项，我们使用的是Schlick-GGX函数，这个函数和前面的GGX很类似，如下所示： 
GSchlickGGX(n,v,α)=n⋅v(n⋅v)(1−k)+k
GSchlickGGX(n,v,α)=n⋅v(n⋅v)(1−k)+k

其中： 
k=α2
k=α2

nn：表面法相向量 
vv：观察或者光照向量 
αα：表面的粗糙程度属性，由表面的Roughness属性计算得来。 

同时由于G项描述了Shadow和Mask两中情况，所以最终的G项为： 
G(n,v,l,α)=GSchlickGGX(n,v,α)⋅GSchlickGGX(n,l,α)
G(n,v,l,α)=GSchlickGGX(n,v,α)⋅GSchlickGGX(n,l,α)

新渲染方程

所以，现在渲染方程变成了： 
Lo=(kd⋅cπ+D⋅F⋅G4⋅(n⋅l)(n⋅v))Li(n⋅wi)
Lo=(kd⋅cπ+D⋅F⋅G4⋅(n⋅l)(n⋅v))Li(n⋅wi)
直接光照处理实现

理论到这里为止了，接下来将以代码的形式给出最终的实现。 

算法的输入

该光照算法，需要如下几个输入参数： 

 

其中， 
Albedo：对于非金属材质，描述的就是基本颜色属性，相当于以前的Diffuse Texture，不过去除了Diffuse Texture中的阴影和高亮，只有纯粹的颜色属性；对于金属材质来说，由于前面说过的，金属材质没有漫反射部分，而金属材质需要一个F0F0参数，所以就用它来描述金属材质的F0F0属性。在前面那张表中也可以看出，金属材质的F0F0参数的各个分量是不一样的，所以刚好可以使用三分量来表示。 

Normal：就是表面的法线信息 

Metallic：表面的金属度，用来表示表面是一个金属还是非金属。从这句话的描述来看，它应该是一个二元的值，要们是金属，要么不是金属。理论上的确是这样的，但是在实际设计的过程中，把这个值设计为0到1范围的一个值。这么做的原因是为了模拟那些金属表面覆盖了一层其他非金属材质的情况，比如金属材质上的锈迹，由于它既不完全是金属，也不完全是非金属，所以使用一个介于他们之间的值来表示更加合理，同时能够模拟两种材质叠加的效果。 

Roughness：物体表面的粗糙属性，从0到1变化，0表示最平滑，1表示最粗糙。 

Fresnel效果计算

由于F项表述了漫反射和镜面反射的比例关系，所以最先计算出这个值来，如下代码计算：

vec3 F0 = mix(vec3(0.04, 0.04, 0.04), albedo, metalic);
vec3 F = calc_frenel(h, v, F0);

第一句就是计算出合适的F0F0值出来，如果是非金属材质，那么就使用0.04来表示；如果是金属材质，就是用albedo里面存放的F0F0属性；如果是介于之间的，就使用金属度进行插值。 

第二句就是简单的带入Fresnel公式，不过有一个注意点。细心的读者可能发现了，这里我带入的公式里面使用的是h(half-way)向量，而不是表面的法线n。这是因为，前面我们在讨论Fresnel公式的时候，都是脱离微表面模型，直接从一个光滑表面的角度去看这个公式的。对于一个光滑表面来说，的确就是使用法线n来计算Fresnel系数。但是在这里，我们使用了微表面模型，而微表面本身并不是一个完全光滑的表面，它是有很多具有一定分布形式的微型光滑表面所构成。对于微型的光滑表面，由于它们都是完美的光滑表面，所以只有那些入射光线经过反射然后进入观察者眼中的表面才能够对宏观表面的外在表现作出贡献，而这些能够做出贡献的微型表面就满足了如下的条件： 
reflect(l,n)=v
reflect(l,n)=v

也就是说，只有具有满足如上条件法线n的微型表面才是有效的。由此我们得出了这个n： 
n=l+v∣∣l+v∣∣=h
n=l+v∣∣l+v∣∣=h

(关于half-way向量的说明，大家可以参考Blin-Phong模型) 

calc_fresnel代码如下：
vec3 calc_frenel(vec3 n, vec3 v, vec3 F0) {
    float ndotv = max(dot(n, v), 0.0);
    return F0 + (vec3(1.0, 1.0, 1.0) - F0) * pow(1.0 - ndotv, 5.0);
}
 

单单看Fresnel效果的计算的话，效果如下所示，你可以根据此图来看看你的计算是否正确： 

 

漫反射比例

由于F项表示的就是镜面反射的比例，那么剩下的就是漫反射的比例，再考虑金属材质没有漫反射部分的情况，得到如下的代码：

vec3 T = vec3(1.0, 1.0, 1.0) - F;
vec3 kD = T * (1.0 - metalic);
 

关于公式中αα的计算方式

前面的D项和G项，都有关于αα值的使用，而在注解中，我提到他们的计算方式需要通过Roughness属性得来，而不是直接使用Roughness。这样做的原因是，我们希望Roughness本身能够保持一个线性的变化，而各函数中αα并不线性变化的，所以对于这两个函数来说，我们需要将Roughness重新映射到αα上去。 

D项

采用Disney的映射方式： 
α=Roughness∗Roughness
α=Roughness∗Roughness

最终得到如下的代码：
float calc_NDF_GGX(vec3 n, vec3 h, float roughness) {
    float a = roughness * roughness;
    float a2 = a * a;
    float ndoth = max(dot(n, h), 0.0);
    float ndoth2 = ndoth * ndoth;
    float t = ndoth2 * (a2 - 1.0) + 1.0;
    float t2 = t * t;
    return a2 / (PI * t2);
}
 

单单看这个函数的结果，根据不同的Roghness值，如下图所示： 

 

G项

同样采用Disney的映射方式: 
α=Roughnessremapping∗Roughnessremapping
α=Roughnessremapping∗Roughnessremapping

同时需要注意，上述使用的是RoughnessremappingRoughnessremapping，而不是直接使用Roughness，所以： 

Roughnessremapping=Roughness+12
Roughnessremapping=Roughness+12

所以最终，G项中的： 
k=(Roughess+1)28
k=(Roughess+1)28
所以关于G项的代码如下所示：

float calc_Geometry_GGX(float costheta, float roughness) {
    float a = roughness;
    float r = a + 1.0;
    float r2 = r * r;
    float k = r2 / 8.0;

    float t = costheta * (1.0 - k) + k;

    return costheta / t;
}

float calc_Geometry_Smith(vec3 n, vec3 v, vec3 l, float roughness) {
    float ndotv = max(dot(n, v), 0.0);
    float ndotl = max(dot(n, l), 0.0);
    float ggx1 = calc_Geometry_GGX(ndotv, roughness);
    float ggx2 = calc_Geometry_GGX(ndotl, roughness);
    return ggx1 * ggx2;
}

 

单独看这个函数的效果如下图所示： 

 

线性空间

PBS强烈的依赖于你的计算空间中各个值是否是线性的，所以需要HDR和Gamma修正。我这里使用了最简单的Tone mapping计算，相关代码如下：

// base tone mapping
color = color / (color + vec3(1.0, 1.0, 1.0));

// gamma correction
color = pow(color, vec3(1.0 / 2.2, 1.0 / 2.2, 1.0 / 2.2));

#version 330

in vec3 vs_Vertex;
in vec3 vs_Normal;

uniform vec3 glb_Albedo;
uniform float glb_Roughness;
uniform float glb_Metalic;
uniform vec3 glb_EyePos;

out vec3 oColor;

const float PI = 3.1415927;

vec3 calc_frenel(vec3 n, vec3 v, vec3 F0) {
    float ndotv = max(dot(n, v), 0.0);
    return F0 + (vec3(1.0, 1.0, 1.0) - F0) * pow(1.0 - ndotv, 5.0);
}

float calc_NDF_GGX(vec3 n, vec3 h, float roughness) {
    float a = roughness * roughness;
    float a2 = a * a;
    float ndoth = max(dot(n, h), 0.0);
    float ndoth2 = ndoth * ndoth;
    float t = ndoth2 * (a2 - 1.0) + 1.0;
    float t2 = t * t;
    return a2 / (PI * t2);
}

float calc_Geometry_GGX(float costheta, float roughness) {
    float a = roughness;
    float r = a + 1.0;
    float r2 = r * r;
    float k = r2 / 8.0;

    float t = costheta * (1.0 - k) + k;

    return costheta / t;
}

float calc_Geometry_Smith(vec3 n, vec3 v, vec3 l, float roughness) {
    float ndotv = max(dot(n, v), 0.0);
    float ndotl = max(dot(n, l), 0.0);
    float ggx1 = calc_Geometry_GGX(ndotv, roughness);
    float ggx2 = calc_Geometry_GGX(ndotl, roughness);
    return ggx1 * ggx2;
}

vec3 calc_lighting_direct(vec3 n, vec3 v, vec3 l, vec3 h, vec3 albedo, float roughness, float metalic, vec3 light) {
    vec3 F0 = mix(vec3(0.04, 0.04, 0.04), albedo, metalic);
    vec3 F = calc_frenel(h, v, F0);

    vec3 T = vec3(1.0, 1.0, 1.0) - F;
    vec3 kD = T * (1.0 - metalic);

    float D = calc_NDF_GGX(n, h, roughness);

    float G = calc_Geometry_Smith(n, v, l, roughness);

    vec3 Diffuse = kD * albedo * vec3(1.0 / PI, 1.0 / PI, 1.0 / PI);
    float t = 4.0 * max(dot(n, v), 0.0) * max(dot(n, l), 0.0) + 0.001;
    vec3 Specular = D * F * G * vec3(1.0 / t, 1.0 / t, 1.0 / t);

    float ndotl = max(dot(n, l), 0.0);
    return (Diffuse + Specular) * light * vec3(ndotl, ndotl, ndotl);
}

void main() {
    vec3 view = glb_EyePos - vs_Vertex;
    view = normalize(view);

    vec3 lightPos = vec3(0.0, 0.0, 200.0);
    vec3 light = lightPos - vs_Vertex;
    light = normalize(light);

    vec3 half = normalize(view + light);

    vec3 color = calc_lighting_direct(vs_Normal, view, light, half, glb_Albedo, glb_Roughness, glb_Metalic, vec3(2.5, 2.5, 2.5));

    // base tone mapping
    color = color / (color + vec3(1.0, 1.0, 1.0));

    // gamma correction
    color = pow(color, vec3(1.0 / 2.2, 1.0 / 2.2, 1.0 / 2.2));

    oColor = color;
}

总结

本文的实例代码放在： 
https://github.com/idovelemon/GraphicsLabtory/tree/master/glbcodebase/graphicslab/glb_pbs 
这里，感兴趣的读者可以自行下载了解。

参考文献
[1] LearnOpenGL 
[2] s2010_physically_based_shading_hoffman_notes 
[3] s2013_pbs_epic_notes 
[4] Ray Tracing From Ground Up 
[5] Physical based Rendering: From Theory to Implementation

————————————————
版权声明：本文为CSDN博主「i_dovelemon」的原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/i_dovelemon/article/details/78945950