[图形学]在半球面上均匀采样和cos加权采样

一、简介

本文介绍了如何在半球表面上进行半球面上的均匀采样、半球面上的cos加权采样。
在介绍半球面上的cos加权采样时，介绍了 在球面上采样（方法一）和在投影圆内采样（方法二）两种方法。
并且给出了相关公式推导和python代码实现。

二、在半球上采样

0.预备知识

1).球面坐标系与笛卡尔坐标系

在半球面上采样时，常使用球面坐标系。先采样球面坐标系下的坐标参数$(\varphi ,\theta )$，然后转换为笛卡尔坐标系参数$(x,y,z)$。两个坐标系下的参数示意图如下：

笛卡尔坐标系示意图：

球面坐标系示意图：

从球面坐标系转为笛卡尔坐标系的公式如下：
$$\begin{gathered} x=cos(\varphi )\times sin(\theta ) \\ y=cos(\theta ) \\ z=sin(\varphi )\times sin(\theta )\text{\hspace{1em}(1)} \end{gathered}$$
在单位球面坐标系上的立体角微元等于：
$$\mathrm{d}\omega =\frac{\mathrm{d}\mathrm{A}}{{\mathrm{r}}^2}=\mathrm{d}\mathrm{A}\text{\hspace{1em}(2)}$$
在单位球面上的面积微元$\mathrm{d}\mathrm{A}$等于：
$$\mathrm{d}\mathrm{A}={\mathrm{r}}^2\times \sin (\theta )\mathrm{d}\theta \mathrm{d}\varphi =\sin (\theta )\mathrm{d}\theta \mathrm{d}\varphi \text{\hspace{1em}(3)}$$

2).逆变换采样

逆变换采样是伪随机数采样的一种基本方法，在已知任意概率分布函数(probability distribution function, PDF)的累计分布函数(cumulative distribution function, CDF)时，可以从该PDF中生成随机样本。本文只简单介绍逆变换采样的步骤和示例，对其推导过程和原理不做过多介绍。

使用逆变换采样方法的步骤如下：

• (1).根据已知的PDF函数$f(x)=y$，计算对应的CDF函数 $F(X)$；
• (2).根据 CDF函数$F(X)$ 求出其反函数 $F^{-1}(Y)$；
• (3).在[0,1]上均匀采样随机变量$u$；
• (4).那么$F^{-1}(u)$即为满足PDF函数$f(x)=y$的随机变量。

示例：

已知PDF函数为$f(x)=\frac{1}{2\sqrt{x}}\exp (-\sqrt{x}),x\ge 0$，求采样得到满足$f(x)$分布的随机变量。

• (1).PDF为$f(x)=\frac{1}{2\sqrt{x}}\exp (-\sqrt{x}),x\ge 0$，那么CDF为$F(X)=1-exp(-\sqrt{X})$；
• (2).CDF对应的逆函数为$F^{-1}(Y)=(log(1-Y){)}^2$；
• (3).在[0,1]上均匀采样随机变量$u$；
• (4).那么随机变量$x^{\prime }=F^{-1}(u)=(log(1-u){)}^2$即满足目标PDF函数；

1.在半球面上均匀采样 uniform sample

对于半球面上的均匀采样，PDF应与微元面积成正比，那么对于$\varphi$和$\theta$的累计密度函数应该为：
$$\begin{gathered} CDF(\varphi )=\frac{{\int }_0^{\varphi }{\int }_0^{\pi /2}sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }{{\int }_0^{2\pi }{\int }_0^{\pi /2}sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }=\frac{{\int }_0^{\varphi }{\int }_0^{\pi /2}sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }{2\pi }=\frac{\varphi }{2\pi } \\ CDF(\theta )=\frac{{\int }_0^{2\pi }{\int }_0^{\theta }sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }{{\int }_0^{2\pi }{\int }_0^{\pi /2}sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }=\frac{{\int }_0^{2\pi }{\int }_0^{\theta }sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }{2\pi }=1-cos(\theta )\text{\hspace{1em}(4)} \end{gathered}$$
对应的逆函数分别为：
$$\begin{gathered} CD{F}^{-1}({\varphi }^{\prime })=2\pi \times {\varphi }^{\prime } \\ CD{F}^{-1}({\theta }^{\prime })=arccos(1-{\theta }^{\prime })\text{\hspace{1em}(5)} \end{gathered}$$
之后在[0,1]之间均匀采样随机变量$u1,u2$，然后令：
$$\begin{gathered} \varphi =CD{F}^{-1}(u1)=2\pi \times u1 \\ \theta =CD{F}^{-1}(u2)=arcos(1-u2)\text{\hspace{1em}(6)} \end{gathered}$$
根据上式计算得到的$(\varphi ,\theta )$就满足在半球面上按照面积均匀分布。

下面是采样python代码和结果（需要注意代码中笛卡尔坐标系的Y轴和Z轴做了交换）：

代码：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def uniform_sample_hemisphere(n_samples):
    points = []
    for _ in range(n_samples):
        u1 = np.random.rand()
        u2 = np.random.rand()
        
        phi = 2 * np.pi * u1   # 计算方位角
        theta = np.arccos(1.0-u2)  # 计算极角

        # 转换到笛卡尔坐标系
        x = np.sin(theta) * np.cos(phi)
        y = np.sin(theta) * np.sin(phi)
        z = np.cos(theta)

        points.append((x, y, z))
    return np.array(points)


# 生成样本
n_samples = 5000
samples = uniform_sample_hemisphere(n_samples)

# 3D 可视化
fig = plt.figure(figsize=(12, 6))

# 绘制 3D 采样
ax1 = fig.add_subplot(121, projection='3d')
ax1.scatter(samples[:, 0], samples[:, 1], samples[:, 2], alpha=0.6, s=1)
ax1.set_xlabel('X axis')
ax1.set_ylabel('Y axis')
ax1.set_zlabel('Z axis')
ax1.set_title('Uniform Sampling on Hemisphere')

# 绘制 x-y 平面映射
ax2 = fig.add_subplot(122)
ax2.scatter(samples[:, 0], samples[:, 1], alpha=0.6, s=1)
ax2.set_xlabel('X axis')
ax2.set_ylabel('Y axis')
ax2.set_title('Projection onto XY Plane')
ax2.axis('equal')

plt.tight_layout()
plt.show()

结果：

2.在半球面上的cos加权采样：方法一

对于半球面上的cos加权采样，PDF应与$微元面积\times cos(\theta )$成正比，即在$\theta$越小的区域采样概率越大，在$\theta$越大的区域概率越小，那么对于$\varphi$和$\theta$的累计密度函数应该为：
$$\begin{gathered} CDF(\varphi )=\frac{{\int }_0^{\varphi }{\int }_0^{\pi /2}cos(\theta )sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }{{\int }_0^{2\pi }{\int }_0^{\pi /2}cos(\theta )sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }=\frac{{\int }_0^{\varphi }{\int }_0^{\pi /2}cos(\theta )sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }{\pi }=\frac{\varphi /2}{\pi }=\frac{\varphi }{2\pi } \\ CDF(\theta )=\frac{{\int }_0^{2\pi }{\int }_0^{\theta }cos(\theta )sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }{{\int }_0^{2\pi }{\int }_0^{\pi /2}cos(\theta )sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }=\frac{{\int }_0^{2\pi }{\int }_0^{\theta }cos(\theta )sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }{\pi }=\frac{\pi \times si{n}^2(\theta )}{\pi }=si{n}^2(\theta )=1-co{s}^2(\theta )\text{\hspace{1em}(7)} \end{gathered}$$
对应的逆函数分别为：
$$\begin{gathered} CD{F}^{-1}({\varphi }^{\prime })=2\pi \times {\varphi }^{\prime } \\ CD{F}^{-1}({\theta }^{\prime })=arccos(\sqrt{1-{\theta }^{\prime }})\text{\hspace{1em}(8)} \end{gathered}$$
之后在[0,1]之间均匀采样随机变量$u1,u2$，然后令：
$$\begin{gathered} \varphi =CD{F}^{-1}(u1)=2\pi \times u1 \\ \theta =CD{F}^{-1}(u2)=arcos(\sqrt{1-u2})\text{\hspace{1em}(9)} \end{gathered}$$
根据上式计算得到的$(\varphi ,\theta )$就满足在半球面上按照$微元面积\times cos(\theta )$加权分布。

下面是采样python代码和结果（需要注意代码中笛卡尔坐标系的Y轴和Z轴做了交换）：

代码：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def cos_weighted_sample_hemisphere(n_samples):
    points = []
    for _ in range(n_samples):
        u1 = np.random.rand()
        u2 = np.random.rand()
        phi = 2 * np.pi * u1   # 计算方位角
        theta = np.arccos(np.sqrt(1 - u2))  # 计算极角

        # 转换到笛卡尔坐标系
        x = np.sin(theta) * np.cos(phi)
        y = np.sin(theta) * np.sin(phi)
        z = np.cos(theta)

        points.append((x, y, z))
    return np.array(points)


# 生成样本
n_samples = 5000
samples = cos_weighted_sample_hemisphere(n_samples)

# 3D 可视化
fig = plt.figure(figsize=(12, 6))

# 绘制 3D 采样
ax1 = fig.add_subplot(121, projection='3d')
ax1.scatter(samples[:, 0], samples[:, 1], samples[:, 2], alpha=0.6, s=1)
ax1.set_xlabel('X axis')
ax1.set_ylabel('Y axis')
ax1.set_zlabel('Z axis')
ax1.set_title('Cos-Weighted Sampling on Hemisphere')

# 绘制 x-y 平面映射
ax2 = fig.add_subplot(122)
ax2.scatter(samples[:, 0], samples[:, 1], alpha=0.6, s=1)
ax2.set_xlabel('X axis')
ax2.set_ylabel('Y axis')
ax2.set_title('Projection onto XY Plane')
ax2.axis('equal')

plt.tight_layout()
plt.show()

结果：

3.在半球面上的cos加权采样：方法二

另一种对于半球面上的cos加权采样的方法如下：
直接在X-Z平面的圆内均匀采样点，然后将采样的点映射到上方的球面上，如下如图所示：

在采样时可以使用接受拒绝采样算法，在一个的大小为1*1的正方形内均匀采样，然后只保留距离圆心小于等于1的采样点。
或者使用前面逆变换采样方法，分别采样半径$r$和角度$\varphi$。其中采样角度$\varphi$的过程与前面相同。采样半径$r$时，先计算$r$的CDF。因为在圆面内面积微元$\mathrm{d}\mathrm{A}$为：
$$\mathrm{d}\mathrm{A}=\mathrm{r}\ast \mathrm{d}\mathrm{r}\mathrm{d}\varphi \text{\hspace{1em}(10)}$$
在圆平面内均匀采样一点，那么采样点对应的半径$r$的$CDF(r)$等于：
$$CDF(r)=\frac{{\int }_0^{2\pi }{\int }_0^{r}r\mathrm{d}\mathrm{r}\mathrm{d}\varphi }{{\int }_0^{2\pi }{\int }_0^{1}r\mathrm{d}\mathrm{r}\mathrm{d}\varphi }=\frac{2\pi \ast r^2/2}{\pi }=r^2\text{\hspace{1em}(11)}$$

$$CD{F}^{-1}(u)=\sqrt{u}\text{\hspace{1em}(12)}$$
因此可以在[0,1]上均匀采样随机变量$u$，那么根据公式(12)，令$r=\sqrt{u}$。根据此方法得到的采样点$(r,\varphi )$在圆平面内即为均匀分布的。

证明：

接下来将证明本方法与在半球面上的cos加权采样：方法一的等效性。
我们将目标放在X-Z平面上的圆平面内，在该圆面内可以使用极坐标$(r,\varphi )$表示任意一点。
因为：$CDF(\theta )=PDF(\Theta \le \theta )$
同时在极坐标中：$r=sin(\theta )$，即$R=sin(\Theta )$
那么：$PDF(\Theta \le \theta )=PDF(R\le sin(\theta ))$
又因为：$r=sin(\theta )$
因此：$PDF(R\le sin(\theta ))=PDF(R\le r)=CDF(r)=r^2=si{n}^2(\theta )$
即：$CDF(\theta )=si{n}^2(\theta )$
即映射点$(\theta ,\varphi )$中的坐标值$\theta$，满足$CDF(\theta )=si{n}^2(\theta )$，而这正是公式(7)中我们想要找到一个随机变量$\theta$。
同理可得，使用该方法得到的映射点坐标$\varphi$也满足公式(7)中$CDF(\varphi )=\frac{\varphi }{2\pi }$

下面是采样python代码和结果（需要注意代码中笛卡尔坐标系的Y轴和Z轴做了交换）：

代码：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 生成 cos 加权采样的函数
def cos_weighted_sample(n_samples):
    points = []

    for _ in range(n_samples):
        u1 = np.random.rand()
        u2 = np.random.rand()

        r = np.sqrt(u1)  # 均匀采样半径
        phi = 2 * np.pi * u2   # 计算方位角  # 均匀采样角度
        x = r * np.cos(phi)
        y = r * np.sin(phi)
        z = np.sqrt(1 - x**2 - y**2)
        if (x*x+y*y <= 1.0):  # 只保留圆内的点
            points.append((x, y, z))
    return np.array(points)

# 生成样本
n_samples = 5000
samples = cos_weighted_sample(n_samples)

# 3D 可视化
fig = plt.figure(figsize=(12, 6))

# 绘制 3D 采样
ax1 = fig.add_subplot(121, projection='3d')
ax1.scatter(samples[:, 0], samples[:, 1], samples[:, 2], alpha=0.6, s=1)
ax1.set_xlabel('X axis')
ax1.set_ylabel('Y axis')
ax1.set_zlabel('Z axis')
ax1.set_title('Cos-Weighted Sampling on Hemisphere')

# 绘制 x-y 平面映射
ax2 = fig.add_subplot(122)
ax2.scatter(samples[:, 0], samples[:, 1], alpha=0.6, s=1)
ax2.set_xlabel('X axis')
ax2.set_ylabel('Y axis')
ax2.set_title('Projection onto XY Plane')
ax2.axis('equal')

plt.tight_layout()
plt.show()

结果：

三、参考

[1].Sampling the hemisphere
[2].Cosine-Weighted采样算法