表面均匀采样

球体

球体表面上的点 $(x,y,z)$ 可以通过半径 $R$ ，方位角 $\phi\in [0,2\pi]$ ，天顶角 $\theta\in [0,\pi]$ 来定义。

x = \sin{\theta} \cos{\phi} \\ y = \sin{\theta} \sin{\phi} \\ z = \cos{\theta}

如果只是均匀随机生成方位角 $\phi$ 和天顶角 $\theta$ ，在球体表面所采样到的点，实际上并不是均匀分布，而是在两极附近更密集的分布。这是因为球体表面的微分面积 $dA = r^2\sin{\theta}d{\theta} d\phi$ 不仅同时受到微分方位角和微分天顶角的影响，还受到了天顶角正弦值 $\sin{\theta}$ 的影响。当天顶角 $\theta$ 贴近极轴即 $\theta= \pi$ 时， $\sin{ \theta} = 1$ ，面积最大。

由微分立体角的定义 $d\omega = \sin \theta d \theta d \phi$ 可以得出球体的表面积。

\int_{\Theta} d \omega=\int_{0}^{2 \pi} \int_{0}^{\pi} \sin \theta d \theta d \phi = 4\pi

因此在球体表面均匀采样入射方向 $\omega_i$ 时，其概率密度函数为面积的倒数，是一个常数。由此可以推出在一个微分面积 $dA$ 内均匀采样的概率密度函数。

pdf(\omega_i) = \frac{1}{4\pi} \\ P_A = \frac{1}{4\pi}dA

通过二维独立连续随机变量方位角 $\phi$ 和天顶角 $\theta$ 来实现在球体表面的均匀采样，需要计算出用来参数化二者的概率密度函数 $pdf(\theta,\phi)$ ，那么首先求出对应的边缘密度函数。

f(\theta)=\int_{0}^{2\pi} pdf(\theta, \phi) d \phi=\frac{\sin \theta}{2} \\ f(\phi)=\int_{0}^{\pi} pdf(\theta, \phi) d \theta = \frac{1}{2\pi}

根据边缘密度函数求出累积分布函数 $F_{X}(x)$ 。

F_{\theta}(\theta)=\int_{0}^{\theta} f(\hat{\theta}) d \hat{\theta} = \frac{1-\cos \theta}{2}

F_{\phi}(\phi)=\int_{0}^{\phi} f(\hat{\phi}) d \hat{\phi}=\frac{\phi}{2 \pi}

对累积分布函数做逆变换取样，得到累积分布函数的反函数 $F_X^{-1}(x)$

F_{\theta}^{-1}(x)=\arccos (1-2 x)

F_{\phi}^{-1}(x)=2 \pi x

对于连续随机变量方位角 $\phi$ 和天顶角 $\theta$ ，如果要获得球体表面上均匀分布的点，只需要给定均匀分布变量 $\xi_1, \xi_2 \sim Uniform[0,1]$ 。

\theta =F_{\theta}^{-1}(\xi_1) = \arccos(1-2\xi_1) \\ \phi =F_{\phi}^{-1}(\xi_2) = 2\pi\xi_2

转为三维坐标可得

x = \sqrt{2\xi_1} \cos{(2\pi \xi_2)} \\ y = \sqrt{2\xi_1} \sin{(2\pi \xi_2)} \\ z = \sqrt{1-2\xi_1}

三角形

取两个标准正态分布变量 $u,v \sim Uniform[-1,1]$ 。利用三角形重心坐标，我们可以通过两个随机变量，使得重心坐标 $\alpha , \beta, \gamma$ 线性组合的结果落在三角形内，且三者分布是均匀的，即可实现在三角形内均匀采样。

α = 1 - \sqrt{u}, β = \sqrt{u} \cdot v, γ = \sqrt{u} \cdot (1 - \sqrt{v})

设 $s = \sqrt{u}$ ，则可以进一步简化公式。

α = 1 - s, β = s \cdot v, γ = 1 - α - β

矩形

取两个标准正态分布变量 $u,v \sim Uniform[-1,1]$ 。从矩形 $ABCD$ 中选取相邻两边 $AB,AD$ 。

P = u \cdot AB + v \cdot AD

拒绝采样法

拒绝采样法是一个有效但不太优雅的算法。例如在单位球体内均匀采样时，随机生成一个标准立方体范围的点，计算该点到球心的距离，大于球体半径则重新随机生成点。

再例如在一个二维三角形内均匀采样，首先确定三角形的包围盒范围，即从三角形三个顶点中选取最大最小 $x$ 坐标和 $y$ 坐标。

max = (x_{max}, y_{max}) \\ min = (x_{min}, y_{min})

取两个标准正态分布变量 $u,v \sim Uniform[-1,1]$ ，随机生成一个包围盒内的点。由最大坐标和最小坐标求出在包围盒内的均匀取样点 $P$ 。再通过重心坐标或者叉乘法等方法判断该点是否三角形内，不在则重新随机生成点

P = (x_{min} + u \cdot x_{max}, y_{min} + v \cdot y_{max})

重要性采样

蒙特卡洛积分方法的路径追踪在半球面上进行采样。已知球面上均匀采样的三维坐标，同理可求得在半球上采样时，给定均匀分布变量 $\xi_1, \xi_2 \sim Uniform[0,1]$ ，求出满足其累积分布函数的随机变量 $\theta, \phi$ ，从而得到采样入射方向的三维坐标 $x,y,z$

\theta = \arccos(1-\xi_1)

\phi= 2\pi\xi_2

x = \sqrt{\xi_1} \cos{(2\pi \xi_2)}

y = \sqrt{\xi_1} \sin{(2\pi \xi_2)}

z = \sqrt{1-\xi_1}

均匀采样

若在半球面上进行均匀采样，则概率密度函数 $pdf(x)$ 应该是各处相同的常数 $C$ 。将积分域转换到半球面上，即可求出常数 $C$ 。

\int_{\Omega^+} C d x=1

\int_{0}^{2 \pi} \int_{0}^{\frac{\pi}{2}} C\sin \theta d \theta d \phi = 1

C=\frac{1}{2\pi}

cosine-weighted 采样

对 $\text{cosine}$ 项进行重要性采样，即 $\text{PDF}$ 与 $\text{cosine}$ 成正比。将积分域转换到半球面上，即可求出系数 $c$ 。

pdf(x) = c \cdot \cos{\theta}

\int_{\Omega^+}c \cdot \cos{\theta}d\omega = 1

\int_{0}^{2 \pi} \int_{0}^{\frac{\pi}{2}} c \cdot \cos{\theta} sin{\theta} d \theta d \phi = 1

c = \frac{1}{\pi}

Cook-Torrance BRDF 重要性采样

漫反射项

通过 $\text{Cook-Torrance BRDF}$ 中漫反射的 $\text{BRDF}$ ，我们可以计算得出其渲染方程和蒙特卡洛积分方法下的渲染方程。

f_r(p, \omega_o,\omega_i) = \frac{1}{\pi}

L_o(p,\omega_o) = \frac{\cos \theta_{\mathrm{i}}}{\pi} L_{\mathrm{i}}\left(\mathrm{p}, \omega_{\mathrm{i}}\right)

L_o(p,\omega_o) ≈ \frac{1}{N} \sum_{i=1}^{N}\frac{\cos \theta_{\mathrm{i}}}{\pi} \frac{ L_i(p, \omega_i)}{pdf(\omega_i)}

观察渲染方程，我们可以对漫反射采用 $\text{cosine-weighted}$ 采样以加快渲染方程的收敛。

镜面反射项

镜面反射项中可以排除需要知道采样方向才能计算的菲涅尔项、非连续函数的几何项，因此我们只需要使概率分布函数近似法线分布函数即可。

根据法线分布函数 $D(\omega_h)$ 的物理含义——每单位立体角 $d\omega$ 、每单位面积 $dA$ 上所有法向量为 $\omega_h$ 的微平面的面积，约定微表面中符合法线分布的微表面投影在宏观平面的面积为1，因此我们可以得到其关于法向量为 $\omega_h$ 的概率密度函数。

\int_{\Omega^+}\cos{\theta_h}D(\omega_h)d\omega_h=1

经过推导可得出法线和入射方向二者的换算关系，经过转换得到关于入射方向 $\omega_i$ 的镜面反射项 $pdf(\omega_i)$ 。

\frac{d\omega_h}{d\omega_i} = \frac{1}{4(\omega_i \cdot \omega_h)}

\int_{\Omega^+}\frac{\cos{\theta_h}}{4(\omega_i \cdot \omega_h)}D(\omega_h)d\omega_i=1

pdf(\omega_i) =\frac{\cos{\theta_h}}{4(\omega_i \cdot \omega_h)}D(\omega_h)

但在计算采样方向时我们无需使用 $pdf(\omega_i)$ ，还是使用前面已经计算好的 $pdf(\omega_h)$ ，通过该 PDF 推导出法向量 $\omega_h$ 后，再以 $\omega_h$ 为镜面反射法线计算我们最终需要的入射方向 $\omega_i$ 即可。

如果使用不同的法线分布函数，显然随机采样的方向计算公式也是不同的，积分计算过程较为复杂，这里不一一计算了。

对光源重要性采样

当我们要对光源进行直接采样时，需要将原本渲染方程中对微分立体角的微元，换源成在光源上的面积微元。

d\omega = \sin{\theta}d\theta d\phi

我们知道微分立体角 $d\omega$ 可以认为是一个微小矩形，如果存在矩形光源，则可以根据相似三角形原则，得到单位半球的半径的平方与着色点 $x$ 到矩形光源的距离平方之比，即为 $d\omega$ 的面积与映射到光源上的面积 $dA \cdot \cos{\theta}$ 之比，从而得到新的渲染方程。

\frac{d\omega}{dA \cdot \cos{\theta}} = \frac{1}{\left\| x-x^{\prime} \right\|^2}

L_o ≈ \frac{1}{N} \sum_{i=1}^{N} L_i \cos{\theta_i} \cos{\theta_o} \frac{1}{\left\| x-x^{\prime} \right\|^2} \frac{1}{pdf} dA

对于面光源，在面光源上均匀采样则其概率密度函数是一个常数。对于球体光源，相比面光源较为复杂，在此不做计算。

pdf = \frac{1}{Area}

反走样

SSAA

超采样抗锯齿 $\text{Super-Sampling Anti-aliasing}$ 的原理是提升采样的分辨率。例如先上采样 $2 \times 2$ 分辨率，则需要四倍的帧缓冲，然后再下采样到原来的分辨率。

MSAA

以 $2 \times 2 \text {MSAA}$ 为例，多重采样抗锯齿 $\text{Multi-Sampling Anti-aliasing}$ 将像素分为四个子像素，根据三角形覆盖四个子像素的比例决定该像素颜色。

例如一个三角形占据了四分之一个子像素，则该三角形为该像素贡献了四分之一的自身颜色值。但由于需要进行深度测试，每个子像素也都需要存储一个深度值，因此深度缓冲将变为原来的 $2 \times 2$ 倍。考虑到多个三角形再单个像素区域有可能发生各占据一部分子像素的情况，因此颜色缓冲也将变为原来的 $2 \times 2$ 倍。

三角形滤波

光线追踪的反走样算法，通常分为两种，一种是先均匀采样再加权累加。一种是先按概率密度采样，再累加平均。

后者的方法中，显然最合理的光线采样分布应该是边缘稀疏，越靠近中心越密集。这种分布理论上 $\text{sinc}$ 函数是最优的，其能过滤掉高频信号，保留低频信号。但是完美的 $\text{sinc}$ 函数只存在于理论，实际工程中只能近似模拟，例如三角形滤波。

sinc(x) = sin(x)/x

\int_{-\infty}^{+\infty} sinc(x)=\pi \\ \\ \int_{-\infty}^{+\infty} sinc(\pi x)=1

对于三角形滤波，其在定义域 $[-1,1]$ 中满足积分为 $1$ ，从而可以得出其边缘密度函数 $f(x)$ 。如果要得到一个符合三角形滤波分布的变量，只需要给定均匀分布变量 $x \sim Uniform[0,1]$ 。求出其累积分布函数 $F_{X}(x)$ ，然后做逆变换取样得到其反函数 $F_{X}^{-1}(x)$ ，而反函数 $F_{X}^{-1}(x)$ 就是我们想要的随机变量 $X$ 。

f(x) = \left\{\begin{array}{ll} x+1, &-1≤x≤0\\ 1-x, &0≤x≤1 \end{array}\right.

F_{X}(x) = \int_{-∞}^{x} f(t)dt = \left\{\begin{array}{ll} \frac{1}{2}(x+1)^2, &-1≤x≤0 \\ -\frac{1}{2}(x-1)^2+1, &0≤x≤1 \end{array}\right.

X = F_{X}^{-1}(x) = \left\{\begin{array}{ll} \sqrt{2x} -1, &0≤x≤\frac{1}{2} \\ 1-\sqrt{2-2x}, &\frac{1}{2}≤x≤1 \end{array}\right.

如果将每个像素划分为 $2 \times 2$ 的子像素，分别对每个子像素应用三角形滤波，然后累加后取平均值，即可得到梯形滤波。

纹理采样

采样纹理时需要注意，我们获取了三角形的 $UV$ 坐标和对应纹理贴图大小后，所计算出的纹理坐标是整形，因此小数部分会舍去。而纹理坐标的范围是 $x \sim [0,width],y \sim [0,height]$ 因此需要为纹理坐标各个维度加 $0.5$ 以实现四舍五入，获取距离采样点最近的纹理坐标。

双线性性插值

以上方法显然不够精准，只是获取了最近的纹理坐标，而双线性插值则将周围的纹理坐标一并考虑进来，进行加权累加。

以采样点所在像素为重心，获取其周围四个象限的八个像素。首先依据象限位置决定要对哪个象限的四个像素进行加权累加。

将四个像素的中心之间的范围视为一个新的像素坐标系，则左下角就是 $(0,0)$ ，采样点在坐标系中的坐标 $(s,t)$ 就是各个像素的贡献。先在 $x$ 轴上以 $s$ 对四个像素两两线性插值，最后在 y 轴上以 t 对两个结果进行线性插值。

mipmap

对纹理以相邻四个像素为单位进行加权平均，最终得到一个原图像四分之一大小的纹理，以此类推，最终得到一个金字塔结构的纹理集合，这就是 mipmap。

mipmap 常用于解决纹理瑕疵，例如过远的三角形，却占据很大一块纹理范围。如果根据远近，获取对于层级的mipmap纹理，就可以一定程度上结果这个问题。

三线性插值

$\text{mipmap}$ 显然是不够完美的，各层级之间的衔接过渡可能过于生硬。因此在双线性插值的基础上，计算 $\text{mipmap}$ 的层级不再只取整数，而是允许 $1.3$ 这种层级的存在，我们可以利用 $[0,1]$ 的层级之差对相邻层级再做一次线性插值，也被称为三线性插值。

后处理

mipmap 近似高斯滤波

该技术原理不在此赘述，详情可以阅读另一篇文章

图形学论文阅读 - 使用盒型滤波图像金字塔实时近似大卷积核 - 掘金 (juejin.cn)

渲染中的采样：从概率密度到图像处理