1. 介绍
首先,明确的是,本人在这里不会涉及大数定理、蒙特卡洛积分、重要性采样的数学原理和推导,仅仅是做个介绍。重点还是重要性采样和光照积分的结合。
1.1 快速回顾
蒙特卡洛积分
蒙特卡洛积分,并不是指一种名叫蒙特卡洛的积分,而是采用蒙特卡洛法来估计积分。
参考上图,在
[
a
,
b
]
[a,b]
[a,b]的区间内均匀取足够多的点(N很大),来估算上诉曲线下的面积(假如我们无法求解析解):
F
N
=
∫
a
b
f
(
x
)
d
x
=
b
−
a
N
∑
i
=
1
N
f
(
x
i
)
F_N=\int_a^b{f(x)dx}=\frac{b-a}{N}\sum^N_{i=1}{f(x_i)}
FN=∫abf(x)dx=Nb−ai=1∑Nf(xi)
上诉方程是使用均匀采样的特殊形式,对于下面这种平滑的曲线,确实效果还行,如下图所示:
但如果这个曲线稍微弯曲一点,不给你老老实实呢?误差会急剧增大:
怎么办呢?又不能增加采样数,但仔细看看上图,就会发现——误差大多发生在函数值比较大的地方(陡峭的地方)。如果我们这样采样:
上图有个缺点:就是函数值大的地方,应该分配三个采样,这样更加直观,如下图:
误差明显小了很多。那么表述成文字呢:对于那些值比较大的区域(陡峭的区域),给这些采样更小的权重——以上图为例,就是长方形的宽变小了。转化为数学形式,就是下面:
F
N
=
1
N
∑
i
=
1
N
f
(
x
i
)
p
(
x
i
)
F_N=\frac{1}{N}\sum^N_{i=1}{\frac{f(x_i)}{p(x_i)}}
FN=N1i=1∑Np(xi)f(xi)
重要性采样
随着采样数的增加,使用蒙特卡洛积分估算得到的结果越精确,但高采样数对于游戏而言是无法接受的。那么有什么方法能够在一定的抽样数量基础上来增加准确度,减少方差呢?
一个立马可以想到的思路就是:在将采样尽量分布给那些贡献较大的区域。例如,下图中的圆形部分应该分配更多的采样。
所以以概率密度函数来说,这个函数应该在贡献大的区域也具有较大的值(这说明采样频率高)。本质上,蒙特卡洛积分和重要性采样在数学公式上应该没有区别,只是采样方式不一样:一个是均匀采样(或随机采样);一个是特定采样(集中在贡献大的区域)。
但大部分情况下,我们无法应用目标函数的采样
π
(
x
i
)
\pi(x_i)
π(xi),所以需要引入可以进行应用的、更加简洁的采样分布
p
(
x
i
)
p(x_i)
p(xi)。然后,数学形式,推导每次看完就忘,这里就直接给出了:
如果,我们可以直接在 π ( x i ) 上 采 样 , 那 么 重 要 性 采 样 的 公 式 就 是 之 前 蒙 特 卡 洛 积 分 的 形 式 \pi(x_i)上采样,那么重要性采样的公式就是之前蒙特卡洛积分的形式 π(xi)上采样,那么重要性采样的公式就是之前蒙特卡洛积分的形式
1.2 图形学如何使用呢?
直接看上诉分析,会云里雾里。这里结合图形学谈一下:
-
我们在求解间接光的时候,需要在着色点的半球域进行采样积分。直接均匀采样求解,需要足够多的采样数来保证渲染效果的真实性,这在实时是不行的。
-
Diffuse光由于其低频性,我们或许可以在采样数量降低的情况下,保持渲染质量。但Specular光则不行,因为它的主要贡献都在一个特定的
lobe
里面——而不像diffuse
在整个半球域的贡献是均匀的。 -
参考之前的分析:在保持低采样数的情况下,我们为了更好的拟合
Specular
,应该尽可能采样lobe
区域——而决定这个lobe
的,就是我们所熟悉的BRDF
,这个时候,积分离散形式是这样的:
L o = 1 N ∑ i = 1 N B R D F ⋅ L i ⋅ cos B R D F L_o=\frac{1}{N}\sum_{i=1}^{N}\frac{BRDF \cdot L_i \cdot \cos}{BRDF} Lo=N1i=1∑NBRDFBRDF⋅Li⋅cos -
那为什么,我们见到的都是GGX重要性采样,而不是上诉的BRDF重要性采样呢?原因不在公式的计算式,而是采样
i
:这里的N
个采样都是根据BRDF
造成的lobe
分布来生成的,但我们很难使用复杂的BRDF公式来生成采样方向。所以,上诉公式无法使用的原因是:我们很难生成对应的采样点。 -
而
GGX
法线分布项,不仅是构成Specular lobe
的主要原因,而且还很容易生成对应的采样(相比于BRDF
),所以我们一般都是使用GGX重要性采样。
2. 常用的重要性采样方法
2.1 GGX重要性采样
基本推导
1️⃣根据之前的分析,我们可以很快知道:
f
(
x
i
)
=
L
i
⋅
f
r
⋅
cos
(
n
,
w
i
)
f(x_i)=L_i\cdot f_r \cdot \cos{(n,w_i)}
f(xi)=Li⋅fr⋅cos(n,wi) 就是我们的光照方程,
π
i
(
x
)
=
f
r
∗
cos
(
n
,
L
i
)
\pi_i(x)=f_r*\cos{(n,L_i)}
πi(x)=fr∗cos(n,Li) 是光照情况的原分布 BRDF
,
p
(
x
i
)
p(x_i)
p(xi)就是实际采样分布的概率密度函数。然后我们需要做什么呢?
- 选择GGX NDF作为我们的概率分布函数,使用它来生成采样方向(入射光的采样方向)。
- 计算对应的权重函数。
- 计算光照方程,按权重累加。
2️⃣首先,GGX NDF为:
其中,
θ
=
<
n
,
h
>
\theta=<n,h>
θ=<n,h>。为了密切匹配镜面BRDF的形状,根据分布函数NDF
和
cos
(
n
,
h
)
\cos{(n,h)}
cos(n,h) 建立PDF
(也就是
p
(
x
i
)
p(x_i)
p(xi)):
p
d
f
(
h
)
=
D
(
h
)
∗
cos
θ
h
pdf(h)=D(h)*\cos{\theta_h}
pdf(h)=D(h)∗cosθh
为什么最终的概率分布函数要乘上余弦项?可以看看:https://www.reedbeta.com/blog/hows-the-ndf-really-defined/。
在求解这个
p
d
f
h
pdf_h
pdfh之前,我们需要知道:因为
p
d
f
h
pdf_h
pdfh是在半矢量空间定义的。然而,积分步骤需要对特定视线方向的光向进行积分。因此,我们需要将
p
d
f
h
pdf_h
pdfh从半矢量空间转换成光空间(从
p
d
f
h
pdf_h
pdfh转换成
p
d
f
L
pdf_L
pdfL)。这种PDF的转换是简单的:
p
d
f
L
=
p
d
f
h
⋅
δ
w
h
δ
w
L
=
p
d
f
h
4
cos
(
v
,
h
)
pdf_L=pdf_h \cdot \frac{\delta{w_h}}{\delta{w_L}}=\frac{pdf_h}{4\cos{(v,h)}}
pdfL=pdfh⋅δwLδwh=4cos(v,h)pdfh
原论文是求输出的pdf,而这里是求输入的,所以,分母的i
要变成v
,这是否解释了上面的PDF转换?
因为我们可以用球面坐标
(
ϕ
,
θ
)
(\phi,θ)
(ϕ,θ) 来表示半向量h,所以我们可以把
p
d
f
h
pdf_h
pdfh表示为
p
d
f
ϕ
pdf_{\phi}
pdfϕ 和
p
d
f
θ
pdf_θ
pdfθ 的乘法:
p
d
f
h
=
p
d
f
ϕ
×
p
d
f
θ
pdf_h={pdf}_{\phi}\times {pdf}_{\theta}
pdfh=pdfϕ×pdfθ
p
d
f
h
pdf_h
pdfh不依赖于角度
ϕ
\phi
ϕ,所以我们可以简单地推导出:
p
d
f
ϕ
=
1
2
π
pdf_{\phi}=\frac{1}{2\pi}
pdfϕ=2π1
所以:
p
d
f
θ
=
D
(
h
)
×
cos
(
θ
h
)
p
d
f
ϕ
=
2
×
α
2
×
cos
(
θ
h
)
(
cos
(
θ
h
)
2
×
(
α
2
−
1
)
+
1
)
2
pdf_{\theta}=\frac{D(h)\times \cos{(\theta_h)}}{pdf_{\phi}}=\frac{2\times \alpha^2\times \cos{(\theta_h)}}{(\cos{(\theta_h)}^2\times (\alpha^2-1)+1)^2}
pdfθ=pdfϕD(h)×cos(θh)=(cos(θh)2×(α2−1)+1)22×α2×cos(θh)
生成采样
3️⃣但我们不能直接使用pdf
来生成采样方向,而是要使用对应的反转累积分布函数CDF_Inv
来生成——这里的数学含义,我是一直无法理解透彻。
参考上图,进行理解。
所以,我们直接给出结果:
c
d
f
ϕ
=
∫
0
ϕ
1
2
π
d
x
=
1
2
π
ϕ
c
d
f
θ
=
∫
0
q
2
×
α
2
×
x
(
x
2
×
(
α
2
−
1
)
+
1
)
2
d
x
=
1
−
q
2
1
+
q
2
(
a
2
−
1
)
cdf_{\phi}=\int^{\phi}_0\frac{1}{2\pi}dx=\frac{1}{2\pi}\phi \\ cdf_{\theta}=\int^q_0 \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)}
cdfϕ=∫0ϕ2π1dx=2π1ϕcdfθ=∫0q(x2×(α2−1)+1)22×α2×xdx=1+q2(a2−1)1−q2
其中,
q
=
c
o
s
(
θ
h
)
q=cos(\theta_h)
q=cos(θh)。现在反转CDF函数,以产生从均匀值
ε
1
,
ε
2
ε_1,ε_2
ε1,ε2,到角度
ϕ
,
θ
\phi,\theta
ϕ,θ 的映射:
ϕ
=
ε
1
×
2
π
θ
=
arccos
1
−
ε
2
1
+
ε
2
(
a
2
−
1
)
\phi=\varepsilon_1 \times 2\pi \\ \theta = \arccos{\sqrt{\frac{1-\varepsilon_2}{1+\varepsilon_2(a^2-1)}}}
ϕ=ε1×2πθ=arccos1+ε2(a2−1)1−ε2
来自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);
}
权重和累加光照
我们在生成采样的函数里面顺便计算了
p
d
f
h
pdf_h
pdfh,然后在主函数里面要将产生的半向量h
转化成入射光:
vec4 re = ImportanceSampleGGX(Xi, N, roughness);
vec3 H = re.xyz;
vec3 L = normalize(2.0 * dot(V, H) * H - V);
此时,我们也需要,将 p d f h pdf_h pdfh 转换成 p d f L pdf_L pdfL ,这个转换公式之前已经给了,所以直接计算:
// pdf_L
float pdf = re.z / (4.0 * VoH);
而真正的权重,不仅和我们预设的分布
p
d
f
L
pdf_L
pdfL有关,还和原分布
π
(
w
i
)
=
f
r
⋅
cos
(
n
,
w
i
)
\pi(w_i)=f_r\cdot \cos{(n,w_i)}
π(wi)=fr⋅cos(n,wi)有关:
w
e
i
g
h
t
i
=
f
r
⋅
cos
(
n
,
w
i
)
p
d
f
L
=
F
⋅
G
⋅
cos
(
v
,
h
)
cos
(
n
,
v
)
⋅
cos
(
n
,
h
)
weight_i=\frac{f_r\cdot \cos{(n,w_i)}}{pdf_L}=\frac{F\cdot G\cdot \cos{(v,h)}}{\cos{(n,v)}\cdot\cos{(n,h)}}
weighti=pdfLfr⋅cos(n,wi)=cos(n,v)⋅cos(n,h)F⋅G⋅cos(v,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的预卷积环境贴图中,使用的权重是 N o L NoL NoL,如何解释?Todo。
Mip map的选择
2.2 Blinn重要性采样算法
Blinn的NDF定义如下:
同样的,实际的概率密度函数需要乘上一个余弦:
为了方便计算,将其转换到球面上来,我们可以得到其在球面坐标下的概率密度函数:
由于这个概率密度函数有2个参数,我们需要分别计算边缘概率密度函数:
计算累计概率函数CDF:
现在反转CDF函数:
权重计算和GGX
类似。
2.3 基于亮度的环境贴图重要性采样
Todo.