【离线渲染】【参与性介质渲染】Volume Rendering（一）

我们讨论光线传播方程的时候一般不会考虑上参与性介质，而是假定光子在真空中传播。但是现实世界中充满了各种介质，如空气、水、云等，这些介质在光线传播过程中会参与光线的交互，因此也称参与性介质。

当光在介质中传播时，其最终的辐射亮度会根据一下四种交互发生变化：1.吸收 absorbtion 2.辐射 emission 3.向外散射 out-scattering 4.向内散射 in-scattering

其中向外散射和吸收造成原有入射辐射亮度的损失，向内散射和自放光使得辐射亮度增加，在下文我们会分别讨论这四种交互，而在最后Volume Renderingt Equation中1/3方式会合在一起，2/4方式会合在一起。

光子与粒子发生交互时，它要么被粒子吸收，要么光子被散射到另一个方向。两种时间发生的相对概率分别表示为absorption coefficient(${\sigma }_a$)和scattering coefficient(${\sigma }_s$)。

Absorbtion

${\sigma }_a:$每单位长度被吸收的概率
单位：1/distance
上面的式子可以解释为一柱光线穿过单位长度（ds）的参与性介质后被吸收的能量 $dL(p,w)$为这束光线$L(p,w)$乘上吸收系数。

上式可转化为: $\frac{1}{L(p,w)}dL(p,w)=-{\sigma }_{a}ds(\ast 1)$

对两边进行积分：$ln(L(p,w)){|}_p^{p+sw}=-{\int }_0^s{\sigma }_a(p+s^{{}^{\prime }}w,w)d{s}^{{}^{\prime }}(\ast 2)$

=>$ln(p+sw,w)-ln(p,w)=-{\int }_0^s{\sigma }_a(p+s^{{}^{\prime }}w,w)d{s}^{{}^{\prime }}(\ast 3)$

=>$L(p+sw,w)=e^{-\tau (s)}L(p,w)=Tr(s)L(p,w)(\ast 4)$

由式子4我们可以得到Tr，即光束透射比，它描述的是光束在介质中两点之间沿光线传播后，保留下来的光子数量的分量，这里暂时忽略了散射。

$e^{-\tau (s)}=Tr(s)(\ast 5)$

对比式子3和4我们可以得到$\tau$，即optical depth。

$\tau (s)=-{\int }_0^s{\sigma }_a(p+s^{{}^{\prime }}w,w)d{s}^{{}^{\prime }}(\ast 6)$

我们把式子4中第一三项单独拿出来：

$L(p+sw,w)=Tr(s)L(p,w)(\ast 7)$

Out-Scatter

上面的式子是Out-Scatter的表达式，可以解释为一柱光线穿过单位长度（ds）的参与性介质后向外散射的能量 $dL(p,w)$为这束光线$L(p,w)$乘上散射系数。

吸收和散射描述的均是光子在单位长度距离传播过程中的衰减行为，我们把这两个系数合在一起，并将其命为衰减系数${\sigma }_t$。

${\sigma }_t={\sigma }_a+{\sigma }_s(\ast 8)$

由7式我们再得出一个变量,它也被称作反射率Albedo

$\rho =\frac{{\sigma }_s}{{\sigma }_t}(\ast 9)$

那么上面的5、6、7式就应做出相应修改：

$\tau (s)=-{\int }_0^s{\sigma }_t(p+s^{{}^{\prime }}w,w)d{s}^{{}^{\prime }}(\ast 10)$

在均匀介质中${\sigma }_t$是一个常数（均匀介质与非均匀介质在下文会介绍），对常数进行积分，我们得到：

$\tau (s)=-{\sigma }_{t}s(\ast 11)$

$L(p+sw,w)=Tr(s)L(p,w)(\ast 7)$

现在Volume Rendring Equation的第一项就有了，即7式。

解释一下7式的意义：7式描述了辐射亮度在介质中的传播中怎样被减少了。

11式介绍了如何在均匀介质中估计$\tau (s)$，现介绍一种在非均匀介质中估计$\tau (s)$的方法，RayMarching方法，其他方法在后文介绍。

$\tau (s)=-{\int }_0^s{\sigma }_t(p+s^{{}^{\prime }}w,w)d{s}^{{}^{\prime }}(\ast 10)$

现通过数值方法估计$\tau (s)$:

$\tau (s)\approx \frac{s}{N}{\Sigma }_i^N{\sigma }_t(x_i)$

$x_i=x+\frac{i+0.5}{N}w$

Emission

$dL=Ledt(\ast 12)$即Emission贡献给L的radiance。

In-Scattering

$S(p,w)={\sigma }_s(p){\int }_{S^2}p(p,w^{{}^{\prime }}\to w)L(p,w^{{}^{\prime }})d{w}^{{}^{\prime }}(\ast 13)$

=>$S(p,w)={\sigma }_s(p)L_i(x\to w)$

上面的式子可以解释为In-Scattering对$L_o$的贡献$S(p,w)$为点p处来自所有方向的光照最终散射到w方向上的总和${\int }_{S^2}p(p,w^{{}^{\prime }}\to w)L_i(p,w^{{}^{\prime }})d{w}^{{}^{\prime }}$乘上散射系数${\sigma }_s$

关于13式中的p函数（phase function）会在下文介绍。

Volume Rendering Equation

把12式和13式结合起来我们就得到了Ls

$L_s(p,w)=L_e(p,w)+{\sigma }_s(p,w){\int }_{S^2}p(p,w^{{}^{\prime }}\to w)L_i(p,w^{{}^{\prime }})d{w}^{{}^{\prime }}(\ast 14)$

在对从p到p+tw这段路径进行积分就得到了Volume Rendring Equation的第二项：

${\int }_0^{t}Tr(p^{{}^{\prime }}\to p)Ls(p^{{}^{\prime }},-w)d{t}^{{}^{\prime }}(\ast 15)$

15式表示从p到p+tw这一段上所有点$p^{{}^{\prime }}$来自emission和In-Scattering的radiance的积分和。

再把7式和15式结合起来就得到了Volume Rendering Equation：

$L_o(p,w)=Tr(p0\to p)L_i(p0,w)+{\int }_0^{t}Tr(p^{{}^{\prime }}\to p)Ls(p^{{}^{\prime }},-w)d{t}^{{}^{\prime }}(\ast 15)$

第一项来自7式，第二项来自15式，其中p0是来自surface上的点，而p’是来自Medium上的点。

Phase Function相位函数

$S(p,w)={\sigma }_s(p){\int }_{S^2}p(p,w^{{}^{\prime }}\to w)L(p,w^{{}^{\prime }})d{w}^{{}^{\prime }}(\ast 13)$

13式中的p被称作phase function，它描述光子在介质中某个位置处散射的方向分布。

相位函数的概念和表面的双向散射分布函数BRDF是类似的，不用的是相位函数用于介质内部，因此适用于整个球面空间的所有方向。此外BRDF返回的是一个完整的散射函数的光照值，而相位函数仅仅返回一个相位的值，它必须乘以散射系数${\sigma }_s$才能表示散射的光照结果。

我们重新定义散射函数为：
$\rho (x,w^{{}^{\prime }}\leftrightarrow w)=\{\begin{array}{c}{\rho }_s(x,w^{{}^{\prime }}\leftrightarrow w)x在表面上\\ p(x,w^{{}^{\prime }}\leftrightarrow w){\sigma }_s(x)x在介质中\end{array}(\ast 16)$

相位函数可以是各向同性的，它将光照均匀地散射到各个方向，对所有方向的值为一个常数$\frac{1}{4\pi }$。

广泛使用的一个各向异性相位函数是Henyey-Greenstein相位函数，pHG函数仅与散射的角度$\theta$有关。

$p_{HG}(x,\theta )=\frac{1-g^2}{4\pi (1+g^2-2gcos\theta {)}^{1.5}}(\ast 17)$

这里的唯一参数$g\in [-1,1]$用于决定散射的相对朝向：当g<0时相位呈后向分布，当g>0时相位呈前项分布，当g=0时变为各向同性分布。

PBRT中的phase function

inline Float PhaseHG(Float cosTheta, Float g) {
    Float denom = 1 + g * g + 2 * g * cosTheta;
    return Inv4Pi * (1 - g * g) / (denom * std::sqrt(denom));
}

Homogeneous Medium均匀介质

在非介质渲染中，光线只与物体表面发生碰撞，所以我们采样的时候我们很容易计算出采样到该碰撞点的pdf概率。但是在介质中的传播，光子可能与光线上任意粒子发生交互。

所以我们的目标变成获得一个光子在介质中传输的距离值，使得光子在与任意粒子发生交互之前可以在该距离内自由传播，我们称这样的一个路径称为自由路径。

由于光子在介质中的传播过程中与粒子发生交互的几率取决于前面介绍的衰减系数${\sigma }_t$，而透射比$Tr$实际上给出了一条光线在介质中两点之间无障碍传输的概率，因此我们可以根据一条光线上Tr的分布来对自由路径进行采样，以产生一个光子。

由5式和7式我们可以得出：

$Tr(p\to p^{{}^{\prime }})=e^{-{\sigma }_{t}d}(\ast 18)$

用逆变换方法从概率分布获取采样点

设$p(t)=c{e}^{-{\sigma }_{t}t}$

=>${\int }_0^{\infty }c{e}^{-{\sigma }_{t}t}dt=-\frac{c}{{\sigma }_t}{e}^{-{\sigma }_{t}t}{|}_0^{\infty }=\frac{c}{{\sigma }_t}=1$

=>$c={\sigma }_t$

=>$p(t)={\sigma }_t{e}^{-{\sigma }_{t}t}(\ast 19)$

$P(t)={\int }_0^t{\sigma }_t{e}^{-{\sigma }_t{t}^{{}^{\prime }}}d{t}^{{}^{\prime }}=1-e^{-{\sigma }_{t}t}(\ast 20)$

=>$P^{-1}(t)=-\frac{ln(1-t)}{{\sigma }_t}$

对t进行采样的采样点$t=-\frac{ln(1-rand)}{{\sigma }_t}rand\in [0,1](\ast 21)$

${\sigma }_t$一般随波长变化而变化，所以在PBRT中Sample函数随既选择了一个波长（RGB）的${\sigma }_t$进行采样。

int channel = std::min((int)(sampler.Get1D() * Spectrum::nSamples),
                           Spectrum::nSamples - 1);
    Float dist = -std::log(1 - sampler.Get1D()) / sigma_t[channel];

此文${\sigma }_t^i$表示不同波长的${\sigma }_t$

单个波长$p^i(t)={\sigma }_t^i{e}^{-{\sigma }_t^{i}t}(\ast 22)$

所有波长的平均采样密度(RGB波长表示中n为3)：

$p_t(t)=\frac{1}{n}{\Sigma }_{i=1}^n{\sigma }_t^i{e}^{-{\sigma }_t^{i}t}(\ast 23)$

$p_{surface}=1-{\int }_0^{tmax}{p}_t(t)dt=1-{\int }_0^{tmax}\frac{1}{n}{\Sigma }_{i=1}^n{\sigma }_t^i{e}^{-{\sigma }_t^{i}t}dt=1-\frac{1}{n}{\Sigma }_{i=1}^n\left\{-e^{-{\sigma }_t^{i}t}{|}_0^{tmax}\right\}=\frac{1}{n}{\Sigma }_{i=1}^n{e}^{-{\sigma }_t^{i}tmax}$

即$p_{surface}=\frac{1}{n}{\Sigma }_{i=1}^n{e}^{-{\sigma }_t^{i}tmax}(\ast 24)$

碰撞点要么在surface上要么在Medium上所以$p_{surface}+p_{medium}=1$

surface在【0，tmax】之外，medium在【0，tmax】之内$p_t(t)$是整个空间内采样，所以${\int }_0^{tmax}{p}_t(t)dt$是medium上的采样，$1-{\int }_0^{tmax}{p}_t(t)dt$是surface上的采样。

由24式得：
${\beta }_{surface}=\frac{Tr(p\to p+tw)}{p_{surface}}=\frac{e^{-{\sigma }_{t}t}}{\frac{1}{n}{\Sigma }_{i=1}^n{e}^{-{\sigma }_t^{i}tmax}}(\ast 25)$

由23式得：
${\beta }_{medium}=\frac{{\sigma }_s(p+tw)Tr(p\to p+tw)}{p_t(t)}=\frac{{\sigma }_s{e}^{-{\sigma }_{t}t}}{\frac{1}{n}{\Sigma }_{i=1}^n{\sigma }_t^i{e}^{-{\sigma }_t^{i}tmax}}(\ast 26)$

现在我们概述一下PBRT中Sample函数的流程：
1.从RGB三个波长中随机选择一个波长作为${\sigma }_t$的波长。
2.获取采样点t
3.如果t<tmax，那么就对介质进行采样
3.1生成一个mediumIntersection以返回给调用sample的函数。
4.如果如果t>tmax,则对物体表面进行采样。
5.计算Tr
6.根据对介质还是对表面采样计算不同的pdf，下面几行代码对应25、26两式。

    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);

Spectrum HomogeneousMedium::Sample(const Ray &ray, Sampler &sampler,
                                   MemoryArena &arena,
                                   MediumInteraction *mi) const {
	//求beta，即损失率
    ProfilePhase _(Prof::MediumSample);
    // Sample a channel and distance along the ray
    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));

    // Compute the transmittance and sampling density
    Spectrum Tr = Exp(-sigma_t * std::min(t, MaxFloat) * ray.d.Length());

    // Return weighting factor for scattering from homogeneous medium
    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;
    if (pdf == 0) {
        CHECK(Tr.IsBlack());
        pdf = 1;
    }
    return sampledMedium ? (Tr * sigma_s / pdf) : (Tr / pdf);
}

Heterogeneous Medium

透射比（Tr）计算和自由路径采样（Sample）是参与介质渲染中的两大基础问题。

透射比描述了一束光子在给定距离的传输后剩余下来的能量。

自由路径采样尝试计算一个光子能够在介质中自由传播的距离，也就是光子在第一次与介质中的粒子发生交互之前自由传播的距离。

DeltaTracking Sample（）

参考了此篇博客

PBRT中使用的方法为DeltaTracking，它通过使异质介质同质化，然后使用取舍方法对自由路径进行无偏估计，通俗来说就是把空的地方加上虚拟粒子（图中红色粒子），使他们都变成“均匀介质”。

在均匀介质中我们近似估计光子与介质中的粒子发生交互的距离为$t=-\frac{ln(1-rand)}{{\sigma }_t}rand\in [0,1](\ast 21)$。那么，
1.在将非均匀介质变成均匀介质后我们如何估计${\sigma }_t^{{}^{\prime }}$？
2.我们如何判断光子碰到了真实粒子？

在PBRT中，我们将整个介质划分成NXxNYxNZ个grid，在每个grid中存储了这个grid的密度，假设每个小正方体中最多可以装64个粒子，如果小正方体中装了32个粒子，则对应的密度为0.5。

在前文我们知道${\sigma }_t$描述的是光束穿过一段介质的衰减行为，而对于非均匀介质，从不同方向穿过衰减也肯定是不同的。为了估计转化成“均匀介质”后的${\sigma }_t$，PBRT先找出所有grid里最大的密度，然后乘上衰减系数作为转化后的衰减系数。即：
${\sigma }_{new}=maxDensity\ast {\sigma }_t={\sigma }_t/invMaxDensity$

21式变成了$t=-\frac{ln(1-rand)}{{\sigma }_{new}}rand\in [0,1](\ast 27)$

我们通过一个随机数来判断光子是否碰到了真实粒子。首先根据采样点p计算该点的粒子密度，对于该点粒子密度的计算我们通过一个三线性插值函数计算出来，也就是PBRT中的Density函数。如果$rand<\frac{density(p)}{maxDensity}(\ast 28)$，那么碰到的粒子就是真实粒子

算法流程：
1.用一个bounding box将整个Medium包起来，光线撞击box得到两个点t_min、t_max。

    const Bounds3f b(Point3f(0, 0, 0), Point3f(1, 1, 1));
    Float tMin, tMax;
    if (!b.IntersectP(ray, &tMin, &tMax)) return Spectrum(1.f);

2.根据27式采样一个采样点

t -= std::log(1 - sampler.Get1D()) * invMaxDensity / sigma_t;

2.1如果该采样点$t>tmax$，那么采样点已经离开了Medium，该采样点位于Surface上，返回${\beta }_{surface}=1$

if (t >= tMax) break;

 return Spectrum(1.f);

2.2根据28式，如果碰到了真实粒子，记录交点并返回${\beta }_{medium}={\sigma }_s/{\sigma }_t$

        if (Density(ray(t)) * invMaxDensity > sampler.Get1D()) {
            PhaseFunction *phase = ARENA_ALLOC(arena, HenyeyGreenstein)(g);
            *mi = MediumInteraction(rWorld(t), -rWorld.d, rWorld.time, this,
                                    phase);
            return sigma_s / sigma_t;
        }

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;//beta medium
        }
    }
    return Spectrum(1.f);//beta surface
}

Tr投射比估计/比例追踪

在增量跟踪算法中，一条自由路径的传输总是按照每个距离处的真实粒子的比例来概率性地决定其是否需要继续传输或者被中止，因此每个路径采样对投射比的贡献要么是0（被中止），要么是1（继续传输）。

上述过程可以看做是一个基于俄罗斯轮盘的随机行走的过程，只不过其用的中止概率为0或1.与此相反，我们可以将中止概率设置为局部真实粒子的比例。

比例追踪在每一步距离采样中，它对路径分配一个与该虚拟粒子比例相等的权重系数，并使每条路径持续随机行走直到满足超出指定的输出长度，最终每条路径的贡献值被所有与虚拟粒子发生碰撞的联合概率加权。

$Tr(d)={\prod }_{i=1}^k(1-\frac{density(xi)}{maxDensity})$

 Tr *= 1 - std::max((Float)0, density * invMaxDensity);

这里被乘数表示每个碰撞发生的位置x处虚拟粒子相对于所有粒子的比例。

比例最终的优点是它很好地利用了路径跟踪过程中的信息而不是简单地使用一份二元判断。相比增量追踪，它与真实的透射比曲线的形状更加接近。

Spectrum GridDensityMedium::Tr(const Ray &rWorld, Sampler &sampler) const {
   
    ProfilePhase _(Prof::MediumTr);    
    ++nTrCalls;
   
    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);

    // Perform ratio tracking to estimate the transmittance value
    Float Tr = 1, t = tMin;
    while (true) {
        ++nTrSteps;
        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);
}

参考：PBRT，全局光照技术。