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

体散射处理过程

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

1. 吸取：减小光能，并将其转化为别的能量，例如热量。
2. 发光：由光子发射光能至环境中。
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)\)二者相集合能够获得如下两个变量：

1. albedo：
   \(\rho=\frac{\sigma_s}{\sigma_t}\)
   他描述了发生散射现象的几率
2. 平均自由路径: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);
}