[OpenGL]使用OpenGL实现基于物理的渲染模型PBR（上）

一、简介

本文介绍了 基于物理的渲染(Physically Based Rendering, PBR) 的基本概念，实现流程和简单的代码实现。
在本文中只实现了基于点光源的PBR，在之后的文章中会继续介绍 基于图像光源（IBL）的PBR。
按照本文代码实现后，可以实现以下效果：

二、PBR介绍

1. 什么是 PBR

根据 ChatGPT的介绍，PBR（Physically Based Rendering，基于物理的渲染） 是一类基于物理原理的渲染方法，旨在更准确地模拟现实世界中的光照与材质交互。PBR 依赖物理公式来描述光的反射、折射、散射等特性，使渲染结果更真实，并且在不同光照条件下具有更一致的表现。

可以简单的将 PBR 视作与 Phong、Blinn-Phong 等模型等类似，但是更加复杂也更加符合物理世界（只是更加符合，并不是完全和物理世界相同）的着色模型。如果你能实现一个基于 Phong 或者 Blinn-Phong 的渲染模型，就能实现基于 PBR模型 的渲染模型。

属于 PBR 的渲染模型需要满足以下三个条件：

• 基于微平面(Microfacet)的表面模型。
• 能量守恒。
• 应用基于物理的BRDF（Bidirectional Reflectance Distribution Function)）。

接下来本文主要介绍 PBR 中的基于，Lambert+Cook-Torrance 的 PBR 模型。

2. 基于 Lambert和Cook-Torrance的PBR

在基于 Lambert和Cook-Torrance 的PBR 模型中，渲染方程中的 BRDF 可以分为 漫反射项 和 镜面反射项两部分，如下所示：
$$Lo(p,wo)=\int fr(p,wi,wo)\ast Li(p,wi)\ast n\ast wi\mathrm{dwi}$$
$$fr(p,wo,wi)=kd\ast f_{lambert}+ks\ast f_{Cook-Torrance}$$
其中，$f_{lambert}$ 是入射光在材质内部经过折射后产生的漫反射，$f_{Cook-Torrance}$是入射光在材质表面直接反射产生的镜面反射。$kd$是着色点的折射率，$kd=1-ks$，$ks$是着色点的反射率，可以使用下面的菲涅尔项公式计算（即下面的$F$）。

2.1. 漫反射项(Lambert)

在 Lambert 模型中，漫反射项可以使用以下公式表示：
$$f_{lambert}=\frac{c}{\pi }$$

2.2. 镜面反射项(Cook-Torrance)

在 Cook-Torrance 模型中，镜面反射项可以使用以下公式表示：
$$f_{Cook-Torrance}=\frac{D\ast F\ast G}{4\ast \left(wo\ast n\right)\left(wi\ast n\right)}$$

镜面反射项主要 分为 D，F，G三项。其中：

• D 为法向分布函数，Normal Distribution Function，微表面中法线朝向半程向量 ℎ 的表面片段所占的比例。使用以下公式计算：
  $$D=ND{F}_{GGX-TR}\left(n,h,\alpha \right)=\frac{{\alpha }^2}{\pi {\left({\left(n\ast h\right)}^2\left({\alpha }^2-1\right)+1\right)}^2}$$
  其中 $\alpha$为着色点粗糙度，$n$为着色点法向，$h$为半程向量。
• F 为菲涅尔项，Fresnel Rquation，描述了光在表面上的反射率随视角变化的情况，即上面的$ks$。可以使用以下公式近似表示：
  $$F=F_{Schlick}\left(h,v,F_0\right)=F_0+\left(1-F_0\right){\left(1-\left(h\ast v\right)\right)}^5$$
  其中 $F_0$ 为着色点的基础反射率，$h$为半程向量，$v$为视线向量。
• G 为几何函数，Geometry Function，计算微表面之间的遮挡效应。G 项比价复杂，计算公式如下：
  $$G=G_{Schlick-GXX}\left(n,v,l,k\right)=G_{sub}\left(n,v,k\right)\ast G_{sub}\left(n,l,k\right)$$
  其中：
  $$\begin{gathered} k=(\alpha +1{)}^2/8，针对直接光照 \\ k={\alpha }^2/2，针对IBL光照 \end{gathered}$$
  $$\begin{gathered} G_{sub}\left(n,v,k\right)=G_{Schlick-GGX}\left(n,v,k\right)=\frac{n\ast v}{\left(n\ast v\right)\left(1-k\right)+k} \\ {G}_{sub}\left(n,l,k\right)=G_{Schlick-GGX}\left(n,l,k\right)=\frac{n\ast l}{\left(n\ast l\right)\left(1-k\right)+k} \end{gathered}$$
  $\alpha$为着色点粗糙度，$n$为着色点法向，$l$为光线向量，$v$为视线向量。

2.3. 总结

Cook-Torrance PBR 模型树状图如下所示：

为了实现 PBR 渲染，我们需有给出 着色的的各种材质属性，包括：反照率（albedo）、法向（normal）、金属度（metallic）、粗糙度（roughness）和环境光遮蔽系数（ambient occlision）。他们的含义如下：

• 反照率：反照率(Albedo)纹理为每一个金属的纹素(Texel)（纹理像素）指定表面颜色或者基础反射率。它只包含表面的颜色（或者折射吸收系数）。

• 法线：法线(normal)纹理可以逐片段的指定独特的法线，来为表面制造出起伏不平的假象。

• 金属度：金属(Metallic)纹理逐个纹素的指定该纹素是不是金属质地的。在 PBR 中规定金属只有 镜面反射项，没有漫反射项。

• 粗糙度：粗糙度(Roughness)贴图可以以纹素为单位指定某个表面有多粗糙。采样得来的粗糙度数值会影响一个表面的微平面统计学上的取向度。一个比较粗糙的表面会得到更宽阔更模糊的镜面反射（高光），而一个比较光滑的表面则会得到集中而清晰的镜面反射。

• 环境光遮蔽：环境光遮蔽(Ambient Occlusion)贴图或者说AO贴图为表面和周围潜在的几何图形指定了一个额外的阴影因子。比如如果我们有一个砖块表面，反照率纹理上的砖块裂缝部分应该没有任何阴影信息。然而AO贴图则会把那些光线较难逃逸出来的暗色边缘指定出来。

三、代码实现

在下面代码中实现的场景中存在两个点光源，n个具有相同反照率纹理、法向纹理、金属度纹理、粗糙度纹理和AO纹理的小球，基于 前面介绍的 PBR 渲染模型进行渲染。部分代码的讲解如下：

1. 顶点着色器

pbrShader.vert：

#version 330 core
layout (location = 0) in vec3 aPos;         //顶点位置
layout (location = 1) in vec3 aNormal;      //顶点法向
layout (location = 2) in vec2 aTexCoords;   //顶点纹理坐标

out vec2 TexCoords;
out vec3 WorldPos;
out vec3 Normal;

uniform mat4 projection;
uniform mat4 view;
uniform mat4 model;

void main()
{
    TexCoords = aTexCoords;
    WorldPos = vec3(model * vec4(aPos, 1.0));
    Normal = normalize(transpose(inverse(mat3(model))) * aNormal);
    gl_Position =  projection * view * vec4(WorldPos, 1.0);
}

2. 片段着色器

pbrShader.frag：

#version 330 core
out vec4 FragColor;
in vec2 TexCoords;
in vec3 WorldPos;
in vec3 Normal;

// material parameters
uniform sampler2D albedoMap; 	// Albedo贴图
uniform sampler2D normalMap; 	// 法向贴图
uniform sampler2D metallicMap;	// 金属度贴图
uniform sampler2D roughnessMap;	//粗糙度贴图
uniform sampler2D aoMap;		// AO贴图

// lights
// 两个 点光源的位置 和 颜色
uniform vec3 lightPositions[2];
uniform vec3 lightColors[2];

uniform vec3 camPos;

const float PI = 3.14159265359;

// 根据 法向贴图 得到 片段的法向
vec3 getNormalFromMap()
{
    vec3 tangentNormal = texture(normalMap, TexCoords).xyz * 2.0 - 1.0;
    vec3 Q1  = dFdx(WorldPos); // Q1 是屏幕 X 方向世界坐标系的变化向量
    vec3 Q2  = dFdy(WorldPos); // Q2 是屏幕 Y 方向世界坐标系的变化向量
    vec2 st1 = dFdx(TexCoords);// st1 是屏幕 X 方向纹理坐标的变化向量
    vec2 st2 = dFdy(TexCoords);// st2 是屏幕 Y 方法纹理坐标的变化向量
    // 假设 切线向量为 T，副切线向量为 B
    // 那么应该有:
    /*
    Q1 = T * st1.s + B * st1.t
    Q2 = T * st2.s + B * st2.t
    那么可以得到:
    [Q1 Q2] = [T B] [st1.x st2.x] = [T B] M
                    [st1.y st2.y]
    那么有：
    [T B] = [Q1 Q2] M^-1
    根据二维矩阵的性质得到 M^-1
    整理公式即可得到：
    T = (Q1*st2.t - Q2*st1.t) / det(M)
    那么：
    T = normalized (Q1*st2.t - Q2*st1.t)
    得到 T 即可根据 cross(N,T) 得到 B
    */
    vec3 N   = normalize(Normal);
    vec3 T  = normalize(Q1*st2.t - Q2*st1.t);
    vec3 B  = -normalize(cross(N, T));
    mat3 TBN = mat3(T, B, N);
    return normalize(TBN * tangentNormal);
}
// ----------------------------------------------------------------------------
// DFG 中的 D 项
float DistributionGGX(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 nom   = a2;
    float denom = (NdotH2 * (a2 - 1.0) + 1.0);
    denom = PI * denom * denom;

    return nom / denom;
}
// ----------------------------------------------------------------------------
// DFG 中 G 项的分量 Gsub
float GeometrySchlickGGX(float NdotV, float roughness)
{
    float r = (roughness + 1.0);
    float k = (r*r) / 8.0;

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

    return nom / denom;
}
// ----------------------------------------------------------------------------
// DFG 中的 G 项
float GeometrySmith(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 ggx2 = GeometrySchlickGGX(NdotV, roughness);
    float ggx1 = GeometrySchlickGGX(NdotL, roughness);

    return ggx1 * ggx2;
}
// ----------------------------------------------------------------------------
// DFG 中的 F 项
vec3 fresnelSchlick(float cosTheta, vec3 F0)
{
    return F0 + (1.0 - F0) * pow(clamp(1.0 - cosTheta, 0.0, 1.0), 5.0);
}
// ----------------------------------------------------------------------------
void main()
{		
    // gamma 校正
    vec3 albedo     = pow(texture(albedoMap, TexCoords).rgb, vec3(2.2));
    // 根据 金属度贴图 得到金属度
    float metallic  = texture(metallicMap, TexCoords).r;
    // 根据 粗糙度 得到粗糙度
    float roughness = texture(roughnessMap, TexCoords).r;
    // 根据 AO贴图 得到 AO
    float ao        = texture(aoMap, TexCoords).r;
    // 根据 法向贴图 得到法向
    vec3 N = getNormalFromMap();
    // 根据相机位置计算 V 向量
    vec3 V = normalize(camPos - WorldPos);

    // calculate reflectance at normal incidence; if dia-electric (like plastic) use F0 
    // of 0.04 and if it's a metal, use the albedo color as F0 (metallic workflow)    
    // 计算 基础反射率 F0
    vec3 F0 = vec3(0.04); 
    F0 = mix(F0, albedo, metallic);

    // reflectance equation
    vec3 Lo = vec3(0.0);
    // 假设存在 2 个点光源
    for(int i = 0; i < 2; ++i) 
    {
        // calculate per-light radiance
        // 光源向量 L
        vec3 L = normalize(lightPositions[i] - WorldPos);
        // 半程向量 H
        vec3 H = normalize(V + L);

        float distance = length(lightPositions[i] - WorldPos);
        // radiance 衰减(与距离相关)
        float attenuation = 1.0 / (distance * distance);
        vec3 radiance = lightColors[i] * attenuation;
        

        // Cook-Torrance BRDF
        // Cook-Torrance BRDF = fr
        // fr = fd + fs
        //    = kD*c/pi + kS*(D*F*G)/(4*(wo*n)*(wi*n))
        //    = kD*c/pi + kS*(D*G)/(4*(wo*n)*(wi*n))
        //    = (1-F)*c/pi + F*(D*G)/(4*(wo*n)*(wi*n))
        //    = (1-F)*c/pi + (D*F*G)/(4*(wo*n)*(wi*n))

        // 计算 fs

        // D 项
        // 当 主法向为 N, 理想半程法向为 H, 粗糙度为 roughness 时
        // 实际 半程法向 等于 H 的概率
        float NDF = DistributionGGX(N, H, roughness);   

        // G 项
        // 当 主法向为 N, 相机向量为 V, 光源向量为 L, 粗糙度为 roughness 时
        // 整体 反射光线未被 阴影遮挡的概率 
        float G   = GeometrySmith(N, V, L, roughness);      
        
        // F 项
        // 反射的概率, (1-F 为折射的概率)
        vec3 F    = fresnelSchlick(max(dot(H, V), 0.0), F0);
        
        
        vec3 numerator    = NDF * G * F; 
        float denominator = 4.0 * max(dot(N, V), 0.0) * max(dot(N, L), 0.0) + 0.0001; // + 0.0001 to prevent divide by zero
        // ks*(D*F*G)/(4*(wo*n)*(wi*n)) = 
        // F *(D*G)/(4*(wo*n)*(wi*n))
        vec3 specular = numerator / denominator;
        
        // kS is equal to Fresnel
        // 反射率 ks
        vec3 kS = F;
        // for energy conservation, the diffuse and specular light can't
        // be above 1.0 (unless the surface emits light); to preserve this
        // relationship the diffuse component (kD) should equal 1.0 - kS.
        // 折射(漫反射)的比例
        // 折射率 kd
        vec3 kD = vec3(1.0) - kS;
        // multiply kD by the inverse metalness such that only non-metals 
        // have diffuse lighting, or a linear blend if partly metal (pure metals
        // have no diffuse light).
        // 只有 非金属有漫反射项, 金属没有漫反射项
        kD *= 1.0 - metallic;	  

        // scale light by NdotL
        // cos(N,L)
        float NdotL = max(dot(N, L), 0.0);        

        // add to outgoing radiance Lo
        // Lo = 漫反射项 + 镜面反射项
        //    = (漫反射 + 镜面反射) * radiance * cos(N,L)
        Lo += (kD * albedo / PI + specular) * radiance * NdotL;  // note that we already multiplied the BRDF by the Fresnel (kS) so we won't multiply by kS again
    }   
    
    // ambient lighting (note that the next IBL tutorial will replace 
    // this ambient lighting with environment lighting).
    // vec3 ambient = vec3(0.03) * albedo * ao;
    // vec3 ambient = vec3(0.00) * albedo * ao;
    // vec3 color = ambient + Lo;
    vec3 color = Lo;

    // HDR tonemapping
    color = color / (color + vec3(1.0));
    // gamma correct
    color = pow(color, vec3(1.0/2.2)); 

    FragColor = vec4(color, 1.0);
}

3. 全部代码及模型文件

使用OpenGL实现PBR的全部代码以及模型文件可以在OpenGL使用OpenGL实现基于物理的渲染模型PBR（上） 中下载。
下载源代码后使用以下命令编译运行：

mkdir build
cd build
cmake ..
make
./OpenGL_PBR

渲染结果如下：

四、参考

[1].LearnOpenGL-PBR