在跟着learning openGL学习的过程中，我被基于物理渲染的优秀效果震惊到了。在阅读了一些大佬的文章并且跟着写了一遍代码后，我决定做一些总结，并且希望他能保持足够的基础。第一次写文章难免有表述不清楚的地方，参考资料在最低端标出，可以去看看他们的。
1 基础知识
1.1 辐射度量学
既然想要达成基于物理的渲染效果，那基于物理的光照计算是必不可少的。在基于物理渲染中，辐射度量学的作用就是量化光线的种种属性。
其中，基于物理渲染中用得比较多的是以下几个：辐照度，辐射度，辐射率，立体角。
辐照度代表能量接收方接收了多少能量，是被辐射体上某一点附近某一微元面积上接收的辐射通量。
辐射度代表辐射体上某一点附近某一微元面积上辐射的总辐射通量，正好与辐照度相对应。如果这个微元面积无限小，就可以将其当成一个点光源。
辐射率代表在某一指定方向上单位立体角内，单位投影面积上发出的辐射通量。注意这里使用的是投影面积，如果面积无限小，立体角无限小，辐射率就可以理解为一个点光源发出的一条光线。
其中，辐照度和辐射率的区别在于，辐射度指一个点接收到的半球上的所有光线，而辐射率指一个点接收到的半球上的一条光线。
立体角：首先假设一个圆面上极小的面积dA，dA的面积可以看作由微元积分面的宽和高的积，然后通过面积除以半径的平方求得立体角。
注释：关于如何求得面积dA的宽和高：
在球面坐标系下ΔS就是球面上一块微小的面积，当Δφ和Δθ足够小的时候，ΔS的两边p和q可以看作以O和O'为圆心的圆的微小弧长，两个圆互相垂直。如果两个圆的半径分别为r和a:
1.2 反射方程
反射方程指描述了给定入射光和反射规律的情况下，从p点沿着方向w0射出的出射光。
在上式中，fr指得就是双向反射分布函数（Bidirectional Reflectance Distribution Function, BRDF），BRDF是用于描述反射规律的函数，同样描述反射规律的函数还有BSSDF，BTDF等用来描述散射，投射的函数，但我在这里不做说明（因为我也不会）。
这里的双向反射分布函数，双向指得就是入射方向和出射方向，BRDF描述了入射方向输入的能量，有多少被传递到出射方向上。
2. 基于物理的渲染
现在，我们终于站到了基于物理渲染的起跑线上。所谓基于物理的渲染（Physically Based Rendering，PBR）指使用基于物理原理和微平面理论（microfacet theory）建模的着色/光照模型，以及使用从现实中测量的表面参数来准确表示真实世界材质的渲染理念。
而在本篇文章中，我们根据接受的光照的不同，将PBR渲染分为直接光照和间接光照。而二者都包括漫反射和镜面反射。
注：之后的公式中，V（v）为视角方向向量，N（n）为整体法线向量，L（l）为光线入射向量，H（h）为V以及L的半角向量，如有其他特殊参数会进行相应说明
2.1 直接光照
在PBR渲染中，常用的物理模型为Cook-Torrance 模型，这也是我们这篇文章所实现的模型。其中 f lambert（fd）是漫反射项， f cook-torrance（fs）是镜面反射项，系数 kd 和 ks是为了保证了能量守恒而引入的项。
接下来，我将详细讲解各个选项。
2.1.1 漫反射 fd
漫反射项非常简单。由于在漫反射中，我们假定了所有地方射入的光都是均匀的且大小一致，所以最后得出的值就是从albedo贴图上取得的值除以PI。为什么要除以PI呢？其实是出于能量守恒的考虑。证明如下。
根据上述证明，我们很轻松地就得出了fd必须要小于1/PI。
2.1.2 镜面反射项 fs
镜面反射项是PBR中比较复杂的地方，实际上，他是对于微表面模型（Microfacet）的一个诠释。
那什么是微表面模型呢？微表面模型实际上是将所有含有一定粗糙程度的三角形面片，看作为由无数个法线方向不同的小的镜面组成的。
事实上，通常情况上讲，我们一个三角形面片只有一个面法线，也就是说三角形面片上的所有点都共用该法线。这导致了一种尴尬的情况：如果我们认定反射全部都是镜面的，那么所有的面片都可以看作是完全光滑的，因为他上面所有点的反射方向在入射方向确定的情况下都是唯一确定的。所以，为了能够实现不同粗糙度的效果，我们引入了微表面模型，让程序在计算光照时，考虑光线和每一个微小平面的作用，来实现不同粗糙度下的不同效果。
2.1.2.1 菲涅尔系数 F
菲涅尔项的作用是描述在不同观察角度下会有多少能量被反射。当入射方向接近掠射角度（grazing angle）的时候，光线是被反射的最多的,也就是当你的入射方向与法线几乎垂直时候,反射的辐射率（radiance）是最多的。
下面是对于绝缘体反射率与角度的关系:
而对于导体，有些值会导致反常情况：
菲涅尔项的推导时要考虑光线的S极化和P极化效果，因为要考虑不同介质的差别，所以公式比较复杂。
而在实际运算中，我们不可能采用如此复杂的公式，因为运算效率的关系，我们一般会采用如下近似 Schlick’s approximation：
当然，我们也可以选择Epic提供的拟合版本，进一步提升效率。
具体可以参见下方链接：
Real Shading in Unreal Engine 4
在上述两个公式中，都有一个相同的参数：F0。F0指得是当入射角度和法线方向垂直时的反射率，也称为基础反射率。
对于非金属物体，它的基础反射率是一个比值，而对于金属来说它的基础反射率是带颜色的RGB。如果物体介于金属和非金属之间，基础反射率就是物体颜色(金属)和基础反射率(非金属)的插值，因此我们需要增加一个新的参数金属度(Metalness)来表示物体的金属程度。
求得F0的方式也很简单，代码如下（怎么使用代码块？）：
vec3 F0 = vec3(0.04);
F0 = mix(F0, surfaceColor.rgb, metalness);//metalness为金属度
2.1.2.2 法线分布函数 D
法线分布函数决定了微表面的法线分布情况。当朝向比较分散的时候会得到比较扩散的结果,如果朝向比较集中指向时认为是光滑的，更集中时会呈现出类似镜面的效果。
在描述法线分布的NDF函数中，我们有多种模型可以选择，包括Beckmann模型，GGX模型，等等等等。在这里我只介绍Trowbridge-Reitz GGX模型，也是目前使用最多的一种模型。他的特点是高光周边会呈现出一种类似光晕的效果，过渡比较自然。
GGXTR的计算式如下。
其中，α指得是roughness（粗糙度）的平方。当然，也有的地方（比如寒霜引擎）使用的是roughness的四次方。
//TO DO：对于NDF函数原理的推导
2.1.2.3 几何函数 G
几何函数的作用是解决微表面的自遮挡问题。因为需要借助法线分布函数来进行剔除，所以两者经常是配套使用的。比如我这里使用的是GGXTR的NDF函数，那我的法线分布函数则需要使用由GGX函数推导的Schlick-GGX函数。
Schlick的式子如下：
其中α为粗糙度。同时，因为最后要计算入射时被遮挡的光线和出射时被遮挡的光线，所以要进行两次计算并相乘以得到最后的结果。
2.1.2.4 镜面反射公式的推导
具体可参考：参考资料
2.1.3 能量守恒
//TO DO 等我正式写了相关代码后再更新
2.1.4 相关glsl代码
float DistributionGGX(vec3 N, vec3 H, float roughness) //GGXRT NDF函数
{
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;
}
float GeometrySchlickGGX(float NdotV, float roughness)//Schlick-GGX
{
float r = (roughness + 1.0);
float k = (r*r) / 8.0;
float nom = NdotV;
float denom = NdotV * (1.0 - k) + k;
return nom / denom;
}
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;
}
vec3 fresnelSchlick(float cosTheta, vec3 F0)//菲涅尔项
{
return F0 + (1.0 - F0) * pow(clamp(1.0 - cosTheta, 0.0, 1.0), 5.0);
}
//省去参数的PBR直接光照计算
{
float NDF = DistributionGGX(N, H, rough);
float G = GeometrySmith(N, V, L, rough);
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;
vec3 specular = numerator / denominator;
vec3 kS = F;
vec3 kD = vec3(1.0) - kS;
kD *= 1.0 - metal;
float NdotL = max(dot(N, L), 0.0);
color += (kD * albedo / PI + specular) * radiance * NdotL;
}
2.2 基于图像渲染（IBL）
IBL是基于物理渲染的真实感的重要来源，是对环境光照的一种处理方案。对于大部分情况来说，环境光来自于天空盒，也就是cube-map贴图。因此，IBL的重点就在于如何从图像中获取光照信息。
首先还是老样子，分成漫反射和镜面反射两部分。对于基于图像的反射，我们还是使用Cook-Torrance光照模型。
2.2.1 漫反射积分
原理与直接光相同。我们只需要对cubemap上求卷积（计算朝向 N 的半球 中每个方向 wi的总平均辐射率），然后将结果保存在cubemap上。
那么，我们该如何对这个环境贴图求得对应卷积呢？我们将对环境贴图上的每一个纹素在其对应半球上生成固定数目，均匀分布的采样向量。同时，由于求解有关于立体角的积分的过程中，我们需要用到立体角，但它难以处理。所以我们我们使用球坐标 θ 和 ϕ 来代替立体角。
注：为什么我们要采样呢？因为在计算机中，我们是无法计算连续的积分的，因此我们采用采样的方法。采样的向量越多，就越接近真实积分的结果。为此，我们也演化出了非常多的采样方法。
在采用球坐标 θ 和 ϕ 后，以及明确了采样的思路后，我们的公式就变成了这样：
glsl采样积分代码如下：
#version 430 core
out vec4 Fragcolor;
in vec3 localPos;
uniform samplerCube environmentMap;
const float PI= 3.14159265359;
void main(){
//我们使用世界空间向量localPos作为表面的法线，因为在cubemap中，它将被插值到每个面，并形成一个单位球体
vec3 N = normalize(localPos);
vec3 irradiance = vec3(0.0);
//建立基向量
vec3 up = vec3(0.0, 1.0, 0.0);
vec3 right = cross(up, N);
up = cross(N, right);
float delta = 0.025;
float nSamples = 0.0;
//我们采用等距采样的方式来对环境贴图进行卷积
for(float phi = 0.0; phi < 2.0 * PI; phi +=delta){
for(float theta = 0.0; theta < 0.5 * PI; theta +=delta){
//从球面空间到笛卡尔坐标系
vec3 tangentSample = vec3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
//从切线空间到世界空间
vec3 sampleDir = tangentSample.x * right + tangentSample.y * up + tangentSample.z * N;
irradiance += texture(environmentMap, sampleDir).rgb * cos(theta) * sin(theta);
nSamples++;
}
}
irradiance = PI * irradiance * (1.0 / float(nSamples));
Fragcolor = vec4(irradiance, 1.0);
}
glsl采样漫反射贴图代码如下：
vec3 irradiance = texture(irradianceMap, N).rgb;//之前将朝向N的半球的平均辐射率保存在cubemap的N方向
2.2.2 镜面反射积分
在进行公式的相关内容之前，让我们先科普几个相关知识，以便于我们进行公式的推导：
1.近似公式：
这个近似公式在图形学中终极常用，并且只在两种情况下可以保持其的近似性。
1. g(x)的积分域很小
2. g(x)的值域变化很小
2.蒙特卡洛积分
蒙特卡罗积分建立在大数定律的基础上，并采用相同的方法来求解积分。不为所有可能的（理论上是无限的）样本值 x 求解积分，而是简单地从总体中随机挑选样本 N生成采样值并求平均。随着 N 的增加，我们的结果会越来越接近积分的精确结果。
其中 pdf 代表概率密度函数 (probability density function)，它的含义是特定样本在整个样本集上发生的概率。
镜面反射相对于漫反射更为复杂，因为它对于精确度的要求更高。对于镜面反射部分，我们通常希望通过各种方式获得一个效果更好的，性能更优的近似结果，因此诞生了例如球谐光照之类的种种方法。在这里我要介绍的是Epic所使用的split sum方法，它通过一系列近似获得较优的性能sh以及较好的效果。
先来看看公式：
是不是很复杂？对于镜面反射光照而言，我们不可能通过一个cubemap来解决问题，因为它涉及到了N ，w0，此外还有 BRDF 中的粗糙度和 F0（金属度和 albedo共同决定 ）。所以，我们必须要想出一个办法来简化这个过程。
因此，我们使用split sum 近似：
其中W为n*wi。
这个公式让人一眼蒙蔽。
我们分开来看，首先，左式是对于镜面反射方程进行蒙特卡洛积分。右边则是通过近似公式，将它拆分为光照部分，和BRDF部分。
2.2.2.1 光照部分
首先我们得用到有关于重要性采样的知识。
重要性采样指得是什么呢？
镜面反射的强度跟视线方向，宏观表面的法线以及微表面的粗糙度相关，所以不能像漫反射一样，在半球上均匀采样，这会导致其很难快速的收敛，如下图：
从中我们可以看出，对于镜面反射，只有在镜面波瓣的部分的入射光线才对镜面反射有贡献，所以在采样过程中，我们需要采用重要性采样，来针对性的进行采样。而在IBL中，我们采用GGX进行重要性采样，来采样波瓣部分。
公式推导：
我们得到的公式如下：
对于这个公式，我们还需要进一步简化。
对于BRDF来说，大概会在反射方向
取得值。同时，
我们假设从不同的方向入射，波瓣的形状变化不大。如此可得：
由此，得到公式结果。当然，由于我们没有关心视角的实际方向，我们在掠射角（grazing angles）看表面时没法得到拖长的反射。
2.2.2.2 BRDF部分
对于brdf部分，我们同样要进行一次近似。右侧的积分为
其中含有的变量为w0，wi，n，F0，粗糙度。其中，由于GGX是各向同性，所以我们只需要知道N和w0的夹角。现在需要的变量为θ，n，F0。
所以，我们可以将F0移出表达式。然后用一张图片来预计算θ和粗糙度，然后在实际运算中通过uv读取对应的点的brdf值。
公式计算如下：
至于scale和bias，最后可以求出来，然后放在一张图片中，其中R通道储存scale，G通道储存bias。在这里我就不做求解的公式推导，因为这个图哪里都能找到。
2.2.2.3 预计算环境贴图
在求得这两部分如何求解后，我们就可以进行环境贴图的预计算了。
在预计算环境贴图上，我们需要考虑不同粗糙度时，卷积结果的不同。因为我们不可能对roughness做太细的区分（由于性能关系），所以我们首先要做的，就是根据不同的roughness，来生成一定阶层的mipmap。通常，我们会生成5到6个mipmap。
低差异序列：
低差异序列是一种基于蒙特卡洛积分的采样思路，它能够生成更加均匀的采样向量，Hammersley 序列是基于 Van Der Corput 序列，该序列是把十进制数字的二进制表示镜像翻转到小数点右边而得。
通过将蒙特卡洛采样与低差异序列相结合，并使用重要性采样偏置样本向量的方法，我们可以获得很高的收敛速度。因为我们求解的速度更快，所以要达到足够的近似度，我们所需要的样本更少。
接下来，我们会通过实际代码来演示我之前所写的所有内容。
glsl预滤波贴图：
#version 430 core
out vec4 Fragcolor;
in vec3 localPos;
uniform samplerCube environmentMap;
uniform float roughness;
const float PI= 3.14159265359;
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;
}
// Van Der Corput 序列的巧妙生成方法 其等效于
//float VanDerCorpus(uintn,uintbase){floatinvBase =1.0/float(base);floatdenom =1.0;floatresult=0.0;for(uinti =0u; i //<32u; ++i) {if(n >0u) { denom =mod(float(n),2.0);result+= denom * invBase; invBase = invBase /2.0; n //=uint(float(n) /2.0); } }returnresult;}
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
}
//基于Van Der Corput 序列获得的Hammersley序列 在切线空间中获取采样向量
vec2 Hammersley(uint i, uint N)
{
return vec2(float(i)/float(N), RadicalInverse_VdC(i));
}
//重要性采样
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);
// 从球面坐标转到笛卡尔坐标系
vec3 H;
H.x = cos(phi) * sinTheta;
H.y = sin(phi) * sinTheta;
H.z = cosTheta;
// 从切线坐标系转化到世界坐标系
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);
}
void main()
{
vec3 N = normalize(localPos);
vec3 R = N;
vec3 V = R;
const uint SAMPLE_COUNT = 1024u;
vec3 prefilteredColor = vec3(0.0);
float totalWeight = 0.0;
for(uint i = 0u; i < SAMPLE_COUNT; ++i)
{
//通过重要性采样，生成一个朝向R方向附近的采样向量
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)
{
// 基于roughness来采样不同的环境贴图的mipmap 可以很好的去除粗糙度较高的预滤波贴图上的亮点
float D = DistributionGGX(N, H, roughness);
float NdotH = max(dot(N, H), 0.0);
float HdotV = max(dot(H, V), 0.0);
float pdf = D * NdotH / (4.0 * HdotV) + 0.0001;
float resolution = 2048.0; // resolution of source cubemap (per face)
float saTexel = 4.0 * PI / (6.0 * resolution * resolution);
float saSample = 1.0 / (float(SAMPLE_COUNT) * pdf + 0.0001);
float mipLevel = roughness == 0.0 ? 0.0 : 0.5 * log2(saSample / saTexel);
prefilteredColor += textureLod(environmentMap, L, mipLevel).rgb * NdotL;
totalWeight += NdotL;
}
}
prefilteredColor = prefilteredColor / totalWeight;
Fragcolor = vec4(prefilteredColor, 1.0);
}
glsl采样BRDF和预滤波贴图：
const float MAX_REFLECTION_LOD =5.0;
vec3 prefilteredColor = textureLod(specularMap, R, roughness * MAX_REFLECTION_LOD).rgb;
vec2 brdf = texture(brdfLUT, vec2(max(dot(N, V), 0.0), rough)).rg;
最终效果：
