PBRT阅读：第十五章 蒙特卡罗积分 2: 提高效率 （第15.1-15.5节）

http://www.opengpu.org/forum.php?mod=viewthread&tid=7363&fromuid=10107

第15章 蒙特卡罗积分II：提高效率

蒙特卡罗光线追踪算法里的方差以图像噪声的形式体现出来。如何竭力地减小方差竟成了大多数优化蒙特卡罗算法的基本任务。回忆一下蒙特卡罗算法的收敛速度，这意味着为了对方差减半，就需要四倍的采样。由于估算例程的运行时间跟采样数目成比例关系，所以减少方差的花销可以是很高的。本章将介绍可以改进蒙特卡罗算法的效率而并不增加采样的理论和技术。

一个估计量公式F的效率可以定义为：

ε[F] = 1 / V[F]T[F]

其中V[F]是它的方差，T[F], T[F]是求值所用的运行时间。根据这个度量方法，如果估计量F1跟另一个估计量F2相比，花费更少的时间而产生相同的方差，或者在相同时间内生成更小的方差，则称它更有效率。其中一项技术被认为是改进渲染效率的最为有效的技术，即重要性采样(importance sampling)。回忆一下，在蒙特卡罗估计量的实现中，我们有一定的自由来选择采样分布函数p(x)。如果我们选择跟被积函数很相像的采样分布，那么就可以减小方差。这项技术之所以被称为重要性采样，是因为我们趋向于在函数定义域的“重要”部分进行采样，而这些部分有相对大的函数值。

本章的前半部分讨论重要性采样和几项其它改进蒙特卡罗算法效率的技术。本章的后半部分讨论根据BSDF、光源、体积散射来生成采样的技术，而这些技术可用做重要性采样中的采样分布。第16，17章将使用这些采样方法。

15.1 俄罗斯轮盘法和分割法

俄罗斯轮盘法和分割法是两项可以改进蒙特卡罗算法效率的相关方法，它们的做法是要增加对结果有较大贡献值的采样的可能性。有些求值很耗时的采样却对最终结果的贡献值很小，俄罗斯轮盘法要解决就是这个问题。我们希望避免对这些采样进行求值却仍可以计算出正确的估计量。

最为一个例子，我们考虑对直接光照积分的估算问题，该积分给出了场景中从光源直接照射到某点而产生的反射辐射亮度：

假定我们从某个分布p(ω)取N=2个采样来计算估计量：

其中累加式中绝大部分的计算开销来自于追踪起始于点p(看光源是否被遮挡)的阴影光线。

对于那些被积函数值明显为0的那些方向（例如，对于fr(p, ωo, ωi)为0的那些方向） 而言，我们可以跳过追踪阴影光线的工作，因为是否去追踪并不影响最终的计算结果。使用俄罗斯轮盘法可以让我们也跳过那些被积函数值很低但不一定为0的那些方向的光线追踪，而仍然大致可以计算出正确的结果。例如，在那些fr(p, ωo, ωi)很小或者ωi靠近水平面使得|cosθi|很小的方向上，我们可以避免进行光线追踪。当然，我们不能完全忽略这些采样，否则估算公式就会低估结果值。

为了使用俄罗斯轮盘法，我们先要选定某个终止概率q。我们几乎可以用任意的方式来选择这个概率值；例如，我们可以根据一个选定的采样来估算出被积函数的值，并以此为基础确定一个概率值，当被积函数值减小时再对之增加。在使用概率q时，就不对被积函数求值，而是用一个常数c代替（一般用c=0）。在使用概率1-q时，就要对被积函数求值，但要乘以一个权值1/(1-q)，用来弥补那些被跳过的采样：

F' = (F-qc)/(1-q), ξ > q

F' = c, otherwise

要注意的是它的期望值跟原估计量的期望值相同：

E[F'] = (1-q) ((E[F] - qc)/(1-q)) + qc = E[F]

俄罗斯轮盘法绝不能减少方差。实际上，除非c = F，否则它总是加大方差的。然而，如果概率选择得当的话，就可以跳过那些对最终结果贡献很小的采样，从而提高了效率。

其中一个陷阱是，如果俄罗斯轮盘法选择了一个很差的权值，就会大大地增加方差。想象一下我们用0.99做为终止概率对所有的相机光线使用俄罗斯轮盘法，我们即只追踪其中的1%的相机光线，而用1/0.1 = 100的权值。在严格的数学意义下结果仍然“正确”，但最终的视觉效果无比糟糕：大多数像素是黑色的，只有很少明亮的像素。

15.1.1 分割法

俄罗斯轮盘法可以减少对次要的采样的求值，而分割法却以增加采样数的方式来提供效率。我们再考虑一下计算直接光照下的反射辐射亮度的问题。在忽略像素滤波的情况下，这个问题可以写成一个像素的面积A下以及方向球面δ2下的双重积分，其中Ld(x,y, ω)是图像坐标(x,y)处由方向ω上入射辐射亮度产生的出射辐射亮度：

估算该积分的一个自然方式就是生成N个采样，每个采样包含一个(x,y)图像位置和一个朝向光源的方向ω，然后使用蒙特卡罗估计值公式。如果场景中有许多光源，或者有一个投射软阴影的面光源，为了使结果有一个可以接受的方差水平，就需要数十个或数百个采样。不幸地是，每个采样需要两条要被追踪的光线：第一条用来计算从图像平面(x,y)处看到的第一个可见表面，第二条是沿方向ω到光源的阴影光线。

这个方法的问题是，如果使用N=100个采样来估算这个积分，那么就要追踪100条相机光线和100条阴影光线，而且100条相机光线对于像素反走样运算是过多了，因此对最终结果的方差的减少并没有什么贡献。分割法为了解决这个问题，将这个方法做了正规化处理：即在某些维上做了采样以后，对每一个采样再在其它维上做多个采样。

利用分割法，这个积分的估计值公式可以写成（使用N个图像采样，每个图像采样使用M个光源采样）：

这样一来，我们只使用5个图像采样，但对每个图像采样使用20个光源采样，就可以只追踪105条光线（而不是200条），我们仍旧使用100个面光源采样来得到高质量的软阴影效果。

14.2 精细的采样布置

为了更好地扑捉到“重要的”被积函数的特征，有一类经典且有效的减少方差的技术是基于对采样的精细布置上的。这类技术是对重要性采样等技术的有益补充，并广泛地应用于pbrt中。实际上，第7章中的Sampler的一项任务就是为此目的而生成具有良好分布的采样。

14.2.1 分层采样

我们在第7.3节介绍了分层采样，现在我们有了使用这项技术的工具。分层采样将积分定义域Λ分成n个互不重叠的区域Λ1, Λ2,..., Λn。每个区域称为一个层，它们的并集完全覆盖了Λ。

为了从Λ中进行采样，我们将按照每个层中的密度pi来在Λi中抽取ni个采样。像素周围的区域被分成一个k x k网格，我们利用分层采样的方法在每个网格单元抽取一个采样。这样做要好于随机地抽取k2个采样，因为采样位置不会拥挤到一起。下面我们将表明为什么这种技术会减小方差。

在一个层Λi中，蒙特卡罗估计量为：

其中Xi,j是从密度pi中抽取的第j个采样。总体的估算值为F = Σi=1,nvi Fi，而vi是层i的体积所占的比例(0 < vi <= 1)。

被积函数在层i的真值为：

μi = E [f(Xi,j)] = 1/vi ∫Λi f(x)dx

故得在这个层中的方差：

σi2 = 1 / vi ∫Λi (f(x) - μi )dx

因此，由于层中的采样有ni个，故层的估计量的方差为σi2 / ni。这表明总估计量的方差为：

如果我们做这样一个比较合理的假定：采样个数ni正比于体积vi，则有 ni = vi N，那么总体估计量的方差公式为：

V [ FN] = 1/N Σvi σi2

为了将这个结果跟没有分层的方差相比较，我们注意到选择一个没有分层的采样等价于这样的做法：即按照由体积vi定义的离散的概率分布选择一个分层I，然后在层Λi中选择一个随机采样X。在这种意义下，X是在I的条件下选择的，所以我们用条件概率公式有：

V[F] = ExViF + VxEiF

= 1/N [Σvi σi2 + 1/N Σvi (μi - Q)]

这里Q是f在整个定义域上的平均值。对这个结果的推导见Veach(1997)。

这个表达式有两个值得注意的地方。首先，第二项累加和必须是非负的，因为方差总是非负的。第二，它显示出分层采样绝不会增加方差。实际上，除非第二项的和为0，分层法总是可以减小方差的。只有当f在每个层有相同的平均值时，该式才为0。事实上，为了使分层采样效果达到最佳，我们希望第二项的和达到最大化，这样就使得各层上的平均值尽可能地不均等。如果我们不知道f的任何情况，就需要更“紧”(compact)的分层。如果分层很宽，就会有更大的方差，就会有更靠近真正平均值Q的μi。

分层采样的主要缺点是它也受制于“维数灾难”(curse of dimensionality)。在D维空间中每个维上做S个分层，就会需要SD个采样，随着维数的增加，就很快地变为不可行了。幸运地是，我们常常可以只对某些维做独立的分层，然后随机地将来自不同维的采样联系起来，正如我们在第7.3节中所做的那样。我们应当这样选择需要分层的维：这些维在被积函数值的效果上要尽可能地高度相关联。例如，对于第15.1.1节上的直接光照的例子，对像素位置(x,y)和光线方向(θ,Φ)分层要比对(x, θ)和(y,
Φ)要高效得多。

另一个解决“维数灾难”的方法是前面介绍过的拉丁超立方体采样。不幸地是，拉丁超立方体采样在减少方差方面并不如分层采样那样高效。然而，拉丁超立方体采样的效果不比均匀随机采样差，而且经常要好上许多。

15.2.2 准蒙特卡罗算法

在第7.4节介绍过的低差异采样计算是一类被称为“准蒙特卡罗”的算法的基础。准蒙特卡罗算法的核心部分就是将标准蒙特卡罗算法中所用的随机数用低差异点集来代替，而这些低差异点是由精细设计的确定性算法生成的。

该方法的优点是，对于许多积分问题，准蒙特卡罗算法比基于标准蒙特卡罗的方法有更快的渐近收敛速度。许多用于常规蒙特卡罗算法的技术也可以很好地应用到准随机采样点上，但有些技术并非如此（例如俄罗斯轮盘法和采样拒绝法）。虽然渐近收敛速度通常无法适用于不连续的被积函数（因为收敛速率依赖于被积函数的光滑性），但准蒙特卡罗算法对这些被积函数而言仍然比常规蒙特卡罗算法要好一些，详见Keller 2001。

在pbrt中，我们忽略了这两种方法的差别，并在第7章中对Sampler做了适应性上的改动。这会引起细微错误的可能性，但只要Integrator使用这两类采样都适用的技术（如重要性采样或其它技术），就不会产生问题。

15.2.3 失真的采样和变形

当用分层采样或低差异采样来处理诸如在光源采样等问题时，pbrt生成[0,1]x[0,1]上的一组采样(ξ1,ξ2)，然后用第14.4和14.5节中的变换算法将这些采样映射为光源上的点。我们期望在这个过程中[0,1]x[0,1]上的采样被变换后仍然保持分层性质---换句话说，邻近的采样在光源的表面上应该还是保持邻近关系，而彼此远离的采样在光源上应该还是彼此远离的。如果映射无法保留住这个性质，那么分层的好处就被丧失了。

例如，这可以解释为什么Shirley的“正方形到圆”的映射要好于直接映射，因为直接映射在远离圆心的地方更有宽松些的分层。这也可以解释为什么(0,2)序列法在实际应用中要比分层法好上许多：因为通过定义域之间的变换后，它们比分层法相比更容易保持住良好的分布特性。

15.3 偏差

另一个减少方差的方法是将偏差引入计算中：有时我们知道计算一个估计值并不能跟期望值相等，利用这一点却可以达到减小方差的目的。如果一个估计值的期望值等于正确结果，则称它是无偏差的。否则，它的偏差量等于：

β = E [F] - F

Kalos和Whitlock（1986）给出了下面的一个利用偏差的例子。考虑一下如何估算一组在区间[0,1]之间的随机变量的平均值。我们可以用下面的估计值：

1/N Σi=1,nXi

或者用有偏差的估计值：

1/2 max(X1, X2, ..., XN)

第一个估计值实际上是无偏差的，但方差是O(1/N-1)。第二个估计值的期望值是：

0.5 N / (N+1) ≠ 0.5

所以它是有偏差的，虽然它的方差是O(N-2)，所以要好得多。对于很大的N值而言，我们可选用第二个估计值。

在第7.6节中的像素重构方法也可以视为一个有偏差的估计值。作为一个蒙特卡罗估算问题，我们可以估算下面的公式：

I (x,y) = ∫∫ f(x-x', y-y') L(x',y') dx' dy'

其中I(x,y)是最终的像素值，f(x,y)是像素滤波函数（假定它已经被归一化，即积分为1），L(x,y)是图像辐射亮度函数。

假定我们在图像平面上进行均匀采样，所有的采样有相同的概率密度（记为pc)，那么该方程的蒙特卡罗无偏差估计值为：

I(x,y) ≈ 1/Npc Σi=1,N f(x-xi,y-yi)L(xi,yi)

这个结果跟前面我们所用的像素滤波函数不同，它在第7章：

I(x,y) = Σi f(x-xi,y-yi)L(xi,yi) / Σi f(x-xi,y-yi)

然而带偏差的估计值在实际应用中成为首选，是因为其结果有更小的方差。例如，如果所有的辐射亮度值L(xi,yi)的值为1，那么带偏差的估计值在图像重构时将所有像素值设为1，这正是我们需要的性质。而无偏差估计值总是将像素设为不是0的值，因为Σi f(x-xi,y-yi)一般并不等于pc，由于滤波函数的方差而得到不同的值。甚至对于更复杂的图像而言，由无偏差估计值引起的方差跟公式7.3的偏差相比更令人无法接受。

15.4 重要性采样

重要性采样是一种减小方差的技术，它利用了这个事实：当分布函数p(x)跟函数f(x)很相似时，下列的蒙特卡罗估计值就收敛更快：

FN = 1/N Σi=1,N f(Xi)/p(Xi)

它的基本思想是，将采样工作集中到被积函数值相对比较高的地方，那么就可以更有效率地计算出精确的估算值。

例如，假定我们要对第5章的散射方程5.6进行求值。考虑一下估算积分值的过程：如果我们随机地采样到一个几乎跟表面法向量垂直的方向，那么它的余弦项几乎为0。那么对BSDF的计算以及在采样位置上追踪光线以求得入射辐射亮度的工作就白白浪费掉了，因为最终所得到的贡献值微乎其微。所以，如果我们能够减少到球面赤道附近采样的可能性，就会得到更满意的结果。更一般地，如果采样的分布跟被积函数的其它要素相匹配（例如BSDF，入射照明分布，等等），那么也会相应地提供效率。

只要随机变量所用的概率分布在形状上类似于被积函数，就可以减少方差。我们在这里不做严格的证明，而是做一个非正式的直观的讨论。假定我们对某个积分∫f(x)dx用蒙特卡罗技术进行求值，因为我们可以自由地选择一个采样分布，我们就可以考虑一个跟f(x)呈比例的一个分布，p(x) = c f(x)。很显然地，归一化使得下式成立：

c = 1 / ∫f(x)dx

为了找到这样的PDF，我们需要知道积分值，而这恰恰是我们要计算的东西。然而，对这个例子而言，我们“可以”从这个分布上采样，每个估算值将有下列估计：

f(Xi)/p(Xi) = 1/c = ∫f(x)dx

因为c是一个常数，所以每个估算值有相同的值，那么方差就是0。当然这有点儿搞笑，如果我们能够直接对f积分就不用劳烦蒙特卡罗了。但是，如果能找到跟f(x)很形似的密度函数p(x)，确实可以减少方差。

实际上，如果选择了一个很差的分布，重要性分布可以增加方差。为了理解选择有很差匹配性的分布而带来的影响，我们考虑下面的分布：

p(x) = 99.01, x ∈[0, .01)

0.01, x ∈[.01, 1)

用它来估算下面的积分：

f(x) = .01, x ∈[0, .01)

1.01, x ∈[.01, 1)

我们通过构造使得该积分的值为1，然而用p(x)来做蒙特卡罗积分却得到很糟糕的结果：几乎所有的采样来自区间[0, .01)，这时估计值大约等于f(x)/p(x) ≈ 0.0001。如果不存在这个区间范围之外的采样，那么结果就极其不精确，几乎比正确值小10000倍。如果有些采样取自于区间[.01,1)，那么结果会更糟。虽然这极少发生，但如果发生的话，那么我们就会有相对很高的被积函数值，和相对很低的PDF值： f(x)/p(x) = 101。这就需要大量的采样来平衡这些极端的值，才可以将方差减小到估算值和实际值（1）比较相近的程度。

幸运地是，对于图形学中的许多积分问题，找到重要性采样的良好采样分布并不太困难。例如，在许多情况下，被积函数是多个函数的乘积。构造出跟整个乘积相似的PDF是很难甚至是不可能的，但只要找到跟其中一个因子相似的PDF也是很有用的。在接下来的两章中的光传输算法中，这是一个非常常用的策略，其中我们要的估算积分是光照、可见性、散射和余弦项的乘积。

15.4.1 多重重要性采样

蒙特卡罗方法提供了估算形如∫f(x)dx形式的积分的工具。然而，我们通常要处理两个或多个函数乘积的积分：∫f(x)g(x)dx。如果我们有了f(x)的重要性采样策略，也有了g(x)的策略，我们将使用哪一个呢？一般地说，我们很难将这样两个采样策略组合出一个跟f(x)g(x)呈正比且容易用来采样的概率密度函数。正如我们在讨论重要性采样时所显示的那样，如果所选的采样分布很糟糕，那么其效果被一个均匀分布要差得多。

例如，考虑一下直接光照积分的求值问题，其形式如下：

如果我们对这个积分求值时按照基于Ld 或fr的分布进行重要性采样，那么这两个分布中的一个常常性能不佳。

考虑一下一个被面光源照明的接近镜面效果的BRDF，这里我们使用Ld的分布进行采样。因为BRDF几乎就是一面镜子，除了镜面方向附近的地方，其它所有的方向ωi的被积函数值接近为0。这意味着几乎所有的来自Ld分布的方向采样有零贡献值，那么方差就相当地高。更糟糕地是，如果增大光源面积，就要生成更多的方向采样，PDF的值就相应减小，对于采样方向上BRDF不为0的那些地方，就会出现相当大的被积函数值除以很小的PDF值的情况。在这种特殊情况下，使用BRDF分布进行采样效果要好得多，而对于漫反射或光泽反射的BRDF和小光源而言，使用BRDF分布采样要比使用光分布进行采样有更高的方差。一个显然的方案是从各个分布中采样，在将两个估计值取平均，不幸地是，这种方法没有什么改进。因为在这种情况下方差是累加的，这种方法就无能为力了---一旦方差已经掺入到估计值中，即使它很小，也很难用加另一个估计值的方法摆脱它。

多重重要性采样(Multiple importance sampling, MIS)用了一种简单而又容易实现的方法解决了这些问题。其基本思想是，当我们估算一个积分时，我们应该从多个采样分布中采样，并且期望至少其中一个分布可以比较好地匹配被积函数，即使我们并不知道到底哪一个更好。MIS提供了一种可以对不同分布中的采样进行加权的方法，这样可以摆脱由于被积函数值跟采样密度不匹配而产生的很大的方差尖峰。这个方法甚至鼓励使用那些为特殊情况而设置的特殊采样例程，因为在这些特殊情况发生时，一般而言确实可以减小方差而又没有什么开销。

如果两个采样分布pf和pg被用来估算∫f(x)g(x)dx，MIS的新蒙特卡罗估计值为：

其中nf是来自分布pf的采样个数，ng是来自分布pg的采样个数，而wf和wg是特殊的加权函数，使得该估计值的期望值等于f(x)g(x)的积分值。

这里的加权函数考虑到了生成采样Xi和Yj的所有不同的方式，而不仅仅是我们实际上要用的特定的那个方式。其中一个很好的选择是使用均衡启发法(balance heuristic):

均衡启发法是一种可以证明的用采样加权来减小方差的好方法。

我们考虑一下该项在这样一种情形下的效果：在一个pf(X)值相对较小的点上从pf分布上抽取采样X。假定pf跟f(x)有相似的形状，那么f(X)的值也相对较小。我们假定g(X)有相对较大的值。那么标准的估计值f(X)g(X)/pf(X)就有很大的值，因为pf(X)很小，这样就会出现前面提到的方差问题。

然而，在使用均衡启发法的情况下，X的贡献值是：

只要pg分布跟g(x)比较匹配，那么分子就不会太小（因为有ngPg(X)项的作用），就会避免很大的方差尖峰，即使X的采样来自跟被积函数匹配很差的分布。

这里我们给出两个分布pf和pg这一特殊情形下的均衡启发法实现。在pbrt中，我们不会用到更一般情况下的多重分布。

<Monte Carlo Inline Functions> =

inline float BalanceHeuristic(int nf, float fPdf, int ng, float gPdf) {

return (nf * fPdf) / (nf * fPdf + ng * gPdf);

}

在实际应用中，使用幂函数启发法可以进一步地减小方差。对于指数β，幂函数启发法的公式为：

β = 2就是一个很好的值，我们将β = 2 硬写到代码实现中：

<Monte Carlo Inline Functions> +=

inline float PowerHeuristic(int nf, float fPdf, int ng, float gPdf) {

float f = nf * fPdf, g = ng * gPdf;

return (f*f) / (f*f + g*g);

}

15.5 对反射函数采样

因为重要性采样是一种非常有效的减小方差的技术，所以值得我们去努力寻找图形学中跟常见的被积函数相似的采样分布。在本节中，我们介绍pbrt中那些跟BSDF模型相似的分布，并对其对应的采样策略进行推演并加以实现。

BxDF::Sample_f()函数从跟相应的散射函数相似的分布中随机地进行方向采样。在第9.2节中，我们用这个函数来求来自理想镜面反射表面的反射光和透射光；在本节中，我们将表明这个采样过程是前一章所介绍的采样技术的一种特殊情况。

BxDF::Sample_f()在区间[0,1]x[0,1]中取两个采样值，并将它们用于基于变换的采样算法。由于它们是该函数的参数，这个调用该函数的例程就可以使用分层或低差异采样算法生成这两个采样，以保证这些方向有具有良好的分布。BxDF::Sample_f()所使用的生成方向的算法并不需要这些采样值（例如，利用拒绝法采样的时候），只需忽略这些参数就可以了。

该函数将采样得到的方向放入*wi，将p(ωi)的返回值放到*pdf。被选定的方向上的BSDF值被放入到一个Spectrum对象中并返回。被返回的PDF值是按照半球上的立体角还计算的，出射方向ωo和被采样到的入射方向ωi都是在标准反射坐标系下定义的。

该函数的缺省实现是在带余弦项加权的单位半球分布上采样的。对于所有那些不用delta分布来表示的BRDF，该函数都可以给出正确的结果。

<BxDF Method Definitions> +=

Spectrum BxDF::Sample_f(const Vector &wo, Vector * wi,

float u1, float u2, float *pdf) const {

<Cosine-sample the hemisphere, flipping the direction if necessary>

*pdf = Pdf(wo, *wi);

return f(wo, *wi);

}

我们需要解释一下一个跟反射坐标系中的方向的朝向相关的细节：由Malley法返回的法向量总是位于以(0,0,1)为球心的半球之内。如果方向ωo在相对的另一个半球之内，我们就必须反转ωi的方向使之跟ωo位于同一个半球之内。只是因为pbrt并不自动将法向量反转到跟ωo方向相同的一侧。

<Cosine-sample the hemisphere, flipping the direction if necessary> =

*wi = CosineSampleHemisphere(u1, u2);

if (wo.z < 0.) wi->z *= -1.f;

BxDF::Sample_f()返回所它选定的方向上的PDF值，而BxDF::Pdf()函数返回任意给定方向上的PDF值。这个函数对于多重重要性采样非常有用，因为我们必须能够在来自其它分布的方向采样上求PDF值。所以任何覆盖BxDF::Sample_f()函数的BxD子类也必须要覆盖BxDF::Pdf()函数，使得两者都返回相互一致的结果。

为了求带余弦加权的采样函数所使用的PDF值（我们前面提到过p(ω) = cosθ/π），首先要检查ωi和ωo是否在表面的同一侧；如果不在同一侧的话，其采样概率就是0。否则，该函数计算[n . ωi]，由于n = (0,0,1)，所以可以将其简化为求ωi的z坐标的绝对值。其中一个易错的地方是：ωo和ωi的次序是不可随意交换的。对于这里的余弦加权的分布而言，一般情况下p(ωi)不等于p(ωo)。在调用这个函数时一定要注意使用正确的变量次序。

<BSDF Inline Functions> +=

inline bool SameHemisphere(const Vector &w, const Vector &wp) {

return w.z * wp.z > 0.f;

}

<BxDF Method Definitions> +=

float BxDF::Pdf(const Vector &wo, const Vector &wi) const {

return SameHemisphere(wo, wi) ? fabsf(wi.z) * INV_PI : 0.f;

}

这个采样函数对Lambertian漫反射BRDF和Oren-Nayar模型都很不错，所以就不重写这些类的相应函数了。

15.5.1 对Blinn微平面分布采样

对基于微平面分布函数的BRDF采样要更困难一些。对于这些模型而言，BRDF是D,G,F三项的乘积，然后在除以两个余弦项（见第9章的公式9.8）。要找到跟这个完整的模型相匹配的概率分布是不现实的，所以我们要推导出一个单独的采样方法，用它从微平面分布函数D所描述的分布上采样。对于重要性采样而言这是一个很有效的采样策略，因为大多数的变化来自D项。

因此，所有MicrofacetDistribution的实现类必须实现从它们的分布上进行采样的函数，还有计算其PDF值的函数，它们的调用格式跟相应的BxDF函数相同，只是Sample_f()函数要返回void。

<MicrofacetDistribution Interface> +=

virtual void Sample_f(const Vector &wo, Vector *wi,

float u1 , float u2, float *pdf) const = 0;

virtual float Pdf(const Vector &wo, const Vector &wi) const = 0;

Microfacet BxDF将采样和PDF请求传递给相应的分布函数：

<BxDF Method Definitions> +=

Spectrum Microfacet::Sample_f(const Vector &wo, Vector *wi,

float u1 , float u2, float *pdf) const {

distribution->Sample_f(wo, wi, u1, u2, pdf);

if (!SameHemisphere(wo, *wi)) return Spectrum(0.f);

return f(wo, *wi);

}

<BxDF Method Definitions> +=

float Microfacet::Pdf(const Vector &wo, const Vector &wi) const {

if (!SameHemisphere(wo, wi)) return 0.f;

return distribution->Pdf(wo, wi);

}

Blinn的微平面分布函数也可以等价地写为D(cosθh) = (n+2) (cosθh)n，其中cosθh = | n . ωh|，D(ωh) = ((n+2)/2π) (cosθh)n 。为了用这个分布进行采样，我们先要按照这个分布对半角方向ωh进行采样，然后用它来将ωo按方向ωh做镜像反射从而得到ωi。

<BxDF Method Definitions> +=

void Blinn::Sample_f(const Vector &wo, Vector *wi, float u1, float u2,

float *pdf) const {

<Compute sampled half-angle vector ωh for Blinn distribution>

<Compute incident direction by reflecting about ωh>

<Compute PDF for ωi from Blinn distribution>

*pdf = blinn_pdf;

}

因为Φ的值不对D产生影响，这个分布的PDF，ph(θ,Φ)可以分解程ph(θ)和ph(Φ)。ph(Φ)是一个带1/2π的常量，则Φ的值的值可以通过下面公式来采样：

Φ = 2πξ2

我们总可以事先将微平面分布函数归一化，它们毕竟是分布函数。所以，PDF是：

ph(cosθh) = (n + 2) cosnθh

这是一个关于cosθh的幂函数分布，所以可以用一个均匀随机变量生成采样：

cosθh = (ξ1)1/n+1

由于pbrt将法向量变换到反射坐标系中的(0,0,1)，所以我们几乎总可以从球面坐标直接求出方向。唯一要处理的一个细节是：如果ωo在法向量的相反的一侧，就要反转半角向量。

<Compute sampled half-angle vector ωh for Blinn distribution> =

float costheta = powf(u1, 1.f / (exponent+1));

float sintheta = sqrtf(max(0.f, 1.f - costheta*costheta));

float phi = u2 * 2.f * M_PI;

Vector wh = SphericalDirection(sintheta, costheta, phi);

if (!SameHemisphere(wo, wh)) wh = -wh;

接下来我们只需使用向量镜像公式（如图）：

<Compute incident direction by reflecting about ωh> =

*wi = -wo + 2.f * Dot(wo, wh) * wh;

在计算采样方向上的PDF值时，我们要处理一个重要的细节。微平面分布给出的是半角向量(half-angle vector)附近的法向量分布，但辐射积分是关于入射向量(incoming vector)的，所以我们必须将半角PDF转换为入射角的PDF。换句话讲，我们必须用第14.4节所介绍的技术将关于ωh的密度函数转换为关于ωi的密度函数。为此，我们需要根据dωh / dωi的变化来做相应的调整。

我们通过如图的几何关系来给出两个分布的关系。

考虑朝向ωo的球面坐标系：立体角微分dωi和dωh分别是sinθidθidφi sinθhdθhdφh。所以有：

由于ωi是有ωo做ωh的镜像而得到的，故而有θi = 2θh。在加上Φi = Φh， 就可以得到所想要的转换因子：

<Compute PDF for ωi from Blinn distribution> =

float blinn_pdf = ((exponent + 1.f) * powf(costheta, exponent)) /

(2.f * M_PI * 4.f * Dot(wo, wh));

if (Dot(wo, wh) <= 0.f) blinn_pdf = 0.f;

有了以上这些工作，就可以很容易地给出Blinn::Pdf()的实现了，这里我们重用了前面计算PDF的代码：

<BxDF Method Definitions> +=

float Blinn::Pdf(const Vector &wo, const Vector &wi) const {

Vector wh = Normalize(wo + wi);

float costheta = AbsCosTheta(wh);

<Compute PDF for ωi from Blinn distribution>

return blinn_pdf;

}

15.5.2 对各向异性的微平面模型采样

Ashikhmin和Shirley（2000，2002）给出了对他们所定义的各向异性的微平面模型进行采样的公式。我们这里只陈述它们的结果而不做推导。有兴趣的读者可以直接参考他们的原始文献。

<BxDF Method Definitions> +=

void Anisotropic::Sample_f(const Vector &wo, Vector *wi,

float u1, float u2, float *pdf) const {

<Sample from first quadrant and remap to hemisphere to sample ωh>

<Compute incident direction by reflecting about ωh>

<Compute PDF for ωi from anisotropic distribution>

*pdf = anisotropic_pdf;

}

Anisotropic::sampleFirstQuadrant()中的采样函数在单位半球中的第一象限进行方向采样，也就是说，其球面角范围 (θ,Φ) 是[0, π/2] x [0, π/2]。我们只考虑这一部分的采样，并以此为基础来实现所有方向上的采样，这样做比较容易推导出这个分布的采样函数。这里我们用ξ1来选择一个象限，再将其映射到[0,1]之间来生成来自第一象限的采样，然后将采样映射到所需要的象限。我们要仔细地处理ξ1和Φ采样值的再映射以保证分层(stratification)效果：对于从一个象限到另一个象限的很小的ξ1变化，所生成的半角方向也应该只产生微小的变化。

<Sample from first quadrant and remap to hemisphere to sample ωh> =

float phi, costheta;

if (u1 < .25f) {

sampleFirstQuadrant(4.f * u1, u2, &phi, &costheta);

} else if (u1 < .5f) {

u1 = 4.f * (.5f - u1);

sampleFirstQuadrant(u1, u2, &phi, &costheta);

phi = M_PI - phi;

} else if (u1 < .75f) {

u1 = 4.f * (u1 - .5f);

sampleFirstQuadrant(u1, u2, &phi, &costheta);

phi += M_PI;

} else {

u1 = 4.f * (1.f - u1);

sampleFirstQuadrant(u1, u2, &phi, &costheta);

phi = 2.f * M_PI - phi;

}

float sintheta = sqrtf(max(0.f, 1.f - costheta*costheta));

Vector wh = SphericalDirection(sintheta, costheta, phi);

if (!SameHemisphere(wo, wh)) wh = -wh;

在半球的第一象限的半角向量采样可以用下面的分布函数：

注意，如果ex = ey并且分布是各向同性的，这里的分布即是前面Blinn分布所用到的采样技术。

<BxDF Method Definitions> +=

void Anisotropic::sampleFirstQuadrant(float u1, float u2,

float *phi, float *costheta) const {

if (ex == ey)

*phi = M_PI * u1 * 0.5f;

else

*phi = atanf(sqrtf((ex+1.f) / (ey+1.f)) *

tanf(M_PI * u1 * 0.5f));

float cosphi = cosf(*phi), sinphi = sinf(*phi);

*costheta = powf(u2, 1.f/(ex * cosphi * cosphi +

ey * sinphi * sinphi + 1));

}

跟前面的Blinn模型相似地，用这个函数所得到的采样方向的PDF值是归一化的半球上方向分布上的值，如前所述的Blinn::Pdf()那样，我们仍需将半角分布转换到入射角分布。

<BxDF Method Definitions> +=

float Anisotropic::Pdf(const Vector &wo, const Vector &wi) const {

Vector wh = Normalize(wo + wi);

<Compute PDF for ωi from anisotropic distribution>

return anisotropic_pdf;

}

<Compute PDF for ωi from anisotropic distribution> =

float costhetah = AbsCosTheta(wh);

float ds = 1.f - costhetah * costhetah;

float anisotropic_pdf = 0.f;

if (ds > 0.f && Dot(wo, wh) > 0.f) {

float e = (ex * wh.x * wh.x + ey * wh.y * wh.y) / ds;

float d = sqrtf((ex+1.f) * (ey+1.f)) * INV_TWOPI *

powf(costhetah, e);

anisotropic_pdf = d / (4.f * Dot(wo, wh));

}

15.5.3 对FRESNELBLEND采样

FresnelBlend是漫反射项和光泽反射项的一种混合型反射。有一个对该BSDF采样的很直接的方法就是，既对余弦加权的分布采样，也对微平面分布采样。这里的实现方法根据ξ1是小于0.5还是大于0.5来概率均等地选择这两个模型之一进行采样。做好这个决定之后，该实现将ξ1重映射到整个[0,1]区间。（否则，举例来说，用于余弦加权分布采样的ξ1总小于0.5了）。以这种方式对采样ξ1做双重用途的应用会稍微降低（ξ1，ξ2）的分层质量，但跟随机地使用一个未分层的随机变量来随机地选择采样方法相比起来，这根本就不是什么缺点了。

<BxDF Method Definitions> +=

Spectrum FresnelBlend::Sample_f(const Vector &wo, Vector *wi,

float u1, float u2, float *pdf) const {

if (u1 < .5) {

u1 = 2.f * u1;

<Cosine-sample the hemisphere, flipping the direction if necessary>

}

else {

u1 = 2.f * (u1 - .5f);

distribution->Sample_f(wo, wi, u1, u2, pdf);

if (!SameHemisphere(wo, *wi)) return Spectrum(0.f);

}

*pdf = Pdf(wo, *wi);

return f(wo, *wi);

}

这个采样策略的PDF函数很简单，只是两种PDF的平均而已。

<BxDF Method Definitions> +=

float FresnelBlend::Pdf(const Vector &wo, const Vector &wi) const {

if (!SameHemisphere(wo, wi)) return 0.f;

return .5f * (AbsCosTheta(wi) * INV_PI + distribution->Pdf(wo, wi));

}

15.5.4 镜面反射和透射

用于定义镜面反射BRDF和镜面透射BTDF的Dirac delta分布实际上可以很好地适应这里的采样框架，但在使用这些采样函数和PDF函数时，要需要注意一些约定。

回想一下Dirac delta分布定义为：

那么(x) = 0 对于所有 x != 0

并且∫-∞，+∞ δ(x) dx = 1

这就是一个概率密度函数，其中如果x不为0，则概率为0。生成这个分布的采样是很简单的，因为只有一个可能的值。如果以这个角度看待问题，那么我们就可以很自然地将SpectularReflection和SpectularTransmission的Sample_f()的实现放入蒙特卡罗采样框架之中。

确定PDF要返回哪一个值却不那么容易。严格地讲，delta分布并不是一个真正的函数，它必须定义为另一个函数的极限（例如这个函数可以是一个单位面积的宽度为0的方盒函数）。如果这样看问题的话，那么δ(0)就趋向于无穷大。当然，如果PDF返回一个无穷大或极大的值就不会得到正确的渲染结果。

回忆一下，如果BSDF有delta组成部分，那么这些组成部分的相应的fr函数也会有delta组成部分，在第9章的Sample_f()中，我们没有提到这一点。对于这样的散射方程的蒙特卡罗估计值可以写成：

其中ρhd(ωi)是半球-方向反射率函数(hemispherical-directional reflectance)，ω是理想镜面反射或透射方向。

因为PDF p (ωi)也有一个delta项，p(ωi) = δ(ω - ωi)，则两个delta分布就相互抵消了，则估计值就成了：

ρhd(ωo) Li(p, ω)

这正是Whitted积分器所计算出的值。

因此，当使用Sample_f()进行采样时，这里的实现用常数1作为镜面反射和透射所用的PDF值，这里约定对于镜面BxDF而言，当估算积分值时，PDF值内所隐含的delta分布会跟BSDF内所隐含的delta分布相互抵消。因此，Pdf()函数对所有的方向都返回0，因为另一个采样方法从delta随机地找到方向的概率为0。

<SpecularReflection Public Methods> +=

float Pdf(const Vector &wo, const Vector &wi) const {

return 0.;

}

<SpecularTransmission Public Methods> +=

float Pdf(const Vector &wo, const Vector &wi) const {

return 0.;

}

这个约定有个可能的陷阱：当使用多重重要性采样来计算权值时，那些带隐含delta分布的PDF值不可以跟普通的PDF值混在一起。在实际应用中这不是一个问题，因为如果被积函数有delta分布，就没有理由使用多重重要性采样。下一章中的光传输例程会有相应的逻辑来避免这个错误。

15.5.5 应用： 估算反射率

作为一个应用的例子，我们来看一下如何使用这些BxDF采样例程估算第9.1.1节中任意BRDF的反射率积分。例如，对于半球-方向反射率函数，有：

调用该函数的例程可选择提供采样个数和采样位置。如果没有提供的话，就使用拉丁超立方体采样模式。

BxDF::rho()函数利用它的采样函数以及重要性采样为任意的BxDF估算其反射率的蒙特卡罗估算值。

BxDF Method Definitions> +=

Spectrum BxDF::rho(const Vector &w, int nSamples,

float *samples) const {

if (!samples) {

samples = (float *)alloca(2 * nSamples * sizeof(float));

LatinHypercube(samples, nSamples, 2);

}

Spectrum r = 0.;

for (int i = 0 ; i < nSamples; ++i) {

<Estimate one term of phd>

}

return r / (M_PI * nSamples);

}

实际上我们可以很容易求出这个估计值，因为我们只需对反射率函数进行采样，求得相应的值，除以PDF值，就轻易地求出下列公式的每一项：

BxDF的Sample_f()函数返回所有的ωi，p(ωi)，和fr(p, ω, ωi)。唯一有点麻烦的地方是必须判断p(ωi) = 0的情况以防止被0除的异常。

<Estimate one term of phd>

Vector wi;

float pdf = 0.f;

Spectrum f = Sample_f(w, &wi, samples[2*i], samples[2*i+1], &pdf);

if (pdf > 0.) r += f * fabsf(wi.z) / pdf;

类似地我们可以估算半球-半球反射率函数：

对于下面的每项，要生成两个向量ω', ω''的采样：

<BxDF Method Definitions> +=

Spectrum BxDF::rho(int nSamples, float *samples) const {

if (!samples) {

samples = (float *)alloca(4 * nSamples * sizeof(float));

LatinHypercube(samples, nSamples, 4);

}

Spectrum r = 0.;

for (int i = 0 ; i < nSamples; ++i) {

<Estimate one term of p>

}

return r / (M_PI*nSamples);

}

这里的实现在半球上做均匀采样生成ω'的采样。第二个方向采样是用BxDF::Sample_f()函数生成。

<Estimate one term of Phh> =

Vector wo, wi;

wo = UniformSampleHemisphere(samples[4*i], samples [4*i + 1] ) ;

float pdf_o = INV_TW0PI, pdf_i = 0.f;

Spectrum f = Sample_f(wo, &wi, samples[4*i+2], samples[4*i+3], &pdf_i);

if (pdf_i > 0.)

r += f * fabsf(wi.z * wo.z) / (pdf_o * pdf_i);

15.5.6 对BSDF采样

有了对每个BxDF的采样函数之后，我们就可以定义BSDF类的采样函数了。当SurfaceIntegrator要按照BSDF的分布进行采样是，就要调用该函数；而它再对相应的BxDF::Sample_f()函数进行调用来生成采样。BSDF存放了一个或多个BxDF的指针，我们要用到的密度函数是每个密度函数的累加。

p(ω) = 1/N Σi=1,N p(ωi)

BSDF::Sample_f()函数使用了三个随机变量来驱动这个采样过程。第一个用于对密度函数的选择，另外两个被传送到所选定的采样函数。

<BSDF Method Definitions> +=

Spectrum BSDF::Sample_f(const Vector &woW, Vector *wiW,

float u1, float u2, float u3, float *pdf,

BxDFType flags, BxDFType *sampledType) const {

<Choose which BxDF to sample>

<Sample chosen BxDF>

<Compute overall PDF with all matching BxDFs>

<Compute value of BSDF for sampled direction>

}

该函数先要确定使用哪一个BxDF采样函数。由于调用者必须传入跟所选BxDF相匹配的标志，所以实际过程就很复杂。因而，这里只用到了可用的采样密度函数的一个子集。

<Choose which BxDF to sample> =

int matchingComps = NumComponents(fIags);

if (matchingComps == 0) {

*pdf = 0.f;

return Spectrum(0.f);

}

int which = min(Floor2Int(u3 * matchingComps), matchingComps-1);

BxDF *bxdf = NULL;

int count = which;

for (int i = 0; i < nBxDFs; ++i)

if (bxdfs[ i ]->MatchesFlags(flags))

if (count-- == 0) {

bxdf = bxdfs;

break;

}

一旦选定了适当的BxDF后，就会调用BxDF::Sample_f()。回忆一下这些函数所用到和要返回的向量都是在BxDF的局部坐标系中定义的，所以先要将传入的向量变换到BxDF的坐标系中，所要返回的向量也要被变换到世界坐标系中。

<Sample chosen BxDF> =

Vector wi;

Vector wo = WorldToLocal(woW);

*pdf = 0.f;

Spectrum f = bxdf->Sample_f(wo, &wi, u1, u2, pdf);

if (*pdf == 0.f) return 0.f;

if (sampledType) *sampledType = bxdf->type;

*wiW = LocalToWorld(wi);

为了计算采样方向ωi的实际PDF，我们要将所有用到的BxDF的PDF做平均。因为*pdf已经有了采样所用的分布的PDF值，我们只需加上其它分布的贡献值。如果所选的BxDF是理想镜面反射，就要跳过这一步，因为PDF有隐含的delta分布（否则，加上其它PDF值就产生不正确的结果）。

<Compute overall PDF with all matching BxDFs> =

if (!(bxdf->type & BSDF_SPECULAR) && matchingComps > 1) {

for (int i = 0; i < nBxDFs; ++i) {

if (bxdfs[ i ] != bxdf &&

bxdfs [ i ]->MatchesFlags(flags))

*pdf += bxdfs[ i ]->Pdf(wo, wi);

}

}

if (matchingComps > 1) *pdf /= matchingComps;

有了方向采样以后，该函数需要计算方向对(ωo, ωi)的BSDF值，所有相关的BSDF组成都要对这个方向对进行计算，除非被采样的方向来自于镜面反射，这时就要用Sample_f()早些时候计算出的值。

<Compute value of BSDF for sampled direction> =

if (!(bxdf->type & BSDF_SPECULAR)) {

f = 0.;

if (Dot(*wiW, ng) * Dot(woW, ng) > 0) // ignore BTDFs

flags = BxDFType(flags & ~BSDF_TRANSMISSION);

else / / ignore BRDFs

flags = BxDFType(flags & ~BSDF_REFLECTION);

for (int i = 0; i < nBxDFs; ++i)

if (bxdfs[ i ]->MatchesFlags(flags))

f += bxdfs[ i ]->f(wo, wi);

}

return f;

BSDF::Pdf()要做类似的计算，对这一组BxDF做循环并调用相应的Pdf()函数来求任意方向采样的PDF值。

<BSDF Public Methods> +=

float Pdf(const Vector &wo, const Vector &wi,

BxDFType flags = BSDF_ALL) const;