一、简介

本文介绍了如何在半球表面上进行半球面均匀采样、半球面cos加权采样采样。
给出了相关公式推导和python代码实现。

二、在半球上采样

0.预备知识

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

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

笛卡尔坐标系示意图：

球面坐标系示意图：

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

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 ) \end{gathered}$$
对应的逆函数分别为：
$$\begin{gathered} CD{F}^{-1}({\varphi }^{\prime })=2\pi \times {\varphi }^{\prime } \\ CD{F}^{-1}({\theta }^{\prime })=arccos(1-{\theta }^{\prime }) \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) \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 ) \end{gathered}$$
对应的逆函数分别为：
$$\begin{gathered} CD{F}^{-1}({\varphi }^{\prime })=2\pi \times {\varphi }^{\prime } \\ CD{F}^{-1}({\theta }^{\prime })=arccos(\sqrt{1-{\theta }^{\prime }}) \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}) \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 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(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 = 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()

结果

三、参考

[1].Sampling the hemisphere