均匀采样问题总结

摘要

均匀采样问题是图形学中常遇到的一类问题，这类问题的描述如下：

如何用[0,1]之间的均匀采样获得在xxx上的均匀采样？

其中，xxx可以是圆盘、球面、球体等。我们可以利用的只有[0,1]之间的均匀采样，这个可以通过程序的伪随机实现。

总的来说，这类问题有着类似的思路和解决方案，都可以利用概率论知识一步步推导出来。本文将介绍相关方法，并给出几个常用的示例。

理论

概率论知识回顾

注：如果基础比较扎实，可以直接跳到求解思路

累积分布函数 随机变量$X$的累积分布函数（cumulative distribution function，缩写为CDF）记为$P(x)$：

$$P_x(x)=Pr(X\le x)$$

其中，$Pr(..)$表示某事件的概率。

概率密度函数 随机变量$X$的概率密度函数（probability density function，缩写为PDF）记为$p(x)$，满足：

$$P_x(x)={\int }_{-\infty }^x{p}_x(t)dt$$

反过来，有：

$$p_x(x)=\frac{d}{dx}{P}_x(x)$$

CDF和PDF有一些非常实用的性质，例如

$${\int }_{-\infty }^{\infty }{p}_x(x)dx=1$$

即PDF在整个定义域上积分为1，其物理意义是所有可能的事件的概率和为1。

所谓的“均匀采样”，也就是说PDF得是常量。但是这个PDF必须是采样对象的，将采样对象的PDF列出来是解题的第一步。

随机变量之间的转换

假设随记变量$Y$和随机变量$X$满足：$Y_i=y(X_i)$，那么：

$$Pr(Y\le y(x))=Pr(X\le x)$$

因此：

$$\begin{array}{rl}{P}_y(y)=P_y(y(x))& =P_x(x)\\ p_y(y)\frac{dy}{dx}& =p_x(x)\\ p_y(y)& ={\left(\frac{dy}{dx}\right)}^{-1}{p}_x(x)\end{array}$$

如果随机变量$X$和$Y$是多维的，令$Y=T(X)$，那么$\frac{dy}{dx}$就变成了$T$的雅克比矩阵$J_T(x)$的行列式：

$$p_y(y)=\frac{p_x(x)}{|J_T(x)|}$$

雅克比矩阵定义为：

$$J_T(x)=\left(\begin{array}{ccc}\partial T_1/\partial x_1& \cdots & \partial T_1/\partial x_n\\ ⋮& \ddots & ⋮\\ \partial T_n/\partial x_1& \cdots & \partial T_n/\partial x_n\end{array}\right)$$

例如，对于对于极坐标：

$$\begin{array}{rl}x& =r\cos \theta \\ y& =r\sin \theta \end{array}$$

它的雅克比矩阵为：

$$J_T(x)=\left(\begin{array}{cc}\partial x/\partial r& \partial x/\partial \theta \\ \partial y/\partial r& \partial y/\partial \theta \end{array}\right)=\left(\begin{array}{cc}\cos \theta & -r\sin \theta \\ \sin \theta & r\cos \theta \end{array}\right)$$

那么它的行列式为$|J_T(x)|=r{\cos }^2\theta +r{\sin }^2\theta =r$。因此

$$p(r,\theta )=rp(x,y)$$

同理，对于求坐标：

$$\begin{array}{rl}x& =r\sin \theta \cos \varphi \\ y& =r\sin \theta \sin \varphi \\ z& =r\cos \theta \end{array}$$

它的雅克比矩阵的行列式$|J_T(x)|=r^2\sin \theta$，因此

$$p(r,\theta ,\varphi )=r^2\sin \theta p(x,y,z)$$

边缘概率

边缘概率相当于对多元概率分布函数进行降维，只考虑其中一维时的概率分布函数：

$$p(x)=\int p(x,y)dy$$

条件概率

假设$X$和$Y$为两个随机事件，而且$Pr(Y)>0$，则事件$Y$发生的条件下事件$X$发生的条件概率记为$Pr(X|Y)$，表示为：

$$Pr(X|Y)=\frac{Pr(X\cap Y)}{Pr(Y)}$$

他们的概率密度函数满足：

$$p(x|y)=\frac{p(x,y)}{p(y)}$$

立体角

另外需要补充的一点是，有时候我们会用立体角（solid angle）来计算光线，立体角$\omega$的微分形式是：

$$d\omega =sin\theta d\theta d\varphi$$

关于这个公式是怎么来的，可以参考[1]，在此不再赘述。

求解思路

通常来说，对于多变量的均匀采样问题，我们要求的是一个映射关系，能够把[0,1]上的均匀分布${\xi }_i$映射到要求的空间，并且满足均匀分布。解题思路是首先要用PDF来描述题目要求，并求出所有变量的联合PDF。然后先求出第一个变量的边缘PDF，再根据这个边缘概率分布求其他变量的条件PDF。最后用这些PDF求CDF，并令这些CDF分别等于${\xi }_i$，以求得各个变量与给定均匀分布的关系。

1. 用PDF来描述题目要求，求出所有变量的联合PDF
2. 求出第一个变量的边缘PDF
3. 根据已知的边缘PDF或条件PDF，依次求其余变量的条件PDF
4. 求所有变量的CDf，并令其分别等于${\xi }_i$
5. 求出各个变量和${\xi }_i$的关系

只说思路显得有些空洞，接下来看一下对一些常见问题的求解细节。

实践

单位圆盘均匀采样

问题：如何用[0,1]之间的均匀采样获得在单位圆盘上的均匀采样？

由于我们要在圆盘上面均匀采样，即$p(A)$是常数，其中A表示圆盘的面积。而$dA=dxdy$，所以$p(x,y)$是常数。考虑到PDF在定义域的积分为1，再加上半径$r=1$，那么求出

$$p(x,y)=\frac{1}{\pi }$$

对于圆盘采样来说，用极坐标系会更方便一些，因此我们将$p(x,y)$改写为$p(r,\theta )$：

$$p(r,\theta )=rp(x,y)=\frac{r}{\pi }$$

有了联合PDF，接下来先求边缘PDF：

$$p(r)={\int }_0^{2\pi }p(r,\theta )d\theta =2r$$

以及条件PDF

$$p(\theta |r)=\frac{p(r,\theta )}{p(r)}=\frac{1}{2\pi }$$

然后求他们的CDF：

$$\begin{array}{rl}P(r)& =r^2\\ P(\theta )& =\frac{\theta }{2\pi }\end{array}$$

令CDF等于均匀分布${\xi }_1$和${\xi }_2$：

$$\begin{array}{rl}r& =\sqrt{{\xi }_1}\\ \theta & =2\pi {\xi }_2\end{array}$$

这就是获得圆盘上均匀采样的方法。也可以将其转换成直角坐标系。

单位半球面均匀采样

问题：如何用[0,1]之间的均匀采样获得在单位半球面上的均匀采样？

对半球面的均匀采样同样是面积，但是与圆盘采样不同的是，半球面的面积可以用立体角$\omega$表示，即$p(\omega )$是常量，求得：

$$p(\omega )=\frac{1}{2\pi }$$

而同样的，我们得把$\omega$的PDF转换成球坐标基底r、$\theta$、$\varphi$的PDF，考虑到$d\omega =\sin d\theta d\varphi$，有：

$$p(\theta ,\varphi )=p(\omega )\sin \theta =\frac{\sin \theta }{2\pi }$$

有了PDF，接下来就球边缘PDF和条件PDF就可以了：

$$\begin{array}{rl}p(\theta )& ={\int }_0^{2\pi }p(\theta ,\varphi )d\varphi =\sin \theta \\ p(\varphi |\theta )& =\frac{p(\theta ,\varphi )}{p(\theta )}=\frac{1}{2\pi }\end{array}$$

以及CDF：

$$\begin{array}{rl}P(\theta )& ={\int }_0^{\theta }\sin {\theta }^{\prime }d{\theta }^{\prime }=1-\cos \theta \\ P(\varphi )& ={\int }_0^{\varphi }\frac{1}{2\pi }d{\varphi }^{\prime }=\frac{\varphi }{2\pi }\end{array}$$

求得：

$$\begin{array}{rl}\theta & =\arccos (1-{\xi }_1)\\ \varphi & =2\pi {\xi }_2\end{array}$$

如果把半球面改为球面，思路基本一致，只不过有几个数值会改变：

$$\begin{array}{rl}p(\omega )& =\frac{1}{4\pi }\\ p(\theta ,\varphi )& =p(\omega )\sin \theta =\frac{\sin \theta }{4\pi }\\ p(\theta )& ={\int }_0^{2\pi }p(\theta ,\varphi )d\varphi =\frac{\sin \theta }{2}\\ P(\theta )& =\frac{1-\cos \theta }{2}\\ \theta & =\arccos (1-2{\xi }_1)\end{array}$$

单位半球面余弦采样

问题：如何用[0,1]之间的均匀采样获得在单位球面上的余弦采样？

余弦采样指的是采样点在半球面分布的概率$p(\omega )$与角度$\theta$的余弦$\cos \theta$正相关，可以求出：

$$\begin{array}{rl}\int p(\omega )& =1\\ {\int }_0^{2\pi }{\int }_0^{\frac{\pi }{2}}c\cos \theta \sin \theta d\theta d\varphi & =1\\ c& =\frac{1}{\pi }\\ p(\theta ,\varphi )& =\frac{1}{\pi }\cos \theta \sin \theta \end{array}$$

接下来的求解思路与前几个例子一致，先求边缘PDF和条件PDF：

$$\begin{array}{rl}p(\theta )& ={\int }_0^{2\pi }p(\theta ,\varphi )d\varphi =2\cos \theta \sin \theta =\sin 2\theta \\ p(\varphi |\theta )& =\frac{p(\theta ,\varphi )}{p(\theta )}=\frac{1}{2\pi }\end{array}$$

以及CDF：

$$\begin{array}{rl}P(\theta )& ={\int }_0^{\theta }2\cos {\theta }^{\prime }\sin {\theta }^{\prime }d{\theta }^{\prime }=\frac{1-\cos 2\theta }{2}\\ P(\varphi )& ={\int }_0^{\varphi }\frac{1}{2\pi }d{\varphi }^{\prime }=\frac{\varphi }{2\pi }\end{array}$$

令CDF等于均匀分布${\xi }_1$和${\xi }_2$：

$$\begin{array}{rl}\theta & =\arccos (\sqrt{1-{\xi }_1})\\ \varphi & =2\pi {\xi }_2\end{array}$$

这个最终结果有点耐人寻味，如果我们定义$r=\sin \theta$，那么有：

$$r=\sqrt{{\xi }_1}$$

此时可以发现，$r$、$\varphi$的取值跟前面圆盘均匀采样是一致的。这个结果相当于，我们只需要进行一个圆盘均匀采样，然后把平面上的采样点映射到半球面上，得到的就是半球面余弦采样。这个小技巧可以用在具体的实现中，避免反三角函数运算。

球体均匀采样

问题：如何用[0,1]之间的均匀采样获得在球体上的均匀采样？

对半径为$R$的球体的均匀采样是体积采样，因此单位体素内的样本数量是常数：

$$p(x,y,z)=c$$

结合PDF的积分为1，求得：

$$p(x,y,z)=\frac{3}{4\pi R^3}$$

把直角坐标系转换为球坐标系：

$$\begin{array}{rl}p(r,\theta ,\varphi )=|J_T|p(x,y,z)& =r^2\sin \theta \cdot \frac{3}{4\pi R^3}\\ & =\frac{3\sin \theta r^2}{4\pi R^3}\end{array}$$

然后求边缘PDF和条件PDF：

$$\begin{array}{rl}p(r)& ={\int }_0^{2\pi }{\int }_0^{\pi }p(r,\theta ,\varphi )d\theta d\varphi =\frac{3r^2}{R^3}\\ p(\theta |r)& =\frac{{\int }_0^{2\pi }p(r,\theta ,\varphi )d\varphi }{p(r)}=\frac{\sin \theta }{2}\\ p(\varphi |r,\theta )& =\frac{p(r,\theta ,\varphi )}{p(\theta |r)p(r)}=\frac{1}{2\pi }\end{array}$$

以及CDF：

$$\begin{array}{rl}P(r)& =\frac{r^3}{R^3}\\ P(\theta )& =\frac{1-\cos \theta }{2}\\ P(\varphi )& =\frac{\varphi }{2\pi }\end{array}$$

令CDF等于单位均匀分布${\xi }_1$、${\xi }_2$、${\xi }_3$：

$$\begin{array}{rl}r& =R\cdot \sqrt[3]{{\xi }_1}\\ \theta & =\arccos (1-2{\xi }_2)\\ \varphi & =2\pi {\xi }_3\end{array}$$

这个结果与球面均匀采样和圆盘均匀采样有一定的相似性。

三角形均匀采样

问题：如何用[0,1]之间的均匀采样获得在任意三角形内的均匀采样？

对三角形的均匀采样看起来比较容易，但是推导起来反而最麻烦。当然，也有一种简单的方法，比如对三角形的矩形凸包进行采样，然后将三角形外部的点映射到三角形内部（沿一条边的中点做对称），这个方法也可以获得对三角形的均匀采样，但是这个方法有几个问题

1. 需要构建一个凸包，这个操作对平面三角形还不麻烦，但是对三维空间的三角形就相对繁琐；
2. 需要对一半的点做映射；

我们希望避免这些操作，仍旧用前面提到的方法进行分析。

首先对问题进行细化。假设任意三角形$\Delta ABC$，为了方便起见，我们将坐标轴原点设置在$A$所在的位置，并令$x$轴和$AC$重合，那么对于三角形内的任意点$P(x,y)$，它的坐标$(x,y)$可以写为：

$$\begin{array}{rl}(x,y)& =u\overrightarrow{AB}+v\overrightarrow{AC}\\ x& =u||AB\cos A||+v||AC||\\ y& =u||AB\sin A||\end{array}$$

其中，$u$和$v$满足$1-u-v>0$。

根据坐标之间的关系，求出雅克比矩阵$J_T$：

$$\begin{array}{rl}{J}_T& =\frac{\partial (x,y)}{\partial (u,v)}=\left(\begin{array}{cc}||AB\cos A||& ||AC||\\ ||AB\sin A||& 0\end{array}\right)\\ |J_T|& =||AB\cdot AC\sin A||\end{array}$$

再回到采样问题上。对面积的均匀采样，意味着定义域内$p(x,y)$是常数，即

$$p(x,y)=\frac{1}{S}=\frac{2}{||AB\cdot AC\sin A||}$$

接下来，把$p(x,y)$改写为$p(u,v)$的形式

$$\begin{array}{rl}p(u,v)& =p(x,y)\ast |J_T|=2\end{array}$$

这个结果倒是有些意外地简洁。然后求边缘PDF和条件PDF：

$$\begin{array}{rlr}p(u)& ={\int }_0^{1-u}p(u,v)dv=2(1-u)p(v|u)& =\frac{p(u,v)}{p(u)}=\frac{1}{1-u}\end{array}$$

这里有一个细节，在求$p(u)$的时候，我们需要对$p(u,v)$在$v$上做积分，此时$v$的取值范围是$[0,1-u]$，而不是$[0,1]$。

然后求CDF：

$$\begin{array}{rl}P(u)& ={\int }_0^{u}p(u^{\prime })d{u}^{\prime }=2u-u^2\\ P(v)& ={\int }_0^{v}p(v^{\prime }|u)d{v}^{\prime }=\frac{v}{1-u}\end{array}$$

令CDF等于单位均匀分布${\xi }_1$、${\xi }_2$：

$$\begin{array}{rl}u& =1-{\sqrt{\xi }}_1\\ v& ={\xi }_2{\sqrt{\xi }}_1\end{array}$$

这样，利用${\xi }_1$、${\xi }_2$求得了$u$、$v$之后，再求出$x$和$y$，就能获得一个采样点了。

$$P=(1-{\sqrt{\xi }}_1)A+({\xi }_2{\sqrt{\xi }}_1)B+({\sqrt{\xi }}_1(1-{\xi }_2))C$$

注意到这个公式与$(x,y)$无关，因此可以简单地推广到三维甚至高维。

参考资料

1. 立体角(Solid Angle)详解
2. Physically based rendering, 3rd edition.
3. StackExchange上对这个问题的讨论：q18686，但是这个讨论里面的答案都是从特例证明这个结论的正确性，个人觉得证明不够完善。