3.4、漫反射项的重要性采样计算

3.4.1、漫反射项的二重积分形式极其近似预积分计算

  进一步来说，因为对漫反射项积分来说 ${\kappa }_d、c、\pi$ 等参数对于指定的一点 $\vec{p}$ 来说是常数，可以提到积分外面，所以其中第一个积分可以进一步推导为：
$$\begin{gathered} {\mathrm{L}}_{o_d}(\vec{p},\overrightarrow{{\omega }_o})={\kappa }_d\frac{c}{\pi }{\int }_{\Omega }{\mathrm{L}}_i(\vec{p},\overrightarrow{{\omega }_i})\vec{n}\cdot \overrightarrow{{\omega }_i}d\overrightarrow{{\omega }_i} \\ \because d\omega =\sin (\theta )d\theta d\varphi ,\vec{n}\cdot \overrightarrow{{\omega }_i}=\cos (\theta ) \\ ={\kappa }_d\frac{c}{\pi }\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\frac{\pi }{2}}{\int }}{\mathrm{L}}_i(\vec{p},\overrightarrow{{\omega }_i})\cos (\theta )\sin (\theta )d\theta d\varphi \end{gathered}$$
  上式中，角度 $\theta$ 就是入射光线与点 $\vec{p}$ 处法线 $\vec{n}$ 的夹角也就是入射角。同时在点 $\vec{p}$ 处作直角坐标系如下图所示，其中正 $z$ 方向与点 $\vec{p}$ 处法线 $\vec{n}$ 方向一致，$x$ 轴与 $y$ 轴在与点 $\vec{p}$ 处法线 $\vec{n}$ 垂直的平面内。

  接着取入射光 ${\omega }_i$ 与点 $\vec{p}$ 处法线 $\vec{n}$ 所成平面与 $x$ 轴间的夹角为 $\varphi$ （称之为方位角），然后计算得到入射立体角微分 $d\omega =\sin (\theta )d\omega d\varphi$ 最终就得到了上式中漫反射项的最终可计算的二重积分形式。

  需要注意的是，我们现在使用的模型是针对不透明物体的表面的反射，所以积分的区域是点 $\vec{p}$ 处，朝向法线方向的单位正半球，所以方位角 $\varphi$ 的范围是 $[0,2\pi ]$，对应的天顶角 $\theta$ 范围就是 $[0,\frac{\pi }{2}]$ 。

  这个积分项在 IBL 方法中称之为“漫反射项”。注意到其中积分部分已经完全与 $\overrightarrow{{\omega }_o}$ 也即视方向无关了，所以其实这一项是可以从每一帧循环渲染计算中完全提取出来的。当然提取的前提是被照射的物体与其环境在某段固定的时间中基本不变，或者说针对一个相对固定的场景来说。

  具体来讲，对于一个固定点 $\vec{p}$ 来说，整个式子又可以被拆分成两部分去计算，一部分是针对点 $\vec{p}$ 的常数项 ${\kappa }_d\frac{c}{\pi }$ ，这个可以被延迟到渲染循环中的物体表面 PS 中去计算，因为参数 $c$ 就是物体表面点的漫反射率，也就是我们通常所说的物体表面的颜色，而 ${\kappa }_d$ 就是控制整体能量守恒的系数，需要先计算点 $\vec{p}$ 处的 ${\kappa }_s$ 也就是菲涅尔系数后再确定，所以整个常数项实际是与物体表面相关联的，因此就放到最终渲染循环中去。

  而第二部分积分项，其实本质上也是与物体表面的点 $\vec{p}$ 相关联的，但如果也因此而放到最终渲染循环中的话，就会导致计算量过大，而无法实现实时的 PBR 渲染。因此在这里第一个关于质量和性能的折中处理，就是我们总是假设这个积分中的点 $\vec{p}$ 就是物体的几何中心点所在的位置，然后将整个积分的计算从渲染循环中提取出来，变成一个全局预处理的计算过程，而这个预积分计算的结果就被称之为“预积分辐照度贴图”，因为整个预计算部分正好就是针对所有方向的入射光也就是入射辐照度进行积分。

  实际应用时，这需要假设我们可以提前预知点 $\vec{p}$ 的位置。对于固定的场景中的固定的物体位置来说，这本来就是已知量。这就使得预积分辐照度贴图最终完全可以在场景设计确定的阶段就计算完毕，并存储起来，在真正运行渲染时进行采样使用即可。

  在本章示例中我们只是用了最简单的 Cubemap 纹理（或者叫2D纹理数组 Texture2DArray）的方式，对这个辐照图贴图进行了存储，这样的存储方式还是比较耗费空间的。更高级的可以使用球谐函数的方式简化存储。

  在示例代码运行后的矩形框中，紧跟解算后的6副环境映射纹理之后的6副纹理就展示了预积分贴图的六个面：

  如果场景中的物体是经常变动位置时，而且大多数情况下，我们渲染的核心很多时候就是3D动画的场景，这时点 $\vec{p}$ 的位置就会经常变动，就使的这个预积分计算跟渲染循环，或者说跟具体的渲染帧发生了关联，因此刚才说的假设点 $\vec{p}$ 是某物体中心点的近似预积分方法就可能失效了。

  此时可以在场景中的某些按照一定规律固定排列的点 ${\vec{p}}_n$ 来作为预积分计算的点，然后在渲染循环中，根据物体位置离哪个点 ${\vec{p}}_n$ 近，就取那个点处的预积分贴图来进行漫反射项的近似计算。这样就可以在场景设计定型的时候，预计算一堆关于点 ${\vec{p}}_n$ 的漫反射辐照度预积分贴图。而这些点 ${\vec{p}}_n$ 就被称为 Probe （探针）。更进一步的，在探针技术的基础上，可以使用球谐函数等来简化预积分辐照度贴图计算结果的存储和使用。从而变成一个一次性计算（或称之为预计算）！这些话题后续教程再说。

3.4.2、漫反射辐照度积分项的直接积分计算

  至此其实漫反射项的二重积分形式已经可以直接编码进行计算了，在示例代码 Shader 文件： GRSD3D12Sample/GRS_IBL_Diffuse_Irradiance_Convolution_With_Integration_PS.hlsl 中就展示了直接二重积分计算的方法：

float4 PSMain(ST_GRS_HLSL_PS_INPUT pin):SV_Target
{
	float3 N = normalize(pin.m_v4WPos.xyz);

	float3 irradiance = float3(0.0f, 0.0f, 0.0f);

	float3 up = float3(0.0f, 1.0f, 0.0f);
	float3 right = normalize(cross(up, N));
	up = normalize(cross(N, right));
	N = normalize(cross(right,up));

	float nrSamples = 0.0f;

	float deltaPhi = (2.0f * PI) / 180.0f;
	float deltaTheta = (0.5f * PI) / 90.0f;

	// 积分运算，直接翻译自公式
	for (float phi = 0.0f; phi < 2.0f * PI; phi += deltaPhi)
	{
		for (float theta = 0.0f; theta < 0.5f * PI; theta += deltaTheta)
		{
			// 球坐标转换到笛卡尔坐标（切空间）
			float3 tangentSample = float3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));

			// 切空间转换到世界坐标空间
			float3 sampleVec = tangentSample.x * right + tangentSample.y * up + tangentSample.z * N;

			//float3 sampleVec = right;
			irradiance += (g_texHDREnvCubemap.Sample(g_sapLinear, sampleVec).xyz * cos(theta) * sin(theta));

			nrSamples += 1.0f;
		}
	}

	irradiance = PI * irradiance * (1.0 / float(2 * nrSamples));
	return float4(irradiance, 1.0f);
}

  上面代码中变量 irradiance 代表最终点 $\vec{p}$ 处的所有入射光引起的入射辐照度，所以这个过程最终计算出来的预积分纹理也被称为“辐照度贴图”。

  代码中的 deltaPhi 即 $d\varphi$，根据代码中的表达式，可以看出它的步长（step）是 $\frac{2\pi }{180}$弧度，也就是总共计算180次、deltaTheta 即 $d\theta$ ，对应的其步长是 $\frac{\frac{\pi }{2}}{90}$弧度，总共计算 90 次。最终这使得双重循环黎曼和计算对于一个采样点 $\vec{p}$ 来说总共需要 $180\times 90=16200$ 次！然而这还只是一个很粗糙的颗粒度，这从其步长设置就可以看出来。并且实际的效果也是惨不忍睹的。

  而总体上我们总共有 $6\times 32\times 32=6144$ 个点，这只是低分辨率下的辐照度贴图映射的大小。在示例代码开头处我们用宏定义指定了这个大小：

#define GRS_RTA_IBL_DIFFUSE_MAP_SIZE_W  32
#define GRS_RTA_IBL_DIFFUSE_MAP_SIZE_H  32

  接着代码中主要就是用点 $\vec{p}$ 的方向用作其自身的法线，并假设向上的方向为正 y 轴方向，然后用施密特正交法计算得到了点 $\vec{p}$ 处局部切空间的三个正交坐标轴，也就是得到了切空间转换到世界坐标空间的变换矩阵，这也是一个非常重要和基础的切空间计算的简单技巧。然后用这个切空间的变换矩阵，将点 $\vec{p}$ 处局部坐标系中的立体角微分向量变换到了世界坐标中，并作为 3D 采样（Sample）操作的向量对环境映射进行采样。

  最终再一次提醒，在没有完全把握的情况下，请千万不要使用和运行GRSD3D12Sample/GRS_IBL_Diffuse_Irradiance_Convolution_With_Integration_PS.hlsl 这个 Shader 脚本，以免因计算量过大造成显卡损伤。放在本章示例中，仅仅是为了说明预计算辐照度贴图计算过程的基本原理。或者说是为了说明下基本的数值积分代码如何编写而已。

3.4.3、漫反射辐照度积分项的蒙特卡洛积分重要性采样计算

  了解完了基本的预积分辐照度贴图的计算过程之后，希望你也看了 3D数学系列之——再谈特卡洛积分和重要性采样 并明白了其中介绍的重要性采样方法。索性的是，蒙特卡洛方法可以使我们极大的减小预积分辐照度贴图的计算量，并且还能保持一定的精度。

  根据文中的公式：
$$\mathrm{p}(x)=\frac{|f(x)|}{\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x}$$
  我们继续使用瞪眼法来观察漫反射项中的积分部分：
$$\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\frac{\pi }{2}}{\int }}{\mathrm{L}}_i(\vec{p},\overrightarrow{{\omega }_i})\cos (\theta )\sin (\theta )d\theta d\varphi$$
  虽然这是个二重积分，但是根据蒙特卡洛积分的特性，这也是可以用蒙特卡洛积分来处理的。只是需要我们选取合适的概率密度函数 $\mathrm{p}$ 。通过观察我们可以发现，使用如下的概率密度函数，就可以进行重要性采样计算：
$$\mathrm{p}(\theta ,\varphi )=c\cos (\theta )\sin (\varphi )$$
  首先要注意的是，因为被积函数是关于 $\theta ,\varphi$ 的二元函数（点 $\vec{p}$ 参数相对固定对于整个积分来说相当于常数，而入射光向量 ${\vec{\omega }}_i$ 及立体角微分都是可以用球坐标系下的 $\theta ,\varphi$ 表达的），所以我们设计的概率密度函数也是关于相同自变量的二元函数。

  另外在函数表达式中只是提取了三角函数的因子，这主要是因为一般函数表达式中如果乘入了三角函数因子尤其是正弦或余弦函数时，函数的形状就主要会受到他们的显著影响，所以我们选择如上的概率密度函数，作为我们进一步简化积分计算的切入点。

  根据重要性采样的方法步骤，第一步就是要验证概率密度函数的“归一性”，从而还可以顺道计算出式中常数项 $c$ 的值，这可以根据概率密度函数的性质要求来计算：
$$\begin{gathered} \mathrm{p}(\theta ,\varphi )=\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\frac{\pi }{2}}{\int }}c\cos (\theta )\sin (\theta )d\theta d\varphi =1 \\ \Rightarrow c\underset{0}{\overset{2\pi }{\int }}\mathrm{d}\varphi \underset{0}{\overset{\frac{\pi }{2}}{\int }}\cos (\theta )\sin (\theta )\mathrm{d}\theta =1 \\ \because \mathrm{d}\sin (\theta )=\cos (\theta )\mathrm{d}\theta \\ c\underset{0}{\overset{2\pi }{\int }}\mathrm{d}\varphi \underset{0}{\overset{\frac{\pi }{2}}{\int }}\sin (\theta )\mathrm{d}\sin (\theta )=1 \\ \Rightarrow c\left[\varphi {]}_0^{2\pi }\right[\frac{1}{2}{\sin }^2(\theta ){]}_0^{\frac{\pi }{2}}=1 \\ \Rightarrow c=\frac{1}{\pi } \\ \Rightarrow \mathrm{p}(\theta ,\varphi )=\frac{1}{\pi }\cos (\theta )\sin (\theta ) \end{gathered}$$

  因为这个概率密度函数是二元函数，而我们产生随机数需要一个变量一个变量的产生，这就需要进一步针对每个随机变量计算其边缘概率密度，从而对随机变量进行分离，此时推算如下：
$$\begin{gathered} \because \mathrm{p}(\theta ,\varphi )=\frac{1}{\pi }\cos (\theta )\sin (\theta ) \\ \therefore \mathrm{p}(\theta ,\varphi )=\mathrm{p}(\theta )\times \mathrm{p}(\varphi ) \end{gathered}$$

  上面过程说明我们要使用的随机变量 $\theta ,\varphi$ 是相互独立的，其实从其定义也可以看出，这两个角度变量是没有关联性的。所以就可以继续推导每一个独立随机变量的概率密度函数：
$$\begin{gathered} \mathrm{p}(\varphi )=\underset{0}{\overset{\frac{\pi }{2}}{\int }}\mathrm{p}(\theta ,\varphi )\mathrm{d}\theta \\ =\underset{0}{\overset{\frac{\pi }{2}}{\int }}\frac{1}{\pi }\cos (\theta )\sin (\theta )\mathrm{d}\theta \\ =[\frac{1}{2\pi }{\sin }^2(\theta ){]}_0^{\frac{\pi }{2}} \\ =\frac{1}{2\pi } \end{gathered}$$

$$\begin{gathered} \because \mathrm{p}(\theta ,\varphi )=\mathrm{p}(\theta )\times \mathrm{p}(\varphi ) \\ \frac{1}{\pi }\cos (\theta )\sin (\theta )=\mathrm{p}(\theta )\times \frac{1}{2\pi } \\ \Rightarrow \mathrm{p}(\theta )=2\cos (\theta )\sin (\theta ) \end{gathered}$$

  至此，我们就分析推导清楚了我们需要的随机变量的概率密度函数，接着按照重要性采样的步骤，就需要计算每个随机变量对应的概率分布函数 CDF :
$$\begin{gathered} \mathrm{CDF}(\theta )=\underset{0}{\overset{\theta }{\int }}2\cos (x)\sin (x)\mathrm{d}x \\ =[2\times \frac{1}{2}{\sin }^2(x){]}_0^{\theta } \\ ={\sin }^2(\theta ) \\ \mathrm{CDF}(\varphi )=\underset{0}{\overset{\varphi }{\int }}\frac{1}{2\pi }\mathrm{d}x \\ =\frac{1}{2\pi }\varphi \end{gathered}$$
  接着继续计算 CDF 的反函数如下：
$$\begin{gathered} 令\mu ={\sin }^2(\theta ),\nu =\frac{1}{2\pi }\varphi 则： \\ \theta =\arcsin (\sqrt{\mu }) \\ \varphi =2\pi \nu \end{gathered}$$
   最后我们只需要产生两个均匀分布的伪随机序列 $\mu ,\nu$ 就可以利用蒙特卡洛积分公式计算辐照度积分了。那么让我们将所有推导的结果反代回蒙特卡洛积分先做一下简化：
$$\begin{gathered} \because \underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx \frac{1}{N}\sum _{n=1}^N\frac{f(x_n)}{\mathrm{p}(x_n)} \\ {\mathrm{L}}_{o_d}(\vec{p},\overrightarrow{{\omega }_o})={\kappa }_d\frac{c}{\pi }\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\frac{\pi }{2}}{\int }}{\mathrm{L}}_i(\vec{p},\overrightarrow{{\omega }_i})\cos (\theta )\sin (\theta )d\theta d\varphi \\ \because \mathrm{p}(\theta ,\varphi )=\frac{1}{\pi }\cos (\theta )\sin (\theta ) \\ \therefore {\mathrm{L}}_{o_d}(\vec{p},\overrightarrow{{\omega }_o})\approx {\kappa }_d\frac{c}{\pi }\frac{1}{N}\sum _{n=1}^N\frac{{\mathrm{L}}_i(\vec{p},\overrightarrow{{\omega }_i})\cos (\theta )\sin (\theta )}{\frac{1}{\pi }\cos (\theta )\sin (\theta )} \\ ={\kappa }_{d}c\frac{1}{N}\sum _{n=1}^N{\mathrm{L}}_i(\vec{p},\overrightarrow{{\omega }_i}) \end{gathered}$$
  从上面的最终推导过程我们可以看出，首先在形式上比原先的二重积分形式简单了很多，并且计算过程直接变成了 1 维求和计算，同时还自动消去了一个除以 $\pi$ 的计算；其次求和部分只剩下入射光函数 ${\mathrm{L}}_i(\vec{p},\overrightarrow{{\omega }_i})$ ，这在 IBL 中仅仅只是一个 3D Sample 操作，最终整个求和式就变成了针对环境映射的一个采样平均值。

  具体的在本章示例代码的 Shader GRSD3D12Sample/GRS_IBL_Diffuse_Irradiance_Convolution_With_Monte_Carlo_PS.hlsl 中实现如下：

float4 PSMain(ST_GRS_HLSL_PS_INPUT stPSInput) :SV_Target
{
    float3 N = normalize(stPSInput.m_v4WPos.xyz);

    float3 irradiance = float3(0.0f,0.0f,0.0f);

    // 计算切空间
    float3 Up = float3(0.0f, 1.0f, 0.0f);
    float fdot = dot(Up, N);
    if ((1.0 - fdot) > 1e-6)
    {
        Up = float3(1.0f, 0.0f, 0.0f);
    }

    float3 Right = normalize(cross(Up, N));
    Up = normalize(cross(N, Right));
    N = normalize(cross(Right, Up));


    float theta = 0.0f;
    float phi = 0.0f;
    float2 Xi;
    float3 L;

    const uint SAMPLE_NUM = 4096u;
    // 蒙特卡洛积分循环
    for (uint i = 0u; i < SAMPLE_NUM; i++)
    {
        // 产生均匀分布随机数
        Xi = Hammersley(i, SAMPLE_NUM);

        // 根据 CDF 反函数计算指定分布的随机变量
        phi = TWO_PI * Xi.x;
        theta = asin(sqrt(Xi.y));

        // 球坐标转换到笛卡尔坐标（切空间）
        L = float3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
        // 切空间转换到世界坐标空间
        L = L.x * Right + L.y * Up + L.z * N;

        // 采样求和统计入射辐照度
        irradiance += g_texHDREnvCubemap.Sample(g_sapLinear, L).xyz;
    }
    // 取平均得到最终结果
    irradiance *= 1.0 / float(SAMPLE_NUM);

    return float4(irradiance,1.0f);
}

  代码中的过程及注释已经很清晰的说明了整个计算过程，与我们做的数学推导过程是一致的。所以最终只要理解了整个数学推算的过程，最终的代码就变得很好理解了。

  但其中还是需要注意几个细节：首先，代码中生成均匀分布随机数的方法使用了流行的 Hammersley 算法，关于这个算法，回头在另外篇章的教程中单独介绍，本教程中就不在赘述了，现在只需要知道它能够产生高质量的 2D 均匀分布伪随机序列即可，并且范围在$[0,1]$ 之间。

  其次，整个蒙特卡洛积分循环的次数设置成了 4096 次，其实大概 2000多次的时候效果就已经很好了。但即就是使用4096次循环，也比之前的二重积分双循环运行16200次的直译算法效果好的很多，并且计算次数直接节省了到了四分之一以上。假如要使用直接数值积分算法效果达到重要性采样方法的级别，估计循环次数还需要再增进至少4倍以上，当然为了防止显卡损伤，所以还是不要去做这样的尝试了（主要是因为循环中嵌套着耗时并且费硬件的 Sample 操作）。

  最后，代码中为了做教程的目的，所以对 $\theta$ 角的变量做了反三角函数的运算，其实实际工程代码中，不需要做这个操作，直接将开方后的随机值作为 $\theta$ 角的正弦值即可，余弦值根据三角公式即可算得，这样就节省了反复计算三角函数的开销。而大多数实际的代码中确实是这样的做的，比如 UE 引擎的 Shader 代码中。

  额外需要提醒的是在著名的 [漫反射辐照 - LearnOpenGL CN](https://learnopengl-cn.github.io/07 PBR/03 IBL/01 Diffuse irradiance/) 教程中，既没有详细讲述整个数学过程，也没有推导重要性采样的过程，最后其漫反射辐照度贴图使用了计算量更大的双重循环来计算，如果不加思索的直接使用其代码，可能会造成显卡过热保护而退出，严重的情形下可能会造成显卡硬伤，所以在没有搞清楚的情况下请不要随意执行。而正确的做法是使用本章教程中所讨论的重要性采样方法来计算漫反射辐照度贴图。