PBRT笔记(10)——体积散射

体散射类型

3个影响参与介质在环境中的辐射度分布的主要因素：

吸收：减少光能，并将其转化为别的能量，例如热量。
发光：由光子发射光能至环境中。
散射：由于粒子碰撞，使得一个方向的辐射度散射至其他方向。

吸收

吸收被描述为一段介质截面区域。每单位距离媒介密度与吸收光能的比被定义为 $\sigma_a$ 。截面区域或随着位置与方向的变化而变化，尽管如此，通常还是会定义为位置函数。它通常也被定义为一个光谱变化量， $\sigma_a$ 与距离成反比关系 $m^{-1}$ 。这意味这个值无需被判断是否处于(0,1)。

沿长度为dt的微分光线的辐射度变化可用微分方程来描述：

$L_o(p,\omega)-L_i(p,-\omega)=dL_o(p,\omega)=-\sigma_a(p,\omega)L_i(p,-\omega)dt$

求解分为方程后得积分方程。如果假设光线于介质内p点处，沿方向ω行进距离d，则全部吸收量为：

$e^{-\int_0^d\sigma_a(p+t\omega,\omega)dt}$

自发光

与吸收相对应，自发光过程是将能量转换成可见光的化学、热或核过程。自发光的辐射度变化微分方程为：

$dL_o(p,\omega)=L_e(p,\omega)dt$

但这是基于入射光不会对发光产生影响的前提下建立的假设。

外散射和衰减

当光线穿过某一介质时，由于粒子碰撞，导致其会往别的方向散射。它减少了出射位置处出射的辐射度，这种现象被称为外散射，与之对应的被称为内散射。

每单位距离外散射现象的发生概率由散射系数 $\sigma_s$ 决定。微分长度dt上的外散射辐射度减少量为：

$dL_o(p,\omega)=-\sigma_s(p,\omega)L_i(p,-\omega)dt$

总共减少的辐射度由吸收量与外散射量求和得到。这两种效果综合被称为衰减或者消光。
$\sigma_t(p,\omega)=\sigma_a(p,\omega)+\sigma_s(p,\omega)$ 两者相集合可以得到以下两个变量：

albedo： $\rho=\frac{\sigma_s}{\sigma_t}$ 他描述了发生散射现象的概率
平均自由路径:1/σt即发生散射现象前，光线在介质中行进的平均路径。

对于给定的衰减系数σt，描述整个衰减的微分方程为：

$\frac{dL_o(p,\omega)}{dt}=-\sigma_t(p,\omega)L_i(p,-\omega)$

求解之后可得：

$T_r(p\rightarrow p')=e^{-\int^d_0 \sigma_t(p+t\omega,\omega)dt}$

d表示p到p’的距离，Tr表示透射率。
如图11.8所示：

$T_r(p\rightarrow p'')=T_r(p\rightarrow p')+T_r(p'\rightarrow p'')$

Tr中的逆置属性τ被记作光学厚度：

$\tau(p\rightarrow p')=\int^d_0 \sigma_t(p+t\omega,-\omega)dt$

在均一介质中，σt是一个常量，根据Beer定律， $T_r(p\rightarrow p')=e^{-\sigma_td}$

内散射

内散射是因为其他方向上发生的散射现象，从而增加了出射方向的辐射度的现象。

在忽略介质粒子间互相作用的情况下，phase函数p(ω, ω’)描述了在这一点辐射度的角度分布情况。它用于体积模拟有点类似BSDF。phase函数有一个归一化条件，对于所有ω必须遵守：
$\int_{g2}p(\omega,\omega’)d\omega ‘=1$ 这就意味着phase函数定义特定方向的散射的分布概率。

根据内散射，内散射每单位距离增加的的辐射度为：

$dL_o(p,\omega)=L_s(p,\omega)dt$

他包含了自发光与内散射现象：

$L_s(p,\omega)=L_e(p,\omega)+\sigma_s(p,\omega)\int_{g2}L_i(p,\omega_i)d\omega_i$

相位函数

目前有许多相位函数，一般分为参数化模型(拥有少量变量的近似拟合函数)和散射辐射度解析模型。

大多数天然介质，其相位函数为角度在ωo与ωi之间的一维函数，这些被经常写作p(cos θ)。
因为在旋转时，不会对入射光产生影响，所以被称为各向同性。

除了归一化之外，另一个重要特征是:入射方向与出射方向是可以互换的，且phase函数的值不变。

在由排列成凝聚结构的粒子组成的各向异性介质中，它的相位函数是带两个方向参数的4维函数。它具有较为复杂的互相作用关系。（主要是晶体之类的）

相位函数自身既可以呈现为各项同性也呈现为各项异形的。因为与方向无关，所以可以定义为：

$p(\omega_o,\omega_i)=\frac{1}{4\pi}$

PhaseFunction类位于medium.h中，
p()返回相位函数对于指定方向的值。

Henyey与Greenstein在1941年开发了一种相位函数，现今依然被广泛使用。它的设计目的是为了便于与实测散射数据进行拟合。参数g（被称为非对称变量）用于控制散射分布。PhaseHG()对其进行了实现，其公式为：

$P_{HG}(cos\theta)=\frac{1}{4}\frac{1-g^2}{(1+g^2+2g(cos\theta))^{3/2}}$

这个模型的g值毕竟在(-1,1)范围内，g值为负时代表后向散射，光线主要散射回入射方向，g值为正时代表前向散射。g值越大散射现象就越进阶入射或者出射方向。

HenyeyGreenstein类实现了这个模型。

非对称变量g在这个模型中拥有明确意义，即为相位函数估算值与cos(ω’,ω)乘积的均值。

对于任意相位函数p,g的值可以计算为：

$g=\int_{g^2}p(-\omega,\omega')(\omega\cdot \omega')d\omega '=2\pi\int^\pi_0p(-cos\theta)cos\theta sin\theta d\theta$

所以各向同性介质的相位函数的g值为0。

任意数量的相位函数都可以满足这个方程;光靠一个g值不足以描述散射分布。然而方便地将一个复杂的散射分布转换成一个简单的参数化模型往往更加重要。

复杂的相位函数使用一个非对称变量是不能够很好地描述分布情况的，所以一般会采用相位函数加权和的方式来建模：

$p(\omega,\omega')=\sum_{i=1}^{n}w_i p_i (\omega \rightarrow \omega')$

其中wi的和需要被保持标准化为1。

介质

这段的完全机翻

介质基类实现了多种特性的体积散射属性，在复杂场景中可能有多个介质实例，用于表现不同的散射效果。

一个关键操作是介质类必须计算光束透射率。对应的函数为Tr()，这个方法返回光线起点到终点的透射率的估算值。

Medium-aware积分器负责计算与光线相交的介质几何图元的透射率，具体会使用Tr()与蒙特卡洛方法解积分方程的方式。因此我们假设传递给Tr()的射线是没有被阻挡，且光线完全处于当前介质中。

pbrt会将Medium实例与摄像机、灯光或者几何图元结合使用，用于模拟灯光雾等空间分布效果。

在pbrt中，为了表现场景中一块区域的介质，会使用几何图元来表现，具体的方式是使用MediumInterface类，它存储了外部与内部介质的指针。而不会只存储一种介质。（因为不存在什么介质都没有的情况）当然PBRT允许将测试设为nullptr，即不对光线产生任何影响。

bool GeometricPrimitive::Intersect(const Ray &r,
                                   SurfaceInteraction *isect) const {
    Float tHit;
    if (!shape->Intersect(r, &tHit, isect)) return false;
    r.tMax = tHit;
    isect->primitive = this;
    CHECK_GE(Dot(isect->n, isect->shading.n), 0.);
    // Initialize _SurfaceInteraction::mediumInterface_ after _Shape_
    // intersection
    if (mediumInterface.IsMediumTransition())
        isect->mediumInterface = mediumInterface;
    else
        isect->mediumInterface = MediumInterface(r.medium);
    return true;
}

但是这样有可能会出现不同物体所指定的外部介质不一样的情况。在这种情况下，离开图元朝向相机的光线与离开相机朝向图元的光线将被视为处于不同的介质中。因此光传输算法将无法计算一致的结果。pbrt认为用户能够在场景中指定一致的介质配置，并且不值得为此增加代码的复杂性。（毕竟是渲染核心，指定参数这种操作还是编辑器来操作比较好）

对于场景空间中的区域,我们将实现Scene::IntersectTr()方法,如果相交成功成功，它返回第一个相交的会发生光散射现象的SurfaceInteraction的指针。不然就返回false。

bool Scene::IntersectTr(Ray ray, Sampler &sampler, SurfaceInteraction *isect,
                        Spectrum *Tr) const {
    *Tr = Spectrum(1.f);
    while (true) {
        bool hitSurface = Intersect(ray, isect);
        // Accumulate beam transmittance for ray segment
        if (ray.medium) *Tr *= ray.medium->Tr(ray, sampler);

        // Initialize next ray segment or terminate transmittance computation
        if (!hitSurface) return false;
        if (isect->primitive->GetMaterial() != nullptr) return true;
        ray = isect->SpawnRay(ray.d);
    }
}

均匀介质

均匀介质是最为简单的介质，它代表了一个具有恒定σa和σs的空间，它使用了Henyey–Greenstein相位函数来表现这个介质中发生的散射现象，在这个相位函数中具有恒定的g值。类名为HomogeneousMedium。

因为整个媒介σt是恒定的,所以根据Beer定律可以用来计算沿光线的透射率。这里需要注意浮点数运算误差的问题。使用浮点数的正无穷来初始化Ray::tMax，它确保了光线足够可以取得最大范围内的相交信息。然而，对于Ray::tMax使用Infinity的情况下，在应用Beer定律时会产生了一个小问题。原则上，我们只需要计算射线参数t范围，乘以射线方向的长度,然后乘以σt，但是正无穷乘以0会得到NaN，具体发生在当光线通过一个吸收率为零的介质无限远时，上面的代码将尝试执行0 *∞的乘法运算，最终得到NaN，而不是0的预期透射率。所以这里需要做容错处理。

Spectrum HomogeneousMedium::Tr(const Ray &ray, Sampler &sampler) const {
    ProfilePhase _(Prof::MediumTr);
    return Exp(-sigma_t * std::min(ray.tMax * ray.d.Length(), MaxFloat));
}

3D网格

GridDensityMedium类讲介质的密度存储在3D网格中，类似于ImageTexture类用2D网格采样来表示图像的方式。GridDensityMedium通过对这些样本进行插值以计算他们之间的位置采样点的介质密度。

GridDensityMedium::Tr()会调用Density()，利用提供的样本重构给定点的体积密度函数。在这个过程中它将会用WorldToMedium()把世界坐标转换成局部坐标。σa和σs将会按照对应的位置进行插值运算。

首先调用Density()以周围样本坐标与距离插值计算点的坐标。（具体可以参照7.1.7节，说到底就是个双线性插值）

Float GridDensityMedium::Density(const Point3f &p) const {
    // Compute voxel coordinates and offsets for _p_
    Point3f pSamples(p.x * nx - .5f, p.y * ny - .5f, p.z * nz - .5f);
    Point3i pi = (Point3i)Floor(pSamples);
    Vector3f d = pSamples - (Point3f)pi;

    // Trilinearly interpolate density values to compute local density
    Float d00 = Lerp(d.x, D(pi), D(pi + Vector3i(1, 0, 0)));
    Float d10 = Lerp(d.x, D(pi + Vector3i(0, 1, 0)), D(pi + Vector3i(1, 1, 0)));
    Float d01 = Lerp(d.x, D(pi + Vector3i(0, 0, 1)), D(pi + Vector3i(1, 0, 1)));
    Float d11 = Lerp(d.x, D(pi + Vector3i(0, 1, 1)), D(pi + Vector3i(1, 1, 1)));
    Float d0 = Lerp(d.y, d00, d10);
    Float d1 = Lerp(d.y, d01, d11);
    return Lerp(d.z, d0, d1);
}

为了插值计算样本周围的点

D()返回给定整数坐标的样本密度，它的唯一任务是处理超出边界的样本位置，并为给定的样本计算适当的数组偏移量。与MIPMaps不同，在区域外密度讲返回为0。

BSSRDF

双向散射表面反射分布函数：

$L_o(p_o,\omega_o)=\int_A\int_{H^2(n)} S(p_o,\omega_o,p_i,\omega_i)L_i(p_i,\omega_i) \left|cos\theta_i \right|d\omega_idA$

次表面光线传输使用11.1和11.2中介绍的体积散射过程与15.1中介绍的体积光传输方程来描述。BSSRDF中的S是对边界上给定的一对点和方向之间发生的散射现象进行建模的一种概括表示。

类名为BSSRDF，位于core/bssrdf.h和bssrdf.cpp。
BSSDF必须实现传递当前相交表面以及散射介质折射率给基类的构造函数。因此，这里假设在这个介质中的折射率是不变的，这也是BSSRDF模型中广泛使用的假设。

BSSRDF实现必须提供的关键方法是估算八维分布函数S()。

为了支持一般形状图元，我们将引入一个更简单、可分离的SeparableBSSRDF。

SeparableBSSRDF(const SurfaceInteraction &po, Float eta,
const Material *material, TransportMode mode)
: BSSRDF(po, eta), ns(po.shading.n), ss(Normalize(po.shading.dpdu)),
ts(Cross(ns, ss)), material(material), mode(mode) { }
const Normal3f ns;
const Vector3f ss, ts;
const Material *material;
const TransportMode mode;

SeparableBSSRDF接口将将BSSRDF转换为三个独立的组成部分的可分离形式(一个空间和两个方向):

$S(p_o,\omega_o,p_i,\omega_i)\approx(1-F_r(cos\theta_o))S_p(p_o,p_i)S_\omega(\omega_i)$

Spectrum S(const SurfaceInteraction &pi, const Vector3f &wi) {
    Float Ft = 1 - FrDielectric(Dot(po.wo, po.shading.n), 1, eta);
    return Ft * Sp(pi) * Sw(wi);
}

对于上面给出的可分离变量的表达式，由次表面散射引起的照度的积分可简化为：

$L_o(p_o,\omega_o)=\int_A\int_{H^2(n)}S(p_o,\omega_o,p_i,\omega_i)L_i(p_i,\omega_i)\left|cos\theta_i\right|d\omega_idA(p_i)=(1-F_r(cos\theta_o))\int_AS_p(p_o,p_i)\int_{H^2(n)}S_\omega(\omega_i)L_i(p_i,\omega_i)\left|cos\theta_i\right|d\omega_idA(p_i)$

我们定义Sω(ωi)作为扩展版本的菲涅耳透光率

$S_\omega(\omega_i)=\frac{1-Fr(cos\theta_i)}{c\pi}$

归一化因子c是精心挑选的，所以Sw比上cos权重半球：

$\int_{H^2}S_\omega(\omega)cos\theta d\omega=1$

换句话说：

$c=\int^{2\pi}_0\int^{\frac{\pi}{2}}_{0}\frac{1-Fr(\eta,cos\theta)}{\pi}sin\theta\cos\theta d\theta d\phi=1-2\int^{\frac{\pi}{2}}_{0}F_r(\eta,cos\theta)sin\theta cos\theta d\theta$

这个积分称为菲涅尔反射函数的第一时刻。第i个菲涅耳时刻的一般定义是：

$F_r,i(\eta)=\int^{\frac{\pi}{2}}_{0}F_r(\eta,cos\theta)sin\theta cos^i \theta d\theta$

pbrt提供了FresnelMoment1()和FresnelMoment2()，用于拟合来计算相应的菲尼尔时刻。使用FresnelMoment1()可以很方便地实现SeparableBSSRDF::Sw()

Spectrum Sw(const Vector3f &w) const {
    Float c = 1 - 2 * FresnelMoment1(1 / eta);
    return (1 - FrDielectric(CosTheta(w), 1, eta)) / (c * Pi);
}

将空间和方向参数解耦大大减少了S的维数,但没有解决支持一般形状实现较为困难的问题。这里我们将引入第二个近似，它假设曲面不仅是局部平面的，而且影响BSSRDF值的是点之间的距离，而不是它们的实际位置。这将Sp简化为一个只涉及两点po和pi的距离的函数Sr:

$S_p(p_o,p_i)\approx S_r(\left|\right|p_o-p_i\left|\right|)$

Spectrum Sp(const SurfaceInteraction &pi) const {
    return Sr(Distance(po.p, pi.p));
}

SeparableBSSRDF在平面形状的情况下效果较好，随着形状逐渐偏离平面，误差将会越来越大。

TABULATED BSSRDF

类名为TabulatedBSSRDF，是SeparableBSSRDF的子类。TabulatedBSSRDF 实现了父类的接口，它提供表化的BSSRDF的表现方法。它可以控制广泛的次表面参数，可以用于模拟真实世界中的测算的BSSRDFs。

TabulatedBSSRDF使用与FourierBSDF 反射模型相同类型的自适应样条插值方法。

次表面散射材质

SubsurfaceMaterial，定义于materials/subsurface.h与materials/subsurface.cpp。KdSubsurfaceMaterial,定义于materials/kdsubsurface.h与materials/kdsubsurface.cpp。两者的区别在于如何设定散射属性。

const Float scale;
std::shared_ptr<Texture<Spectrum>> Kr, Kt, sigma_a, sigma_s;
std::shared_ptr<Texture<Float>> uRoughness, vRoughness;
std::shared_ptr<Texture<Float>> bumpMap;
const Float eta;
const bool remapRoughness;
BSSRDFTable table;

ComputeScatteringFunctions()使用贴图计算指定坐标的散射属性值，使用成员变量scale来缩放吸收值与散射系数。

void SubsurfaceMaterial::ComputeScatteringFunctions(SurfaceInteraction *si, MemoryArena &arena, TransportMode mode,bool allowMultipleLobes) const {
    <Perform bump mapping with bumpMap, if present >
    <Initialize BSDF for SubsurfaceMaterial>
    Spectrum sig_a = scale * sigma_a->Evaluate(*si).Clamp();
    Spectrum sig_s = scale * sigma_s->Evaluate(*si).Clamp();
    si->bssrdf = ARENA_ALLOC(arena, TabulatedBSSRDF)(*si, this, mode, eta, sig_a, sig_s, table);
}