【图形学】光照模型-从最简单开始到BRDF实现

最新推荐文章于 2022-09-26 11:22:27 发布

程序员菜鸟

最新推荐文章于 2022-09-26 11:22:27 发布

阅读量6.3k

点赞数 4

分类专栏： OpenGL 文章标签： html5 html css 1024程序员节

本文链接：https://blog.csdn.net/qq_40947718/article/details/120837634

版权

OpenGL 专栏收录该内容

4 篇文章 2 订阅

订阅专栏

不积跬步无以至千里，想要了解复杂的光照模型，就要从最简单的入手，看看他们是如何一步步“进化”成完善的样子。每天300字，看看能否在一个月能完结这个主题吧~

首先从基本光照模型开始：

1.Lambert模型（漫反射）

Lambert = kd * Ia （环境光）

其中Ia表示环境光强度，kd表示反射系数，两者成绩就是环境光和材质交互反射后得到的光强度。

LambertMinor = kd *Im * Cos(θ)

其中kd 和 Im依然分别是入射光强度和反射系数，θ是入射光方向与顶点法线的夹角，称入射角，这里考虑了光照和平面的夹角来模拟反射光强度

最终综合两个公式Lambert光照模型即为
ILambert = ILambertdiff + ILambertMinor = kd * Ia + kd * Im * Cos(θ)

2.Phong模型（镜面反射）

Phong模型认为镜面板设的光强度与光线和视线的夹角相关。

$I_{spec }= K_s * I_i *dot(V,R)^{N_s}$

其中 $K_s$ 是镜面反射系数， $I_i$ 是光强度，乘以视角和反射光线的点乘的N次幂， $N_s$ 是高光指数

3.Blinn-Phong光照模型（修正镜面光）

$I_{spec }= K_s * I_i *dot(N,H)^{N_s}$

Blinn-Phong是Phong的改良版本，区别就在于把视线和反射光线的点乘换成了半程向量和法向量的点乘，半程向量H是光入射方向L和视线方向V的中间向量。

$H = (L + V) / |L + V|$

接下来就是直接进阶，基于物理的渲染，PBR

说起PBR，我们首先要了解它的理论是什么，理论分为三个部分
1 基于微平面的表面模型

2 能量守恒

3 基于物理的BRDF

接下来我们来逐一了解这些理论，首先是微平面理论：

微平面理论认为一个平面是由若干个微小的平面组成的，这些微小平面的朝向随着整个平面的粗糙程度的不同，分布也会不同，一般来说，一个平面越是粗糙，那么微平面分布就越混乱。

我们在渲染时用到的PS是逐像素渲染的，一个像素所对应的平面可能仍然由若干个微平面组成，这些微平面的朝向整齐程度，也会影响整个平面的反射能力，微平面的朝向越趋于整个平面的法线，那么就越光滑，反射能力越强。

这时候我们可以借用之前的向量来“计算”每一个微平面的朝向。那就是H向量，试想：对于视线在反射方向上，此时整个平面的H向量就代表着这个平面的法线，微观化讨论这个问题：对于每一个微平面的入射光线和观察向量得到的H向量，就代表了这个平面的法线，当然这是近似得到的，因为观察向量不可能正好在这个微平面的反射方向上。

然后是能量守恒理论：

这个比较简单，属于高中知识，简单来说就是出射光线的能量不能超过入射光线的能量，还有一部分能量进入了物体发生了折射，但是折射的光线并不完全转换成了物体的热能，有一部分在物体内兜兜转转又反射出了物体表面，沿着随机的方向开始发散。这部分折射的能量就构成了物体的漫反射的颜色，也就是最开始的Lambert模型的漫反射（只和入射光有光，和角度还有观察方向无关）

对于金属来说，没有漫反射的颜色，所有折射光都会被吸收，所以我们可以简单的把物体进入我们视线中的能量分为两部分：折射，反射，两者相加占比为1（理想情况啦~

接下来我们就要引入一个反射率方程，也就是渲染方程的一个精细版本，在讲解这个复杂的方程之前，先说一些能量学的概念。

辐射率：单位方向上光线的强度。

辐射通量：一个光源所输出的能量。

立体角：如下图：

辐射强度：单位立体角所包含的辐射通量，如上图，球心所在位置是光源的位置。

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

这里在learnopengl中的解释我看的不是很明白，感觉和我的理解正好相反，以下是我从别的地方找来的解释：辐射率概念：每单位面积下每单位立体角的辐射通量密度。

这里的光是入射光，定义为 $L=\frac{d \Phi}{d \omega d A^{\perp}}$

辐射通量对角度取微分，之后再对面积取微分，这里为什么要对A⊥取微分而不是A呢，因为！！！！因为这里少了两个关键字！！投影，我找了若干篇文章才提到，辐射率也就是辐射亮度，其实是单位投影面积单位立体角的辐射通量，这也就解释了，为什么角度越大，强度越大，因为角度越大，对于dA来说， L是固定的，角度越大，单位面积上包含的L就越多，因为角度越大，dA⊥投影到dA的面积就越小，这样对整个dA取辐射通量得到的值自然是越大的，所以辐射率最终的解释是：单位投影面积下每立体角的辐射强度为中（此处为阿拉伯符号）的辐射通量。

在实际计算光照的时候，我们要考虑到每个方向辐射能量的影响，也就是一个以平面法线为轴所环绕形成的半球领域

要计算这个领域里的所有辐射亮度，就要对这个半球领域进行二重积分，我们也可以进行离散计算，把积分变成一个等步长值的求和取均值，这种方法就是黎曼和。

接下来，前面的概念都介绍完了，该介绍重头戏了，双向反射分布函数，它的作用是基于表面材质属性来对入射辐射率进行缩放或者加权，通俗理解就是，基于入射光的角度，材质粗糙度，金属度，来得到反射的辐射亮度。

BRDF函数的输入是：入射光，观察方向，平面法线和一个平面粗糙度，就可以得到这束光线对于某个平面最终反射出来的光线所作的贡献程度。前面说到的blinn-Phone也是一种BRDF，但它比较简单，没有能量守恒，现有的几乎所有实时渲染管线使用的都是一种Cook-Torrance BRDF模型（这一点笔者没有考察过，完全是抄过来的

Cook-Torrance BRDF模型分为两个部分，镜面部分和漫反射部分。

$f_r = k_df_{lambert}+k_sf_{cook-torrance}$

左边是漫反射，右边是镜面反射部分。继续拆解公式

$f_{lambert} = \frac{c}{\pi }$

这里的c是表面颜色，Π是为了标准化（后面会讲）漫反射里的和法线点乘被挪到了积分方程里了。

$f_{cook-torrance} = \frac{DFG}{4(\omega_o*n)(\omega_i*n)))}$

这里包含三个函数：D，F，G。分母还是一个标准化因子

D正态分布函数：估算在受到表面粗糙度的影响下，取向方向与中间向量一致的微平面的数量。这是用来估算微平面的主要函数。（微观法线和宏观法线相同的平面的比例

F几何函数：描述了微平面自成阴影的属性。当一个平面相对比较粗糙的时候，平面表面上的微平面有可能挡住其他的微平面从而减少表面所反射的光线。

G菲涅尔方程：在不同的表面角下所反射的光线所占的比例。

下面我们详细介绍下每个函数：

正态分布函数：从统计学上近似的表示了微观平面法线和宏观法线一致的比例。输入为粗糙度和宏观平面的法线，宏观平面的半向量。公式如下图：

$NDF_{GGXTR}(n,h,\alpha ) = \frac{\alpha^2}{\pi((n*h)^2)(\alpha^2-1)+1)^2}$

我们可以得到不同粗糙度下的效果图：

当粗糙度很低时，与中间向量取向一致的微平面会高度集中在一个范围，随着粗糙度的增加，这个范围越来越大，并且中央最大值也开始减小。怎么理解这个函数的表现呢：

我们先看一下正态分布函数的曲线图：

正太分布函数沿着x轴的积分和为1。粗糙度alpha相当于公式中的sigma，粗糙度越小，这个函数越”尖锐“，中间最大值越大。这个时候把n，h的点乘作为输入，当n = h，也就是我们的观察方向刚好和出射方向相同时，这个平面应该获得最大概率值，看公式，n*h越大，由于alpha平方-1 是负值，分母反而越小，结果反而越大，是不是和我们观察结果对应上了呢。

记下来是几何函数：也是从统计学上近似求微平面相互遮蔽的概率。我们采用的公式是Schlick-GGX:

$G_{SchilickGGX}(n,v,k) = \frac{n*v}{(n*v)(1-k)+k}$

公式里的k就是粗糙度基于直接光照还是间接光照的重映射：

$k_{direct} = \frac{(\alpha+1)^2}{8}$

$k_{IBL} = \frac{\alpha^2}{2}$

为了更好的考量这个公式，我们需要把观察方向和光线入射方向都考虑进去，我们做两次函数计算，并把结果都作为因子相乘。

$G_{SchilickGGX}(n,v,l,k) = G_{SchilickGGX}(n,l,k) * G_{SchilickGGX}(n,v,k)$

几何遮蔽函数是值域为[0.0,1.0]的乘数。效果图如下

菲涅尔方程：描述的是被反射的光线与被折射光线的比例。

当垂直观察的时候，任何物体或者材质表面都有一个基础反射率(Base Reflectivity)，但是如果以一定的角度往平面上看的时候所有反光都会变得明显起来。你可以自己尝试一下，用垂直的视角观察你自己的木制/金属桌面，此时一定只有最基本的反射性。但是如果你从近乎90度（译注：应该是指和法线的夹角）的角度观察的话反光就会变得明显的多。如果从理想的90度视角观察，所有的平面理论上来说都能完全的反射光线。这种现象因菲涅尔而闻名，并体现在了菲涅尔方程之中。

菲涅尔方程还存在一些细微的问题。其中一个问题是Fresnel-Schlick近似仅仅对电介质或者说非金属表面有定义。对于导体(Conductor)表面（金属），使用它们的折射指数计算基础折射率并不能得出正确的结果，这样我们就需要使用一种不同的菲涅尔方程来对导体表面进行计算。由于这样很不方便，所以我们预先计算出平面对于法向入射（）的反应（处于0度角，好像直接看向表面一样）然后基于相应观察角的Fresnel-Schlick近似对这个值进行插值，用这种方法来进行进一步的估算。这样我们就能对金属和非金属材质使用同一个公式了。

$Fresnel-Schlick = F_0+(1-F_0)(1-cos(\theta ))^5$

了解完这些概念我们就可以继续往下走，看PBR光照到底是怎么样产生的啦~

//坚持是一件很难的事，但在这纷纷扰扰的世界里，除了坚持，有什么能让你我和平庸区分呢~

前面的内容，我们大搞搞懂了方程里的每一个符号是用来干什么的，它有什么含义，接下来我们要把这个方程具体运用到编码上去，看看究竟能实现怎样的效果。

这是方程最终的样子。

先看一张简单的图，对于p点而言，它所有的直接入射光都来自这个点光源wi方向的辐射通量，我们可以获得光源的辐射通量，入射方向，法线方向，夹角，以及衰减，这不就是我们之前几个简单模型的漫反射的计算方法嘛，现在它改名叫辐射率了~

如果我们知道场景里所有直接光源的位置，我们只需要把他们的辐射率累加起来然后除以项数就能得到平均值了。当然对于每一个辐射率，我们都希望能用上BRDF项。

先从菲涅尔反射入手，因为这个公式最为简单，

vec3 fresnelSchlick(float cosTheta, vec3 F0)
{
    return F0 + (1.0 - F0) * pow(1.0 - cosTheta, 5.0);
}  
vec3 F0 = vec3(0.04); 
F0      = mix(F0, albedo, metallic);
vec3 F  = fresnelSchlick(max(dot(H, V), 0.0), F0);

F0我们取（0.04，0.04，0.04），这对大多数非金属材质都适用，最终我们可以用F0和金属度来插值得到最终的F0。

接下来是正态分布函数D和G，直接上代码：

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;
}

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;
}
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;
}

有了这些函数，我们就能最终得到BRDF的计算结果。

vec3 nominator    = NDF * G * F;
float denominator = 4.0 * max(dot(N, V), 0.0) * max(dot(N, L), 0.0) + 0.001;  //分母中加了0.001防止除0
vec3 specular     = nominator / denominator;

现在我们终于可以计算每个光源在反射率方程中的贡献值了！因为菲涅尔方程直接给出kS，我们可以使用F表示镜面反射在所有打在物体表面上的光线的贡献。从kS我们很容易计算折射的比值kD：

const float PI = 3.14159265359;

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

L0就是最终出射光下的辐射率，其实是我们直接求和（最简单的积分）的结果，最后加上环境光项，就得到了这个片段最终的颜色。

现在我们已经可以把BRDF简单的运用到我们的PS里了，现在抛开直接光照，我们把整个环境都当成一个光源，这样可以更好的模拟环境的全局光照和氛围。那这和刚才的直接光照计算有什么区别呢，再看一下反射方程

我们需要对半球领域内的每一个wi，都能拿到对应的入射辐射度，并且实时地把积分计算出来，这就需要很多预计算。

首先我们可以把反射方程拆开成两个部分，一个是漫反射部分，一个是镜面反射部分。

要计算漫反射部分，我们可以创建一个新的贴图，这个贴图的采样结果，就是全局环境光在这个方向上的入射辐射度，暂时把这个贴图叫做预计算积分贴图，怎么得到这个贴图呢，首先我们知道原始的环境光贴图，通过对原来的环境光贴图进行采样，可以得到这个方向上的单一的环境光的入射辐射度，那么我们可以对半球领域内的原始环境光贴图进行离散采样，并且取平均值，最终在据计算贴图中存储采样结果，预计算积分贴图又叫做辐照度图。接下来就是比较复杂的离散采样过程了：

在半球领域内采样，主要是对两个立体角方向采样，

航向角的采样范围是0-2pi，倾斜角的采样范围是0-1/2 pi。最终积分方程变成上图，为什么加了一个sin(theta) 是因为，倾斜角越高，积分面积就越小，贡献度就越小，所以需要平衡一下权重。

在两层循环内，我们获取一个球面坐标并将它们转换为 3D 直角坐标向量，将向量从切线空间转换为世界空间，并使用此向量直接采样 HDR 环境贴图。我们将每个采样结果加到 irradiance，最后除以采样的总数，得到平均采样辐照度。请注意，我们将采样的颜色值乘以系数 cos(θ) ，因为较大角度的光较弱，而系数 sin(θ) 则用于权衡较高半球区域的较小采样区域的贡献度。

（这一段是笔者复制来的，因为讲的很好

这样我们就得到了一个预处理过后的辐照度图，对于每一个片段（或者说每一个方向的采样向量），它所对应的半球领域都是不同的，最终存储在采样方向所对应的贴图位置上。现在IBL的漫反射部分就从卷积计算，转换成了获取辐照度图的一次采样计算了。

还有一个数字需要计算就是菲涅尔系数，之前在讨论的时候，我们考虑用的是微观平面的h向量，但是具体到一个片段上，由于入射光来自四面八方，没法得到一个h向量，我们就用法线和视线夹角临时替代一下（也能大概反映视线关系，我是这么理解的），但是这样就没法体现粗糙度了，所以可以魔改一下菲尼尔方程，把粗糙度也考虑进去

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

大概全局光照的漫反射就是长这样了，或许看起来效果还是不那么突出，是因为还没有加上镜面反射，表面还没办法真正的反射出来环境的倒影。

---------------------------

最近遇上了一些小问题，很是苦恼，希望能尽快过去吧。

接下来我们继续看反射方程的镜面部分：

和辐照图相同的原因，我们也没法实时的计算这个卷积，但是和辐照图不一样，它的积分部分值以来与wi，但是这一次，它的积分部分来以来wo，我们没办法用两个方向向量去采样贴图，相当于先采样一张辐照度图，然后再根据wo计算得到若干个不同的二级辐照度图，这是不可能的事情。不要着急，我们先把积分拆开

积分的第一部分，叫做预滤波环境贴图，它有点类似于辐照度图，是预先计算的环境卷积贴图，（不明白为什么一定要叫卷积，明明就是一个积分而已）但是我们还需要考虑粗糙度，因为越粗糙，我们的采样向量越分散。（采样向量就是N，法线）

接下来至关重要的一步，我们要消元啦，

我们使用 Cook-Torrance BRDF 的正态分布函数(NDF)生成采样向量及其散射强度，该函数将法线和视角方向作为输入。由于我们在卷积环境贴图时事先不知道视角方向，因此 Epic Games 假设视角方向——也就是镜面反射方向——总是等于输出采样方向ωoωo，以作进一步近似。翻译成代码如下：

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

这里让每个视角方向都等于微观平面上的光线出射方向，那么正态分布函数就恒等于1，整个NDR公式相当于只剩下了D，R。

再看方程的第二部分，就是BRDF的积分，首先wo已经是已知的了，对于任意给定的wi和N还有粗糙度，我们都可以直接计算出BRDF值，把计算好的粗糙度和入射角的组合响应结果存在一张2D查找纹理上，我们称为BRDF积分贴图，（这个贴图只用到了R和G来保存值），这里的R和G代表什么我们后面再说。

先看几个“简单”（看了好几周才明白）的概念。

1 蒙特卡洛积分，是一个求解连续积分近似的方法，这里直接拷贝了别人的图片来说明

相当于以1/N为步长取值求和最后取平均。

2 重要度采样，用来优化蒙特卡洛的均值采样的。

对于被积分的函数，每个采样值可能所占的比重不同，简单说就是，被积函数的自变量，xi，可能不是均匀分布在[a,b]区间的，对于那些采样概率高的点，我们应该给更多的权重，才能更加贴近真实情况，（其实这里=后面的积分表达式不准确，因为x是属于[a,b]的，但不是均匀分布在[a,b]上的，没法通过数据积分公式直接得到，只是借用积分公式来表达准确结果。

那怎么判断xi的分布呢，我们需要用到f(x)本身的概率密度函数来构造采样点的分布情况，简单说概率密度函数的定积分结果就是自变量在某个区间的概率。（密度积分即为质量~）

把积分区间无穷地缩小，那么对于被积函数fx来说，p(xi)就是自变量取值为xi的概率，F(x）相当于是若干个f(xi）值的叠加，p(xi）不仅是自变量取值为xi的概率，也是值域为f(xi)的概率，也就是说，p(xi)是f(xi)占最终积分结果的比重

有上面这个公式，我们就可以根据一个单一的函数值求得最终的积分结果，不过误差比较大，对此我们再采用蒙特卡洛积分则有

这就是重要性采样，我们取概率密度最高的N个点进行采样，这结果肯定比之前均匀采样得到的采样结果要准确的多。对于原始的蒙特卡洛积分的均匀分布采样，pdf为

//我又来啦，虽迟但到

采样分布映射

现在我们知道需要采样N个pdf较大的点了，问题是，我们怎么知道哪些点pdf较大呢？

有一种方法，叫反函数法，我们可以通过均匀采样的N个点，通过累计密度函数cdf的反函数，来映射出新的N个采样点，这N个新的点就是我们需要的结果。

cdf是pdf从样本区间下限到当前采样点区间的概率值：

我们在cdf的y轴上均匀选取N个点，对于任意的i，有

取反函数就是

那么x1,x2。。。就是我们要采样的点。

怎么理解这件事呢，

你买面值不超过1元的刮刮乐中奖的概率不超过0.1，买超过1元但不超过2元的中奖概率是不超过0.3，买超过2元但不超过3元的中奖概率不超过0.6。

那pdf是这样的：

p([0, 1]) = [0.0, 0.1] — 记pdf_a

p((1, 2]) = (0.1, 0.3] — 记pdf_b

p((2, 3]) = (0.3, 0.6] — 记pdf_c

cdf就是

P([0, 1]) = [0.0, 0.1]

P((1, 2]) = (0.1, 0.4]

p((2, 3]) = (0.4, 1.0]

我们对y均匀的取值，最终得到的x，肯定是落在概率较大的区间的数量较多，

超过2元的 -> 1.0, 0.9, 0.8, 0.7, 0.6, 0.5

超过1元不超过2元的 -> 0.4, 0.3, 0.2

不超过1元的 -> 0.1, 0.0

最终我们得到的xi 肯定是pdf越大的地方，采样点越密集/

直观上也不难理解，因为pdf可以理解为cdf的导数，反应到函数图像上也就是cdf增长的速度，自然在增长速度快的地方，其取值区间越大，在均匀采样的情况下，落在区间内的点自然越多。

我们需要怎么去实现呢：首先BRDF的pdf函数如下：

物理意义是法线分布函数D在法线n方向上的投影，也就是在宏观平面上的投影，抄个图：

也许读者在其他地方，或者在本文后续的部分，发现pdf函数多了个分母，变成了下面这个样子，请不要惊讶。

这是因为(1.10)(1.10)是以dhi→dhi→为微元的，而在有些积分方程的计算里(后面就会碰到)，我们是以dli→dli→为微元的，所以要乘以一下才能将积分微元转换过来。

最终我们还是回到一开始的蒙特卡洛积分公式上去，BRDF的积分就可以转换成对应的函数除以PDF的形式

好了公式的问题解决了，我们最终可以通过若干离散的点求得BRDF在二维平面上的积分了。

那么我们怎么去在二维球面上均匀采样呢，有一个很简单的方法叫做Hammersley sampling，这里就不过多赘述了，感兴趣的可以去查查。

需要注意的是，反函数法是一个一维的算法，我们需要的是二维的重要性采样，那怎么办呢这里略去若干字，可以采用边缘分布降维打击，最终可以求得我们通过均匀采样得到的Pdf的较大值的极角坐标值。

有了pdf较大的坐标值，我们就可以得到N个h向量去求解我们上面说到的L光照公式了公式了。

时隔两周，由于工作的繁忙（呸由于本人的懒惰，这篇文章现在才得以继续

让我看看上次讲到哪里了。。。上回书说到，我们利用反函数法得到了若干个PDF较大的h向量，来求解BRDF公式。

再次回顾来看我们的镜面IBL公式

第一部分是类似于漫反射的辐射照度图采样，第二部分我们用蒙特卡洛积分计算得到了一个近似值，输入是法线和出射角的夹角以及粗糙度，输出结果是一个BRDF积分贴图(专业一点叫LUT Look Up Texture) 具体看下shader中的代码实现

vec2 IntegrateBRDF(float NdotV, float roughness)
{
    vec3 V;
    V.x = sqrt(1.0 - NdotV*NdotV);
    V.y = 0.0;
    V.z = NdotV;

    float A = 0.0;
    float B = 0.0;

    vec3 N = vec3(0.0, 0.0, 1.0);

    const uint SAMPLE_COUNT = 1024u;
    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(L.z, 0.0);
        float NdotH = max(H.z, 0.0);
        float VdotH = max(dot(V, H), 0.0);

        if(NdotL > 0.0)
        {
            float G = GeometrySmith(N, V, L, roughness);
            float G_Vis = (G * VdotH) / (NdotH * NdotV);
            float Fc = pow(1.0 - VdotH, 5.0);

            A += (1.0 - Fc) * G_Vis;
            B += Fc * G_Vis;
        }
    }
    A /= float(SAMPLE_COUNT);
    B /= float(SAMPLE_COUNT);
    return vec2(A, B);
}
// ----------------------------------------------------------------------------
void main() 
{
    vec2 integratedBRDF = IntegrateBRDF(TexCoords.x, TexCoords.y);
    FragColor = integratedBRDF;
}

看着是不是一脸懵逼，感觉和之前讲的蒙特卡洛积分积分好像一点关系没有，其实不是，我们细细品一下~

首先，我们用用Hammersley产生随机向量，再得到重要性采样的h向量，根据h得到输入光线向量接下来就是我们的蒙特卡洛积分积分部分了，看仔细了，根据蒙特卡洛积分有

展开fs和ps得到

上面就是纯BRDF的积分结果，懒得打字了直接截图下面内容

用蒙特卡洛积分表示上述公式就是

Fo* scale + bias

shader代码中输出的就是F0的系数和偏移值，我们最终只需要去BRDF积分贴图中采用即可。大功告成。最终的镜面IBL就是两个部分，一个是辐照度图的积分，一个是BRDF的采用结果，两者相乘就是镜面IBL啦~

最后奉上镜面IBL的fs部分：

其中留了一个小坑~ 最终结果没有加上直接光照的积分结果~ 看起来除了四个光点，其他的IBL效果都差不多

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

// material parameters
uniform vec3 albedo;
uniform float metallic;
uniform float roughness;
uniform float ao;

// IBL
uniform samplerCube irradianceMap;
uniform samplerCube prefilterMap;
uniform sampler2D brdfLUT;

// lights
uniform vec3 lightPositions[4];
uniform vec3 lightColors[4];

uniform vec3 camPos;

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;
}
// ----------------------------------------------------------------------------
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;
}
// ----------------------------------------------------------------------------
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);
}
// ----------------------------------------------------------------------------
vec3 fresnelSchlickRoughness(float cosTheta, vec3 F0, float roughness)
{
    return F0 + (max(vec3(1.0 - roughness), F0) - F0) * pow(clamp(1.0 - cosTheta, 0.0, 1.0), 5.0);
}   
// ----------------------------------------------------------------------------
void main()
{		
    vec3 N = Normal;
    vec3 V = normalize(camPos - WorldPos);
    vec3 R = reflect(-V, N); 

    // 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)    
    vec3 F0 = vec3(0.04); 
    F0 = mix(F0, albedo, metallic);

    // reflectance equation
    vec3 Lo = vec3(0.0);
    for(int i = 0; i < 4; ++i) 
    {
        // calculate per-light radiance
        vec3 L = normalize(lightPositions[i] - WorldPos);
        vec3 H = normalize(V + L);
        float distance = length(lightPositions[i] - WorldPos);
        float attenuation = 1.0 / (distance * distance);
        vec3 radiance = lightColors[i] * attenuation;

        // Cook-Torrance BRDF
        float NDF = DistributionGGX(N, H, roughness);   
        float G   = GeometrySmith(N, V, L, roughness);    
        vec3 F    = fresnelSchlick(max(dot(H, V), 0.0), F0);        
        
        vec3 numerator    = NDF * G * F;
        float denominator = 4 * max(dot(N, V), 0.0) * max(dot(N, L), 0.0) + 0.0001; // + 0.0001 to prevent divide by zero
        vec3 specular = numerator / denominator;
        
         // kS is equal to Fresnel
        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.
        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
        float NdotL = max(dot(N, L), 0.0);        

        // add to outgoing radiance Lo
        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 (we now use IBL as the ambient term)
    vec3 F = fresnelSchlickRoughness(max(dot(N, V), 0.0), F0, roughness);
    
    vec3 kS = F;
    vec3 kD = 1.0 - kS;
    kD *= 1.0 - metallic;	  
    
    vec3 irradiance = texture(irradianceMap, N).rgb;
    vec3 diffuse      = irradiance * albedo;
    
    // sample both the pre-filter map and the BRDF lut and combine them together as per the Split-Sum approximation to get the IBL specular part.
    const float MAX_REFLECTION_LOD = 4.0;
    vec3 prefilteredColor = textureLod(prefilterMap, R,  roughness * MAX_REFLECTION_LOD).rgb;    
    vec2 brdf  = texture(brdfLUT, vec2(max(dot(N, V), 0.0), roughness)).rg;
    vec3 specular = prefilteredColor * (F * brdf.x + brdf.y);

    vec3 ambient = (kD * diffuse + specular) * ao;
    
    vec3 color = ambient;

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

    FragColor = vec4(color , 1.0);
}

BRDF至此就告一段落了。再见！