渲染方程背景知识
本来想写BSDF的推导，但是发现基础知识很多，正好整理一下作为复习。部分内容参考Physically Based Rendering, Third Edition，非常推荐阅读。
概率论入门
概率密度函数(PDF: probability distribution function)
以一维为例，如果随机变量\(X\)是连续的，就需要使用概率密度函数\(p(x)\)来描述它。\(P(x)\)是累积分布函数。
\[ P(x)=\int_{-\infty}^{x} p(x)dx. \]
因为所有事件发生概率之和为1，所以如果\(x\in[a,b]\)，
\[ P(b)=1. \]
期望与方差(Expected Values and Variance)
如果有一个\(x\)的分布\(p(x)\)和一个函数\(f\)，对应的期望\(E_p[f(x)]\)的计算公式为：
\[ E_p[f(x)]=\int_{D} f(x)p(x)dx. \]
举个二维的例子，向一个半径为\(r_0\)的靶子投飞镖，得分计算方式如下：
\[ f(x,y)=\sqrt{x^2+y^2} \]
如果投手投出的飞镖是均匀分布在该靶子上，那么
\[ p(x,y)=\frac{1}{\pi r_0^2} \]
期望则是：
\[ E_p[f(x)]=\int_{D} \frac{\sqrt{x^2+y^2}}{\pi r_0^2}dxdy, \]
\(D\)为半径为\(r\)的圆内部。
该例子用极坐标变换计算：
\[ x=r\cos\theta, y=r\sin\theta.\\ \frac{d(x,y)}{d(r,\theta)}=\begin{vmatrix}\frac{\partial x}{\partial r}&\frac{\partial x}{\partial \theta}\\\frac{\partial y}{\partial r}&\frac{\partial y}{\partial \theta}\end{vmatrix}=\begin{vmatrix}\cos\theta& -r\sin\theta\\\sin \theta & r \cos \theta\end{vmatrix}=r\\ \begin{align*} \int_{D} \frac{\sqrt{x^2+y^2}}{\pi r^2}dxdy &=\int_{D} \frac{r}{\pi r_0^2}|\frac{d(x,y)}{d(r,\theta)}|drd\theta\\ &=\int_{0}^{2\pi}\int_{0}^{r_0} \frac{r^2}{\pi r_0^2}drd\theta\\ &=\frac{2r_0}{3} \end{align*} \]
也就是说，按照上述规则投手最终得分为\(\frac{2r_0}{3}\)分。
方差描述了随机变量的值与期望之间差别的偏离程度，方差越大，说明该随机变量的值与期望有较大差距。比如某地降雨量的期望是100mm，但是产生的值只为0mm或者200mm且等概率出现，那么此随机变量的方差与每次降雨量产生的值都为100mm的方差会大很多，后者为0。
计算公式为:
\[ \begin{align*} V[f(x)]&=E[(f(x)-E[f(x)])^2]\\ &= E[f(x)^2-2f(x)E[f(x)]+{(E[f(x)])}^2]\\ &= E[f(x)^2]-2E[f(x)E[f(x)]]+E[{(E[f(x)])}^2]\\ &= E[f(x)^2]-2E[f(x)]E[f(x)]+{(E[f(x)])}^2\\ &= E[f(x)^2]-{E[f(x)]}^2 \end{align*} \]
因为\(E[f(x)]\)是一个常数，因此他的期望就是自己，可以由此化简。
蒙特卡洛估计(The Monte Carlo Estimator)
想要计算
\[ F(x)=\int_a^bf(x)dx， \]
给定一个均匀分布的随机变量\(X_i\in[a,b]\)，可以通过下面的公式估计\(F(x)\)：
\[ F_N=\frac{b-a}{N}\sum_{i=1}^Nf(X_i)， \]
可以证明\(E[F_N]\)=F(x)： 首先，由均匀分布可知：
\[ p(x)=\frac{1}{b-a} \]
因此，
$$ \[\begin{align*} E[F_N]&=E[\frac{b-a}{N}\sum_{i=1}^Nf(X_i)]\\ &=\frac{b-a}{N}\sum_{i=1}^NE[f(X_i)]\\ &=\frac{b-a}{N}\sum_{i=1}^N\int_a^bf(x)p(x)dx\\ &=\frac{1}{N}\sum_{i=1}^N\int_a^bf(x)dx\\ &=\int_a^bf(x)dx \end{align*}\] $$
如果随机变量服从分布\(p(x)\)，依然可以通过下面的公式估计\(F(x)\)：
\[ F_N=\frac{1}{N}\sum_{i=1}^N\frac{f(X_i)}{p(X_i)}， \]
证明如下：
$$ \[\begin{align*} E[F_N]&=E[\frac{1}{N}\sum_{i=1}^N\frac{f(X_i)}{p(X_i)}]\\ &=\frac{1}{N}\sum_{i=1}^NE[\frac{f(X_i)}{p(X_i)}]\\ &=\frac{1}{N}\sum_{i=1}^N\int_a^b\frac{f(x)}{p(x)}p(x)dx\\ &=\frac{1}{N}\sum_{i=1}^N\int_a^bf(x)dx\\ &=\int_a^bf(x)dx \end{align*}\] $$
该估计方式可以扩展到高维，在高维的情况下，使用该方法是能大幅度提高计算效率的。要注意的是，该方法的误差为\(O(\sqrt{N})\)，也就是说采样数为10000的估计结果只比采样数为100精确一个小数点。 渲染方程实际上无法直接解，要通过离散化计算结果。最常用的方法是蒙特卡洛积分(Monte Carlo Integration)。
采样随机变换(Sampling Random Variables)：逆变换采样(The Inversion Method)
动机
实际使用随机变量时，我们希望能够在给定的分布中采样随机变量，这样就可以用更快收敛的方式用蒙特卡洛积分来估计积分。
步骤
- 计算累积分布函数CDF:\(P(x)=\int_0^xp(x')dx'\)。
- 求CDF的反函数(inverse) \(P^{-1}(x)\)。
- 使用均匀分布的随机变量采样得到\(\xi\)。
- \(X_i=P^{-1}(\xi)\)。
以cosine-weighted采样为例，已知，在球坐标系中任意球面的极小面积为：
\[ dA = (r\sin\theta\,d\varphi)(r d\theta)=r^2(\sin\theta\,d\theta\,d\varphi) \]
立体角是投影面积与球半径平方值的比，立体角与球面坐标的微分有如下转换关系:
\[ d\omega=\sin\theta d\theta d\phi. \]
易得半个球面的曲面积分为：
\[ \begin{align*} \int_H d\omega&= \oiint_H \sin\theta d\theta d\phi\\ &=\int_0^{2\pi}d\phi\int_0^{\frac{\pi}{2}}\sin\theta d\theta \\ &=2\pi. \end{align*} \]
我们希望采样的点满足\(p(\omega)=C cos\theta\)。代入积分得
\[ \begin{align*} \int_H p(\omega)d\omega&= \oiint_H C cos\theta\sin\theta d\theta d\phi\\ &=C\int_0^{2\pi}d\phi\int_0^{\frac{\pi}{2}}cos\theta\sin\theta d\theta \\ &=C\pi \end{align*} \]
由此可知
$$ \[\begin{align*} p(\omega)&= \frac{cos\theta}{\pi}\\ p(\theta,\phi)&=\frac{sin\theta cos\theta}{\pi} \end{align*}\] $$
对两个变量分别积分可得：
\[ \begin{align*} p(\theta)&=\int_{0}^{2\pi}\frac{sin\theta cos\theta}{\pi}d\phi=\sin2\theta\\ P(\theta)&=\int_{0}^{\theta}p(\theta)d\theta\\&=\frac{1-\cos2\theta}{2}\\ &=\sin^2\theta\\ p(\phi)&=\int_{0}^{\frac{\pi}{2}}\frac{sin\theta cos\theta}{\pi}d\theta=\frac{1}{2\pi}\\ P(\phi)&=\int_{0}^{\phi}p(\phi)d\phi=\frac{\phi}{2\pi}\\ \end{align*} \]
计算反函数：
\[ \begin{align*} P^{-1}(\theta)&=\arcsin{\sqrt{\theta}}\\ P^{-1}(\phi)&=2\pi \phi\\ \end{align*} \]
采样步骤为：
- 在\([-1,1]\)上均匀采样\(x\)，计算\(\arcsin{\sqrt{x}}\)得到\(\theta\)
- 在\([-1,1]\)上均匀采样\(y\)，计算\(2\pi y\)得到\(\phi\)
The Rejection Method
对于无法求反函数的\(f(x)\)的CDF，可以用The Rejection Method： 循环：
- 从已知的\(p(x)\)采样X
- 如果\(xi<f(X)/c\cdot p(X)\)
- 返回X