PBRT阅读：第十七章 光传输 II : 体渲染

http://www.opengpu.org/forum.php?mod=viewthread&tid=7368&fromuid=10107

第17章 光传输II : 体渲染

SurfaceIntegrator是场景几何、材质、光源、于光传输方程相关的复杂算法、对场景中辐射亮度的分布的确定等等的交会点，同样地，VolumeIntegrator负责将参与介质的效果加入到这个进程之中，并确定它对辐射亮度的分布的影响。这一章主要介绍描述参与介质沿光线方向改变辐射亮度的传输方程
，并描述VolumeIntegrator的接口极其一些简单的实现。前一章中关于表面上的光传输算法的许多方法可以被扩展开来，对参与介质做相应的处理。

17.1 传输方程

传输方程是控制介质对光进行吸收、放射、散射等行为的基本方程。它解释了在第12章中介绍的所有的体散射过程--吸收、放射、内散射、外散射，并给出了描述辐射亮度在环境中分布的方程。光传输方程实际上是传输方程的一个特殊形式，只是缺少了参与介质，并只处理表面散射的情况。

传输方程的最基本的形式是一个积分-微分方程，它描述了沿一条光束上的辐射亮度在空间一个点上的变化情况。它可以变换成一个描述源自光线上无限个点的参与介质效果的纯积分方程。我们可以用很直接的方式得到该方程，即从沿一个光束增加能量的过程（放射，内散射）的效应减去同一个光束上减少能量的过程（吸收，外散射）的效应即可。回忆一下第12.1.4节的源项，它给出了给定点p在特定方向ω上由放射和内散射引起的辐射亮度变化：

源项负责说明所有使光线辐射亮度增加的过程。
衰减系数σt(p, ω)
负责说明在某点上所有使光线辐射亮度减少的过程，包括吸收和外散射。描述其效应的微分方程是：

我们将这两种效应加在一起，就得到沿光线上点p'的总体上的辐射亮度微分变化，从而得到传输方程的积分-微分方程：

在适当的边界条件下，我们可以将这个方程转换为一个纯积分方程。例如，如果我们假定场景中没有表面，那么光线就没有阻挡，故有无限的长度，传输积分方程就是：

其中p'= p + tω。这个方程的意义还是很直观的：它是说从给定方向到达一个点的辐射亮度等于从该点出发的光线上所有点对该点的累加辐射亮度。沿光线上每点上到光线原点的辐射亮度增加值要按从原点到该点的总光束透光率的比例减少。如图：

在更一般的情况下，如果场景中有反射或放射表面，光线长度就不是无限的，被光线碰到的表面会对其辐射亮度产生影响，使得该点上出离该表面的辐射亮度增加，并致使交点后面的光线上的点对光线原点的辐射亮度没有什么贡献值。如果光线(p,ω)在距离t处跟表面相交于点p0，那么传输积分方程是（17.1）：

其中p0= p + tω是表面上的点，而p' = p + tω是沿光线上的点，如图：

这个方程描述了沿光线的辐射亮度的两个贡献值：第一项是沿光线从表面上反射回来的辐射亮度，由Lo项给出，它包括表面上的放射辐射亮度和反射辐射亮度。这个辐射亮度可以被参与介质所衰减。第二项代表了由于体积散射和放射所产生的增加出来的辐射亮度，但到光线和表面的交点为止，过了这个交点对辐射亮度就没有什么影响了。

为了简明起见，这里不再过多地讨论传输方程。然而，就像光传输方程可以写成不同路径的累加和并引入重要性函数，我们也可以对传输方程做如此处理。这里只列出几个VolumeIntegrator的实现，其它诸如路径追踪、双向路径追踪、光子映射等等用于表面积分的算法也可以应用于体积分。

17.2 体积分器接口

VolumeIntegrator接口继承于Interator，将其中的Preprocess()、RequestSample()和Li()函数包括进来。体积分器的前两个函数的用法跟表面积分器的用法相同。Li()函数跟表面积分器的版本相似，都返回沿给定光线的辐射亮度，而体积分器应该假定光线已经跟场景中几何体相交，如果光线确实跟某个表面相交，Ray::maxt要设置成交点。这样，体积分器应该只计算参数范围[mint,maxt]中的体积散射效果。

VolumeIntegrator接口还增加了一个函数，Transmittance()，它负责计算从Ray::mint到Ray::maxt的光线上的光束透光率。

<Volume Scattering Declarations> +=
classVolumeIntegrator : public Integrator {
public:
virtualSpectrum Transmittance(const Scene *scene,
constRay &ray, const Sample *sample,

float*alpha) const = 0;
};

有了这个背景知识之后，我们就可以完全理解Scene::Li()函数了。它是对方程17.1的一个直接的实现。表面积分器计算在光线交点处的出射辐射亮度Lo，忽略了到光线原点的衰减效应。体积分器的Transmittance()函数计算到表面上点的光束透光率Tr，它的Li()函数给出参与介质所引起的沿光线的辐射亮度。LoTr跟来自参与介质的辐射亮度增加量的和给出了光线原点出的总辐射亮度值。

17.3 只有放射的积分器

最简单的体积分器可能就是忽略了内散射而只考虑放射和衰减的积分器了。因为忽略了内散射，在源项中的球面积分就不见了，简化后的传输方程为：

EmissionIntegrator用蒙特卡罗积分来解这个方程。

<Emissionlntegrator Declarations> =
classEmissionIntegrator : public VolumeIntegrator {
public:
<EmissionIntegratorPublic Methods>
private:
<EmissionIntegratorPrivate Data>
};

EmissionIntegrator的Transmittance()和Li()函数都要沿光线做一维积分求值。这里没有用固定数目的采样，而是根据光线在体积区域中的跑过的距离来决定采样数，距离越长，则采样越多。这个方法从直观意义上讲是很值得的：光线在介质的区间越长，所要求的精确度就越高，就需要越多的采样来捕捉沿光线的光学性质上的变化。采样个数由一个用户给定的参数（步长值）来间接地确定。光线被分成给定长度的线段，在每个线段中取一个采样。

<EmissionIntegrator Public Methods> =
EmissionIntegrator(floatss) { stepSize = ss; }
<EmissionIntegratorPrivate Data> =
floatstepSize;

Transmittance()和Li()函数各需一个一维采样值来做相应的积分求值计算。

<EmissionIntegrator Method Definitions> =
voidEmissionIntegrator::RequestSamples(Sample *sample,
constScene *scene)
tauSampleOffset =sample->Add1D(1);
scatterSampleOffset= sample->Add1D(1);
}

<EmissionIntegrator Private Data> +=
inttauSampleOffset, scatterSampleOffset;

Transmittance()函数还是很直接了当的。VolumeRegion的Tau()函数负责计算从光线起点到终点之间的光学厚度τ。这里积分器的工作只是选择一个步长大小，并将一个采样值传给函数，并返回e-τ。如果Tau()函数可以用解析法计算τ，就忽略这些额外的值。

这个函数应用了这样一个事实：Sample值只对相机光线是非NULL的，它被用来为阴影和间接光线增加步长大小，这样就降低了计算需求（和精度）。对于这些光线而言，精度的降低通常不会产生什么影响。

<EmissionIntegrator Method Definitions> +=
SpectrumEmissionIntegrator::Transmittance(const Scene *scene,
constRay &ray, const Sample *sample,
float*alpha) const {
if(!scene->volumeRegion) return Spectrum(1.f);
floatstep = sample ? stepSize : 4.f * stepSize;
floatoffset = sample ? sample->oneD[tauSampleOffset][0] :
RandomFloat();
Spectrumtau = scene->volumeRegion->Tau(ray, step, offset);
returnExp(-tau);
}

Li()函数负责对方程17.2的和的第二项求值。如果光线在t=t0进入体区域，那么从光线起点一直到点t0都不会有放射或衰减现象发生，Li()函数可以考虑在t0到t1进行积分，其中t1是光线出离体区域或碰上一个表面的最小参数偏置量。这个积分的估计值如下：

该值可以通过在t0到t1的光线上均匀选择采样点pi，然后求下列值：

其中p(pi)= 1 /(t1-t0)。在估计值公式中的Lve项可以直接用相应的volumeRegion函数来求得，用来求Tr的τ也可以直接求得（均质型天气或指数型天气），或通过15.7节中的蒙特卡罗积分求得。

为了做这个计算，Li()的实现先求积分的t的范围，并对t0,t1分别初始化：

<EmissionIntegrator Method Definitions> +=
SpectrumEmissionlntegrator::Li(const Scene *scene,
constRayDifferential &ray, const Sample *sample,
float*alpha) const {
VolumeRegion*vr = scene->volumeRegion;
floatt0,t1;
if(!vr || !vr->IntersectP(ray, &t0, &t1)) return O.f;
<Doemission-only volume integration in vr>
}

这里使用了两个另外两项技术。第一，就像第15.7节中VolumeRegion::Tau()函数使用均匀步长的采样点，这里的实现基于类似的理由也使用均匀的步长。第二，如果将点pi按照离光线原点p的距离排好序，就可以高效地对光束透光率Tr求值。然后，我们可以利用Tr的连乘的性质来逐步地从前一个点的值来计算当前点的值：

因为Tr(pj-> pi-1)所覆盖的距离比Tr(pj-> p)所覆盖的要短，如果使用蒙特卡罗求值，就可以使用更少的采样。

<Do emission-only volume integration in vr> =
SpectrumLv(0.);
<Preparefor volume integration stepping>
for(int i = 0; i < nSamples; ++i, t0 += step) {
<Advanceto sample at t0 and update T>
<Computeemission-only source term at p>
}
*T =Tr;
returnLv * step;

<Prepare for volume integration stepping> =
int N=Ceil2Int((t1-t0) / stepSize);
floatstep = (t1 - t0) / N;
SpectrumTr(1.f);
Point p= ray(t0), pPrev;
Vectorw = -ray.d;
if(sample)
t0+= sample->oneD[scatterSampleOffset][0] * step;
else
t-0+= RandomFloat() * step;

为了求当前点上的总透光率，需要求前前一个点到当前点之间的透光率，在乘上从光线原点到前一个点的透光率。

<Advance to sample at t0 and update T> =
pPrev =p;
p =ray(t0);
SpectrumstepTau = vr->Tau(Ray(pPrev, p - pPrev, 0, 1),
.5f* stepSize, RandomFloat());
Tr *=Exp(-stepTau);
<Possiblyterminate ray marching if transmittance is small>

在一个很厚的介质中，当光线穿过足够长的距离，透光率可能会变得很低。为了减少花在源项求值（它可能只有很小的对光线原点的辐射亮度的贡献值）上的时间，当透光率非常小时，光线的分步计算会随机地被俄罗斯轮盘法中断。

<Possibly terminate ray marching if transmittance issmall> =
if(Tr.y() < 1e-3) {
constfloat continueProb = .5f;
if(RandomFloat() > continueProb) break;
Tr/= continueProb;
}

有了前面的工作之后，计算该点上的源项就非常简单了：

<Computeemission-only source term at p> =
Lv +=Tr * vr->Lve(p, w);

17.4 单散射积分器

SingleScattering积分器除了考虑沿光线每点上的放射之外，来考虑由直接照明所引起的入射辐射亮度，但忽略了多方向散射引起的入射辐射亮度。这样，Li()函数要估算下列积分：

其中Ld只包括来自直接照明的辐射亮度。如图：

这个辐射亮度可能会被场景中的几何体阻挡住，也可能被光源和点p'之间的参与介质所衰减。虽然这种效果要比EmissionIntegrator所耗的计算时间要多很多，但它可以给出令人吃惊的光幕效果。

<SingleScattering Declarations> =
classSingleScattering : public VolumeIntegrator {
public:
<SingleScatteringPublic Methods>
private:
<SingleScatteringPrivate Data>
};
这个积分器的Li()函数使用了跟EmissionIntegrator相同的光线前进方法来求解传输方程。有一个不同之处是，该积分器在进入对采样位置的循环之前，先要为光源采样计算采样值。因为在调用Li()函数之前我们无法知道到底需要多少个采样（因为采样数取决于积分求值所用的光线线段长度），所以就不可能令Sampler生成采样并通过Sample传入Li()。因此，我们需要生成一组三维的拉丁超级立方体采样。第一个维用于在每个散射点上到底选择那个光源去采样，另外两个维用于在光源上采样。

<Compute sample patterns for single scatteringsamples> =
float*samp = (float *)alloca(3 * N * sizeof(float));
LatinHypercube(samp,N, 3);
intsampOffset = 0;

跟EmissionIntegrator相比另一个不同之处是对源项求值的方式。在沿光线的每个采样点上，下面的片段计算在点p上的源项的单散射近似值。它跟前面的<Computeemission-only source term at p>有相同的作用。我们用跟EmissionIntegrator相同的方法将体积放射包括进来，然后求在该点上的σs值，用适当的拉丁超级立方体采样选择要进行采样的光源，即它在该点上的散射贡献值。因为要在沿光线多个点上计算源项，我们在每个点上只对一个光源采样，然后将它的贡献值乘以光源的个数，这跟直接光照积分器中的“对一个光源采样”的策略相同。

<Compute single scattering source term at p> =
Lv +=Tr * vr->Lve(p, w);
Spectrumss = vr->sigma_s(p, w);
if(!ss.Black()) {
intnLights = scene->lights.size();
intlightNum = min(Floor2Int(samp[samp0ffset] * nLights),
nLights-1);
Light*light = scene->lights[lightNum];
<Addcontribution oflight due to scattering at p>
}
sampOffset+= 3;

在估算点p上的直接照明贡献值时，需要估算下列积分：

这里的实现并不同时对相函数和光源采样并使用重要性采样，而是在光源上选择采样位置，并直接计算估计值。对于那些并不是极端各向异性的介质而言，这个方法的效果还是很好的。

<Add contribution of light due to scattering at p> =
floatpdf;
VisibilityTestervis;
Vectorwi;
floatu1 = samp[samp0ffset+l], u2 = samp[samp0ffset+2];
SpectrumL = light->Sample_L(p, u1, u2, &wi, &pdf, &vis);
if(!L.B1ack() && pdf > O.f && vis.Unoccluded(scene)) {
SpectrumLd = L * vis.Transmittance(scene);
Lv+= Tr * ss * vr->p(p, w, -wi) * Ld * float(nLights) / pdf;
}