重要性采样

1. 介绍

首先，明确的是，本人在这里不会涉及大数定理、蒙特卡洛积分、重要性采样的数学原理和推导，仅仅是做个介绍。重点还是重要性采样和光照积分的结合。

1.1 快速回顾

蒙特卡洛积分

蒙特卡洛积分，并不是指一种名叫蒙特卡洛的积分，而是采用蒙特卡洛法来估计积分。

参考上图，在$[a,b]$的区间内均匀取足够多的点（N很大），来估算上诉曲线下的面积（假如我们无法求解析解）：
$$F_N={\int }_a^{b}f(x)dx=\frac{b-a}{N}\sum _{i=1}^{N}f(x_i)$$

上诉方程是使用均匀采样的特殊形式，对于下面这种平滑的曲线，确实效果还行，如下图所示：

但如果这个曲线稍微弯曲一点，不给你老老实实呢？误差会急剧增大：

怎么办呢？又不能增加采样数，但仔细看看上图，就会发现——误差大多发生在函数值比较大的地方（陡峭的地方）。如果我们这样采样：

上图有个缺点：就是函数值大的地方，应该分配三个采样，这样更加直观，如下图：

误差明显小了很多。那么表述成文字呢：对于那些值比较大的区域（陡峭的区域），给这些采样更小的权重——以上图为例，就是长方形的宽变小了。转化为数学形式，就是下面：
$$F_N=\frac{1}{N}\sum _{i=1}^N\frac{f(x_i)}{p(x_i)}$$

重要性采样

随着采样数的增加，使用蒙特卡洛积分估算得到的结果越精确，但高采样数对于游戏而言是无法接受的。那么有什么方法能够在一定的抽样数量基础上来增加准确度，减少方差呢？

一个立马可以想到的思路就是：在将采样尽量分布给那些贡献较大的区域。例如，下图中的圆形部分应该分配更多的采样。

所以以概率密度函数来说，这个函数应该在贡献大的区域也具有较大的值（这说明采样频率高）。本质上，蒙特卡洛积分和重要性采样在数学公式上应该没有区别，只是采样方式不一样：一个是均匀采样（或随机采样）；一个是特定采样（集中在贡献大的区域）。

但大部分情况下，我们无法应用目标函数的采样 $\pi (x_i)$，所以需要引入可以进行应用的、更加简洁的采样分布 $p(x_i)$。然后，数学形式，推导每次看完就忘，这里就直接给出了：

如果，我们可以直接在$\pi (x_i)上采样，那么重要性采样的公式就是之前蒙特卡洛积分的形式$

1.2 图形学如何使用呢？

直接看上诉分析，会云里雾里。这里结合图形学谈一下：

• 我们在求解间接光的时候，需要在着色点的半球域进行采样积分。直接均匀采样求解，需要足够多的采样数来保证渲染效果的真实性，这在实时是不行的。

• Diffuse光由于其低频性，我们或许可以在采样数量降低的情况下，保持渲染质量。但Specular光则不行，因为它的主要贡献都在一个特定的lobe里面——而不像diffuse在整个半球域的贡献是均匀的。

• 参考之前的分析：在保持低采样数的情况下，我们为了更好的拟合Specular，应该尽可能采样lobe区域——而决定这个lobe的，就是我们所熟悉的BRDF，这个时候，积分离散形式是这样的：
  $$L_o=\frac{1}{N}\sum _{i=1}^N\frac{BRDF\cdot L_i\cdot \cos }{BRDF}$$

• 那为什么，我们见到的都是GGX重要性采样，而不是上诉的BRDF重要性采样呢？原因不在公式的计算式，而是采样i：这里的N个采样都是根据BRDF造成的lobe分布来生成的，但我们很难使用复杂的BRDF公式来生成采样方向。所以，上诉公式无法使用的原因是：我们很难生成对应的采样点。

• 而GGX法线分布项，不仅是构成Specular lobe的主要原因，而且还很容易生成对应的采样（相比于BRDF），所以我们一般都是使用GGX重要性采样。

2. 常用的重要性采样方法

2.1 GGX重要性采样

基本推导

1️⃣根据之前的分析，我们可以很快知道：$f(x_i)=L_i\cdot f_r\cdot \cos (n,w_i)$ 就是我们的光照方程，${\pi }_i(x)=f_r\ast \cos (n,L_i)$ 是光照情况的原分布 BRDF，$p(x_i)$就是实际采样分布的概率密度函数。然后我们需要做什么呢？

1. 选择GGX NDF作为我们的概率分布函数，使用它来生成采样方向（入射光的采样方向）。
2. 计算对应的权重函数。
3. 计算光照方程，按权重累加。

2️⃣首先，GGX NDF为：

其中，$\theta =<n,h>$。为了密切匹配镜面BRDF的形状，根据分布函数NDF和 $\cos (n,h)$ 建立PDF（也就是$p(x_i)$）：
$$pdf(h)=D(h)\ast \cos {\theta }_h$$

为什么最终的概率分布函数要乘上余弦项？可以看看：https://www.reedbeta.com/blog/hows-the-ndf-really-defined/。

在求解这个$pd{f}_h$之前，我们需要知道：因为$pd{f}_h$是在半矢量空间定义的。然而，积分步骤需要对特定视线方向的光向进行积分。因此，我们需要将$pd{f}_h$从半矢量空间转换成光空间（从$pd{f}_h$转换成$pd{f}_L$）。这种PDF的转换是简单的：
$$pd{f}_L=pd{f}_h\cdot \frac{\delta w_h}{\delta w_L}=\frac{pd{f}_h}{4\cos (v,h)}$$

原论文是求输出的pdf，而这里是求输入的，所以，分母的i要变成v，这是否解释了上面的PDF转换？

因为我们可以用球面坐标 $(\varphi ,\theta )$ 来表示半向量h，所以我们可以把$pd{f}_h$表示为 $pd{f}_{\varphi }$ 和 $pd{f}_{\theta }$​ 的乘法：
$$pd{f}_h={pdf}_{\varphi }\times {pdf}_{\theta }$$

$pd{f}_h$不依赖于角度$\varphi$，所以我们可以简单地推导出：
$$pd{f}_{\varphi }=\frac{1}{2\pi }$$
所以：
$$pd{f}_{\theta }=\frac{D(h)\times \cos ({\theta }_h)}{pd{f}_{\varphi }}=\frac{2\times {\alpha }^2\times \cos ({\theta }_h)}{(\cos {({\theta }_h)}^2\times ({\alpha }^2-1)+1{)}^2}$$

生成采样

3️⃣但我们不能直接使用pdf来生成采样方向，而是要使用对应的反转累积分布函数CDF_Inv来生成——这里的数学含义，我是一直无法理解透彻。

参考上图，进行理解。

所以，我们直接给出结果：
$$\begin{gathered} cd{f}_{\varphi }={\int }_0^{\varphi }\frac{1}{2\pi }dx=\frac{1}{2\pi }\varphi \\ cd{f}_{\theta }={\int }_0^q\frac{2\times {\alpha }^2\times x}{(x^2\times ({\alpha }^2-1)+1{)}^2}dx=\frac{1-q^2}{1+q^2(a^2-1)} \end{gathered}$$
其中，$q=cos({\theta }_h)$。现在反转CDF函数，以产生从均匀值${\epsilon }_1,{\epsilon }_2$，到角度 $\varphi ,\theta$​​​ 的映射：
$$\begin{gathered} \varphi ={\epsilon }_1\times 2\pi \\ \theta =\arccos \sqrt{\frac{1-{\epsilon }_2}{1+{\epsilon }_2(a^2-1)}} \end{gathered}$$

来自LearnOpenGL的源码就是复现上面的公式，并做了必要的坐标变换（球坐标到笛卡尔坐标，法线空间到世界空间）：

vec2 Hammersley(uint i, uint N)
{
    return vec2(float(i)/float(N), RadicalInverse_VdC(i));
}  
...
const uint SAMPLE_COUNT = 4096u;
for(uint i = 0u; i < SAMPLE_COUNT; ++i)
{
    vec2 Xi = Hammersley(i, SAMPLE_COUNT);   
}
...
vec4 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);

    // from spherical coordinates to cartesian coordinates
    vec3 H;
    H.x = cos(phi) * sinTheta;
    H.y = sin(phi) * sinTheta;
    H.z = cosTheta;
	
	// 计算pdf_h : D * cos
	float d = (cosTheta * a - cosTheta) * cosTheta + 1;
	float D = a / (PI * d * d);
	float PDF = D * cosTheta;

    // from tangent-space vector to world-space sample vector
    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 vec4(normalize(sampleVec), PDF);
}

权重和累加光照

我们在生成采样的函数里面顺便计算了$pd{f}_h$，然后在主函数里面要将产生的半向量h转化成入射光:

vec4 re  = ImportanceSampleGGX(Xi, N, roughness);
vec3 H = re.xyz;
vec3 L  = normalize(2.0 * dot(V, H) * H - V);

此时，我们也需要，将 $pd{f}_h$ 转换成 $pd{f}_L$ ，这个转换公式之前已经给了，所以直接计算：

// pdf_L
float pdf = re.z / (4.0 * VoH);

而真正的权重，不仅和我们预设的分布 $pd{f}_L$有关，还和原分布 $\pi (w_i)=f_r\cdot \cos (n,w_i)$有关：
$$weigh{t}_i=\frac{f_r\cdot \cos (n,w_i)}{pd{f}_L}=\frac{F\cdot G\cdot \cos (v,h)}{\cos (n,v)\cdot \cos (n,h)}$$
UE4的源码就是这样做的：

if( NoL > 0 )
{
	float3 SampleColor = AmbientCubemap.SampleLevel( AmbientCubemapSampler, L, 0 ).rgb;

	float Vis = Vis_SmithJointApprox( Pow4(Roughness), NoV, NoL );
	float Fc = pow( 1 - VoH, 5 );
	float3 F = (1 - Fc) * SpecularColor + Fc;

	// Incident light = SampleColor * NoL
	// Microfacet specular = D*G*F / (4*NoL*NoV) = D*Vis*F
	// pdf = D * NoH / (4 * VoH)
	SpecularLighting += SampleColor * F * ( NoL * Vis * (4 * VoH / NoH) );
}
return SpecularLighting / NumSamples;

关于NoL权重的解释

在LearnOpenGL的预卷积环境贴图中，使用的权重是 $NoL$，如何解释？Todo。

Mip map的选择

2.2 Blinn重要性采样算法

来自知乎博客：图形程序的点点滴滴

Blinn的NDF定义如下：

同样的，实际的概率密度函数需要乘上一个余弦：

为了方便计算，将其转换到球面上来，我们可以得到其在球面坐标下的概率密度函数：

由于这个概率密度函数有2个参数，我们需要分别计算边缘概率密度函数：

计算累计概率函数CDF：

现在反转CDF函数：

权重计算和GGX类似。

2.3 基于亮度的环境贴图重要性采样

Todo.