PBRT笔记(14)——光线传播2：体积渲染

传输公式

传输方程是控制光线在吸收、发射和散射辐射的介质中的行为的基本方程。它解释了第11章中描述的所有体积散射过程——吸收、发射和内、外散射。并给出了一个描述环境中辐射分布的方程。光传输方程实际上是传输方程的一个特例，由于忽略参与介质而进行简化，并专门用于从表面散射。

在方程的基本形式中，传输方程式是一个描述光束在空间中某一点上的辐射量如何变化的积分微分方程。它可以转化为一个纯积分方程，用以描述描述了沿射线无穷多个坐标点的参与介质的效果。

回忆一下11.1.4章节的辐射源项Ls，它表示p点处ω方向上由于从其他介质中自发光和内散射产生的光线而导致的辐射量的变化。

辐射源项包含了所有光线辐射量增加的过程。其中衰减系数σt(p, ω)包含了所有光线辐射量减少的过程，包括了吸收与外散射。描述其效果的微分方程为：

当计算光线上一点p’的辐射量变化微分时，将以上两种不同的效应结合在一起，从而得到传输方程的积分微分形式：

当采用适当的边界条件时，这个方程可以转换为纯积分方程。例如：如果我们假定场景中没有反射与自发光表面，那么光线就不会被遮挡并且光的长度是无限，那么传输积分方程为：

其中p’=p+tω的意思为：从既定方向到达某一点的辐射量，即该点与光线上所有点的辐射量累加结果。其辐射量将根据光束总的透射率减少。

通常来说，如果场景中存在反射或者自发光表面，光线的长度就不会无限长，并且光线与该表面相交后会影响自身的辐射量，增加表面处该点的出射辐射量以及阻止
光线中超出交点的点对源点的辐射量产生影响。如果光线在表面的p0处相交，其距离为t，则传输积分公式为：

其中p0=p+tω为表面上的点，p’=p+t’ω为光线上的点。

这个方程描述了沿着光线发生的两种对辐射量产生影响的现象。首先，Lo项代表了从表面反射回来的辐射量，它包括了表面自发光与反射回来的辐射量。这种辐射量可能会因为参与介质而减弱，光线源点到p0点之间的透射率解释了这点。第二项解释了由于体积散射和自发光使得沿射线辐射量增加，但只到射线与表面交点为止;超过这一点的点不再影响光线的辐射量。

采样空间散射

在讨论参与介质中光散射效应模型算法之前，我们首先定义一些用于从参与介质与空间位置无关介质的光速透射率相关分布中采样的功能块。

Medium接口定了Sample(),它接受一个世界空间的ray(p,ω)以及MediumInteraction对象指针。输入ray通常将会与场景多边形进行相交测试，因此该方法不应该在光线上超出tmax值的点上采样介质。在不失一般性的前提下，下面的讨论是在总是有一个tmax小于正无穷的曲面的假设下。

这种方法的目的是对传递方程的积分形式进行采样：

其中p0= p + tmax ω是表面上的点。我们讲忽略介质的自发光并假设介质属性为定向常数，此时：

两种情况可能会发生：如果Sample()采样的交点没有在给定光线的区间[0,tMax]中，那么表面相关项Tr与Lo就为估算值。如果它采样了一个点，那么第二个积分项将会是估算项，并且提供的MediumInteraction对象应该相应的初始化。

假设pt(t)表示每单位距离的概率p+tω产生交互作用的位置。由于不采样介质相互作用的可能性，这个函数一般不积分为1，我们将psurf定义为对表面项采样的相关离散概率:

有了这些定义，我们可以指定Sample()的语义，它不同于之前遇到的散射函数所用的技术，比如BSDF::Sample_f()，因为它不向调用者提供关于函数值和采样位置处的PDF信息。一般来说并不需要这种信息，一些介质模型(特别是不均匀的介质)在计算psurf的比值时，会使用更有效的采样方案。

当选择表面项时，方法应该返回一个权重值：

它符合采样第一个和。注意，出射辐射量的值Lo不包括在βsurf内。调用者有责任去解释这个项。在介质中时，方法返回：

它对应了除了内散射光积分（公式15.6）外所有介质相关项都会被采样，必须单独处理。

散射系数与透射率允许光谱产生变化，因此这个方法返回一个光谱值权重因子，用以来更新路径的表面出射权重β或介质散射事件。

均匀介质

HomogeneousMedium类的这个方法实现得较为简单。唯一的复杂性源自需要处理随波长变化的衰减系数。

在13.3.1章节，我们推导了定义域为[0,∞)的指数分布的采样方法。对于$f(t)=e^{-\sigma_tt}$,可以有：

其中pdf为：

然而衰减系数σt通常随波长变化而变化。它在介质中多点采样是不可取的，因此首先采用均匀采样来选择光谱通道i;相应的标量σi t值用于样本距离的分布。

采用15.9中的方法。所得到的采样密度为各个$p^i_t$的均值：

在t = tmax处对表面相互作用进行采样的概率（离散）是在t=0和t=tmax之间产生中间散射事件的补充。计算出的概率等于所有n个光谱通道的平均透过率:

根据15.11的实现提取出的样本；如果采样距离在射线与图元相交之前(如果有)，则通过初始化MediumInteraction来记录介质散射事件。否则，忽略介质中的采样点，积分器应使用相应的表面相互作用作为下一个路径顶点。这种采样方法符合自然规律:产生介质相互作用而不是表面相互作用的概率正好等于1减去选定波长的光束透过率。因此，给定光学薄介质(或短射线范围)，更经常使用表面相互作用，而对于厚介质(或长射线)，更可能取样介质相互作用。

int channel = std::min((int)(sampler.Get1D() * Spectrum::nSamples),
Spectrum::nSamples - 1);
Float dist = -std::log(1 - sampler.Get1D()) / sigma_t[channel];
Float t = std::min(dist * ray.d.Length(), ray.tMax);
bool sampledMedium = t < ray.tMax;
if (sampledMedium)
*mi = MediumInteraction(ray(t), -ray.d, ray.time, this,
ARENA_ALLOC(arena, HenyeyGreenstein)(g));

在这两种情况下，利用比尔定律(11.3)可以很容易地计算光束透射率Tr，就像在HomogeneousMedium::Tr()。

Spectrum Tr = Exp(-sigma_t * std::min(t, MaxFloat) * ray.d.Length());

最后,该方法使用方程(15.11)或(15.12)计算样本密度并返回采样权重βsurf与βmed。

Spectrum density = sampledMedium ? (sigma_t * Tr) : Tr;
Float pdf = 0;
for (int i = 0; i < Spectrum::nSamples; ++i)
    pdf += density[i];
pdf *= 1 / (Float)Spectrum::nSamples;
return sampledMedium ? (Tr * sigma_s / pdf) : (Tr / pdf);

不均匀的介质

对于GridDensityMedium，需要额外的工作来处理介质的异质性。当空间变化可以分解为均匀区域(如分段常数体素)时，常规追踪技术将标准均匀介质技术分别应用于各个体素;这种方法的一个缺点是，当体素量较多时，它会消耗过多性能。由于网格密度介质依赖于线性插值，因此不能采用这种方法。

其他技术是建立在直接从具有空间变化的衰减公式(15.10)中均匀采样pdf。

其中σt(t)=σt(p+tω)估算射线上距离为t的衰减值。最常用的方法是重要抽样方程(15.13)，称为ray marching。该方法通过将范围[0,tmax]划分为若干子区间，对每个区间内的积分进行数值近似计算，最后对该离散表示形式进行反演，从而反演近似的累积分布。然而这种方法会产生系统性的统计偏差，这意味着使用ray marching的积分器通常不会收敛到正确的结果(即使每个像素使用无限多个样本)。

因此，我们更喜欢Woodcock等人(1965)提出的另一种无偏方法，该方法最初是为了模拟原子反应堆中中子的体积散射而开发的。这种技术被称为delta追踪，并且在衰减系数σt是单色的情况下更容易实现。PBRT的实现包括一个断言测试(这里没有显示)，以验证情况确实如此。
注意,散射和吸收系数仍然允许波长不同——不过,它们的和σt=σs+σa必须统一。

delta追踪可以解释为在介质中填充额外粒子（虚拟），直到它的衰减系数处处为常数。利用基础指数公式(15.9)可以很容易地对均匀介质进行采样。然而，无论何时发生与粒子的交互，仍然有必要确定它所涉及的是“真实”粒子还是“虚拟”粒子(在这种情况下交互被忽略)。woodcock等人的见解是：根据基于“真实粒子”的局部分数随机做出判断，因此这也引出了样本匹配方程(15.13)分布。

下面是GridDensityMedium构造函数的一部分，它会预计算整个介质最大密度因子的倒数，在接下来的delta追踪中这将十分有用。

<Precompute values for Monte Carlo sampling of GridDensityMedium> ≡ 690
sigma_t = (sigma_a + sigma_s)[0];
Float maxDensity = 0;
for (int i = 0; i < nx * ny * nz; ++i)
maxDensity = std::max(maxDensity, density[i]);
invMaxDensity = 1 / maxDensity;
<GridDensityMedium Private Data> +≡ 690
Float sigma_t;
Float invMaxDensity;

Sample()方法首先将射线转换为介质坐标系，并对射线方向进行归一化;ray.tMax因为归一化的关系会进行适当缩放。

Spectrum GridDensityMedium::Sample(const Ray &rWorld, Sampler &sampler,MemoryArena &arena,MediumInteraction *mi) const {
    ProfilePhase _(Prof::MediumSample);
    Ray ray = WorldToMedium(
        Ray(rWorld.o, Normalize(rWorld.d), rWorld.tMax * rWorld.d.Length()));
    // Compute $[\tmin, \tmax]$ interval of _ray_'s overlap with medium bounds
    const Bounds3f b(Point3f(0, 0, 0), Point3f(1, 1, 1));
    Float tMin, tMax;
    if (!b.IntersectP(ray, &tMin, &tMax)) return Spectrum(1.f);

    // Run delta-tracking iterations to sample a medium interaction
    Float t = tMin;
    while (true) {
        t -= std::log(1 - sampler.Get1D()) * invMaxDensity / sigma_t;
        if (t >= tMax) break;
        if (Density(ray(t)) * invMaxDensity > sampler.Get1D()) {
            // Populate _mi_ with medium interaction information and return
            PhaseFunction *phase = ARENA_ALLOC(arena, HenyeyGreenstein)(g);
            *mi = MediumInteraction(rWorld(t), -rWorld.d, rWorld.time, this,
                                    phase);
            return sigma_s / sigma_t;
        }
    }
    return Spectrum(1.f);
}

接下来实现计算光线与介质边界重叠的参数范围，即单位立方体[0,1]^3。这一步在技术上并不需要较为正确的操作，但通常是一个好主意:减少考虑的射线段的长度可以转化为相对较少的增量跟踪迭代。(适当的减少追踪光线的步数，以节约性能)

假设整个介质中的最大消光值为σt，max。每个delta追踪的迭代次数为i，其在均匀介质中每步执行的标准指数方程为：

其中t0=tmin。上述步骤会不断重复执行，直至符合以下两项停止准则之一:

1. 如果ti>tmax，那么光线就离开了介质，不会产生交互作用，并且Medium::Sample()不会再采样散射事件。另外每次追踪循环终止的概率为σt(ti)/σt,max,也就是真实粒子的局部分数。 ξ2i+1为随机消耗量，每次迭代两个均匀样本。
   <Run delta-tracking iterations to sample a medium interaction>

2. 不采样介质相互作用的概率等于射线[tmin, tmax]区间的透射率;因此在方程(15.7)中为采样权重βsurf返回1.0。
   <Populate mi with medium interaction information and return>

最终我们还必须实现tr()来计算沿射线方向其中一段的透射率。

由于通过介质的概率等于透射率，因此该随机变量具有正确的均值，可用于无偏蒙特卡罗积分。多次调用Tr()并对结果求平均值将会产生越来越精确的透射率估算值，尽管在实践中这样做通常成本太高。

<GridDensityMedium Method Definitions> +≡
Spectrum GridDensityMedium::Tr(const Ray &rWorld,
Sampler &sampler) const {
Ray ray = WorldToMedium(Ray(rWorld.o, Normalize(rWorld.d),
rWorld.tMax * rWorld.d.Length()));
<Compute [tmin, tmax] interval of ray’s overlap with medium bounds 896>
<Perform ratio tracking to estimate the transmittance value 898>
}

<Perform ratio tracking to estimate the transmittance value> ≡ 898
Float Tr = 1, t = tMin;
while (true) {
t -= std::log(1 - sampler.Get1D()) * invMaxDensity / sigma_t;
if (t >= tMax)
break;
Float density = Density(ray(t));
Tr *= 1 - std::max((Float)0, density * invMaxDensity);
}
return Spectrum(Tr);

采样相位函数

能够从相位函数描述的分布中提取样本也很有用——应用包括使用多重重要性采样来计算参与介质中的直接照明，以及对参与介质中的间接照明样本的散射方向进行采样。对于这些应用，PhaseFunction类必须实现Sample_p()。

注意，与BxDF的采样方法不同，Sample_p()不会同时返回相位函数的值及其PDF。相反，pbrt假设相位函数是用完全匹配其分布的PDFs采样的。结合相位函数必须归一化的要求(式(11.4))，一个返回值对两个值进行编码。当只需要PDF的值时，调用PhaseFunction::p()就足够了。

virtual Float Sample_p(const Vector3f &wo, Vector3f *wi,const Point2f &u) const = 0;

Henyey-Greenstein相位函数的pdf为分离的θ与φ组件，其中p(φ)=1/(2π)。主要的任务为采样cosθ。

对于Henyey–Greenstein分布的θ为：

如果g不为0，否则$cos\theta=1-2\xi$返回一个球体坐标方向的均匀采样。

Float cosTheta;
if (std::abs(g) < 1e-3)
    cosTheta = 1 - 2 * u[0];
else {
    Float sqrTerm = (1 - g * g) /(1 - g + 2 * g * u[0]);
    cosTheta = (1 + g * g - sqrTerm * sqrTerm) / (2 * g);
}

将角度坐标(cosθ,φ)转换为ωi方向。

Float sinTheta = std::sqrt(std::max((Float)0,
1 - cosTheta * cosTheta));
Float phi = 2 * Pi * u[1];
Vector3f v1, v2;
CoordinateSystem(wo, &v1, &v2);
*wi = SphericalDirection(sinTheta, cosTheta, phi, v1, v2, -wo);

体积光透射

这些采样构建块让各种各样的在参与介质中光线投射算法的实现成为可能。我们现在可以实现14.3.1节中EstimateDirect()中的内容，该函数处理与参与媒体相关的案例。

首先，在光线被采样后，如果在参与介质中的产生了散射事件，那么久必须为出射方向与入射光方向计算相位函数的值，以及用于多重重要性采样的采样方向的pdf值。因为我们假设相位函数时完美采样的，所以这些值都是相同的。

const MediumInteraction &mi = (const MediumInteraction &)it;
Float p = mi.phase->p(mi.wo, wi);
f = Spectrum(p);
scatteringPdf = p;

方向光计算需要取得一个相位函数分布样本，Sample_p()提供了这种功能；如上所述，它返回的值同时给出了相位函数的值和pdf的值。

const MediumInteraction &mi = (const MediumInteraction &)it;
Float p = mi.phase->Sample_p(mi.wo, &wi, uScattering);
f = Spectrum(p);
scatteringPdf = p;

路径追踪

VolPathIntegrator是一种SamplerIntegrator，它可以用来描述参与介质、次表面中散射与衰减。实现在integrators/volpath.h和integrators/volpath.cpp中。它跟PathIntegrator的一般结构非常相似。因此在这里我们只讨论这两个类之间的区别。

作为SamplerIntegrator，VolPathIntegrator主要实现了Li(),它的实现与PathIntegrator::Li()非常相似，只是在参与介质的相关部分有些改动。

Spectrum VolPathIntegrator::Li(const RayDifferential &r, const Scene &scene,Sampler &sampler,MemoryArena &arena,int depth) const {
    ProfilePhase p(Prof::SamplerIntegratorLi);
    Spectrum L(0.f), beta(1.f);
    RayDifferential ray(r);
    bool specularBounce = false;
    int bounces;
    // Added after book publication: etaScale tracks the accumulated effect
    // of radiance scaling due to rays passing through refractive
    // boundaries (see the derivation on p. 527 of the third edition). We
    // track this value in order to remove it from beta when we apply
    // Russian roulette; this is worthwhile, since it lets us sometimes
    // avoid terminating refracted rays that are about to be refracted back
    // out of a medium and thus have their beta value increased.
    Float etaScale = 1;

    for (bounces = 0;; ++bounces) {
        // Intersect _ray_ with scene and store intersection in _isect_
        SurfaceInteraction isect;
        bool foundIntersection = scene.Intersect(ray, &isect);

        // Sample the participating medium, if present
        MediumInteraction mi;
        if (ray.medium) beta *= ray.medium->Sample(ray, sampler, arena, &mi);
        if (beta.IsBlack()) break;

        // Handle an interaction with a medium or a surface
        if (mi.IsValid()) {
            // Terminate path if ray escaped or _maxDepth_ was reached
            if (bounces >= maxDepth) break;

            ++volumeInteractions;
            // Handle scattering at point in medium for volumetric path tracer
            const Distribution1D *lightDistrib =
                lightDistribution->Lookup(mi.p);
            L += beta * UniformSampleOneLight(mi, scene, arena, sampler, true,
                                              lightDistrib);

            Vector3f wo = -ray.d, wi;
            mi.phase->Sample_p(wo, &wi, sampler.Get2D());
            ray = mi.SpawnRay(wi);
            specularBounce = false;
        } else {
            ++surfaceInteractions;
            // Handle scattering at point on surface for volumetric path tracer

            // Possibly add emitted light at intersection
            if (bounces == 0 || specularBounce) {
                // Add emitted light at path vertex or from the environment
                if (foundIntersection)
                    L += beta * isect.Le(-ray.d);
                else
                    for (const auto &light : scene.infiniteLights)
                        L += beta * light->Le(ray);
            }

            // Terminate path if ray escaped or _maxDepth_ was reached
            if (!foundIntersection || bounces >= maxDepth) break;

            // Compute scattering functions and skip over medium boundaries
            isect.ComputeScatteringFunctions(ray, arena, true);
            if (!isect.bsdf) {
                ray = isect.SpawnRay(ray.d);
                bounces--;
                continue;
            }

            // Sample illumination from lights to find attenuated path
            // contribution
            const Distribution1D *lightDistrib =
                lightDistribution->Lookup(isect.p);
            L += beta * UniformSampleOneLight(isect, scene, arena, sampler,
                                              true, lightDistrib);

            // Sample BSDF to get new path direction
            Vector3f wo = -ray.d, wi;
            Float pdf;
            BxDFType flags;
            Spectrum f = isect.bsdf->Sample_f(wo, &wi, sampler.Get2D(), &pdf,
                                              BSDF_ALL, &flags);
            if (f.IsBlack() || pdf == 0.f) break;
            beta *= f * AbsDot(wi, isect.shading.n) / pdf;
            DCHECK(std::isinf(beta.y()) == false);
            specularBounce = (flags & BSDF_SPECULAR) != 0;
            if ((flags & BSDF_SPECULAR) && (flags & BSDF_TRANSMISSION)) {
                Float eta = isect.bsdf->eta;
                // Update the term that tracks radiance scaling for refraction
                // depending on whether the ray is entering or leaving the
                // medium.
                etaScale *=
                    (Dot(wo, isect.n) > 0) ? (eta * eta) : 1 / (eta * eta);
            }
            ray = isect.SpawnRay(wi);

            // Account for attenuated subsurface scattering, if applicable
            if (isect.bssrdf && (flags & BSDF_TRANSMISSION)) {
                // Importance sample the BSSRDF
                SurfaceInteraction pi;
                Spectrum S = isect.bssrdf->Sample_S(
                    scene, sampler.Get1D(), sampler.Get2D(), arena, &pi, &pdf);
                DCHECK(std::isinf(beta.y()) == false);
                if (S.IsBlack() || pdf == 0) break;
                beta *= S / pdf;

                // Account for the attenuated direct subsurface scattering
                // component
                L += beta *
                     UniformSampleOneLight(pi, scene, arena, sampler, true,
                                           lightDistribution->Lookup(pi.p));

                // Account for the indirect subsurface scattering component
                Spectrum f = pi.bsdf->Sample_f(pi.wo, &wi, sampler.Get2D(),
                                               &pdf, BSDF_ALL, &flags);
                if (f.IsBlack() || pdf == 0) break;
                beta *= f * AbsDot(wi, pi.shading.n) / pdf;
                DCHECK(std::isinf(beta.y()) == false);
                specularBounce = (flags & BSDF_SPECULAR) != 0;
                ray = pi.SpawnRay(wi);
            }
        }

        // Possibly terminate the path with Russian roulette
        // Factor out radiance scaling due to refraction in rrBeta.
        Spectrum rrBeta = beta * etaScale;
        if (rrBeta.MaxComponentValue() < rrThreshold && bounces > 3) {
            Float q = std::max((Float).05, 1 - rrBeta.MaxComponentValue());
            if (sampler.Get1D() < q) break;
            beta /= 1 - q;
            DCHECK(std::isinf(beta.y()) == false);
        }
    }
    ReportValue(pathLength, bounces);
    return L;
}

在每一次散射路径采样中，光线先与场景中的表面进行相交测试，来寻找最近的表面交点(如果有的话)。接下来通过调用Medium::Sample()方法来解析参与介质信息，如果一个在下一个路径上的顶点为交互介质，那么则初始化提供的MediumInteraction对象。在任一情况下，Sample()会返回一个光束透射率与表面、交互介质的采样pdf因子。
<Sample the participating medium, if present>

在散射介质密度非常大的场景中，由于Medium::Sample()通常会生成交互介质，所以在第一次寻找表面交点时回浪费性能。对于这种场景，一种更有效的实现方法是首先采样一个交互介质，在光线与场景图元相交之前更新对应光线的tMax值。反过来锁，由于要测试的光线通常相当短，所以表面相交测试将会更加高效。

取决于是否为该射线在参与介质或表面上的某个点进行交互采样，其中之一控制计算该点直接光照以及下一个采样方向。

由于之前章节定义的代码，UniformSampleOneLight()已经支持估算参与介质中某点的直接光照，所以我们只需要将MediumInteraction对象传递给这个函数就可以了，然后调用Sample_p()就可以得到光线离开介质的方向。

Spectrum SeparableBSSRDF::Sample_S(const Scene &scene, Float u1,const Point2f &u2, MemoryArena &arena,SurfaceInteraction *si, Float *pdf) const {
    ProfilePhase pp(Prof::BSSRDFSampling);
    Spectrum Sp = Sample_Sp(scene, u1, u2, arena, si, pdf);
    if (!Sp.IsBlack()) {
        // Initialize material model at sampled surface interaction
        si->bsdf = ARENA_ALLOC(arena, BSDF)(*si);
        si->bsdf->Add(ARENA_ALLOC(arena, SeparableBSSRDFAdapter)(this));
        si->wo = Vector3f(si->shading.n);
    }
    return Sp;
}

Spectrum f(const Vector3f &wo, const Vector3f &wi) const {
    Spectrum f = bssrdf->Sw(wi);
    // Update BSSRDF transmission term to account for adjoint light
    // transport
    if (bssrdf->mode == TransportMode::Radiance)
        f *= bssrdf->eta * bssrdf->eta;
    return f;
}

Spectrum SeparableBSSRDF::Sample_Sp(const Scene &scene, Float u1,const Point2f &u2, MemoryArena &arena,SurfaceInteraction *pi, Float *pdf) const {
    ProfilePhase pp(Prof::BSSRDFEvaluation);
    // Choose projection axis for BSSRDF sampling
    Vector3f vx, vy, vz;
    if (u1 < .5f) {
        vx = ss;
        vy = ts;
        vz = Vector3f(ns);
        u1 *= 2;
    } else if (u1 < .75f) {
        // Prepare for sampling rays with respect to _ss_
        vx = ts;
        vy = Vector3f(ns);
        vz = ss;
        u1 = (u1 - .5f) * 4;
    } else {
        // Prepare for sampling rays with respect to _ts_
        vx = Vector3f(ns);
        vy = ss;
        vz = ts;
        u1 = (u1 - .75f) * 4;
    }

    // Choose spectral channel for BSSRDF sampling
    int ch = Clamp((int)(u1 * Spectrum::nSamples), 0, Spectrum::nSamples - 1);
    u1 = u1 * Spectrum::nSamples - ch;

    // Sample BSSRDF profile in polar coordinates
    Float r = Sample_Sr(ch, u2[0]);
    if (r < 0) return Spectrum(0.f);
    Float phi = 2 * Pi * u2[1];

    // Compute BSSRDF profile bounds and intersection height
    Float rMax = Sample_Sr(ch, 0.999f);
    if (r >= rMax) return Spectrum(0.f);
    Float l = 2 * std::sqrt(rMax * rMax - r * r);

    // Compute BSSRDF sampling ray segment
    Interaction base;
    base.p =
        po.p + r * (vx * std::cos(phi) + vy * std::sin(phi)) - l * vz * 0.5f;
    base.time = po.time;
    Point3f pTarget = base.p + l * vz;

    // Intersect BSSRDF sampling ray against the scene geometry

    // Declare _IntersectionChain_ and linked list
    struct IntersectionChain {
        SurfaceInteraction si;
        IntersectionChain *next = nullptr;
    };
    IntersectionChain *chain = ARENA_ALLOC(arena, IntersectionChain)();

    // Accumulate chain of intersections along ray
    IntersectionChain *ptr = chain;
    int nFound = 0;
    while (true) {
        Ray r = base.SpawnRayTo(pTarget);
        if (r.d == Vector3f(0, 0, 0) || !scene.Intersect(r, &ptr->si))
            break;

        base = ptr->si;
        // Append admissible intersection to _IntersectionChain_
        if (ptr->si.primitive->GetMaterial() == this->material) {
            IntersectionChain *next = ARENA_ALLOC(arena, IntersectionChain)();
            ptr->next = next;
            ptr = next;
            nFound++;
        }
    }

    // Randomly choose one of several intersections during BSSRDF sampling
    if (nFound == 0) return Spectrum(0.0f);
    int selected = Clamp((int)(u1 * nFound), 0, nFound - 1);
    while (selected-- > 0) chain = chain->next;
    *pi = chain->si;

    // Compute sample PDF and return the spatial BSSRDF term $\Sp$
    *pdf = this->Pdf_Sp(*pi) / nFound;
    return this->Sp(*pi);
}

Float TabulatedBSSRDF::Sample_Sr(int ch, Float u) const {
if (sigma_t[ch] == 0)
    return -1;
return SampleCatmullRom2D(table.nRhoSamples, table.nRadiusSamples,
table.rhoSamples.get(), table.radiusSamples.get(),
table.profile.get(), table.profileCDF.get(),
rho[ch], u) / sigma_t[ch];
}

Float TabulatedBSSRDF::Pdf_Sr(int ch, Float r) const {
    // Convert $r$ into unitless optical radius $r_{\roman{optical}}$
    Float rOptical = r * sigma_t[ch];

    // Compute spline weights to interpolate BSSRDF density on channel _ch_
    int rhoOffset, radiusOffset;
    Float rhoWeights[4], radiusWeights[4];
    if (!CatmullRomWeights(table.nRhoSamples, table.rhoSamples.get(), rho[ch],
                           &rhoOffset, rhoWeights) ||
        !CatmullRomWeights(table.nRadiusSamples, table.radiusSamples.get(),
                           rOptical, &radiusOffset, radiusWeights))
        return 0.f;

    // Return BSSRDF profile density for channel _ch_
    Float sr = 0, rhoEff = 0;
    for (int i = 0; i < 4; ++i) {
        if (rhoWeights[i] == 0) continue;
        rhoEff += table.rhoEff[rhoOffset + i] * rhoWeights[i];
        for (int j = 0; j < 4; ++j) {
            if (radiusWeights[j] == 0) continue;
            sr += table.EvalProfile(rhoOffset + i, radiusOffset + j) *
                  rhoWeights[i] * radiusWeights[j];
        }
    }

    // Cancel marginal PDF factor from tabulated BSSRDF profile
    if (rOptical != 0) sr /= 2 * Pi * rOptical;
    return std::max((Float)0, sr * sigma_t[ch] * sigma_t[ch] / rhoEff);
}

<Account for subsurface scattering, if applicable> ≡ 877
if (isect.bssrdf && (flags & BSDF_TRANSMISSION)) {
<Importance sample the BSSRDF 915>
<Account for the direct subsurface scattering component 915>
<Account for the indirect subsurface scattering component 916>
}

<Importance sample the BSSRDF> ≡ 915
SurfaceInteraction pi;
Spectrum S = isect.bssrdf->Sample_S(scene, sampler.Get1D(),
sampler.Get2D(), arena, &pi, &pdf);
if (S.IsBlack() || pdf == 0)
break;
beta *= S / pdf;

<Account for the direct subsurface scattering component> ≡ 915
L += beta * UniformSampleOneLight(pi, scene, arena, sampler);

<Account for the indirect subsurface scattering component> ≡ 915
Spectrum f = pi.bsdf->Sample_f(pi.wo, &wi, sampler.Get2D(), &pdf,
BSDF_ALL, &flags);
if (f.IsBlack() || pdf == 0)
break;
beta *= f * AbsDot(wi, pi.shading.n) / pdf;
specularBounce = (flags & BSDF_SPECULAR) != 0;
ray = pi.SpawnRay(wi);