PBRT笔记(12)——蒙特卡洛积分

这里还涉及到pdf、方差等概念，推荐去看《全局光照技术:从离线到实时渲染》
积累分布函数 cumulative distribution function (CDF)

蒙特卡洛估算

为了计算式蒙特卡罗估算量，就有必要从选择的概率分布中抽取随机样本。

逆推法

逆推法使用一个或多个均匀的随机变量映射到随机变量的期望分布中。

为了从该分布中获取样本，我们首先计算CDF P(x),在这个函数是连续的情况下，P表示为p的不定积分。

当从该分布中获取样本时，我们可以使用均匀随机数ξ，并根据CDF选取某一可能结果。对此，其中一种方法是利用自身概率选取某一特定结果，如图13.3所示，时间概率投射至垂直轴上，随机变量ξ沿其进行取值。不难发现，这是从正确的分布中得出的——均匀样本击中任何特定条棒的概率恰好等于该条棒的高度。

为了将这种技术运用到到连续分布中，考虑当离散概率趋于无穷时的情形。图13.1中的PDF会变成了一条光滑的曲线，而图13.2中的CDF则会变成了它的积分。尽管函数是连续的，但投影拥有以下的数学解释——逆CDF。使用以上解释，通过随机数ξ，计算该位置的逆CDF值，因此被称为反演法。

准确地说，可以通过以下步骤根据任意PDF获取样本值。

计算CDF， $P(x)=\int^x_0 p(x')dx'$
计算逆CDF， $P^{-1}(x)$
获取均匀分布的随机数ξ
计算样本， $P^{-1}(ξ)$

舍选法

对于，某些函数，可能无法通过对其进行积分从而获得pdf，或者无法通过计算获得逆cdf。舍选法是一种无需进行以上任意一个步骤就能根据函数分布生成样本的技术。本质上它是一种投飞镖法，假设我们想从某个函数f(x)中抽取样本，但是我们有一个pdf p(x)满足f(x)< c p(x),假设我们知道如何从p中获取样本，那么舍选法为：

loop forever:
    sample X from p’s distribution
    if ξ < f (X)/(c p(X)) then
        return X

这个过程反复选取一对随机变量(X,ξ)，如果点(X,ξc p(X))是在f(X)下面,那么接受样本X。否则，它将被拒绝，并再次选取一个新的样本对。这个方法的效率取决于f(x)与cp(x)的接近程度，越相似收敛地越快。而且适用于任意维度的函数。

分布间的转换

在讨论逆推法时，我们曾采用一种使用特定方式转换均匀随机变量的分布，来生成样本。在这里，我们将研究一个更普遍的问题，即我们将一个任意分布的样本集从转换成另一个基于函数f(x)分布的样本集。

多维转换

在通用n维的情况下，类似的推导给出了不同密度之间的相似关系。假设有一个n维随机变量X以及它的密度函数是 $P_x(x)$ ,令Y = T(x),其中T表示为双射关系(既是单射又是满射的映射称为双射)。此处，pdf通过下列式子进行关联：

$P_y(y)=P_y(T(x))=\frac{P_x(x)}{\left| J_T(x)\right|}$

其中JT表示为T的雅可比矩阵矩阵。

$\left[\begin{matrix} \partial T_1/\partial x_1 & ... & \partial T_1/\partial x_n \\ ... & ... & ... \\ \partial T_n/\partial x_1 & ... & \partial T_n/\partial x_n \\ \end{matrix} \right]$

此处Ti通过T(x)=(T1(x),……,Tn(x))加以定义。

极坐标

极坐标转换通过下列式表示：

$x=rcos\theta$ $y=rsin\theta$

根据密度分布函数p(r,θ)获取样本。dang

$J_T=\left[\begin{matrix} \frac{\partial_x}{\partial_r} & \frac{\partial_x}{\partial_\theta} \\ \frac{\partial_y}{\partial_r} & \frac{\partial_y}{\partial_\theta} \end{matrix} \right]=\left[\begin{matrix} cos\theta & -rsin\theta \\ sin\theta & rcos\theta \end{matrix} \right]$

对应 $r(cos^2 \theta+sin^2 \theta)=r$ ,因而有 $p(x,y)=p(r,\theta)/r$ ,采样策略始于笛卡尔坐标系，并将其在极坐标系中转换，对此有： $rp(x,y)=p(r,\theta)$

球体坐标系

使用方向变量来表达球体坐标系：

$x=rsin\theta cos\phi$ $y=rsin\theta sin\phi$ $z=rcos\theta$

这个雅可比变换矩阵的决定式为|JT| = r2 sin θ,所以对应的pdf为：

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

这个变换胃肠重要，因为这样我就可以通过p(x,y,z)来表示球面上的点。回忆一下，立体角为单位的标准球面点集，在球面坐标系下： $d\omega=sin\theta d\theta d\phi$

如果我们定义一个立体角为 $\Omega$ 的密度函数，那么：

$Pr\left\{ \omega \in \Omega\right\}=\int_\Omega p(\omega)d\omega$

因此对于θ与φ的pdf可有以下推导：

$p(\theta,\phi)d\theta d\phi=p(\omega)d \omega$ $p(\theta,\phi)=sin\theta p(\omega)$

具有多维变换的二维采样

假设我们有一个二维联合密度函数p(x,y)，并使用它进行采样(x,y)。有时，多维密度呈可分离状态，并可以表示为一维密度的乘积。例如 $p(x,y)=p_x(x)p_y(y)$

其中随机变量(x,y)可以分别通过Px与Py分别计算对应密度。但相应的有的函数是不可分离的。当给定二维密度函数后，临界概率密度函数p(x)通过对其中一个维度“积分”得到:

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

这可以被认为是单独的p(x)密度函数,更准确地说，这是全部y值上针对特定x的平均密度。条件密度函数p(y|x)是给定某一特定x的情况下y的密度函数：

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

从联合分布中进行二维采样的基本思想是：首先计算临界概率密度函数，分离出一个特定的变量，然后使用标准一维技术从该密度中抽取样本。一旦有了样本，就可以计算给定该值的条件密度函数并使用同样的标准的一维采样技术从该分布中获取样本。

半球均匀采样

考虑在以立体角为准在均匀半球选择方向。均匀分布意味着密度函数为常数，因而可以得知 p(ω) = c。结合密度函数必须在定义域上积分到1的事实，我们可以得到：

$\int_{H^2}p(\omega)d\omega=1\rightarrow c\int_{H^2}d\omega=1 \rightarrow c=\frac{1}{2 \pi}$

以此可知：p(ω)=1/(2π)或者p(θ,φ)=sinθ/(2π)。需要注意的是这个密度函数是可分离的。所以可以使用上面所说的方法。对于θ：

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

计算条件条件密度函数：

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

这里使用逆推法对各个pdf进行采样：

$P(\theta)=\int^{\theta}_{0}sin\theta'd\theta'=1-cos\theta$ $P(\phi | \theta)=\int^{\phi}_{0}\frac{1}{2\pi}d\phi'=\frac{\phi}{2\pi}$

这里使用1-ξ来代替ξ：

$\theta=cos^{-1}\xi_1$$ $$\phi=2\pi \xi_2$

之后转化成笛卡尔坐标系：

$x=sin\theta cos\phi=cos(2\pi \xi_2)\sqrt{1-\xi_1^2}$ $y=sin\theta sin\phi=sin(2\pi \xi_2)\sqrt{1-\xi_1^2}$ $z=cos\theta=\xi_1$

Vector3f UniformSampleHemisphere(const Point2f &u) {
    Float z = u[0];
    Float r = std::sqrt(std::max((Float)0, (Float)1. - z * z));
    Float phi = 2 * Pi * u[1];
    return Vector3f(r * std::cos(phi), r * std::sin(phi), z);
}

圆盘采样

当均匀采样时：p(x,y)=1/π,转化成极坐标后：p(r,θ)=r/π。以此可以计算临界以及条件密度函数：

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

对于r与θ取随机值：

$r=\sqrt {\xi_1}$$ $$\theta=2\pi \xi_2$

Point2f UniformSampleDisk(const Point2f &u) {
    Float r = std::sqrt(u[0]);
    Float theta = 2 * Pi * u[1];
    return Point2f(r * std::cos(theta), r * std::sin(theta));
}

但是这么做会使得圆盘变形，所以这里会做以下映射： $r=x$ $\theta=\frac{y}{x}/\frac{\pi}{4}$

Point2f ConcentricSampleDisk(const Point2f &u) {
    // Map uniform random numbers to $[-1,1]^2$
    Point2f uOffset = 2.f * u - Vector2f(1, 1);

    // Handle degeneracy at the origin
    if (uOffset.x == 0 && uOffset.y == 0) return Point2f(0, 0);

    // Apply concentric mapping to point
    Float theta, r;
    if (std::abs(uOffset.x) > std::abs(uOffset.y)) {
        r = uOffset.x;
        theta = PiOver4 * (uOffset.y / uOffset.x);
    } else {
        r = uOffset.y;
        theta = PiOver2 - PiOver4 * (uOffset.x / uOffset.y);
    }
    return r * Point2f(std::cos(theta), std::sin(theta));
}

余弦加权半球采样

执行规范化操作后：

$\int_{H^2}cp(\omega)d\omega=1$ $\int^{2\pi}_0\int^{\frac{\pi}{2}}_{0}c cos\theta sin\theta d\theta d\phi=1$ $c2\pi\int^{\pi/2}_{0}cos\theta sin\theta d\theta=1$ $c=\frac{1}{\pi}$

因而有：

$p(\theta,\phi)=\frac{1}{\pi}cos\theta sin\theta$

这里可以使用Malley法来生成余弦权重点。即首先在圆盘上生成均匀的点，之后再投影至半球上，进而生成方向向量，最终得到含有余弦分布的方向分布。

Vector3f CosineSampleHemisphere(const Point2f &u) {
    Point2f d = ConcentricSampleDisk(u);
    Float z = std::sqrt(std::max((Float)0, 1 - d.x * d.x - d.y * d.y));
    return Vector3f(d.x, d.y, z);
}

圆锥体采样

无论是基于球面的区域光源还是基于聚光灯的区域光源，在方向锥上均匀地采样射线对两者是很有用的。
这个分布(θ,φ)是可分离的,p(φ) = 1/(2π)，因此我们需要得到一个方法示例θ方向一致的锥方向围绕一个中心方向的最大角,θmax。

$1=c\int^{\theta max}_{0}sin\theta d\theta=c(1-cos\theta_{max})$

通过积分pdf的方式得到cdf，再通过其进行采样：

$cos\theta=(1-\xi)+\xi cos\theta_{max}$

(0,0,1)轴采样

Vector3f UniformSampleCone(const Point2f &u, Float cosThetaMax) {
    Float cosTheta = ((Float)1 - u[0]) + u[0] * cosThetaMax;
    Float sinTheta = std::sqrt((Float)1 - cosTheta * cosTheta);
    Float phi = u[1] * 2 * Pi;
    return Vector3f(std::cos(phi) * sinTheta, std::sin(phi) * sinTheta,cosTheta);
}

Vector3f UniformSampleCone(const Point2f &u, Float cosThetaMax,const Vector3f &x, const Vector3f &y,const Vector3f &z) {
    Float cosTheta = Lerp(u[0], cosThetaMax, 1.f);
    Float sinTheta = std::sqrt((Float)1. - cosTheta * cosTheta);
    Float phi = u[1] * 2 * Pi;
    return std::cos(phi) * sinTheta * x + std::sin(phi) * sinTheta * y +cosTheta * z;
}

三角形采样

我们将用(u, v)来表示这两个质心坐标,因为基于面积的采样，所以pdf p(u,v)一定等于面积的倒数，面积为1/2，所以p(u,v)=2

其临界概率密度函数为：

$p(u)=\int^{1-u}_0 p(u,v)dv=2\int^{1-u}_0dv=2(1-u)$

所以cdf为：

$p(v|u)=\frac{p(u,v)}{p(u)}=\frac{2}{2(1-u)}=\frac{1}{1-u}$

可以分别得到以下积分：

$P(u)=\int^{u}_0p(u')du'=2u-u^2$ $P(v)=\int^v_0p(v'|u)dv'=\frac{v}{1-u}$

得到最终的抽样策略:

$u=1-\sqrt{\xi_1}$ $u=\xi_2\sqrt{\xi_1}$

Point2f UniformSampleTriangle(const Point2f &u) {
    Float su0 = std::sqrt(u[0]);
    return Point2f(1 - su0, u[1] * su0);
}

分段的二位常量分布（贴图）

二维函数f(u,v)通过一个数量为 $n_u\times n_v$ 的 $f[u_i,v_j]$ 来定义。f的积分为 $f[u_i,v_j]$ 值的简单求和：

$I_f=\iint f(u,v)dudv=\frac{1}{n_u n_v} \sum^{n_u-1}_{i=0}\sum^{n_v-1}_{j=0} f[u_i,v_j]$

使用pdf的定义以及f的积分，我们可以得到f的pdf：

$p(u,v)=\frac{f(u,v)}{\iint f(u,v)dudv}=\frac{f[\hat{u},\hat{v}]}{1/(n_un_v)\sum_i \sum_j f[u_i,v_j]}$

回忆一下公式，临界密度函数可以作为 $f[u_i,v_j]$ 值之和计算：

$p(v)=\int p(u,v)du=\frac{(1/n_u)\sum_i f[u_i,\hat{v}]}{I_f}$

因为这个函数只取决于 $\hat{v}$ ,它奔丧是一个一维的分段常数函数，有了p(v)就可以计算p(u|v):

$p(u|v)=\frac{p(u,v)}{p(v)}=\frac{f[\hat{u},\hat{v}]/I_f}{p[\hat{v}]}$

具体实现在Distribution2D类中，构造函数会先计算 $p[\hat{u}|\hat{v}]$ ,之后计算临界密度函数 $p[\hat{v}]$

Distribution2D::Distribution2D(const Float *func, int nu, int nv) {
    pConditionalV.reserve(nv);
    for (int v = 0; v < nv; ++v) {
        // Compute conditional sampling distribution for $\tilde{v}$
        pConditionalV.emplace_back(new Distribution1D(&func[v * nu], nu));
    }
    // Compute marginal sampling distribution $p[\tilde{v}]$
    std::vector<Float> marginalFunc;
    marginalFunc.reserve(nv);
    for (int v = 0; v < nv; ++v)
        marginalFunc.push_back(pConditionalV[v]->funcInt);
    pMarginal.reset(new Distribution1D(&marginalFunc[0], nv));
}

Point2f SampleContinuous(const Point2f &u, Float *pdf) const {
    Float pdfs[2];
    int v;
    Float d1 = pMarginal->SampleContinuous(u[1], &pdfs[1], &v);
    Float d0 = pConditionalV[v]->SampleContinuous(u[0], &pdfs[0]);
    *pdf = pdfs[0] * pdfs[1];
    return Point2f(d0, d1);
}
Float Pdf(const Point2f &p) const {
    int iu = Clamp(int(p[0] * pConditionalV[0]->Count()), 0, pConditionalV[0]->Count() - 1);
    int iv =
        Clamp(int(p[1] * pMarginal->Count()), 0, pMarginal->Count() - 1);
    return pConditionalV[iv]->func[iu] / pMarginal->funcInt;
}

俄罗斯赌博轮盘与划分机制

俄罗斯赌博轮盘与划分机制通过增加各采样的贡献率来提高蒙特卡洛估算的效率。对于俄罗斯赌博轮盘，pbrt这里举了个表面关照的例子：

分析当前函数
寻找可以省略计算的地方：光线被遮挡或是函数积分值较小时（BRDF值较小或者光线与表面的cos值较小）
需要注意的是如果是在较为黑暗的地方，这个方法就需要改良（会出现一些异常的亮点）。

对于划分机制是将原本一个完整的图像函数划分成多个图像采样结果之和，大概是延迟渲染的意思。

采样点分布策略

好的采样点分布会有效的提高效率。
这里推荐去看文刀秋二的文章：
https://zhuanlan.zhihu.com/p/20197323

重要性采样

重要性采样时一种缩小方差的技术：

$F_N=\frac{1}{N}\sum^N_{i=1}\frac{f(X_i)}{p(x_i)}$

如果被积函数f(x)的分布与p(x)越相似，那么蒙特卡洛估算就收敛越快。其基本思想是，通过将工作集中在被积函数值相对较高的地方，可以更有效地计算出准确的估计值。

也就是说如果采样点分布形状与被积函数相似，那方差将会被缩减。