pbrt-v4 Ep. XIV: 光线传播: 体渲染

辐射转移方程

#

LTE忽略了介质, 而辐射转移方程(radiative transfer equation, RTE)会考虑介质, 因此也被叫做体积光线传播方程(volumetric light transport equation). RTE通常为积分微分方程, 因为外散射过程涉及到积分. 令$p’=p+t\omega$, 其余符号见之前章节, 此时RTE形式如下.

$$\begin{equation} \frac{\partial L_o(p’,\omega)}{\partial t} = -\sigma_t(p’,\omega)L_i(p’,-\omega)+\sigma_t(p’,\omega)L_s(p’,\omega) \end{equation}$$

对于$L_i(p,\omega)$, 它的变化是体渲染的逆过程, 等式右侧需要乘上$-1$, 解积分方程可得如下结果.

$$\begin{equation} L_i(p’,\omega)=\frac{-\int_0^{t} T_r(p’’ \to p)\sigma_t(p’’,\omega)L_s(p’’,-\omega)dt’+C}{T_r(p’ \to p)} \end{equation}$$

若第一个散射点位于$t_s$, 该处可以是由表面BSDF表示的散射, 也可以是介质相位方程表示的散射. 由于$L_i(p_s, \omega)=L_o(p_s,-\omega)$, 此时得到如下结果.

$$\begin{equation} L_i(p’,\omega)=\frac{-\int_0^t T_r(p’’ \to p)\sigma_t(p’’,\omega)L_s(p’’,-\omega)dt’+L_o(p_s,-\omega)T_r(p_s \to p)+\int_0^{t_s} T_r(p’’ \to p)\sigma_t(p’’,\omega)L_s(p’’,-\omega)dt’}{T_r(p’ \to p)} \end{equation}$$

根据上式可得$L_i(p,\omega)$, 若光路上没有任何交点, 则无穷远处辐亮度为0, 只剩下积分项.

$$\begin{equation} L_i(p,\omega)=L_o(p_s, -\omega)T_r(p_s \to p)+\int_0^{t_s} \sigma_t(p’,\omega)T_r(p’ \to p)L_s(p’,-\omega)dt \end{equation}$$

空散射扩展

#

Monte Carlo需要根据与$L_i(p,\omega)$中积分项相似的分布来提高采样效率, 异质介质的$\sigma_t(p)$是复杂的分布不利于从分布中采样. 如果在一定范围内添加空散射$\sigma_n(p)$构造同质介质$\sigma_{\text{maj}}(p)$, 即将$\sigma_t(p)$转为常数或分段常数$\sigma_{\text{maj}}(p)$, 可以有效提高渲染效率. 在RTE右侧添加$-\sigma_n(p)L_i(p,\omega)+\sigma_n(p)L_i(p,\omega)$可得如下结果, 此时衰减过程使用$\sigma_{\text{maj}}(p)$, 光线传播过程中在自发光与散射之外也会在任意$\sigma_n(p)$不为0处发生空散射.

$$\begin{equation} \begin{aligned} \frac{\partial L_i(p’,\omega)}{\partial t} &=\sigma_{\text{maj}}(p’)L_i(p’,\omega)-\sigma_{\text{maj}}(p’)L_n(p’,-\omega)\\ &=(\sigma_a(p’,\omega)+\sigma_s(p’,\omega)+\sigma_n(p’,\omega))L_i(p’,\omega)-\sigma_a(p’,\omega)L_e(p,-\omega)-\sigma_s(p’,\omega)\int_\Theta p(p’,\omega_i,-\omega)L_i(p’,\omega_i)d\omega_i-\sigma_n(p’,\omega)L_i(p’,\omega) \end{aligned} \end{equation}$$

此时可以得到相交点为$p_s$时在$p$处的$L_i(p,\omega)$.

$$\begin{equation} L_i(p,\omega)=L_o(p_s,-\omega)T_{\text{maj}}(p_s \to p)+\int_0^{t_s} \sigma_{\text{maj}}(p’)T_{\text{maj}}(p’ \to p)L_n(p’,-\omega)dt \end{equation}$$

辐射转移方程求解

#

添加空散射后重要性抽样的PDF为$\sigma_{\text{maj}}(p)e^{-\int_0^t \sigma_{\text{maj}}(p’) dt’}$, 令选择$L_i(p,\omega)$第一项的概率为$q=T_{\text{maj}}(p_s \to p)$, 其求解方法如下. 由于之前章节所介绍的归一化, 这里$1-q$被抵消, 同时$L_n$可以按照上述方式求解.

$$\begin{equation} \begin{aligned} L_i(p,\omega)\approx \begin{cases} \frac{T_{\text{maj}}(p_s \to p)L_o(p_s,\omega)}{q}=L_o(p_s,-\omega) & \text{with probability}\ q\\ \frac{\int_0^t \sigma_{\text{maj}}(p’)T_{\text{maj}}(p’ \to p)L_n(p’,-\omega)dt}{1-q}\approx\frac{\sigma_{\text{maj}}(p’)T_{\text{maj}}(p’ \to p)L_n(p’,-\omega)}{p(p’)}=L_n(p’,-\omega) & \text{otherwise} \end{cases} \end{aligned} \end{equation}$$

采样主透射率

#

pbrt中Medium是$\sigma_{\text{maj}}$的分段函数, 积分可以按如下方式表示. 此时可以执行递归采样, 首先采样$T_{\text{maj}}^1$, 若得到的$t’_1<t_1$则采样第一项, 采样结果为$\frac{\sigma_{\text{maj}}^1 T_{\text{maj}}^1(p \to p’)f(p’)}{p_1(t’_1)}=f(p’)$, 若$t’_1>t_2$则采样第二项, 采样结果为$\frac{T_{\text{maj}}^1(p \to p_1)\sigma_{\text{maj}}^2 T_{\text{maj}}^2(p \to p’)f(p’)}{P(t’_1 > t_1)p_1(t’_1)}=f(p’)$, 如此一直递归, 可以看出采样结果始终为$f(p’)$.

$$\begin{equation} \begin{aligned} &\int_0^t \sigma_{\text{maj}}(p’)T_{\text{maj}}(p \to p’)f(p’)dt’\\ &=\sigma_{\text{maj}}^1\int_0^{t_1}T_{\text{maj}}^1(p \to p’)f(p’)dt’\\ &+T_{\text{maj}}^1(p \to p_1)\sigma_{\text{maj}}^2 \int_{t_1}^{t_2}T_{\text{maj}}^2(p_1 \to p’)f(p’)dt’\\ &+T_{\text{maj}}^1(p \to p_1)T_{\text{maj}}^2(p_1 \to p_2)\sigma_{\text{maj}}^3 \int_{t_2}^{t_3}T_{\text{maj}}^3(p_2 \to p’)f(p’)dt’\\ &+ \cdots \end{aligned} \end{equation}$$

pbrt通过在SampleT_maj中实现$\sigma_{\text{maj}}$的采样, 每个样本都通过当前$\sigma_{\text{maj}}$对应的指数分布采样获取, 其中T_maj是当前样本与上个样本之间的$T_{\text{maj}}$. 也就是说, 这里第$n$次采样对应$\sigma_{\text{maj}}^n$, 而不是分段函数$\sigma_{\text{maj}}(p)$的每一段对应一个$\sigma_{\text{maj}}^n$, 可以看出在每段内使用指数分布仍然是符合重要性抽样的. 由于比率跟踪等算法需要多个样本, pbrt通过callback不断获取样本并返回是否需要更多样本.

template <typename ConcreteMedium, typename F>
PBRT_CPU_GPU SampledSpectrum SampleT_maj(Ray ray, Float tMax, Float u, RNG &rng,
                                         const SampledWavelengths &lambda, F callback) {
    // Normalize ray direction and update _tMax_ accordingly
    tMax *= Length(ray.d);
    ray.d = Normalize(ray.d);

    // Initialize _MajorantIterator_ for ray majorant sampling
    ConcreteMedium *medium = ray.medium.Cast<ConcreteMedium>();
    typename ConcreteMedium::MajorantIterator iter = medium->SampleRay(ray, tMax, lambda);

    // Generate ray majorant samples until termination
    SampledSpectrum T_maj(1.f);
    bool done = false;
    while (!done) {
        // Get next majorant segment from iterator and sample it
        pstd::optional<RayMajorantSegment> seg = iter.Next();
        if (!seg)
            return T_maj;
        // Handle zero-valued majorant for current segment
        if (seg->sigma_maj[0] == 0) {
            Float dt = seg->tMax - seg->tMin;
            // Handle infinite _dt_ for ray majorant segment
            if (IsInf(dt))
                dt = std::numeric_limits<Float>::max();

            T_maj *= FastExp(-dt * seg->sigma_maj);
            continue;
        }

        // Generate samples along current majorant segment
        Float tMin = seg->tMin;
        while (true) {
            // Try to generate sample along current majorant segment
            Float t = tMin + SampleExponential(u, seg->sigma_maj[0]);
            PBRT_DBG("Sampled t = %f from tMin %f u %f sigma_maj[0] %f\n", t, tMin, u,
                     seg->sigma_maj[0]);
            u = rng.Uniform<Float>();
            if (t < seg->tMax) {
                // Call callback function for sample within segment
                PBRT_DBG("t < seg->tMax\n");
                T_maj *= FastExp(-(t - tMin) * seg->sigma_maj);
                MediumProperties mp = medium->SamplePoint(ray(t), lambda);
                if (!callback(ray(t), mp, seg->sigma_maj, T_maj)) {
                    // Returning out of doubly-nested while loop is not as good perf. wise
                    // on the GPU vs using "done" here.
                    done = true;
                    break;
                }
                T_maj = SampledSpectrum(1.f);
                tMin = t;

            } else {
                // Handle sample past end of majorant segment
                Float dt = seg->tMax - tMin;
                // Handle infinite _dt_ for ray majorant segment
                if (IsInf(dt))
                    dt = std::numeric_limits<Float>::max();

                T_maj *= FastExp(-dt * seg->sigma_maj);
                PBRT_DBG("Past end, added dt %f * maj[0] %f\n", dt, seg->sigma_maj[0]);
                break;
            }
        }
    }
    return SampledSpectrum(1.f);
}

泛化路径空间

#

这里主要参考Light transport on path-space manifolds第三章和A null-scattering path integral formulation of light transport, 证明涉及到泛函分析, 本人水平不够, 这里只讨论证明结果.

令$A$为表面空间, $V$为体积空间, $V_\emptyset$为空散射空间, 此时路径空间定义如下.

$$\begin{equation} P=\bigcup_{n=1}^{\infty}(A \cup V \cup V_\emptyset)^{n} \end{equation}$$

对于长度为$n$的路径, 其微分项如下.

$$\begin{equation} \begin{aligned} d\bar{p}_n&=\prod_{i=1}^n dp_i\\ dp_i&= \begin{cases} dA(p_i) & p_i \in A\\ dV(p_i) & p_i \in V\\ dV_\emptyset(p_i) & p_i \in V_\emptyset \end{cases} \end{aligned} \end{equation}$$

立体角与表面积微分的转换之前介绍过, 由于引入了空间中的介质, 我们需要立体角与体积微分的转换关系. 将体积积分视为在一个不断扩张的球上积分, 球表面积微分与立体角转换关系为$d\omega=\frac{dA}{t^2}$, 通过Fubini定理可以得到体积微分.

$$\begin{equation} \begin{aligned} &\int_\Theta\int_0^\infty f(p’)V(p \leftrightarrow p’) dt d\omega\\ =&\int_0^\infty \frac{1}{t^2}\int_\Theta f(p’)V(p \leftrightarrow p’) dA(p’) dt\\ =&\int_V \frac{V(p \leftrightarrow p’)}{\Vert p - p’ \Vert^2} f(p’) dV(p’) \end{aligned} \end{equation}$$

对于空散射体积空间, 由于不发生散射以改变光路, 空散射顶点必须位于相邻的实散射顶点形成的边上, 此时通过Dirac delta分布表示.

$$\begin{equation} dV_\emptyset(p_i)=d\delta_{p_i^{r-} \leftrightarrow p_i^{r+}}(p_i) \end{equation}$$

可以看出体积微分与表面积微分相比缺失了余弦项, 同时不难看出只有某个位于表面时, 其上的余弦项才会被考虑, 此时可以定义泛化几何方程. 空散射由于不改变光路不需要求解几何方程.

$$\begin{equation} \begin{aligned} \hat{G}(p \leftrightarrow p’)&=V(p \leftrightarrow p’)\frac{D(p,p’)D(p’,p)}{\Vert p-p’ \Vert^2}\\ D(x, y)&= \begin{cases} |n_x \cdot \frac{x-y}{\Vert x-y \Vert}| & p \in A\\ 1 & p \in V \end{cases} \end{aligned} \end{equation}$$

在介质中同样可以发生散射, 这由相位方程和散射率决定, 因此需要定义泛化的BSDF, 其中$\chi+$为Heaviside函数, 参数为正时为1否则为0, 对于空散射这始终为1, 用于限制空散射顶点位于边上.

$$\begin{equation} \begin{aligned} \hat{f}(p_{i+1} \to p_i \to p_{i-1})= \begin{cases} f(p_{i+1} \to p_i \to p_{i-1}) & p_i \in A\\ \sigma_s(p_i)p(p_{i+1} \to p_i \to p_{i-1}) & p_i \in V\\ \sigma_n(p_i)\chi^+((p_i-p_{i-1})(p_{i+1}-p_i)) & p_i \in V_\emptyset \end{cases} \end{aligned} \end{equation}$$

此时可以得到泛化路径通量, 其中$m$为实散射顶点数量, $r_i$为实散射点, 透射率求和从0开始是因为$p_0$已知.

$$\begin{equation} \hat{T}(\bar{p}_n)=\prod_{i=1}^{n-1}\hat{f}(p_{i+1} \to p_i \to p_{i-1})\prod_{i=0}^{n-1}T_r(p_{i+1} \to p_i)\prod_{i=1}^{m-1}\hat{G}(r_{i+1} \leftrightarrow r_i) \end{equation}$$

自发光亮度同样具有泛化形式, 空散射顶点不会发光.

$$\begin{equation} \begin{aligned} \hat{L_e}(p_n \to p_{n-1})= \begin{cases} L_e(p_n \to p_{n-1}) & p_i \in A\\ \sigma_a(p_n)L_e(p_n \to p_{n-1}) & p_i \in V \end{cases} \end{aligned} \end{equation}$$

$p_0$为相机, 所有长度为$n+1$的路径如下.

$$\begin{equation} \hat{P}(\bar{p}_n)=\int_{P_{n}}\hat{L_e}(p_n \to p_{n-1})\hat{T}(\bar{p}_n)d\bar{p}_{n} \end{equation}$$

此时Monte Carlo结果如下, 同时可以定义辐射通量权重$\beta(\bar{p}_n)$.

$$\begin{equation} \begin{aligned} \hat{P}(\bar{p}_n) &=\hat{L_e}(p_n \to p_{n-1})\beta(\bar{p}_n)\\ &=\frac{\hat{L_e}(p_n \to p_{n-1})\hat{T}(\bar{p}_n)}{p(\bar{p}_n)} \end{aligned} \end{equation}$$

求解路径空间积分

#

令$p_\text{maj}(p_{i+1}|p_i,\omega_i)$为从$p_i$在$\omega_i$方向上采样到$p_{i+1}$的概率, $p_e(p_i)$为在$p_i$选择吸收, 实散射或空散射的概率, 概率分别为$\frac{\sigma_{\{a,s,n\}}(p_i)}{\sigma_{\text{maj}}(p_i)}$, $p_\omega(\omega_{i+1}|r_i)$为实散射入射方向为$\omega_{i+1}$的概率. 场景中没有表面散射时Monte Carlo结果如下, 根据之前的结论可以看出部分项可以抵消, 同时由于只有$p_n$为光源, 当且仅当在$p_n$选择吸收.

$$\begin{equation} \begin{aligned} \hat{P}(\bar{p}_n) &=\frac{\hat{T}(\bar{p}_n)\hat{L}_e(p_n \to p_{n-1})}{\prod_{i=0}^{n-1}p_{\text{maj}}(p_{i+1}|p_i,\omega_i)\prod_{i=1}^{n}p_e(p_i)\prod_{i=1}^{m-1}p_\omega(\omega_{i+1}|r_i)\hat{G}(r_i \leftrightarrow r_{i+1})}\\ &=\frac{\hat{L}_e(p_n \to p_{n-1})\prod_{i=1}^{n-1}\hat{f}(p_{i+1} \to p_i \to p_{i-1})}{\prod_{i=0}^{n-1}\sigma_{\text{maj}}(p_{i+1})\prod_{i=1}^{n}p_e(p_i)\prod_{i=1}^{m-1}p_\omega(\omega_{i+1}|r_i)}\\ &=\frac{\hat{L}_e(p_n \to p_{n-1})\prod_{i=1}^{n-1}\hat{f}(p_{i+1} \to p_i \to p_{i-1})}{\sigma_{a}(p_n)\prod_{i=1}^{n-1}\sigma_{\{s,n\}}(p_i)\prod_{i=1}^{m-1}p_\omega(\omega_{i+1}|r_i)} \end{aligned} \end{equation}$$

若场景中没有表面散射, 此时$\hat{f}$和$\hat{L}_e$都是确定的, 可以进一步简化.

$$\begin{equation} \hat{P}(\bar{p}_n)=\prod_{i=1}^{m-1}\frac{p(r_{i+1} \to r_i \to r_{i-1})}{p_\omega(\omega_{i+1}|r_i)}L_e(p_n \to p_{n-1}) \end{equation}$$

体散射积分器

#

简单体积积分器

#

SimpleVolPathIntegrator不支持光源采样与表面散射, 即只用来渲染体积效果. 该类的成员只有maxDepth, 表示路径顶点数.

class SimpleVolPathIntegrator : public RayIntegrator {
  public:
    // SimpleVolPathIntegrator Public Methods
    // ...

  private:
    // SimpleVolPathIntegrator Private Members
    int maxDepth;
};

Li执行路径追踪, 与PathIntegrator的积分过程类似. 由于体数据可能随波长变化, Li只保留一个波长样本以简化计算.

Li首先计算光线求交, 因为场景中仍然会有一些几何物体, 例如用于表示介质分界的几何形状, 或者不进行散射的面积光源. 此时可以获取光线传播的最大距离, 体渲染样本不会超过这个范围.

Li通过SampleT_maj获取体渲染样本并执行差值跟踪, 在回调函数中首先选取介质事件, 根据吸收率, 散射率与空散射率进行概率选择.

// Compute medium event probabilities for interaction
Float pAbsorb = mp.sigma_a[0] / sigma_maj[0];
Float pScatter = mp.sigma_s[0] / sigma_maj[0];
Float pNull = std::max<Float>(0, 1 - pAbsorb - pScatter);

// Randomly sample medium scattering event for delta tracking
int mode = SampleDiscrete({pAbsorb, pScatter, pNull}, uMode);

若为吸收则停止采样, 获取光照并返回结果.

L += beta * mp.Le;
terminated = true;
return false;

若为散射则根据相位方程采样散射方向, 由于是实散射路径深度会增加. 发生散射返回false, 因为路径已经被改变, 不能再执行当前SampleT_maj. 注意到由于可以直接从Henyey-Greenstein中采样, 而pbrt目前只实现了这一种相位函数. 这里ps->p / ps->pdf始终为1.

// Stop path sampling if maximum depth has been reached
if (depth++ >= maxDepth) {
    terminated = true;
    return false;
}

// Sample phase function for medium scattering event
Point2f u{rng.Uniform<Float>(), rng.Uniform<Float>()};
pstd::optional<PhaseFunctionSample> ps =
    mp.phase.Sample_p(-ray.d, u);
if (!ps) {
    terminated = true;
    return false;
}

// Update state for recursive evaluation of $L_\roman{i}$
beta *= ps->p / ps->pdf;
ray.o = p;
ray.d = ps->wi;
scattered = true;
return false;

若为空散射则返回true以继续采样, 同时由于只有空散射会继续执行当前回调函数, uMode需要重新生成.

uMode = rng.Uniform<Float>();
return true;

路径追踪结束后, 若与表面相交则获取面积光源, 否则使用无限光源.

if (terminated)
    return L;
if (scattered)
    continue;
// Add emission to surviving ray
if (si)
    L += beta * si->intr.Le(-ray.d, lambda);
else {
    for (const auto &light : infiniteLights)
        L += beta * light.Le(ray, lambda);
    return L;
}

pbrt会检查是否有发生散射的几何物体, 在SimpleVolPathIntegrator中这是无效的.

BSDF bsdf = si->intr.GetBSDF(ray, lambda, camera, buf, sampler);
if (!bsdf)
    si->intr.SkipIntersection(&ray, si->tHit);
else {
    // Report error if BSDF returns a valid sample
    Float uc = sampler.Get1D();
    Point2f u = sampler.Get2D();
    if (bsdf.Sample_f(-ray.d, uc, u))
        ErrorExit("SimpleVolPathIntegrator doesn't support surface scattering.");
    else
        break;
}

改进采样技术

#

VolPathIntegrator支持表面散射, 多重波长介质光谱属性以及光源重要性抽样. 采样时若样本距离超过表面相交点, 则可以进行表面散射.

彩色介质

#

光谱变化介质属性的实现问题主要在于采样, 如果只采用单个波长做重要性抽样, 会降低其余波长的采样效果. 例如采样$\sigma_{\text{maj}}$时, 其余波长下可能会发生更多的空散射. 但如果为每个波长分别积分, 这在时间开销上又是低效的. 同时, 对于前文的路径追踪求解方法, 这会导致有些项无法抵消, 极有可能增大方差.

多个波长上不同的概率分布可以通过单抽样MIS解决, pbrt使用均匀分布来选择波长概率分布, 由于构造光谱渲染使用的波长时已经完成了概率选择, 这里直接选用第一个波长. 单抽样MIS下平衡启发式是较优的, pbrt直接通过各个分布密度的归一化实现. 单抽样MIS形式如下.

$$\begin{equation} \begin{aligned} f_\lambda(x) &=\frac{p_{\lambda_1}(x)}{\frac{1}{n}\sum_1^n p_{\lambda_i}(x)}\frac{\left[f_{\lambda_1}(x),f_{\lambda_2}(x),\dots,f_{\lambda_n}(x)\right]}{p_{\lambda_1}(x)}\\ &=\frac{\left[f_{\lambda_1}(x),f_{\lambda_2}(x),\dots,f_{\lambda_n}(x)\right]}{\frac{1}{n}\sum_1^n p_{\lambda_i}(x)} \end{aligned} \end{equation}$$

直接光照

#

引入体渲染后, 直接光源与表面间的透射率也需要被考虑进PDF, 因此pbrt对路径追踪与直接光照使用不同的采样技术. 路径追踪采用SimpleVolPathIntegrator使用的差值跟踪, pbrt称之为单向路径采样, PDF用$p_u$表示. 直接光照则采用比率跟踪获取透射率, 称之为光源路径采样, PDF用$p_l$表示.

单向路径采样概率如下, $p_s$为在$p_{n-1}$选择折射到并沿$\omega_n$传播的概率, $p_{\emptyset}$为实散射之间的空散射顶点对应的概率, 最后一项的$T_{\text{maj}}(p_m \to r_n)$代表着所有采样到超过$p_n$的点的累积概率.

$$\begin{equation} \begin{aligned} p_u(\bar{p}_n) &=p_u(\bar{p}_{n-1})p_s(r_{n-1})p_{\emptyset}(r_{n-1},r_n)\\ &=p_u(\bar{p}_{n-1})\frac{\sigma_s(p_{n-1})}{\sigma_{\text{maj}}(p_{n-1})}p_\omega(r_n|r_{n-1}, \omega_n)(\prod_{k=1}^m\frac{\sigma_n(p_k)}{\sigma_{\text{maj}}(p_k)}\sigma_{\text{maj}}(p_k)T_{\text{maj}}(p_{k-1} \to p_k))T_{\text{maj}}(p_m \to r_n)\\ &=p_u(\bar{p}_{n-1})\frac{\sigma_s(p_{n-1})}{\sigma_{\text{maj}}(p_{n-1})}p_\omega(r_n|r_{n-1}, \omega_n)(\prod_{k=1}^m\sigma_n(p_k)T_{\text{maj}}(p_{k-1} \to p_k))T_{\text{maj}}(p_m \to r_n) \end{aligned} \end{equation}$$

光源路径采样的概率与单向路径采样类似, 由于直接光照下散射点与光源之间不会发生折射, 这里不包含选择折射的概率.

$$\begin{equation} \begin{aligned} p_l(\bar{p}_n) &=p_u(\bar{p}_{n-1})p_s(r_{n-1})p_{\emptyset}(r_{n-1},r_n)\\ &=p_u(\bar{p}_{n-1})\frac{\sigma_s(p_{n-1})}{\sigma_{\text{maj}}(p_{n-1})}p_{l,\omega}(\omega_n)(\prod_{k=1}^m\sigma_{\text{maj}}(p_k)T_{\text{maj}}(p_{k-1} \to p_k))T_{\text{maj}}(p_m \to r_n) \end{aligned} \end{equation}$$

与光谱单抽样MIS结合后的结果如下.

$$\begin{equation} \begin{aligned} \hat{P}(\bar{p}_n) &=\omega_u(\bar{p}_n)\frac{\hat{T}(\bar{p}_n)L_e(p_n \to p_{n-1})}{p_{u,\lambda_s}(\bar{p}_n)}+\omega_l(\bar{p}’_n)\frac{\hat{T}(\bar{p}’_n)L_e(p’_n \to p’_{n-1})}{p_{l,\lambda_s}(\bar{p}’_n)} \end{aligned} \end{equation}$$

平衡权重如下.

$$\begin{equation} \begin{aligned} \omega_u(\bar{p}_n)=\frac{p_{u,\lambda_1}(\bar{p}_n)}{\frac{1}{m}\left(\sum_{i=1}^m p_{u,\lambda_i}(\bar{p}_n)+\sum_{i=1}^m p_{l,\lambda_i}(\bar{p}’_n)\right)}\\ \omega_l(\bar{p}_n)=\frac{p_{l,\lambda_1}(\bar{p}_n)}{\frac{1}{m}\left(\sum_{i=1}^m p_{u,\lambda_i}(\bar{p}_n)+\sum_{i=1}^m p_{l,\lambda_i}(\bar{p}’_n)\right)} \end{aligned} \end{equation}$$

改进的体积积分器

#

VolPathIntegrator的成员与PathIntegrator类似, 构造函数只做成员初始化.

class VolPathIntegrator : public RayIntegrator {
  public:
    // VolPathIntegrator Public Methods
    // ...

  private:
    // VolPathIntegrator Private Methods
    // ...

    // VolPathIntegrator Private Members
    int maxDepth;
    LightSampler lightSampler;
    bool regularize;
};

可以看出平衡权重的分子与采样项分母可以抵消, 但这可能会导致$\beta$溢出, 例如BSDF采样结果过大, 分母上的PDF可以有效防止溢出. 同样的平衡权重也需要考虑浮点数精度, 因此定义重缩放路径概率.

$$\begin{equation} \begin{aligned} r_{u,\lambda_i}(\bar{p}_n)=\frac{p_{u,\lambda_i}(\bar{p}_n)}{p_{\text{path}}(\bar{p}_n)}\\ r_{l,\lambda_i}(\bar{p}_n)=\frac{p_{l,\lambda_i}(\bar{p}_n)}{p_{\text{path}}(\bar{p}_n)}\\ \end{aligned} \end{equation}$$

$p_{\text{path}}$是采样当前路径的概率, 在单向路径采样中为$p_{u,\lambda_1}$, 在光源路径采样中为$p_{l,\lambda_1}$, 因此重缩放后的平衡权重如下. 注意这两项中的$p_{\text{path}}$是不同的, 这里只是防止数值误差的手段, 和之前的权重实际上是一致的.

$$\begin{equation} \begin{aligned} \omega_u(\bar{p}_n)=\frac{1}{\frac{1}{m}\left(\sum_{i=1}^m r_{u,\lambda_i}(\bar{p}_n)+\sum_{i=1}^m r_{l,\lambda_i}(\bar{p}’_n)\right)}\\ \omega_l(\bar{p}_n)=\frac{1}{\frac{1}{m}\left(\sum_{i=1}^m r_{u,\lambda_i}(\bar{p}_n)+\sum_{i=1}^m r_{l,\lambda_i}(\bar{p}’_n)\right)} \end{aligned} \end{equation}$$

VolPathIntegrator中所有采样点都会添加自发光, 因为是Monte Carlo, 所以这只是多统计了一条路径来提高统计效率, 不影响最终结果. 令新添加的路径顶点为$p’$, 光照结果如下. 由于MIS, $p_{\text{maj}}$中的$T_{\text{maj}}$不会再被抵消, 同时由于这条路径始终添加自发光$p_e$也始终为1.

$$\begin{equation} \begin{aligned} \hat{P}(\left[\bar{p}_n + p’\right]) &=\beta(\left[\bar{p}_n+p’\right])\sigma_a(p’)L_e(p’ \to p_n)\\ &=\frac{\beta(\bar{p}_n)T_{\text{maj}(p_n \to p’)}}{p_e(p’)p_{\text{maj}}(p’|p_n,\omega)}\sigma_a(p’)L_e(p’ \to p_n) \end{aligned} \end{equation}$$

由于波长的影响这里同样需要做MIS, MIS权重如下.

$$\begin{equation} \omega_e(\left[\bar{p}_n + p’\right])=\frac{1}{\frac{1}{m}\sum_{i=1}^m r_{e,\lambda_i}([\bar{p}_n+p’])} \end{equation}$$

代码实现如下.

if (depth < maxDepth && mp.Le) {
    // Compute $\beta'$ at new path vertex
    Float pdf = sigma_maj[0] * T_maj[0];
    SampledSpectrum betap = beta * T_maj / pdf;

    // Compute rescaled path probability for absorption at path vertex
    SampledSpectrum r_e = r_u * sigma_maj * T_maj / pdf;

    // Update _L_ for medium emission
    if (r_e)
        L += betap * mp.sigma_a * mp.Le / r_e.Average();
}

之后按照与SimpleVolPathIntegrator一样的方法选择散射事件. 由于自发光以及被统计过, 选择吸收时只停止追踪, 不做其它操作. 这里再次强调一下我的理解, Monte Carlo方法选择某条路径的概率是任意的, 只需要符合实际使用的概率就行, 使用不同的PDF只是方差上有区别, 因此每次统计都选择自发光时概率设置为1没有问题, 这里通过散射事件的概率来选择也没有问题.

terminated = true;
return false;

实散射事件中选择当前顶点概率为$\sigma_{\text{maj}}(p’)T_{\text{maj}}(p_n \to p’)$, 选择实散射概率为$\frac{\sigma_s(p’)}{\sigma_{\text{maj}}(p’)}$, 这里$\sigma_{\text{maj}}(p’)$被抵消. 这里$\beta$只更新了与选择顶点相关的部分, 散射部分还没添加, 因此不必在这里纠结.

Float pdf = T_maj[0] * mp.sigma_s[0];
beta *= T_maj * mp.sigma_s / pdf;
r_u *= T_maj * mp.sigma_s / pdf;

直接光照部分后续再介绍, Li开头是计算了距离最近的表面的距离的, 采样不会超过这个距离, 因此这里还没有接触到表面, 需要使用相位函数. 对于直接光照路径, 它的顶点除光源外与单向路径$\bar{p}_n$是相同的, 因此$r_l(\left[\bar{p}_n \to p’\right])$可以直接从$r_u(\bar{p}_n)$中推导出来, 这里除以ps->pdf是因为MIS使用的$p_{\text{path}}$为$p_{u,\lambda_1}$, 分子上的光源路径概率后续再更新. $r_u(\left[\bar{p}_n \to p\right]$不需要更新是因为相位方程与波长无关, 分子分母抵消了.

// Sample direct lighting at volume-scattering event
MediumInteraction intr(p, -ray.d, ray.time, ray.medium,
                       mp.phase);
L += SampleLd(intr, nullptr, lambda, sampler, beta, r_u);

// Sample new direction at real-scattering event
Point2f u = sampler.Get2D();
pstd::optional<PhaseFunctionSample> ps =
    intr.phase.Sample_p(-ray.d, u);
if (!ps || ps->pdf == 0)
    terminated = true;
else {
    // Update ray path state for indirect volume scattering
    beta *= ps->p / ps->pdf;
    r_l = r_u / ps->pdf;
    prevIntrContext = LightSampleContext(intr);
    scattered = true;
    ray.o = p;
    ray.d = ps->wi;
    specularBounce = false;
    anyNonSpecularBounces = true;
}

空散射事件中概率中的$\sigma_{\text{maj}}(p’)$同样被抵消. 光源路径采样认为当前采样是计算与光源之间的透射率, 因此在路径传播过程中不需要选择散射事件, 只需要记录选择顶点的概率, 因此不包含$\frac{\sigma_n(p’)}{\sigma_{\text{maj}}(p’)}$, 分母上的pdf同样是因为MIS使用的$p_{\text{path}}$为$p_{u,\lambda_1}$.

SampledSpectrum sigma_n = ClampZero(sigma_maj - mp.sigma_a - mp.sigma_s);
Float pdf = T_maj[0] * sigma_n[0];
beta *= T_maj * sigma_n / pdf;
if (pdf == 0) beta = SampledSpectrum(0.f);
r_u *= T_maj * sigma_n / pdf;
r_l *= T_maj * sigma_maj / pdf;
return beta && r_u;

退出SampleT_maj后, 如路径大于最大深度, 遇到吸收事件, 散射PDF为0或任意一个波长的$\beta$或$r_u$为0, 则返回光照结果. 若遇到散射事件则继续循环.

if (terminated || !beta || !r_u)
    return L;
if (scattered)
    continue;

剩下的情况是没有发生任何散射事件或只经历了空散射, 这部分的实现在SampleT_maj中, 因为介质迭代器迭代到末尾而返回当前累计的$T_\text{maj}$. 没有散射事件则返回值在各个波长上都为$1$, 不影响这些系数; 若只经历空散射则代表光线沿直线从介质起点传播到终点, 最后一个采样点超过了介质最大距离, 根据之前章节的结论, 这种情况概率可以从CDF中获取, 值为$T_{\text{maj}}$. $\beta$分子乘上$T_{\text{maj}}$是因为路径通量需要统计透射率, $r_u$和$r_l$则是各个波长上超过最大距离的概率.

beta *= T_maj / T_maj[0];
r_u *= T_maj / T_maj[0];
r_l *= T_maj / T_maj[0];

处理与物体表面相交的情况, 此时若已经达到最大深度, 则统计表面自发光后就结束路径追踪.

// Terminate path if maximum depth reached
if (depth++ >= maxDepth)
    return L;

获取直接光照, 同样使用SampleLd. 与PathIntegrator一样, 由于镜面反射概率分布为Dirac delta分布, 无需统计直接光照.

// Sample illumination from lights to find attenuated path contribution
if (IsNonSpecular(bsdf.Flags())) {
    L += SampleLd(isect, &bsdf, lambda, sampler, beta, r_u);
    DCHECK(IsInf(L.y(lambda)) == false);
}
prevIntrContext = LightSampleContext(isect);

根据BSDF生成追踪路径.

// Sample BSDF to get new volumetric path direction
Vector3f wo = isect.wo;
Float u = sampler.Get1D();
pstd::optional<BSDFSample> bs = bsdf.Sample_f(wo, u, sampler.Get2D());
if (!bs)
    break;

更新$\beta$与$r_l$, $r_u$因为BSDF与波长无关, 并不需要更新. $\beta$更新时添加的余弦项可参照泛化路径空间的定义. $r_l$与之前一样只更新分母上的PDF, pdfIsProportional与某些特殊材质相关, 具体见后续章节.

// Update _beta_ and rescaled path probabilities for BSDF scattering
beta *= bs->f * AbsDot(bs->wi, isect.shading.n) / bs->pdf;
if (bs->pdfIsProportional)
    r_l = r_u / bsdf.PDF(wo, bs->wi);
else
    r_l = r_u / bs->pdf;

俄罗斯轮盘部分与PathIntegrator类似, 在其基础上添加了$r_u$相关的系数, $r_u$在各个波长上的均值越小退出概率越小. 猜测这里是因为较小的$r_u$代表当前使用的波长越重要.

SampledSpectrum rrBeta = beta * etaScale / r_u.Average();
Float uRR = sampler.Get1D();
PBRT_DBG("%s\n",
         StringPrintf("etaScale %f -> rrBeta %s", etaScale, rrBeta).c_str());
if (rrBeta.MaxComponentValue() < 1 && depth > 1) {
    Float q = std::max<Float>(0, 1 - rrBeta.MaxComponentValue());
    if (uRR < q)
        break;
    beta /= 1 - q;
}

估计直接光照

#

SampleLd实现直接光照, 与PathIntegrator相比这里需要计算透射率.

散射点可能是表面或体积, 若为表面则和PathIntegrator一样对交点做偏移, 使得反射时位于面外, 折射时位于面内.

LightSampleContext ctx;
if (bsdf) {
    ctx = LightSampleContext(intr.AsSurface());
    // Try to nudge the light sampling position to correct side of the surface
    BxDFFlags flags = bsdf->Flags();
    if (IsReflective(flags) && !IsTransmissive(flags))
        ctx.pi = intr.OffsetRayOrigin(intr.wo);
    else if (IsTransmissive(flags) && !IsReflective(flags))
        ctx.pi = intr.OffsetRayOrigin(-intr.wo);

} else
    ctx = LightSampleContext(intr);

生成随机变量并采样光源, 与PathIntegrator一致.

Float u = sampler.Get1D();
pstd::optional<SampledLight> sampledLight = lightSampler.Sample(ctx, u);
Point2f uLight = sampler.Get2D();
if (!sampledLight)
    return SampledSpectrum(0.f);
Light light = sampledLight->light;
DCHECK(light && sampledLight->p != 0);

// Sample a point on the light source
pstd::optional<LightLiSample> ls = light.SampleLi(ctx, uLight, lambda, true);
if (!ls || !ls->L || ls->pdf == 0)
    return SampledSpectrum(0.f);
Float p_l = sampledLight->p * ls->pdf;

分别处理BSDF与相位方程. BSDF的计算方式与PathIntegrator一致, 注意这里$\hat{f}$把$D(p,p’)$即余弦项包括进去了. 相位方程的$\hat{f}$没有包含$\sigma_s$, 因为已经包括在$\beta$中了.

// Evaluate BSDF or phase function for light sample direction
Float scatterPDF;
SampledSpectrum f_hat;
Vector3f wo = intr.wo, wi = ls->wi;
if (bsdf) {
    // Update _f_hat_ and _scatterPDF_ accounting for the BSDF
    f_hat = bsdf->f(wo, wi) * AbsDot(wi, intr.AsSurface().shading.n);
    scatterPDF = bsdf->PDF(wo, wi);

} else {
    // Update _f_hat_ and _scatterPDF_ accounting for the phase function
    CHECK(intr.IsMediumInteraction());
    PhaseFunction phase = intr.AsMedium().phase;
    f_hat = SampledSpectrum(phase.p(wo, wi));
    scatterPDF = phase.PDF(wo, wi);
}
if (!f_hat)
    return SampledSpectrum(0.f);