WIP：PBR基础知识整理

概述

本文主要用于个人的学习总结。系统地梳理PBR相关知识，并结合工作中使用引擎的实现来进一步加深理解效果和性能的取舍。

什么是PBR？

PBR全程physical based rendering，意即为基于物理的渲染，他通过模拟真实的物理光照来计算颜色。PBR主要通过三个维度体现：

• 基于物理的材质
• 基于物理的光照
• 基于物理适配的摄像机

PBR的前提

为了正确地计算出颜色（能量），需要保证几个前提：

线性工作流

线性空间渲染为光照计算提供了正确的数学运算。在线性空间中，能够还原现实世界方式的光与物质的交互的方式。所以颜色值的计算和颜色操作必须在线性空间中执行。而为了将渲染图像正确地呈现给观看者，需要将图像编码为伽马空间，所以基于物理的渲染会往往会涉及到线性空间和伽马空间之间的相互转换。

Tone Maping

也称色调复制（Tone Reproduction），在图形学中，表示以感知上令人信服的方式将HDR场景的强度值转换为显示强度的过程，也可以理解为将宽范围的光照级别拟合到屏幕有限色域内的过程。通常，由于通过HDR渲染出来的亮度值会超过显示器能够显示最大亮度，所以需要结合色调映射（Tone Mapping），将光照结果从HDR转换为显示器能够正常显示的LDR。

基于真实世界测量的材质参数（Real-World Measurement Based Substance Properties）

PBR的正统材质参数往往都基于真实世界测量。真实世界中的物质可分为三大类：绝缘体（Insulators），半导体（semiconductors）和导体（conductors）。在渲染和游戏领域，我们一般只对其中的两个感兴趣：导体（金属）和绝缘体（电解质，非金属）。菲涅尔反射率代表材质的镜面反射颜色与强度，是真实世界材质的核心测量数值。其中非金属具有非彩色的镜面反射颜色，而金属具有彩色的镜面反射颜色，即非金属的F0是一个float,金属的F0是一个float3。

能量守恒

出射光线的能量永远不能大于入射光线的能量。随着粗糙度的上升镜面反射区域的面积会增加，作为平衡，镜面反射区域的平均亮度则会下降。
这一类守恒可以通过构建一个简单的线性关系来简单地近似

float kS = CalculateSpecularComponent(); //镜面反射部分
float kD = 1.0 - ks;                //折射、漫反射部分

PBR的原理

要了解PRB的计算方法，首先要了解PBR的物理原理。

渲染方程

渲染方程(The Rendering Equation)作为渲染领域中的重要理论，其描述了光能在场景中的流动，是渲染中不可感知方面的最抽象的正式表示。根据光学的物理学原理，渲染方程在理论上给出了一个完美的结果，而各种各样的渲染技术，只是这个理想结果的一个近似。

反射方程 The Reflectance Equation

而在PBR中，只关心反射到我们眼睛中的颜色，因此取渲染方程的后半部分，称为反射方程。舍弃的Le主要代表了光源本身的发光，这个涉及到了全局光照的计算，在此不做赘述

辐射度量学和PBR Radiometry

那么通过反射方程到底在计算什么样的一个物理量？这里就需要引出辐射度量学。辐射度量学是一种用来度量电磁场辐射（包括可见光）的手段。有很多种辐射度量(radiometric quantities)可以用来测量曲面或者某个方向上的光。为了更好的计算PBR，首先我们要使用一套物理正确的计量标准，因此我们需要利用辐射度量学来衡量我们计算的能量。实际转为RGB编码后物理不正确，但在图形学中常说，一个普遍的调侃就是一个东西如果看起来是对的，那么就是对的。（狗头）下文将会探讨几个概念，方便读者理解PBR的度量，以及如何将其在计算机光栅化中进行模拟。

辐射通量

辐射通量表示的是一个光源所输出的能量，以瓦特为单位。光是由多种不同波长的能量所集合而成的，而每种波长则与一种特定的（可见的）颜色相关。因此一个光源所放射出来的能量可以被视作这个光源包含的所有各种波长的一个函数。

辐射通量将会计算这个由不同波长构成的函数的总面积。我们通常不直接使用波长的强度而是使用三原色编码，也就是RGB（或者按通常的称呼：光色）来作为辐射通量表示的简化

立体角 solid angle

立体角用表示，它可以为我们描述投射到单位球体上的一个截面的大小或者面积。投射到这个单位球体上的截面的面积就被称为立体角(Solid Angle)

辐射强度 radiant intensity

辐射强度(Radiant Intensity)表示的是在单位球面上，一个光源向每单位立体角所投送的辐射通量。举例来说，假设一个全向光源向所有方向均匀的辐射能量，辐射强度就能帮我们计算出它在一个单位面积（立体角）内的能量大小：

辐射率 radiance

一个拥有辐射强度的光源在单位面积，单位立体角上的辐射出的总能量。

辐射率是一个非常接近计算机中的光栅化的模型。我们可以把一个像素看作是一个无限小的平面，把立体角也看作无限小，那么我们可以将这个定义建模成，一束光线照射到一个像素上的强度（RGB值）。
辐射率也受到光纤与平面法线夹角的影响。

float cosTheta = dot(lightDir, N);

反射方程实际描述的就是某个无限小的立体角在某个点上的辐射率。他在辐射率中只考虑一个单位立体角的基础上扩展到了需要计算整个半球范围内，所有立体角上所有光源照射到这个平面的辐射强度。这个立体角可以视作入射方向向量。而我们实际看到的颜色是这个点上反射到我们眼睛的辐射率。因此我们最终计算的是一个反射方程。PBR的核心也就是对于这个反射方程的近似求解。他计算了半球上所有光线通过所有单位立体角照射到单位平面上经过反射被我们的眼睛所观察到的辐照度。

BRDF

有了前述的概念之后，参考反射方程的公式，我们知道实际需要计算的就是fr函数，有了fr函数我们就知道光线照射到一个平面上反射到我们的眼睛里具体的比例有多少。有了这个比例，再乘以光照的辐照度即可。 这个fr函数称作BxDF:

• BxDF一般而言是对BRDF、BTDF、BSDF、BSSRDF等几种双向分布函数的一个统一的表示。
• BSDF可以看做BRDF和BTDF更一般的形式，而且BSDF = BRDF + BTDF。
• BSSRDF和BRDF的不同之处在于，BSSRDF可以指定不同的光线入射位置和出射位置。
  而本文将会聚焦BRDF。BRDF也有许多种模型，业界一般漫反射部分采用lambert，镜面反射部分采用Cook-Torrance BRDF

那么BRDF就可以表示为：

其中,Kd和Ks是漫反射和镜面反射分别的比例。需要注意的是这里的漫反射包含了折射和次表面散射，次表面散射是否要纳入考量重点计算取决于观测的尺度。本文暂时不考虑,仅将其近似为漫反射。接下来我们就要分别对BRDF的漫反射和镜面反射部分进行求解。

Diffsue BRDF

采用lambert光照模型。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-VwClWakt-1677993306164)(https://note.youdao.com/yws/res/2802/WEBRESOURCEec36d30d000d891ba9f054b82ba68102)]
π用来对漫反射光进行标准化。这里的计算非常简单：

(Kd * albedo) * rcp(PI) * radiance * NdotL;

Specular BRDF

Microfacet Cook-Torrance BRDF

• D(h) : 法线分布函数 （Normal Distribution Function），描述微面元法线分布的概率，即正确朝向的法线的浓度。即具有正确朝向，能够将来自l的光反射到v的表面点的相对于表面面积的浓度。
• F(l,h) : 菲涅尔方程（Fresnel Equation），描述不同的表面角下表面所反射的光线所占的比率。
• G(l,v,h) : 几何函数（Geometry Function），描述微平面自成阴影的属性，即m = h的未被遮蔽的表面点的百分比。
• 分母 4(n·l)(n·v）：校正因子（correctionfactor），作为微观几何的局部空间和整个宏观表面的局部空间之间变换的微平面量的校正。

attention

• 对于分母中的点积，仅仅避免负值是不够的 ,也必须避免零值。通常通过在常规的clamp或绝对值操作之后添加非常小的正值来完成。
• Microfacet Cook-Torrance BRDF是实践中使用最广泛的模型，实际上也是人们可以想到的最简单的微平面模型。它仅对几何光学系统中的单层微表面上的单个散射进行建模，没有考虑多次散射，分层材质，以及衍射。Microfacet模型，实际上还有很长的路要走。

D(h) Specular D 法线分布函数

微表面理论

微平面理论是将物体表面建模成无数微观尺度上有随机朝向的理想镜面反射的小平面（microfacet）的理论。微观几何（microgeometry）的效果是在表面上的不同点处改变微平面的法线，从而改变反射和折射的光方向。出于着色的目的，通常会用统计方法处理这种微观几何现象，将表面视为具有微观结构法线的随机分布，并将宏观表面视为在每个点处多个方向上反射（和折射）光的集合。在微观尺度上，表面越粗糙，反射越模糊，表面越光滑，反射越集中。

常见模型

Beckmann[1963]
Blinn-Phong[1977]
GGX [2007] / Trowbridge-Reitz[1975]
Generalized-Trowbridge-Reitz(GTR) [2012]
Anisotropic Beckmann[2012]
Anisotropic GGX [2015]

业界： GGX (Trowbridge-Reitz)

float D_GGX_TR(vec3 N, vec3 H, float a)
{
    float a2     = a*a;
    float NdotH  = max(dot(N, H), 0.0);
    float NdotH2 = NdotH*NdotH;

    float nom    = a2;
    float denom  = (NdotH2 * (a2 - 1.0) + 1.0);
    denom        = PI * denom * denom;

    return nom / denom;
}

Specular F 菲涅尔项

基于F0建模的菲涅尔反射

菲涅尔效应表示观察看到的反射光线的量与视角相关的现象，且掠射角度（90度）下反射率最大。万物皆有菲涅尔效应。在宏观层面看到的实际上是微观层面微平面菲涅尔效应的平均值，即影响菲涅尔效应的关键参数在于每个微平面的法向量和入射光线的角度，而不是宏观平面的法向量和入射光线的角度。F0即0度角入射的菲涅尔反射率。任意角度的菲涅尔反射率可由F0和入射角度计算得出。

常见模型

Cook-Torrance [1982]
Schlick [1994]
Gotanta [2014]

业界: Schlick approximation

F0

F0是从0度观察时反射的辐照度比例,它通过折射系数（index of refraction IOR）计算得出。Schlick approximation仅仅对电解质或者非金属表面有定义。对于导体(Conductor)表面（金属），使用它们的折射指数计算基础折射率并不能得出正确的结果，这样我们就需要使用一种不同的菲涅尔方程来对导体表面进行计算。由于这样很不方便，所以我们预计算出平面对于法向入射的结果（，处于0度角，好像直接看向表面一样），然后基于相应观察角的Fresnel-Schlick近似对这个值进行插值，用这种方法来进行进一步的估算。这样我们就能对金属和非金属材质使用同一个公式了。

使用金属度插值计算F0:

float3 f0 = float3(0.04);
F0 = lerp(F0, surfaceColor.rgb, metalness);

Schlick approximation:

half3 fresnelSchlick(float cosTheta, half3 F0)
{
    return F0 + (1.0 - F0) * pow(1.0 - cosTheta, 5.0);
}

Specular G 几何遮挡函数

常见模型

Smith [1967]
Cook-Torrance [1982]
Neumann [1999]
Kelemen [2001]
Implicit [2013]

联合遮蔽阴影函数 （Smith Joint Masking-Shadowing Function）

四种形式：
分离遮蔽阴影型（Separable Masking and Shadowing）
高度相关掩蔽阴影型（Height-Correlated Masking and Shadowing）
方向相关掩蔽阴影型（Direction-Correlated Masking and Shadowing）
高度-方向相关掩蔽阴影型（Height-Direction-Correlated Masking and Shadowing）

分离遮蔽阴影（Separable Masking and Shadowing Function）

将G分为俩个独立的部分：光线方向和视角方向，两者使用相同的分布函数来描述。
结合法线分布函数和smith几何阴影函数，就有了：
Smith-GGX
Smith-Beckmann
Smith-Schlick
Schlick-Beckmann
Schlick-GGX

业界：schlick-GGX

float GeometrySchlickGGX(float NdotV, float k)
{
    float nom   = NdotV;
    float denom = NdotV * (1.0 - k) + k;

    return nom / denom;
}

float GeometrySmith(vec3 N, vec3 V, vec3 L, float k)
{
    float NdotV = max(dot(N, V), 0.0);
    float NdotL = max(dot(N, L), 0.0);
    float ggx1 = GeometrySchlickGGX(NdotV, k);
    float ggx2 = GeometrySchlickGGX(NdotL, k);

    return ggx1 * ggx2;
}

有了以上这些功能函数，那么完整的PBR直接光照计算是：

// cook-torrance brdf
float3 Lo = 0;
float NDF = DistributionGGX(N, H, roughness);       
float G   = GeometrySmith(N, V, L, roughness);
float3 F    = fresnelSchlick(max(dot(H, V), 0.0), F0);

float3 kS = F;
float3 kD = vec3(1.0) - kS;
kD *= 1.0 - metallic;     

float3 nominator    = NDF * G * F;
float denominator = 4.0 * max(dot(N, V), 0.0) * max(dot(N, L), 0.0) + 0.001; 
float3 specular     = nominator / denominator;

float NdotL = max(dot(N, L), 0.0);                
Lo += (kD * albedo / PI + specular) * radiance * NdotL; 

Environment (IBL)

基于图像的光照(Image based lighting, IBL)是一类光照技术的集合。其光源不是如前一节教程中描述的可分解的直接光源，而是将周围环境整体视为一个大光源。IBL 通常使用（取自现实世界或从3D场景生成的）环境立方体贴图 (Cubemap) ，我们可以将立方体贴图的每个像素视为光源，在渲染方程中直接使用它。这种方式可以有效地捕捉环境的全局光照和氛围，使物体更好地融入其环境。
与直接光照不同的是，我们现在需要处理半球上的积分。通过对反射方程拆解，我们可以将IBL的计算分成2个部分，diffuse和specular：

HDR贴图

需要指出，IBL的贴图依赖HDR格式来实现记录超越[0,1]之间值的数据。HDR格式每个通道存储8位，再以alpha通道存放指数。是一种有损压缩，但能提高效率。

等距柱状投影图equiretangular map

这张环境贴图是从球体投影到平面上，以使我们可以轻松地将环境信息存储到一张等距柱状投影图(Equirectangular Map) 中。有一点确实需要说明：水平视角附近分辨率较高，而底部和顶部方向分辨率较低,在大多数情况下，这是一个不错的折衷方案，因为对于几乎所有渲染器来说，大部分有意义的光照和环境信息都在水平视角附近方向。

diffuse IBL

要求解积分，我们只能通过离散求数值解的方法来进行，但对于shader来说，去通过大量采样实时求解半球积分的黎曼和不太现实。因此需要通过预计算的方式来优化。我们可以预先计算出每一个法线朝向半球的黎曼和，将其存储在一个贴图中，例如hdr图中。在运行时，用法线作为方向向量去采样这个hdr图。即可得到当前法线朝向半球上的黎曼和。

// the sample direction equals the hemisphere's orientation 
float3 normal = normalize(localPos);
float3 irradiance = float3(0.0);  
float3 up = float3(0.0, 1.0, 0.0);
float3 right = cross(up, normal);
up = cross(normal, right);

float sampleDelta = 0.025;
float nrSamples = 0.0; 
for(float phi = 0.0; phi < 2.0 * PI; phi += sampleDelta)
{
    for(float theta = 0.0; theta < 0.5 * PI; theta += sampleDelta)
    {
        // spherical to cartesian (in tangent space)
        vec3 tangentSample = vec3(sin(theta) * cos(phi),  sin(theta) * sin(phi), cos(theta));
        // tangent space to world
        vec3 sampleVec = tangentSample.x * right + tangentSample.y * up + tangentSample.z * N; 

        irradiance += texture(environmentMap, sampleVec).rgb * cos(theta) * sin(theta);
        nrSamples++;
    }
}
irradiance = PI * irradiance * (1.0 / float(nrSamples));

其中，采样后乘以sin(theta)是为了模拟天顶角在接近球顶时面积减小的现象。tangent相关的计算是为了将HDR的采样坐标转换到球面坐标。
之后，只需要在PBR计算的shader中采样这个图就可以得到IBL的diffuse部分。
为了引入粗糙度对于IBL diffuse的影响，我们使用一个大佬提出的roughness相关的frensel schlick

vec3 fresnelSchlickRoughness(float cosTheta, vec3 F0, float roughness)
{
    return F0 + (max(vec3(1.0 - roughness), F0) - F0) * pow(1.0 - cosTheta, 5.0);
}   

float3 kS = fresnelSchlickRoughness(max(dot(N, V), 0.0), F0);
float3 kD = 1.0 - kS;
float3 irradiance = texture(irradianceMap, N).rgb;
float3 diffuse    = irradiance * albedo;
float3 ambient    = (kD * diffuse) * ao; 

specular IBL

IBL Specular在正常情况下受两个变量的影响。

• 入射光方向
• 视角方向
  变量过多的积分难以求解， 因此需要使用到split sum approximation

业界: Split Sum Approximation

将镜面反射部分的积分拆解为2部分:

预滤波环境贴图

第一部分是一个关于入射角的半球积分。
计算的过程类似diffuse ibl，不同的是引入了粗糙度来模拟镜面反射的模糊程度。因为镜面反射的计算依赖于法线分布函数。因此需要计算不同粗糙度下的specular ibl。利用mipmap特性，我们不需要计算太多张类似的specular ibl，只需要将粗糙度更大的ibl塞入更高的mipmap等级中。利用mipmap的特性，shader会根据dxdy来插值两张mipmap之间的值，效率很高。而且粗糙度高的贴图有了mipmap的降采样，会自动呈现出更模糊的状态。
与之前求均匀球半球和的方式不同，镜面反射如果采用这种方式会非常模糊，效果不好。因此我们需要将采样的范围限定在镜面反射可能的范围中。也就是反射的镜面波瓣中。

为了实现这个特性，我们需要利用重要性采样来做到。
而我们求解积分的方式也从黎曼和改为了求解蒙特卡洛积分，利用随机采样的概率来近似积分的结果。而为了让蒙特卡洛积分有更快的收敛速度，可以引入一个低差异序列来生成。

float RadicalInverse_VdC(uint bits) 
{
    bits = (bits << 16u) | (bits >> 16u);
    bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
    bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
    bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
    bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
    return float(bits) * 2.3283064365386963e-10; // / 0x100000000
}
// ----------------------------------------------------------------------------
vec2 Hammersley(uint i, uint N)
{
    return vec2(float(i)/float(N), RadicalInverse_VdC(i));
}  

有了低差异序列的样本，我们需要将其映射到镜面波瓣中，可以利用之前提到过的ggx法线分布函数：

vec3 ImportanceSampleGGX(vec2 Xi, vec3 N, float roughness)
{
    float a = roughness*roughness;

    float phi = 2.0 * PI * Xi.x;
    float cosTheta = sqrt((1.0 - Xi.y) / (1.0 + (a*a - 1.0) * Xi.y));
    float sinTheta = sqrt(1.0 - cosTheta*cosTheta);

    // from spherical coordinates to cartesian coordinates
    vec3 H;
    H.x = cos(phi) * sinTheta;
    H.y = sin(phi) * sinTheta;
    H.z = cosTheta;

    // from tangent-space vector to world-space sample vector
    vec3 up        = abs(N.z) < 0.999 ? vec3(0.0, 0.0, 1.0) : vec3(1.0, 0.0, 0.0);
    vec3 tangent   = normalize(cross(up, N));
    vec3 bitangent = cross(N, tangent);

    vec3 sampleVec = tangent * H.x + bitangent * H.y + N * H.z;
    return normalize(sampleVec);
}

最后specular第一部分的计算如下：

vec3 N = normalize(localPos);    
    vec3 R = N;
    vec3 V = R;

    const uint SAMPLE_COUNT = 1024u;
    float totalWeight = 0.0;   
    vec3 prefilteredColor = vec3(0.0);     
    for(uint i = 0u; i < SAMPLE_COUNT; ++i)
    {
        vec2 Xi = Hammersley(i, SAMPLE_COUNT);
        vec3 H  = ImportanceSampleGGX(Xi, N, roughness);
        vec3 L  = normalize(2.0 * dot(V, H) * H - V);

        float NdotL = max(dot(N, L), 0.0);
        if(NdotL > 0.0)
        {
            prefilteredColor += texture(environmentMap, L).rgb * NdotL;
            totalWeight      += NdotL;
        }
    }
    prefilteredColor = prefilteredColor / totalWeight;

    FragColor = vec4(prefilteredColor, 1.0);

需要注意的是，由于计算ibl中不知道视角朝向，因此假设了n=v=r来简化计算。

预计算BRDF

split sum approximation的第二部分是BRDF的求解。由于我们假设了n=v=r，因此参数就只有粗糙度和出射光线。
这里分为两种流派。

2D LUT

UE4在[Real Shading in Unreal Engine 4, 2013]中提出，第二个求和项 ，使用Schlick近似后， F0可以从积分中分出来：

上式留下了两个输入（Roughness 和 cos θv）和两个输出（缩放和向F0的偏差（a scale and bias to F0）），即把上述方程看成是F0 * Scale + Offset的形式。 我们预先计算此函数的结果并将其存储在2D查找纹理（LUT，look-up texture）中。

这张红绿色的贴图，输入roughness、cosθ，输出环境BRDF镜面反射的强度。是关于roughness、cosθ与环境BRDF镜面反射强度的固有映射关系。可以离线预计算。
具体的取出方式为：

即UE4是通过把Fresnel公式的F0提出来，组成F0 * Scale +Offset的方式，再将Scale和Offset的索引存到一张2D LUT上。靠roughness和 NdotV进行查找。

函数拟合

Blinn-Phong:

float3 EnvironmentBRDF( float g, float NoV, float3 rf0 )
{
    float4 t = float4( 1/0.96, 0.475, (0.0275 - 0.25 \* 0.04)/0.96, 0.25 );
    t *= float4( g, g, g, g );
    t += float4( 0, 0, (0.015 - 0.75 * 0.04)/0.96, 0.75 );
    float a0 = t.x * min( t.y, exp2( -9.28 * NoV ) ) + t.z; float a1 = t.w;
    return saturate( a0 + rf0 * ( a1 - a0 ) );
}

GGX:

float3 EnvDFGLazarov( float3 specularColor, float gloss, float ndotv )
{
    float4 p0 = float4( 0.5745, 1.548, -0.02397, 1.301 );
    float4 p1 = float4( 0.5753, -0.2511, -0.02066, 0.4755 );
    float4 t = gloss * p0 + p1;
    float bias = saturate( t.x * min( t.y, exp2( -7.672 * ndotv ) ) + t.z );
    float delta = saturate( t.w );
    float scale = delta - bias;
    bias *= saturate( 50.0 * specularColor.y );
    return specularColor * scale + bias;
}

上式中的specularColor即F0。
EnvironmentBRDF函数的输入参数分别为光泽度gloss，NdotV，F0。和UE4的做法有异曲同工之妙，但COD：Black Ops 2的做法不需要额外的贴图采样，这在进行移动端优化时，是不错的选择。

参考

https://learnopengl-cn.github.io/07%20PBR/01%20Theory/
https://zhuanlan.zhihu.com/p/53086060