[图形学]在半球面上按照微表面模型采样

一、简介

本文介绍了如何根据 pbr 中的法向分布函数（Normal Distribution Function, NDF）在半球上进行采样。

二、微表面模型

在微表面模型理论中，任何宏观上的表面都由无数法向不同的微表面组成。宏观表面的粗糙度越大，各微表面的法向与宏观表面的法向差异越大。宏观表面的粗糙度越小，那么各微表面的法向与宏观表面的法向越一致。如下图所示（左边粗糙度大，右边粗糙度小）：

体现在渲染结果中即粗糙度越大，镜面反射越扩散，渲染结果越粗糙如下图所示：

在 pbr 中常使用 法向分布函数（Normal Distribution Function, NDF）描述微表面与宏观表面法向的一致性。NDF函数定义为：单位立体角内，法向为 $h$ 的微表面面积与 宏观表面面积 A 的比，公式如下：
$$D(h)=\frac{d{A}_h}{A\ast d{\omega }_h}$$
其中，$h$ 为半程向量，$h=\frac{l+v}{||l+v||}$，$A$为宏观表面的面积，$d{A}_h$为宏观表面上法向在立体角 ${\omega }_h$ 内的微表面的面积，函数$D(h)$的单位为 $\frac{1}{sr}$。
常用的GGXTR NDF 公式如下：
$$D(h,n,\alpha )=\frac{{\alpha }^2}{\pi ((n\cdot h{)}^2({\alpha }^2-1)+1{)}^2}$$
其中，$h$ 为半程向量，$n$ 为宏观表面的法向，$\alpha$为粗糙度。

三、在半球上根据微表面模型采样

在半球上根据微表面模型采样$h$，我们需要按照以下步骤进行：
(1). 确定采样方向$h$ 的 pdf 函数，$pdf(h)$；
(2). 将半球面上微分 $d{\omega }_h$ 变换到微分 $d\theta ,d\varphi$；
(3). 计算得到关于 $\theta$ 和 $\varphi$的 pdf 函数，$pdf(\theta ),pdf(\varphi )$；
(4). 计算对应的 cdf 函数，$cdf(\theta ),cdf(\varphi )$；
(5). 计算 cdf 函数的逆函数，$cd{f}^{-1}(\theta ),cd{f}^{-1}(\varphi )$；
(6). 采样 [0,1] 范围内部的均匀分布变量${ϵ}_1,{ϵ}_2$，那么方向 $(\theta =cd{f}^{-1}({ϵ}_1),\varphi =cd{f}^{-1}({ϵ}_2))$ 即为满足微表面模型的方向。

1. 确定 $pdf(h)$

(1). 为什么 D(h) 不是 pdf(h)

一个需要强调的事情是，函数 $D(h)$ 并不是 $pdf(h)$！！！
这是因为 $D(h)$ 描述的是，单位立体角内，法向为 $h$ 的微表面面积与 宏观表面面积 A 的比。
而 $pdf(h)$ 应该为，单位立体角内，法向为 $h$ 的微表面面积与 总微表面面积 Am 的比。
如下图所示，很明显 宏观表面面积A 要小于 总微表面面积Am。

我们通过在半球上积分 $D(h)$ 也可以得到：
$${\int }_{\Omega }D(h)d{\omega }_h={\int }_{\Omega }\frac{d{A}_h}{A\ast d{\omega }_h}d{\omega }_h=\frac{{\int }_{\Omega }d{A}_h}{A}=\frac{A_m}{A}$$
但是 $A$ 小于 $Am$。这说明 $D(h)$ 不能满足$pdf(h)$ 在定义域上积分等于 1 的要求。

(2). pdf(h) 公式是什么

根据 $pdf(h)$ 的定义，pdf(h)是：单位立体角内，法向为 $h$ 的微表面面积与 总微表面面积 Am 的比。
可以得到：
$$pdf(h)=\frac{d{A}_h}{{\int }_{\Omega }d{A}_{h^{\prime }}}=\frac{D(h)\ast A\ast d{\omega }_h}{{\int }_{\Omega }D(h^{\prime })\ast A\ast d{\omega }_{h^{\prime }}}=\frac{D(h)\ast A}{A_m}=\frac{D(h)}{(A_m/A)}$$
因此，一个直接的方案是，根据 $D(h)$，使用数值积分或者计算积分的方法求得 $A_m={\int }_{\Omega }D(h^{\prime })\ast A\ast d{\omega }_{h^{\prime }}$，那么理想的 $pdf(h)$ 即为：
$$pdf(h)=\frac{D(h)}{(A_m/A)}$$
现在我们已经知道应该如根据给定的 $D(h)$计算理想的 $pdf(h)$ 了，只要简单的求解一个数值积分 $A_m={\int }_{\Omega }D(h^{\prime })\ast A\ast d{\omega }_{h^{\prime }}$ 即可。
可是问题是，求解积分 $A_m={\int }_{\Omega }D(h^{\prime })\ast A\ast d{\omega }_{h^{\prime }}$ 有时候并不是一个简单的问题，尤其是 $D(h)$ 还跟粗糙度$\alpha$相关。

(3). 一个近似的 pdf’(h)

根据前面的分析，根据 $D(h)$ 求理想的 $pdf(h)$ 似乎并不是一个可行的方案，我们得重新思考下为什么需要 $pdf(h)$。

这是因为在 重要性采样 中，我们需要根据蒙特卡洛采样方法，利用数值积分方法求解渲染方程：
$$\begin{gathered} L(p,wo)=\int (kd\frac{c}{\pi }+\frac{D\ast G\ast F}{4(n\cdot \omega i)(n\cdot \omega o)})\ast L(p,\omega i)\ast (n\cdot \omega i)d\omega i \\ =\int kd\frac{c}{\pi }\ast L(p,\omega i)\ast (n\cdot \omega i)d\omega i+\int \frac{D\ast G\ast F}{4(n\cdot \omega i)(n\cdot \omega o)}\ast L(p,\omega i)\ast (n\cdot \omega i)d\omega i \end{gathered}$$
我们主要关注的是镜面反射部分，即：
$$Ls(p,wo)=\int \frac{D\ast G\ast F}{4(n\cdot \omega i)(n\cdot \omega o)}\ast L(p,\omega i)\ast (n\cdot \omega i)d\omega i$$

我们想要根据一定的 $pd{f}^{\prime }(h)$，对 $h$ 进行采样，再利用蒙特卡洛积分方法可以得到上述积分公式的值：
$$Ls(p,wo)\approx \frac{1}{N}\sum \frac{\frac{D\ast G\ast F}{4(n\cdot \omega i)(n\cdot \omega o)}\ast L(p,\omega i)\ast (n\cdot \omega i)}{pd{f}^{\prime }(h)}$$
其中 $\omega i$ 可以根据采样得到的 $h$ 和 $\omega o$ 计算得到。
该 $pd{f}^{\prime }(h)$ 需要满足以下条件：

• 🎯积分为1，即 $\int pd{f}^{\prime }(h)d{\omega }_h=1$；
• 🎯$pd{f}^{\prime }(h)$ 最好等于 $D(h)\ast f^{\prime }(h)$。其中 $f^{\prime }(h)$ 是一个简单的函数形式或者定值 $C$，这样可以消去蒙特卡洛积分中分子里的 $D$ 项；

回到我们之间讲的微表面理论中，我们将微表面 $d{A}_h$ 从 $-n$ 方向投影到 宏观表面 上可以得到（乘以$n\cdot h$即为对微表面进行投影）：
$$\int (n\cdot h)d{A}_h=A$$
即：
$$\int D(h)\ast A\ast (n\cdot h)d{\omega }_h=A$$
那么：
$$\int D(h)\ast (n\cdot h)d{\omega }_h=1$$
假如我们令 $pd{f}^{\prime }(h)=D(h)\ast (n\cdot h)$，就正好满足刚刚我们提到的两个条件：

• ✅积分值为1；
• ✅是 $D(h)\ast f(h)$的形式。

因此我们可以按照 $pd{f}^{\prime }(h)=D(h)\ast (n\cdot h)$ 对 $h$ 进行采样（以下为了方便，将$pd{f}^{\prime }(h)$写作$pdf(h)$），这样得到的半程向量 $h$ 既能保证在半球上积分值为1，又能便于在蒙特卡洛积分中根据 $D$ 分布进行重要性采样，降低方差，加快收敛速度。

2. 微分 $d{\omega }_h$ 变换到 $d\theta ,d\varphi$

根据上面的介绍我们已经确定，采样 $h$ 的 $pdf(h)=D(h)\ast (n\cdot h)$。因为在半球上采样方向 $h$ 是通过采样角度球面坐标上的 $\theta ,\varphi$ 得到的，因此在采样时我们需要将半球上的方向向量 $h$ 转为球面坐标系上的坐标 $(\theta ,\varphi )$。
因为：
$$d{\omega }_h=sin(\theta )d\theta d\varphi$$
因此：
$${\int }_{{\Omega }^{+}}pdf(h)d{\omega }_h={\int }_0^{2\pi }{\int }_0^{\pi /2}pdf(\theta ,\varphi )\ast sin(\theta )d\theta d\varphi$$
根据 $pdf(h)=D(h)\ast (n\cdot h)$ 采样 $h$，可以通过根据 $pdf(\theta ,\varphi )\ast sin(\theta )=D(\theta ,\varphi )\ast cos(\theta )\ast sin(\theta )$ 采样 ($\theta ,\varphi$) 得到。

3. 计算边缘概率密度 $pdf(\theta ),pdf(\varphi )$

边缘概率密度 $pdf(\theta )$：
$$\begin{gathered} pdf(\theta )={\int }_0^{2\pi }pdf(\theta ,\varphi )\ast sin(\theta )d\varphi \\ ={\int }_0^{2\pi }D(\theta ,\varphi )\ast cos(\theta )\ast sin(\theta )d\varphi \\ ={\int }_0^{2\pi }\frac{{\alpha }^2}{\pi (co{s}^2(\theta )({\alpha }^2-1)+1{)}^2}\ast cos(\theta )\ast sin(\theta )d\varphi \\ =\frac{2{\alpha }^2\ast cos(\theta )\ast sin(\theta )}{(co{s}^2(\theta )({\alpha }^2-1)+1{)}^2} \end{gathered}$$
边缘概率密度 $pdf(\varphi )$：
$$\begin{gathered} pdf(\varphi )={\int }_0^{\pi /2}pdf(\theta ,\varphi )\ast sin(\theta )d\theta \\ ={\int }_0^{\pi /2}D(\theta ,\varphi )\ast cos(\theta )\ast sin(\theta )d\theta \\ ={\int }_0^{\pi /2}\frac{{\alpha }^2}{\pi (co{s}^2(\theta )({\alpha }^2-1)+1{)}^2}\ast cos(\theta )\ast sin(\theta )d\theta \\ =\frac{{\alpha }^2}{2\pi }{\int }_{\pi /2}^0\frac{d(co{s}^2(\theta ))}{(co{s}^2(\theta )({\alpha }^2-1)+1{)}^2} \\ =\frac{{\alpha }^2}{2\pi }{\int }_0^1\frac{dt}{(({\alpha }^2-1)t+1{)}^2} \\ =\frac{{\alpha }^2}{2\pi ({\alpha }^2-1)}{\int }_1^{{\alpha }^2}\frac{dx}{x^2} \\ =\frac{{\alpha }^2}{2\pi ({\alpha }^2-1)}\ast (\frac{1}{x}){|}_{{\alpha }^2}^1 \\ =\frac{1}{2\pi } \end{gathered}$$

4. 计算 $cdf(\theta ),pdf(\varphi )$

计算 $cdf(\theta )$：
$$\begin{gathered} cdf(\theta )={\int }_0^{\theta }pdf(t)dt \\ ={\int }_0^t\frac{2{\alpha }^2\ast cos(t)\ast sin(t)}{(co{s}^2(t)({\alpha }^2-1)+1{)}^2}dt \\ ={\alpha }^2{\int }_{\theta }^0\frac{dco{s}^{2}t}{(({\alpha }^2-1)co{s}^{2}t+1{)}^2} \end{gathered}$$
令 $x=({\alpha }^2-1)co{s}^{2}t+1$，那么：
$$\begin{gathered} cdf(\theta )={\alpha }^2{\int }_{\theta }^0\frac{dco{s}^{2}t}{(({\alpha }^2-1)co{s}^{2}t+1{)}^2} \\ =\frac{{\alpha }^2}{({\alpha }^2-1)}{\int }_{({\alpha }^2-1)co{s}^2\theta +1}^{{\alpha }^2}\frac{dx}{x^2} \\ =\frac{{\alpha }^2}{({\alpha }^2-1)}\ast (\frac{1}{x}){|}_{{\alpha }^2}^{({\alpha }^2-1)co{s}^2\theta +1} \\ =\frac{{\alpha }^2}{({\alpha }^2-1{)}^{2}co{s}^2\theta +({\alpha }^2-1)}-\frac{1}{{\alpha }^2-1} \end{gathered}$$
计算 $cdf(\varphi )$：
$$\begin{gathered} cdf(\varphi )={\int }_0^{\varphi }pdf(\varphi )d\varphi \\ ={\int }_0^{\varphi }\frac{1}{2\pi }d\varphi \\ =\frac{\varphi }{2\pi } \end{gathered}$$

5. 计算 $cd{f}^{-1}(\theta ),cd{f}^{-1}(\varphi )$

计算 $cd{f}^{-1}(\theta )$：
$$cd{f}^{-1}(\theta )=arccos\sqrt{\frac{1-\theta }{\theta ({\alpha }^2-1)+1}}$$
计算 $cd{f}^{-1}(\varphi )$：
$$cd{f}^{-1}(\varphi )=2\pi \varphi$$

6. 采样均匀随机变量 ${ϵ}_1,{ϵ}_2$ 得到目标方向 $(\theta ,\varphi )$

根据 $cd{f}^{-1}(\theta )$ 和 $cd{f}^{-1}(\varphi )$，我们可以在 [0,1] 范围内随机采样得到随机变量 ${ϵ}_1,{ϵ}_2$。
那么，
$$(\theta =cd{f}^{-1}({ϵ}_1),\varphi =cd{f}^{-1}({ϵ}_2))$$
即为满足概率分布 $pdf(\theta ,\varphi )\ast sin(\theta )=D(\theta ,\varphi )\ast cos(\theta )\ast sin(\theta )$ 的 ($\theta ,\varphi$)，那么对应的 $h$ 也满足我们的目标概率密度分布 $pdf(h)=D(h)\ast (n\cdot h)$。

7. 实现代码

在 UE 中对GGX 法向分布根据 $pdf(h)=D(h)\ast (n\cdot h)$ 的重要性采样实现如下Github-EpicGames-UnrealEngine：

// PDF = D * NoH / (4 * VoH)
float4 ImportanceSampleGGX( float2 E, float a2 )
{
	float Phi = 2 * PI * E.x;
	float CosTheta = sqrt( (1 - E.y) / ( 1 + (a2 - 1) * E.y ) );
	float SinTheta = sqrt( 1 - CosTheta * CosTheta );

	float3 H;
	H.x = SinTheta * cos( Phi );
	H.y = SinTheta * sin( Phi );
	H.z = CosTheta;
	
	float d = ( CosTheta * a2 - CosTheta ) * CosTheta + 1;
	float D = a2 / ( PI*d*d );
	float PDF = D * CosTheta;

	return float4( H, PDF );
}

四、参考引用

[1].GGX重要性采样-博客园
[2].How Is The NDF Really Defined?
[3].Sampling Microfacet BRDF