前言

IBL，是image based lighting的简称，实质上是一些列光照技术的总和，这些技术有个共同的特点，那就是不会像传统方式那样分析光源的具体性质，而是将整个环境（环境image）视作是一个巨大的环绕光源，贴图上的每一个像素都各自代表了一个方向上的光亮度。这听起来与全局光照(GI)非常契合，而且IBL相较于传统的ambient lighting有着更加精确的输入以及更加基于物理的公式描述，因此人们一般会使用IBL技术来处理PBR渲染流程中的GI部分，特别是环境光漫反射和环境光镜面反射这两项。

环境光漫反射部分

这部分的产出反应到美术资源上既是：irradiance map

首先我们关注下面的漫反射项定义

1

其中kd近似认为是菲涅尔折射比率，小d代表diffuse。我们说fd是BRDF的漫反射公式，具体可以参考Cook-Torrence模型中的fd定义，作为BRDF分布方程，其作用在于控制朝任意给定方向的出射光亮度与入射光辉度的比例关系。

因为是漫反射模型，此时fd可以近似为均匀分布，顾写成ρ/π的形式，其中π代表各个角度，而ρ表示物体固有色（albedo），此时fd退化成了一个常量，所以可以轻松从积分项中提出。

kd项也需要先近似成与积分参数无关的量，才能提出积分项:

2

如此kd转换为一个和金属度，粗糙度，F0以及出射角相关的常量。

简化后的漫反射模型如下：

3

可以看到积分项只剩L，n和w三项，其中L是入射角为wi时的光强度，wi就是积分Li的立体角，n则是控制Ω的法向量，它控制了整个积分半球的朝向。

展开积分项（按照天顶角和方位角的方式在半球空间展开），然后使用蒙特卡洛积分转换成采样求和:

4

上式中的cosθ来自于n dot wi的结果（请忽略此前的n，属笔误）；另外我们知道立体角的微分dw = sinθdθdφ，其中θ是天顶角，φ是方位角。而转换式中第二行从积分号到求和的转变就是依据的蒙特卡洛重要性采样，分母p是在给定角度下的概率密度函数，由均匀分布的假设可以获得p = 1 / (2π * 0.5π)。

为何第一重积分中，积分上限时π/2？，这其实和问为什么p = 1 / (2π * 0.5π)中含有0.5π因子是一样的。我想因为这边考察的是半球积分，对于方位角如果定义在360°上的话，那么天顶角只要定义在0°到90°的区间上就能100%覆盖半球区域。

为何要把 （1/π）放在上述公式中？因为这样可以把除以π的操作放在预积分过程中进行，运行时就不用处理这个操作了（采样的Irradiance map里直接包含了除以π的值）

现在使用概率密度函数p是个均匀的函数，对于不同角度的概率处处相等 -> 1 / （2π * 0.5π）
我们也可以不使用这种简单的均匀分布，而用余弦函数替代，此时概率密度p = （n dot wi）/ π，带入这个p再算一遍上式可得：

5

总结下，我们在预积分（预处理）过程中，对环境半球（SkyBox）进行采样和求和，把数据存入IrrandanceMap中（Unity里就是TexCube）；运行时使用物体表面宏观法线n对这个3D的纹理进行采样（可能需要BoxProjection修正），采样的返回值乘上（kd * ρ）就是我们要的漫反射项。

origin(left), GI Diffuse(right)

下面附上Unity shaderlab格式的预积分代码，frag函数计算并返回每个IrrandanceMap上像素点的记录值。

fixed4 frag(v2f i) : SV_Target
{
    float3 irradiance = float3(0, 0, 0);
    float3 normal = normalize(i.vertex);
    float3 up = float3(0, 1, 0);
    float3 right = cross(up, normal);
    up = cross(right, normal);
    int nrSamples = 0;
    for (float phi = 0.0; phi < 2.0 * UNITY_PI; phi += _SamplerDelta) {
        for (float theta = 0.0; theta < UNITY_PI * 0.5; theta +=  _SamplerDelta) {
            // spherical to cartesian (in tangent space)
            float3 tangentSample = float3(sin(theta) * cos(phi), sin(theta)  * sin(phi), cos(theta));
            //tangent to world space
            float3 sampleVec = tangentSample.x * right + tangentSample.y *  up + tangentSample.z * normal;
            irradiance += texCUBE(_Skybox, sampleVec).rgb * sin(theta) *  cos(theta);
            nrSamples++;
        }
    }
    irradiance = UNITY_PI * irradiance / nrSamples;
    return float4(irradiance, 1);
}

环境光镜面反射部分

这部分的产出反应到美术资源上既是：prefilterMap

镜面项从标准的BRDF光照模型开始：

6

其中fs相当于BRDF，影响fs的因素和变量有很多：入射立体角wi，出射立体角wo，法线n，粗糙度，F0，金属度以及albedo等，维度较多，而一个TexCube只有3个维度，不可能一次性全部预计算并保持到这个三维纹理里，而我们要做的是：提取其中的一部分，拆分出来，预计算并存入数据结构。

那么提取那部分呢？从结构功能来看fs，n都与待渲染物体材质有关，wi是积分项，余下的Li与环境相关，与待渲染物体无关。显然我们首先应该预积分的项是与环境关联的Li，因为这项数据在实时计算时开销巨大，且没有很好的简化模型可用。

7

Lc是我们要求的累积入射光强度，将上式的fs替换为F*G*D，并使用D（法线分布函数）作为重要性采样的概率分布来源，可得如下推算：

8

说明：

1. 关于约去的项有：分子上D项被分母上概率密度函数p所带有的D分量约去了; wo dot n项因为上下两个积分都有，且与微表面无关，也被约去了。
2. 最后的约等号意思是，使用单项遮蔽G1近似替代联合遮蔽函数G。
3. 4(wo dot n) (wi dot n) 源自几何遮蔽项G，定义了从宏观到微观的缩放尺度。
4. 概率密度函数p(wi)的展开形式是：D(wh)(wh dot n) / (4wo dot wh)，其本质是关于分子上所有项DFGL或DFG总和的概率密度，之所以选择以D项的输出为主是因为它是概率分布函数p的主要塑型者，而余下的 (wh dot n) / (4wo dot wh)用于近似FG项的影响（个人理解）。
5. 额外的，上式结果中的F项是菲涅尔反射率，可以先默认为成常数1从而进一步简化公式，这主要是因为F项和物体本身以及Li无关，是个可以在实时计算出来的影响因子。

注意上式分子和分母的求和公式中，都各自包含了（wo dot wh）以及（wh dot n）部分，我们不能轻易约去这2个分量，因为它们从属于不同的重要性采样，彼此独立不关联。

但是我们知道在积分过程中，光线的出射方向wo、以及物体宏观法线方向n是与微表面固有属性无关的，如果我们让这2个参数统一成一个变量R，势必可以大大简化等式右边的变量个数。更进一步，我们不妨假设wo和n在数值上也完全一致，让它们都等于R：

9

稍微解释下，摄像机观察物体表面某个点时，该处所呈现的镜面反射图像来自于R方向确定的skybox局部图像。上述假设意味着，我们在积分过程中更加关心沿着物体表面法线方向的天空盒采样结果（因为wo == n）。

先说使用R替代wo和n的好处，那就是（wo dot wh）以及（wh dot n）的计算结果都是1，可以自我约去，给后续计算带来了极大便利。
再说使用R近似的问题：

1. 从理论方面考察：退化了BRDF公式的功能，从原本wi和wo共2个方向入参变成一个wi入参，固定了wo方向为n方向（或者说R方向），也就是说只有摄像机正视物体表面时，计算结果才正确，斜视物体时，物体表面像素的返回值等效于垂直观察物体时的计算结果。
2. 从效果方面考虑：反射变化会更趋向于镜面反射，而忽略了掠射角时形成的拖影效果（参考下图）

grazing specular reflection

不论如何，这种假设使得我们能够将环境光镜面反射预处理数据压缩金一个3D纹理中，并最大程度的保留了数据的本身特性。
来看看最终化解后的Lc项：

10

可以预见，目前为止能控制掠射角附近光强的只有G项了，我们一般采用虚幻所使用的近似权重：

11

来替代G1的功能（吐槽，这不就是Lambert模型嘛）。
回顾下完整的积分公式：

12

我们很费力的把Li从积分项中提出，变成Lc并压缩到预计算纹理中，余下的积分项则没有Li参与（或者认为很等于1），可以在运行时交给BRDF去处理，也可以另起一套预积分纹理来存储。注意最后合成结果时，两部分使用乘法算子混合。
下面给出预计算Lc时使用的着色器代码片段：

fixed4 frag(v2f i) : SV_Target
{
    const float resolution = 512.0;
    float3 N = normalize(i.normal);
    float3 R = N;
    float3 V = R;

    const uint SAMPLE_COUNT = 1024u;
    float totalWeight = 0.0;
    float3 filterColor = float3(0.0, 0.0, 0.0);
    for (uint i = 0u; i < SAMPLE_COUNT; i++) {  //采样求和 -> 积分 
        float2 Xi = Hammersley(i, SAMPLE_COUNT);
        float3 H = ImportanceSampleGGX(Xi, N, _Roughness);  //重要性采样，可以提高主要反射方向的采样率 
        float3 L = normalize(2.0 * dot(V, H) * H - V);  //光线方向，既采样方向
        float nl = max(dot(N, L), 0.0);  //n dot wi 
        if (nl > 0) {
            float mipmap = _Roughness == 0.0 ? 0.0 : CalMip(N, H, V, _Roughness, resolution);  //引入粗糙度和方向向量，计算mip等级
            filterColor += texCUBElod(_Skybox, half4(L, mipmap)).rgb * nl;  //求和公式Lc的分子部分，对Skybox采样的返回值既Li(wi)
            totalWeight += nl;   //求和公式Lc的分母部分
        }
    }
    filterColor /= totalWeight;

    return float4(filterColor, 1);
}

BRDF LUT

还是环境光镜面反射的话题，上面一部分对Lc的提取和计算做了详细说明，那么余下的另一半积分项：

13

也可以自然也可以采用预积分技术提前将结果输出到一张纹理上存储起来，这张纹理就叫LUT。

在谈论如何压缩前，我们先分析下上式中所含独立变量有哪些：代表出射方向的wo，法线n，菲尼尔系数F0（一般由albedo贴图和金属度贴图共同控制）以及粗糙度roughness。对于各向同性的D和G项来说，wo与n本身指向不重要，重要的是它们间的夹角θ，故我们合并wo与n得到最终一共3个影响因子：θ，F0和roughness。
由于F0只影响菲涅尔项F，且通常采用的Shlick Fresnel近似公式比较简单，可以方便提出F0，所以为了进一步压缩预计算内容，我们不妨将常量F0从积分公式中提取出来：

14

简单说明下：

1. fs / F代表的是将F项移除后的BRDF公式
2. Shlick Fresnel近似：F0 + (1-F0)(1-wo dot h)^5
3. 上式通过展开F项，提出了F0因子，并且将原积分式转化为2个子积分项之和：F0*scale + bias

对于上式中scale项（bias类似），我们视法线分布函数D为重要性采样的概率密度函数p(wi)，进行蒙特卡洛积分：

15

可以看到，结果中D项被约去了，只余下了G项以及去除F0的F项余部。该函数与wo与n构成的夹角θ，以及粗糙度roughness（影响G）相关，至于wi与n构成的夹角φ是我们在积分过程中的微元，是确定值。我们使用余弦函数cosθ将坐标轴规范到[0,1]区间，将scale和bias分别存入一张2D纹理的2个通道里去，所得结果就是LUT贴图！

LUT参考：

LUT

其中：

1. 红色通道存放了预计算结果的scale部分。
2. 绿色通道存放了预计算结果的bias部分。
3. 坐标轴按cosθ和roughness展开。

附上参考代码：
做了简单英文注释，结合上面推导公式，相信大家可以看懂。

fixed4 frag (v2f i) : SV_Target{
    float nv = i.uv.x；
    float roughness = i.uv.y；

    float3 V;
    V.x = sqrt(1.0 - nv * nv);  //sqrt(1-(cosθ)^2) = sinθ
    V.y = 0;
    V.z = nv;   //cosθ
    float A = 0.0;  //output scale
    float B = 0.0;  //output bias
    float3 N = float3 (0.0, 0.0, 1.0); //default normal in model space

    const uint SAMPLE_COUNT = 1024u;
    for (uint i = 0u; i < SAMPLE_COUNT; i++) {
        float2 Xi = Hammersley(i, SAMPLE_COUNT);  //sample solid angle as incident ray direction
        float3 H = ImportanceSampleGGX(Xi, N, roughness);  //apply importance sampling to H
        float3 L = normalize(2.0 * dot(V, H) * H - V); //cal light dir according to H and V
        //prepare params
        float nl = max(L.z, 0.0);
        float nh = max(H.z, 0.0);
        float vh = max(dot(V, H), 0.0);
        if (nl > 0.0)  //the incoming light should above the horizon line
        {
            float G = GeometrySmith(nv, nl, roughness);  //cal SmithG
            float G_Vis = (G * vh) / (nl * nv); //trun G to macro space by appling coefs
            float Fc = pow(1.0 - vh, 5.0);      //the core part of Shlick Fresnel
            A += (1.0 - Fc) * G_Vis;            //accumulate results
            B += Fc * G_Vis;
        }
    }
    A /= float(SAMPLE_COUNT);   //divide by N
    B /= float(SAMPLE_COUNT);

    return float4(A, B, 0, 1);
}

合并输出

考察下我们的拼图，它由GI漫反射部分和GI镜面反射部分组成，漫反射我们预积分了Li，范围是法线规定的半球空域，但是还缺少Kd * ρ 这2部分。

对于GI镜面反射，我们只需要知道F0，但是得采样2张纹理，分别是存放了环境入射光积分的prefilter map，以及环境BRDF项的LUT。

我们最后组合这两部分获得输出结果时，往往还会对其乘以环境光遮蔽贴图（AO），用以平衡因重要性采样带来的局部能量守恒问题。
具体处理流程简便起见请直接参考如下代码和注释。

fixed4 frag(v2f i) : SV_Target
{
    ...
    //GI Diffuse
    F0 = lerp(0.001, Albedo, _Metallic);  //用金属度插值出F0
    float3 F_roughness = fresnelSchlickRoughness(nv, F0, roughness);   //perceptualRoughness，考虑了粗糙度的菲涅尔项
    float kd_ibl = (1 - F_roughness) * (1 - _Metallic);   //计算Kd，Kd是一个与金属度，粗糙度，F0以及出射角相关的常量
    float3 irradiance = texCUBE(_IrradianceMap, i.normal).rgb;  //采样由i.normal方向规定的半球空域内的预积分值，代表Li的影响
    float3 indiffuse = kd_ibl * Albedo * irradiance; //Kd * ρ * irradiance 注：1/π已经包含在irradiance内了

    //GI Specular
    //首先是prefilter map
    float3 reflectVector = reflect(-viewDir, i.normal); //反射方向
    float3 prefilter_Specular = texCUBE(_PrefilterMap, reflectVector).rgb;  //采样预积分的 Lc

    //其次是env BRDF
    float2 envBRDF = tex2D(_LUT, float2(lerp(0, 0.99, nv), lerp(0, 0.99, roughness))).rg;  //采样scale + bias

    //合成GI的镜面反射部分  
    float3 inspecular = prefilter_Specular * (envBRDF.r * F0 + envBRDF.g);

    //组合Diffuse 和 Specular
    float3 indirectLight = (indiffuse + inspecular) * _AO;
    ...
}

参考

[1] 深入理解 PBR/基于图像照明 (IBL)
[2] IBL推导及实现
[3] 深入浅出基于物理的渲染二
[4]Diffuse-irradiance