计算机图形学基础-13-采样

《Fundamentals of Computer Graphics》5th（计算机图形学基础/虎书），中文翻译。

第 13 章 Sampling 采样

许多图形应用程序需要对不寻常空间进行“公正”的采样，例如所有可能的直线空间。例如，在一个像素内生成随机边缘，或在像素上随机采样点，其密度根据某个密度函数而变化。本章提供了这些概率操作所需的基础知识。这些技术也将被用于使用蒙特卡罗积分 (Monte Carlo integration) 来数值评估复杂积分（也在本章中介绍）。

13.1 积分

尽管“积分”和“测度”这两个词经常显得令人生畏，但它们与数学中最直观的概念之一相关，并且不应该被害怕。对于我们非常不严格的目的而言，测度只是一个将子集映射到 $\mathbb{R}^{+}$ 的函数，其方式与我们对长度、面积和体积的直觉概念一致。例如，在二维实平面 $\mathbb{R}^{2}$ 上，我们有面积测度 $A$ ，它为平面上一组点分配一个值。请注意，A 只是一个函数，它将平面的部分映射到面积。这意味着 $A$ 的定义域是 $\mathbb{R}^{2}$ 的所有可能子集，我们将其表示为幂集 $\mathcal{P}(\mathbb{R}^{2})$ 。因此，我们可以使用箭头符号来描述 A：

A: \mathcal{P}\left(\mathbb{R}^{2}\right) \rightarrow \mathbb{R}^{+}

将面积测度应用于一个正方形的示例，可以得到边长为 1 的正方形的面积为 1：

A([a,a+1]×[b,b+1]) = 1

其中 $(a, b)$ 是正方形左下角的坐标。请注意，单个点（如 $(3,7)$ ）是 $\mathbb{R}^{2}$ 的有效子集，并且具有零面积： $A((3,7)) = 0$ 。同样，点集 S 在 x 轴上， $S = {(x,y)∈ \mathbb{R}^{2}}$ 并且 $y=0$ ，即 $A(S)=0$ 。这些集合被称为零测度集合。

为了被视为测度，函数必须遵守某些类似于面积的属性。例如，我们有一个函数 $\mu: \mathcal{P}(\mathbb{S}) \rightarrow \mathbb{R}^{+}$ 。对于 $\mu$ 成为测度，以下条件必须成立：

空集的测度为零： $\mu(∅) = 0$ 。
两个不同集合的测度之和等于它们各自测度的总和减去它们的交集的测度。即 $\mu(A \cup B)=\mu(A)+\mu(B)-\mu(A \cap B)$ ，其中 $∪$ 是集合并运算符， $∩$ 是集合交运算符。

当实际计算测度时，通常使用积分。我们可以将积分视为只是符号：

A(S) \equiv \int_{x \in S} d A(\mathbf{x})

您可以非正式地将右侧解读为“取在区域 $S$ 中的所有点 $\text x$ ，并汇总它们关联的微分面积。” 积分通常可以用其他方式表示，包括

\int_{S} d A, \quad \int_{\mathbf{X} \in S} d \mathbf{x}, \quad \int_{\mathbf{X} \in S} d A_{\mathbf{X}}, \quad \int_{\mathbf{X}} d \mathbf{x}

上述所有公式都表示“区域 $S$ 的面积”。我们将使用我们最初使用的那个公式，因为它非常冗长，避免了歧义。要分析地评估这样的积分，通常需要建立一些坐标系，并使用微积分技巧来解方程。但是如果你忘记了这些技巧，也不用担心，因为我们通常需要对积分进行数值近似，而这只需要使用第 13.3 节中描述的简单技术。

给定集合 $\mathbb{S}$ 上的一个测度，我们可以通过用非负函数 $w: \mathbb{S} \rightarrow \mathbb{R}^{+}$ 进行加权来创建一个新的测度。这最好以积分符号表示。例如，我们可以从 $[0,1]^2$ 中的简单面积测度开始：

\int_{\mathbf{x} \in[0,1]^{2}} d A(\mathbf{x})

然后我们可以插入一个半径平方的加权函数来使用“径向加权”测度：

\int_{\mathbf{x} \in[0,1]^{2}}\|\mathbf{x}\|^{2} d A(\mathbf{x})

要分析地评估这个积分，我们可以使用笛卡尔坐标系展开，其中 $dA ≡ dx dy$ ：

\int_{\mathbf{X} \in[0,1]^{2}}\|\mathbf{x}\|^{2} d A(\mathbf{x})=\int_{x=0}^{1} \int_{y=0}^{1}\left(x^{2}+y^{2}\right) d x d y

关键在于，如果你认为 $||x||^2$ 项与 $dA$ 项相乘，那么它们共同形成了一个新的测度，我们可以将该测度称为 $ν$ 。这样，我们可以写出 $ν(S)$ 而不是整个积分。如果你认为这只是一堆符号和簿记，你是正确的。但是这确实使我们能够根据我们的喜好编写紧凑或展开的方程式。

13.1.1 测度和平均值

当计算函数的平均值时，测度真正开始发挥作用。你只能根据特定的测度取平均值，并且希望选择一个在应用或领域中“自然”的测度。一旦选择了测度，函数 $f$ 在区域 $S$ 上相对于测度 $μ$ 的平均值为

\operatorname{average}(f) \equiv \frac{\int_{x \in S} f(\mathbf{x}) d \mu(\mathbf{x})}{\int_{x \in S} d \mu(\mathbf{x})}

例如，对于函数 $f(x,y) = x^2$ ，在 $[0,2]^2$ 上相对于面积测度的平均值为

\operatorname{average}(f) \equiv \frac{\int_{x=0}^{2} \int_{y=0}^{2} x^{2} d x d y}{\int_{x=0}^{2} \int_{y=0}^{2} d x d y}=\frac{4}{3}

这种方法有助于解决看似困难的问题，其中选择测度是棘手的部分。这些问题经常出现在积分几何学中，这是一门研究几何实体（如线条和平面）上的测度的学科。例如，人们可能想知道通过 $[0,1]^2$ 的直线的平均长度。也就是说，根据定义，

\text { Average(length) }=\frac{\int_{\text {lines } L \text { through }[0,1]^{2}} \text { length }(L) d \mu(L)}{\int_{\text {lines } L \text { through }[0,1]^{2}} d \mu(L)}

一旦我们知道这个公式，剩下的就是为应用选择适当的 $μ$ 。这在下一节中处理线条时会涉及。

13.1.2 例子：二维平面上的线条测度

哪个测度 $μ$ 是“自然”的？

如果你将直线参数化为 $y=mx+b$ ，你可以将给定的直线视为“斜截式”空间中的一点 $(m,b)$ 。一个易于使用的测度是 $dm \space db$ ，但这不是一个“好的”测度，因为不是所有相等大小的“束”都具有相同的测度。更准确地说，该测度在坐标系统变换下不具有不变性。例如，如果你考虑通过正方形 $[0,1]^2$ 上的所有直线，则它们的测度与通过旋转 45° 的单位正方形上的直线的测度不同。我们真正想要的是一种“公平”的测度，它不会随着一组直线的旋转或平移而改变。这个想法在图 13.1 和 13.2 中说明了。

图 13.1。这两束直线应该具有相同的测度。它们与 y 轴的交点长度不同，因此使用 $db$ 作为微分测度将是一个很差的选择。

图 13.2。这两束直线应该具有相同的测度。由于它们的斜率变化值不同，因此使用 $dm$ 作为微分测度将是一个很差的选择。

要在直线上开发一种自然的测度，我们首先应该将它们视为对偶空间中的点来思考。这是一个简单的概念：直线 $y=mx+b$ 可以被指定为斜截式空间中的点 $(m,b)$ 。这个概念在图 13.3 中说明了。在 $(ϕ,b)$ 空间中开发测度更为直接。在那个空间里， $b$ 是 $y-$ 截距，而 $ϕ$ 是直线相对于 x 轴的夹角，如图 13.4 所示。在这种情况下，微分测度 $dϕ \space db$ 几乎可以工作，但由于图 13.1 中所示的影响，它并不公平。为了解决常数宽度直线束产生的较大跨度 $b$ 的影响，我们必须添加一个余弦因子：

dμ = cos \space ϕ \space dϕ \space db。

图 13.3。在 $(x,y)$ 空间中，直线 $y=mx+b$ 上的点集也可以用 $(m,b)$ 空间中的单个点表示，因此上方的直线和下方的点表示相同的几何实体：一个二维直线。

图 13.4。在角度截距空间中，我们使用角度 $ϕ ∈ [-π∕2,π∕2)$ 来参数化直线，而不是斜率。

可以证明，这种测度是旋转和平移不变性的唯一测度，除了一个常数。

这种测度可以转换为其他直线参数化的适当测度。例如， $(m,b)$ 空间的适当测度为

d \mu=\frac{d m d b}{\left(1+m^{2}\right)^{\frac{3}{2}}}

对于在 $(u,v)$ 空间中参数化的直线，

ux+vy+ 1 = 0，

适当的测度为

d \mu=\frac{d u d v}{\left(u^{2}+v^{2}\right)^{\frac{3}{2}}}

对于以 $(a,b)$ 参数化的直线，即 $x-$ 截距和 $y-$ 截距，测度为

d \mu=\frac{a b d a d b}{\left(a^{2}+b^{2}\right)^{\frac{3}{2}}}

注意，这些空间中的任何一个都是指定直线的同样有效的方式，哪一个更好取决于具体情况。然而，人们可能想知道是否存在一种坐标系，其中直线集的测度只是对偶空间中的面积。事实上，存在这样的坐标系，它非常简单；它是法线坐标系，它用法线距离和法线相对于 x 轴的夹角指定一条直线（图 13.5）。这些直线的隐式方程为

x \space cosθ + y \space sinθ − p = 0

图 13.5。法线坐标系使用法线距离和角度来指定一条直线。

而且，实际上在该空间中的测度为

dμ = dp \space dθ

我们将在稍后的一节中使用这些测度来选择公平的随机直线。

13.1.3 例子：三维空间中的直线测度

在三维空间中，有许多方法可以参数化直线。也许，最简单的方法是使用它们与特定平面的交点以及它们方向的某些规范。例如，我们可以记录与 xy 平面的交点以及其方向的球面坐标。因此，每条直线都可以被指定为一个 $(x,y,θ,ϕ)$ 四元组。这表明，三维空间中的直线是四维实体；即它们可以用一个四维空间中的点来描述。

一条直线的微分测度不应随 $(x,y)$ 变化，但具有相等横截面的直线束应具有相等的测度。因此，一种公平的微分测度是

dμ = dx \space dy \space sin⁢ θ \space dθ \space dϕ

另一种参数化直线的方法是通过图表两个平行平面的交点。例如，如果直线与平面 $z=0$ 在 $(x=u,y=v)$ 相交，并且与平面 $z=1$ 在 $(x=s,y=t)$ 相交，则可以通过四元组 $(u,v,s,t)$ 来描述该直线。请注意，与以前的参数化方式一样，对于与 $xy$ 平面平行的直线，这种方式也是退化的。尽管这种参数化的微分测度更为复杂，但它可以近似为

dμ≈du⁢ \space dv⁢ \space a⁢ \space ds⁢ \space dt，

对于几乎与 z 轴平行的直线束。这是在基于图像的渲染中经常隐式使用的测度。

对于与球相交的直线集，我们可以使用直线与球相交的两个点的参数化。如果这些点是用球面坐标表示的，则该点可以用四元组 $(θ1,ϕ1,θ2,ϕ2)$ 表示，测度就是与每个点相关的微分面积：

dμ=sinθ1 \space dθ1 \space dϕ1 \space sinθ2 \space dθ2 \space dϕ2。

这意味着在球面上随机选择两个均匀的端点会得到等密度的直线。Mateu Sbert 在他的论文中使用这个发现来计算形状因子（Sbert，1997）。

请注意，有时我们想要参数化定向直线，有时我们不希望端点的顺序影响测度。这是一个簿记细节，在渲染应用程序中特别重要，因为沿着一条直线流动的光量在沿着该直线的两个方向上是不同的。

13.2 连续概率

许多图形算法使用概率构建随机样本以解决积分和平均问题。这是应用连续概率的领域，它与测度理论有基本联系。

13.2.1 一维连续概率密度函数

粗略地说，连续随机变量 $x$ 是一个随机从实数线 $ℝ=(-∞,+∞)$ 中取某个值的标量或向量量。 $x$ 的行为完全由它所取值的分布描述。这些值的分布可以通过与 $x$ 相关的概率密度函数 (pdf)， $p$ ，来定量描述（这种关系表示为 $x$ ~ $p$ ）。x 在某个区间 $[a,b]$ 上取特定值的概率由以下积分给出：

\operatorname{Probability}(x \in[a, b])=\int_{a}^{b} p(x) d x \tag{13.1}

粗略地说，概率密度函数 $p$ 描述了随机变量取某个值的相对可能性；如果 $p(x_1)=6.0$ 而 $p(x2)=3.0$ ，则具有密度 $p$ 的随机变量在“接近” $x_1$ 的值上比在 $x2$ 的值上两倍可能。密度 $p$ 具有两个特征：

p(x)≥0（概率非负） \tag{13.2}

\int_{-\infty}^{+\infty} p(x) d x=1 \quad(\operatorname{Probability}(x \in \mathbb{R})=1) \tag{13.3}

例如，典型的随机变量 $ξ$ 取值在 $[0,1)$ 区间内，概率均匀（这里的“均匀”意味着 $ξ$ 的每个值等可能）。这意味着 $ξ$ 的概率密度函数 $q$ 为

q(\xi)=\left\{\begin{array}{ll} 1 & \text { if } 0 \leq \xi<1 \\ 0 & \text { otherwise } \end{array}\right.

$ξ$ 定义的空间只是区间 $[0,1)$ 。 $ξ$ 取某个值在区间 $[a,b] ∈ [0,1)$ 内的概率为

\operatorname{Probability}(a \leq \xi \leq b)=\int_{a}^{b} 1 d x=b-a .

13.2.2 一维期望值

对于具有基础概率密度函数 $p$ 的一维随机变量的实数函数 f，它将取平均值的期望值称为它的期望值 $E(f(x))$ （有时写作 $Ef(x)$ ）：

E(f(x))=\int f(x) p(x) d x

一维随机变量的期望值可以通过设置 $f(x) = x$ 来计算。期望值具有出人意料且有用的性质：两个随机变量的和的期望值是这些变量的期望值之和：

E(x+y) =E(x) +E(y)

对于随机变量 $x$ 和 $y$ 。因为随机变量的函数本身就是随机变量，所以这种期望的线性性也适用于它们：

E(f(x) +g(y)) =E(f(x)) +E(g(y))

一个显然的问题是，如果进行加法的随机变量是相关的（不相关的变量称为独立），这个性质是否仍然成立。事实上，这个线性性质无论这些变量是否独立都成立！这个求和性质对于大多数蒙特卡罗应用程序至关重要。

13.2.3 多维随机变量

关于随机变量及其期望值的讨论自然地扩展到多维空间。大多数图形问题都将在这样的高维空间中进行讨论。例如，许多光照问题是在半球面上提出的。幸运的是，如果我们在随机变量所占据的空间上定义一个测度 $μ$ ，一切都与一维情况非常相似。假设空间 $S$ 具有相关的测度 $μ$ ；例如， $S$ 是一个球体的表面，而 $μ$ 用于测量面积。我们可以定义一个概率密度函数 $p：S↦ℝ$ ，如果 $x$ 是随机变量且 $x ~ p$ ，则 $x$ 将在某个区域 $S_i ⊂ S$ 上取某个值的概率由积分

\operatorname{Probability}\left(x \in S_{i}\right)=\int_{S_{i}} p(x) d \mu

给出。这里， $Probability (event)$ 是事件为真的概率，因此积分是 x 取某个值在区域 Si 内的概率。

在图形学中， $S$ 经常是一个区域 $(dμ=dA=dxdy)$ 或者一组方向（单位球面上的点： $dμ=dω=sinθdθdϕ$ ）。例如，一个二维随机变量 $α$ 是一个均匀分布在半径为 $R$ 的圆盘上的随机变量。这里，均匀意味着相对于面积是均匀的，例如，一个不好的飞镖手在飞镖板上的命中点分布方式。由于它是均匀分布的，我们知道 $p(α)$ 是某个常数。从圆盘的面积为 $πr^2$ 以及总概率为一的事实中，我们可以推断出

p(\alpha)=\frac{1}{\pi R^{2}} .

这意味着 $α$ 在圆盘的某个子集 $S_1$ 中的概率只是

\operatorname{Probability}\left(\alpha \in S_{1}\right)=\int_{S_{1}} \frac{1}{\pi R^{2}} d A \text {. }

这一切都很抽象。要实际使用这些信息，我们需要以可以评估的形式表示积分。假设 $S_i$ 是靠近中心而不是周边的圆盘部分。如果我们转换为极坐标，则 $α$ 可以表示为 $(r,ϕ)$ 对，而 $S1$ 是 $r < R/2$ 的区域。请注意，仅仅因为 $α$ 是均匀的，并不意味着 $ϕ$ 或 $r$ 一定是均匀的（实际上， $ϕ$ 是均匀的，而 $r$ 不是均匀的）。微分面积 $dA$ 就是 $r \space dr \space dϕ$ 。因此，

\text { Probability }\left(r<\frac{R}{2}\right)=\int_{0}^{2 \pi} \int_{0}^{\frac{R}{2}} \frac{1}{\pi R^{2}} r d r d \phi=0.25 \text {. }

对于多维情况，实数函数的期望值公式同样适用：

E(f(x))=\int_{S} f(x) p(x) d \mu

其中 $x ∈ S$ ， $f : S↦ℝ$ ， $p : S↦ℝ$ 。例如，在单位正方形 $S=[0,1]×[0,1]$ 上且 $p(x,y)=4xy$ 时，对于 $(x,y) ~ p$ ， $x$ 坐标的期望值为

\begin{aligned} E(x) & =\int_{S} f(x, y) p(x, y) d A \\ & =\int_{0}^{1} \int_{0}^{1} 4 x^{2} y d x d y \\ & =\frac{2}{3} . \end{aligned}

请注意，这里 $f(x,y) = x$ 。

13.2.4 方差

一维随机变量的方差 $V(x)$ 的定义是 $x$ 与 $E(x)$ 差值的平方的期望值：

V(x)=E[x-E(x)]^{2}

一些代数处理得到了一个不显然的表达式：

V(x)=E\left(x^{2}\right)-[E(x)]^{2}

表达式 $E([x − E(x)]^2)$ 对于直观地思考方差更有用，而在计算中通常使用等价的代数表达式 $E(x2)−[E(x)]^2$ 。如果这些变量是独立的，则随机变量的和的方差是方差之和。方差的这种求和性质是其在概率模型分析中经常使用的原因之一。方差的平方根称为标准差 σ，它给出了离期望值的平均偏离程度的大致指示。

13.2.5 估计均值

许多问题涉及到独立随机变量 $x_i$ 的和，其中这些变量共享一个公共密度 $p$ 。这些变量被称为独立同分布 (iid) 随机变量。当将总和除以变量数量时，我们得到 $E(x)$ 的估计值：

E(x) \approx \frac{1}{N} \sum_{i=1}^{N} x_{i}

随着 $N$ 的增加，此估计的方差减少。我们希望 $N$ 足够大，以便我们对估计结果有足够的信心。然而，在蒙特卡罗中没有确定的事情；我们只是获得统计上的置信，认为我们的估计是好的。要确保，我们需要 $N = ∞$ 。大数定律 (Law of Large Numbers) 表达了这种置信度：

\text { Probability }\left[E(x)=\lim _{N \rightarrow \infty} \frac{1}{N} \sum_{i=1}^{N} x_{i}\right]=1 \text {. }

13.3 蒙特卡罗积分

本节概述了确定积分的基本蒙特卡罗解决方法。然后，将这些技术直接应用于某些积分问题。本节的所有基本材料也在几篇经典蒙特卡罗文本中涵盖。（请参见本章末尾的“注”部分。）

如前所述，给定一个函数 $f : S↦ℝ$ 和一个随机变量 $x ~ p$ ，我们可以通过求和来近似 $f(x)$ 的期望值：

E(f(x))=\int_{x \in S} f(x) p(x) d \mu \approx \frac{1}{N} \sum_{i=1}^{N} f\left(x_{i}\right) \tag{13.4}

由于期望值可以表示为积分，因此积分也可以用总和来近似。方程式 (13.4) 的形式有点笨拙；我们通常希望近似一个单个函数 $g$ 的积分，而不是 $fp$ 的积分。我们可以通过将 $g = fp$ 作为被积函数来实现这一点：

\int_{x \in S} g(x) d \mu \approx \frac{1}{N} \sum_{i=1}^{N} \frac{g\left(x_{i}\right)}{p\left(x_{i}\right)} \tag{13.5}

要使此公式有效，当 $g$ 不为零时， $p$ 必须是正值。

因此，为了获得良好的估计值，我们希望尽可能多地采样，并且 $g/p$ 应具有较低的方差（ $g$ 和 $p$ 应具有相似的形状）。聪明地选择 $p$ 被称为重要性采样，因为如果 $p$ 在 $g$ 大的区域大，则在重要区域中会有更多的样本。方程式 (13.4) 还显示了蒙特卡罗积分的基本问题：收益递减。由于估计值的方差与 $1/N$ 成比例，标准偏差与 $1/\sqrt{N}$ 成比例。由于估计误差的行为类似于标准偏差，我们需要将 $N$ 四倍才能使误差减半。

另一种减少方差的方法是将积分的定义域 $S$ 划分为几个较小的区域 $S_i$ ，并将积分评估为对 $S_i$ 上积分的总和。这被称为分层采样，即抖动在像素采样中使用的技术（第 4 章）。通常，在每个 $S_i$ 中只采取一个样本（密度为 $p_i$ ），在这种情况下，估计值的方差为

\operatorname{var}\left(\sum_{i=1}^{N} \frac{g\left(x_{i}\right)}{p_{i}\left(x_{i}\right)}\right)=\sum_{i=1}^{N} \operatorname{var}\left(\frac{g\left(x_{i}\right)}{p_{i}\left(x_{i}\right)}\right) . \tag{13.6}

可以证明，如果所有层具有相等的度量，则分层采样的方差永远不会高于未分层采样：

\int_{S_{i}} p(x) d \mu=\frac{1}{N} \int_{S} p(x) d \mu .

图形学中分层采样最常见的例子是像素采样的抖动。

作为蒙特卡罗方法求解积分 $I$ 的示例，将 $g(x)$ 设置为区间 $(0,4)$ 上的 $x$ ：

I=\int_{0}^{4} x d x=8 \tag{13.7}

表 13.1 显示了函数 $p$ 的形状对 $N$ 个样本估计值的方差的影响。请注意，当 $p$ 的形状类似于 g 的形状时，方差会减小。如果 $p = g/I$ ，则方差降至零，但通常不知道 $I$ ，否则我们就不必求助于蒙特卡罗。表 13.1 中展示的一个重要原则是，分层采样通常比重要性采样好得多 (Mitchell，1996)。尽管这种分层对 $I$ 的方差与样本数的立方成反比，但在分层下方差行为的一般结果没有。有些函数对于分层没有好处。其中之一是白噪声函数，在所有区域的方差都是恒定的。另一方面，大多数函数都将从分层采样中受益，因为每个子区域的方差通常小于整个域的方差。

13.3.1 拟蒙特卡罗积分

一种常用的数值积分方法是用拟随机点替换蒙特卡罗积分中的随机点。这些点是确定性的，但在某种意义上是均匀的。例如，在单位正方形 $[0,1]2$ 上，一组 $N$ 个拟随机点应具有以下属性在该正方形内部面积为 $A$ 的区域：

\text{区域中的点数}≈AN

例如，一个格点中的一组规则样本具有此属性。

拟随机点可以提高许多积分应用程序的性能。有时，必须小心确保它们不会引入混叠。特别好的是，在任何调用 $[0,1]^d$ 中的随机或分层点的应用程序中，都可以将 $d$ 维拟随机点替换为无需其他更改。

促使拟蒙特卡罗积分的关键直觉是，在估计被积函数的平均值时，只要采样点“公平”，任何一组采样点都可以使用。

13.4 随机点选择

我们经常希望为分布光线跟踪等应用程序在单位正方形上生成一组随机或伪随机点。有几种方法可以做到这一点，例如抖动。这些方法给出了一组在单位正方形 $[0,1]^2$ 上相当均匀分布的 $N$ 个点：从 $(u_1,v_1)$ 到 $(u_N,v_N)$ 。

有时，我们的采样空间可能不是正方形（例如，圆形透镜）或可能不是均匀的（例如，以像素为中心的滤波器函数）。如果我们可以编写一个数学变换，使其将我们均匀分布的点 $(u_i,v_i)$ 作为输入，并输出一组具有所需密度的点，则好极了。例如，要对相机透镜进行采样，该变换将获取 $(u_i,v_i)$ 并输出 $(r_i,ϕ_i)$ ，以便新的点大致均匀分布在透镜的圆盘上。虽然我们可能会尝试使用变换

ϕ_i=2πu_i, \\ r_i=v_iR,

但它存在严重问题。虽然这些点覆盖了透镜，但它们并不均匀（图 13.6）。在这种情况下，我们需要的是一种将等面积区域映射到等面积区域的变换，它将在正方形上的均匀采样分布映射到新域上的均匀分布。

图 13.6。将水平和垂直维度均匀映射到 $(r,ϕ)$ 的变换不保持相对面积；并非所有的结果区域都相同。

有几种方法可以生成这种非均匀点或非矩形域上的均匀点，下面的章节将回顾最常用的三种方法：函数反演、拒绝采样和 Metropolis 方法。

13.4.1 函数反演

如果密度函数 $f(x)$ 是一维的，并定义在区间 $x∈[x_{min},x_{max}]$ 上，则可以从一组均匀随机数 $ξ_i$ 中生成具有密度 $f$ 的随机数 $α_i$ ，其中 $ξ_i ∈ [0,1]$ 。为此，我们需要累积概率分布函数 $P(x)$ ：

\text { Probability }(\alpha<x)=P(x)=\int_{x_{\min }}^{x} f\left(x^{\prime}\right) d \mu \text {. }

要获取 $α_i$ ，我们只需将 $ξ_i$ 进行变换：

\alpha_{i}=P^{-1}\left(\xi_{i}\right),

其中 $P-1$ 是 $P$ 的逆函数。如果 $P$ 不能解析地求逆，则数值方法就足够了，因为所有有效的概率分布函数都存在逆。

请注意，解析反演函数比应该更令人困惑，这是由于符号表示法造成的。例如，如果我们有函数

y=x^2

对于 $x>0$ ，则逆函数用 $y$ 表示为 $x$ 的函数：

x=\sqrt{y}

当函数可以解析求逆时，它几乎总是那么简单。然而，使用标准符号表示法时比较不透明：

\begin{aligned} f(x) & =x^{2}, \\ f^{-1}(x) & =\sqrt{x} . \end{aligned}

这里， $x$ 只是一个虚拟变量。您可能会发现以下非标准符号表示法更容易理解：

y=x^2, \\ x=\sqrt{y}，

同时牢记它们是彼此的逆函数。

例如，在 $[-1, 1]$ 上选择具有密度

p(x) = {3x^2 \over 2}

的随机点 $x_i$ ，我们看到

P(x) = {x^3+1 \over 2}

和

P^{-1}(x)=\sqrt[3]{2 x-1}

因此，我们可以将一组规范化的随机数（ξ1,...,ξN）“扭曲 (wrap)”为合适分布的数字

\left(x_{1}, \ldots, x_{N}\right)=\left(\sqrt[3]{2 \xi_{1}-1}, \ldots, \sqrt[3]{2 \xi_{N}-1}\right)

当然，这种相同的扭曲函数也可用于将“均匀”的抖动采样转换为具有所需密度的漂亮分布采样。

如果我们有一个随机变量 $α=(α_x,α_y)$ ，其二维密度 $(x,y)$ 定义在 $[x_{min},x_{max}] × [y_{min},y_{max}]$ 上，则需要二维分布函数：

\operatorname{Probability}\left(\alpha_{x}<x \text { and } \alpha_{y}<y\right)=F(x, y)=\int_{y_{\min }}^{y} \int_{x_{\min }}^{x} f\left(x^{\prime}, y^{\prime}\right) d \mu\left(x^{\prime}, y^{\prime}\right) \text {. }

我们首先使用边际分布 $F(x,y_{max})$ 选择 xi，然后根据 $F(x_i,y)∕F(x_i,y_{max})$ 选择 $y_i$ 。如果 $f(x,y)$ 是可分离的（可表示为 $g(x)h(y)$ ），则可以对每个维度使用一维技术。

回到我们之前的例子，假设我们从半径为 $R$ 的圆盘中均匀采样，因此 $p(r,ϕ) = 1∕(πR^2)$ 。二维分布函数是

\operatorname{Probability}\left(r<r_{0} \text { and } \phi<\phi_{0}\right)=F\left(r_{0}, \phi_{0}\right)=\int_{0}^{\phi_{0}} \int_{0}^{r_{0}} \frac{r d r d \phi}{\pi R^{2}}=\frac{\phi r^{2}}{2 \pi R^{2}}

这意味着，规范化的一对 $(ξ1,ξ2)$ 可以变换为圆盘上的均匀随机点：

\begin{aligned} \phi & =2 \pi \xi_{1} \\ r & =R \sqrt{\xi_{2}} \end{aligned}

该映射在图 13.7 中显示。

图 13.7。一个将单位正方形上的等面积区域映射到圆盘上的等面积区域的映射。

为了选择一些逼真渲染应用程序的反射光线方向，我们根据密度

p(\theta, \phi)=\frac{n+1}{2 \pi} \cos ^{n} \theta

选择单位半球上的点，其中 $n$ 是类似于 Phong 的指数， $θ$ 是与表面法线成角度， $θ ∈ [0,π∕2]$ （处于上半球）， $ϕ$ 是方位角（ $ϕ ∈ [0,2π]$ )。累积分布函数是

P(\theta, \phi)=\int_{0}^{\phi} \int_{0}^{\theta} p\left(\theta^{\prime}, \phi^{\prime}\right) \sin \theta^{\prime} d \theta^{\prime} d \phi^{\prime} \tag{13.8}

这里的 $sinθ'$ 项是因为在球体上， $dω = cosθ dθ dϕ$ 。当找到边缘密度 $p$ 时（如预期），它是可分离的，并且我们发现规范化的一对 $(ξ1,ξ2)$ 可以通过以下方式转换为方向：

\begin{array}{l} \theta=\arccos \left(\left(1-\xi_{1}\right)^{\frac{1}{n+1}}\right) \\ \phi=2 \pi \xi_{2} \end{array}

同样，好的一点是，可以轻松地将单位正方形上的一组抖动点转换为具有所需分布的半球上的一组抖动点。请注意，如果将 $n$ 设为 $1$ ，则会得到一种漫反射分布，这通常是必要的。

通常，我们必须将在球面上的点映射到相对于 $uvw$ 基向量的适当方向。为此，我们可以首先将角度转换为单位向量 $\vec{a}$ ：

\mathbf{a}=(\cos \phi \sin \theta, \sin \phi \sin \theta, \cos \theta)

作为效率改进，我们可以避免取反三角函数的三角函数（例如 $cos(arccosθ)$ ）。例如，当 $n = 1$ （漫反射分布）时，向量 $\bold{a}$ 简化为

\mathbf{a}=\left(\cos \left(2 \pi \xi_{1}\right) \sqrt{\xi_{2}}, \sin \left(2 \pi \xi_{1}\right) \sqrt{\xi_{2}}, \sqrt{1-\xi_{2}}\right)

13.4.2 拒绝采样

拒绝方法根据某些简单分布选择点，并拒绝一些在更复杂分布中的点。有几种情况使用拒绝方法，我们通过示例展示其中一些。

假设我们想要在单位圆内均匀随机选取点。我们可以首先选择均匀随机点 $(x,y)∈[-1,1]2$ ，然后拒绝那些在圆外的点。如果函数 $r()$ 返回规范化的随机数，则过程为

如果我们想要一个随机数 $x~p$ ，我们知道 $p:[a,b]↦ℝ$ ，并且对于所有 $x，p(x)<m$ ，则我们可以在矩形 $[a,b]×[0,m]$ 中生成随机点，并取其中 $y<p(x)$ 的点：

这个想法也可以应用于在球面上取随机点。为了选择具有均匀方向分布的随机单位向量，我们首先在单位球中选择一个随机点，然后通过取同一方向上的单位向量将该点视为方向向量：

虽然拒绝方法通常很容易编码，但它很少与分层兼容。因此，它的收敛速度较慢，因此主要用于调试或特别困难的情况。

13.4.3 Metropolis

Metropolis 方法使用随机变异来产生具有所需密度的样本集。这个概念在本章注释中引用的 Metropolis 光传输算法中被广泛使用。假设我们有定义在区域 $S$ 中的随机点 $x_0$ 。此外，假设对于任何点 $x$ ，我们都可以生成随机 $y~p_x$ 的方法。我们使用边缘符号 $px(y) ≡ p(x → y)$ 来表示该密度函数。现在，假设我们让 $x_1$ 成为 $S$ 中的一个随机点，该点由基础密度 $p(x_0 → x_1)$ 选择而成。我们使用密度 $p(x_1 → x_0)$ 生成 $x^2$ ，以此类推。在生成无限多个样本的极限情况下，可以证明这些样本将具有由 p 确定的某些基础密度，而与初始点 $x_0$ 无关。

现在，假设我们希望选择 $p$ ，使得我们收敛到的样本的基础密度与非负函数 $f(x)$ 成比例，其中 $f$ 具有 $S$ 作为定义域。此外，假设我们可以评估 $f$ ，但我们几乎没有关于其性质的其他知识（这种函数在图形学中很常见）。此外，假设我们有能力从 $x_i$ 到 $x_i+1$ 进行“转换”，其基础密度函数为 $t(x_i → x_i+1)$ 。为了增加灵活性，进一步假设我们添加了可能为非零的概率，即 $x_i+1$ 转换为自身，即 $x_i+1 = x_i$ 。我们将其表述为生成一个潜在候选项 $y~t(x_i→y)$ ，并以概率 $a(x_i→y)$ “接受”该候选项（即 $x_i+1=y$ ），以概率 $1−a(x_i→y)$ “拒绝”它（即 $x_i+1=x_i$ ）。请注意，序列 $x_0,x_1,x2,…$ 将是一个随机集合，但样本之间会存在某些相关性。它们仍然适合于蒙特卡罗积分或密度估计，但分析这些估计的方差要困难得多。

现在，假设我们已经给出了一个转换函数 $t(x→y)$ 和一个我们想要模拟其分布的函数 $f(x)$ ，我们可以使用 $a(y→x)$ 来使点以 $f$ 的形状分布吗？或者更精确地说，

\left\{x_{0}, x_{1}, x_{2}, \ldots\right\} \sim \frac{f}{\int_{s} f}

事实证明，这可以通过确保 $x_i$ 在某种强意义下是平稳的来实现。如果您将大量样本点 $x$ 可视化，则希望两个点之间的“流”在每个方向上都相同。如果我们假设与 $x$ 和 $y$ 相邻的点的密度分别与 $f(x)$ 和 $f(y)$ 成比例，则两个方向上的流应该相同：

\begin{array}{l} \text { flow }(x \rightarrow y)=k f(x) t(x \rightarrow y) a(x \rightarrow y), \\ \text { flow }(y \rightarrow x)=k f(y) t(y \rightarrow x) a(y \rightarrow x), \end{array}

其中 k 是某个正常数。设置这两个流恒定会对 a 产生约束：

\frac{a(y \rightarrow x)}{a(x \rightarrow y)}=\frac{f(x) t(x \rightarrow y)}{f(y) t(y \rightarrow x)} .

因此，如果 $a(y → x)$ 或 $a(x → y)$ 中的任何一个已知，则另一个也已知。使它们变大可以提高接受的机会，因此通常的技巧是将它们中较大的一个设为 $1$ 。

在使用 Metropolis 样本生成技术时的一个困难是很难估计在集合点“好”的情况下需要多少点。如果丢弃前 $n$ 个点，则可以加速处理，尽管明智地选择 $n$ 是不易的。

13.4.4 示例：在正方形中选择随机直线

作为设计采样策略的完整过程的示例，考虑找到与单位正方形 $[0,1]^2$ 相交的随机直线。我们希望这个过程是公平的；也就是说，我们希望直线在正方形内均匀分布。直观地看，我们可以看到这个问题有些微妙；以斜角度量的直线“更多”，而水平或垂直方向上的直线较少。这是因为正方形的横截面不是均匀的。

我们的第一个目标是找到一个函数反演方法（如果存在），然后在失败时退回到拒绝或 Metropolis 方法。这是因为我们希望在线空间中有分层样本。我们首先尝试使用正常坐标，因为选择正方形中的随机直线的问题只是在对应于正方形中的直线的 $(r,θ)$ 空间的任何部分中找到均匀随机点的问题。

考虑 $−π∕2 < θ < 0$ 的区域。哪些 $r$ 值对应于击中正方形的直线？对于那些角度， $r < cosθ$ 是所有击中正方形的直线，如图 13.8 所示。在其他四个象限中类似的推理找到了在 $(r,θ)$ 空间中必须采样的区域，如图 13.9 所示。该区域的边界方程 $r_{max}(θ)$ 是

图 13.8。最大距离 $r$ 对应于 $θ∈[−π∕2,0]$ 的击中正方形的直线。因为正方形的边长为 $1$ ，所以 $r=cosθ$ 。

图 13.9。随 $θ$ 变化而变化，击中单位正方形 $[0,1]^2$ 的直线的最大半径。

r_{\max }(\theta)=\left\{\begin{array}{ll} 0 & \text { if } \theta \in\left[-\pi,-\frac{\pi}{2}\right], \\ \cos \theta & \text { if } \theta \in\left[-\frac{\pi}{2}, 0\right], \\ \sqrt{2} \cos \left(\theta-\frac{\pi}{4}\right) & \text { if } \theta \in\left[0, \frac{\pi}{2}\right], \\ \sin \theta & \text { if } \theta \in\left[\frac{\pi}{2}, \pi\right] . \end{array}\right.

因为 $r_{max}(θ)$ 下面的区域是一个简单函数，下界由 $r=0$ 界定，所以我们可以通过首先根据密度函数选择 θ 来对其进行采样：

p(\theta)=\frac{r_{\max }(\theta)}{\int_{-\pi}^{\pi} r_{\max }(\theta) d \theta} .

这里的分母为 $4$ 。现在，我们可以计算累积概率分布函数：

P(\theta)=\left\{\begin{array}{ll} 0 & \text { if } \theta \in\left[-\pi,-\frac{\pi}{2}\right] \\ (1+\sin \theta) / 4 & \text { if } \theta \in\left[-\frac{\pi}{2}, 0\right] \\ \left(1+\frac{\sqrt{2}}{2} \sin \left(\theta-\frac{\pi}{4}\right)\right) / 2 & \text { if } \theta \in\left[0, \frac{\pi}{2}\right], \\ (3-\cos \theta) / 4 & \text { if } \theta \in\left[\frac{\pi}{2}, \pi\right] . \end{array}\right.

我们可以通过将 ξ_1=P(θ) 操作成 θ=g(ξ_1) 的形式来反演这个函数。这得到

\theta=\left\{\begin{array}{ll} \arcsin \left(4 \xi_{1}-1\right) & \text { if } \xi_{1}<\frac{1}{4} \\ \arcsin \left(\frac{\sqrt{2}}{2}\left(2 \xi_{1}-1\right)\right)+\frac{\pi}{4} & \text { if } \xi_{1} \in\left[\frac{1}{4}, \frac{3}{4}\right] \\ \arccos \left(3-4 \xi_{1}\right) & \text { if } \xi_{1}>\frac{3}{4} \end{array}\right.

一旦我们有了 $θ$ ，则 $r$ 简单地是

r=\xi_{2} r_{\max }(\theta)

如前所述，直线有许多参数化方式，每种方式都有一个相关的“公平”度量。我们也可以在这些空间中生成随机直线。例如，在斜率-截距空间中，击中正方形的区域如图 13.10 所示。与正常空间类似的推理，斜率的密度函数为

p(m)=\frac{1+|m|}{4}

图 13.10。包含与单位正方形 $[0,1]^2$ 相交的直线的 $(m,b)$ 空间区域。

相对于差分测量

d \mu=\frac{d m}{\left(1+m^{2}\right)^{\frac{3}{2}}}

这导致了累积分布函数：

P(m)=\left\{\begin{array}{ll} \frac{1}{4}+\frac{m+1}{4 \sqrt{1+m^{2}}} & \text { if } m<0 \\ \frac{3}{4}+\frac{m-1}{4 \sqrt{1+m^{2}}} & \text { if } m \geq 0 \end{array}\right.

可以通过解两个二次方程来反演这些分布函数。给定使用 $ξ_1$ 生成的 $m$ ，则有

b=\left\{\begin{array}{ll} (1-m) \xi_{2} & \text { if } \xi<\frac{1}{2} \\ -m+(1+m) \xi_{2} & \text { otherwise } \end{array}\right.

这并不比使用正常坐标更好；它只是另一种方式。

常见问题解答

本章讨论了概率而不是统计学。两者有何不同？

概率是研究事件发生的可能性。统计学推断出随机变量的大但有限的总体特征。在这个意义上，统计学可以被视为应用概率的一种特定类型。

Metropolis 抽样和 Metropolis 光传输算法是相同的吗？

不是。Metropolis 光传输（Veach & Guibas，1997）算法在其过程中使用 Metropolis 抽样，但它专门用于渲染，并且还有其他步骤。

注释

几何概率的经典参考书是《Geometric Probability》（Solomon，1978）。在正方形中选择随机边的另一种方法在《Random - Edge Discrepancy of Supersampling Patterns》（Dobkin＆Mitchell，1993）中给出。有关图形学中准蒙特卡罗方法的更多信息，请参见《Efficient Multidimensional Sampling》（Kollig＆Keller，2002）。关于 Monte Carlo 方法的三本经典且易读的书是《Monte Carlo Methods》（Hammersley＆Handscomb，1964）、《Monte Carlo Methods,Basics》（Kalos＆Whitlock，1986）和《The Monte Carlo Method》（Sobel、Stone、＆Messer，1975）。

练习题

在单位立方体 $(x,y,z) ∈ [0,1]^3$ 中，函数 $xyz$ 的平均值是多少？
在单位半径的圆盘上， $r$ 的平均值是多少： $(r,ϕ) ∈ [0,1] × [0,2π)$ ？
证明将规范随机点 $(ξ_1,ξ_2)$ 均匀映射到任意三角形的重心坐标的映射是： $β=1− \sqrt{1−ξ_1}$ ， $γ=(1-u)ξ_2$ 。
单位正方形中一条线段的平均长度是多少？通过在单位正方形内生成一千万条随机线段并计算其长度的平均值来验证您的答案。
单位立方体中一条线段的平均长度是多少？通过在单位立方体内生成一千万条随机线段并计算其长度的平均值来验证您的答案。
根据方差的定义，证明 $V(x) = E(x^2) - [E(x)]^2$ 。

目录