这个系列是Physically Based Rendering: From Theory To Implementation的读书笔记，本节主要介绍PBRT中的体散射。
目前我们只考虑了光线在真空中进行传输的情况。对于现实世界中的物体，由于参与介质(participating media)的存在相机接收到的能量会比物体实际反射的能量要少一些，这种现象在包含云、烟雾等介质的场景中尤为明显。因此我们需要体散射(volume scattering)来描述光线与参与介质的相互作用。
Volume Scattering Processes
光线在参与介质中主要有三种不同的行为：
- 吸收(absorption)：光线的能量会被介质吸收并转化成其它形式的能量，如热能等；
- 发射(emission)：介质自身会发射出能量；
- 散射(scattering)：某个方向上光线的能量会散射到其它方向上；
如果介质的这些性质在空间中是一致的则称为同质(homogeneous)，否则称为是异质(inhomogeneous)的。
Absorption
对于吸收的情况，出射和入射光线的能量(radiance)满足方程：
\[L_o(\text{p}, \omega) - L_i(\text{p}, -\omega) = d L_o(\text{p}, \omega) = -\sigma_a(\text{p}, \omega) L_i(\text{p}, -\omega) dt\]其中$\sigma_a$为光线在介质中沿$\omega$方向传输单位长度被吸收能量的概率密度。$\sigma_a$的单位是m-1，它是一个正数但不要求在区间(0, 1)上。
求解吸收能量的微分方程可以得到光线从$\text{p}$点出发沿$\omega$方向传输$d$距离后剩余能量的比例为：
\[e^{-\int_0^d \sigma_a(\text{p} + \omega t, \omega) dt}\]Emission
对于介质发射能量的情况，假设单位距离增加的能量为$L_e(\text{p}, \omega)$且与入射能量$L_i(\text{p}, -\omega)$无关，可以得到微分方程：
\[d L_o(\text{p}, \omega) = L_e(\text{p}, \omega) dt\]
Out-Scattering and Attenuation
对于散射的现象有两种情况：其一是当前方向上的能量散射到其它方向，称为外散射(out-scattering)；另一种情况是其它方向上的能量散射到当前方向上，称为内散射(in-scattering)。
外散射的情况与吸收类似，我们可以引入一个系数$\sigma_s$来描述散射的能量损失，这样光线传输能量的衰减(attenuation or extinction)为：
\[\sigma_t(\text{p}, \omega) = \sigma_a(\text{p}, \omega) + \sigma_s(\text{p}, \omega)\]
与$\sigma_t$相关的还有两个量：首先是反照率(albedo)，它表示散射占能量总体损失的比例，记为$\rho$：
\[\rho = \frac{\sigma_s}{\sigma_t}\]显然$\rho$介于0和1之间，它描述了发生散射的概率。另一项是平均自由程(mean free path)，$1/\sigma_t$，它表示光线和介质中的粒子发生碰撞的平均路径长度。
利用衰减系数$\sigma_t$我们可以得到光线能量损失的微分方程：
\[\frac{d L_o(\text{p}, \omega)}{d t} = -\sigma_t L_i(\text{p}, -\omega)\]求解后可以得到光线在两点上的透射(transmittance)系数：
\[T_r(\text{p} \rightarrow \text{p}') = e^{-\int_0^d \sigma_t(\text{p} + t \omega, \omega) dt}\]其中$d = \Vert \text{p} - \text{p}’ \Vert$为$\text{p}$和$\text{p}’$两点之间的距离；$\omega$为$\text{p}$指向$\text{p}’$的方向向量。注意到$T_r$是一个(0, 1)区间上的量，从$\text{p}$点发射的光线经过介质后在$\text{p}’$处接收到的能量为：
\[T_r(\text{p} \rightarrow \text{p}') \ L_o(\text{p}, \omega)\]
透射系数$T_r$有很多重要的性质：
- 真空中$T_r = 1$；
- 对于大多数介质衰减系数是对称的$\sigma_t(\omega) = \sigma_t(-\omega)$，因此$T_r(\text{p} \rightarrow \text{p}’) = T_r(\text{p}’ \rightarrow \text{p})$；
- 对于任意介质$T_r$是累乘的$T_r(\text{p} \rightarrow \text{p}’’) = T_r(\text{p} \rightarrow \text{p}’) \cdot T_r(\text{p}’ \rightarrow \text{p}’’)$；
$T_r$的指数部分取正号后称为光学深度(optical depth)，一般用$\tau$来表示：
\[\tau (\text{p} \rightarrow \text{p}') = \int_0^d \sigma_t(\text{p} + t \omega, \omega) \ dt\]对于同质的介质，$\sigma_t$为常数，此时光学深度即为衰减系数$\sigma_t$与距离$d$的乘积，对应的透射系数$T_r$为：
\[T_r(\text{p} \rightarrow \text{p}') = e^{-\sigma_t d}\]In-scattering
内散射的情况与发射类似，它们都会使观测方向上的能量增加。为了考虑不同方向上内散射的贡献，我们需要引入相位函数(phase function)$p(\omega, \omega’)$。相位函数类似于BSDF，不过它们的归一化条件不同：
\[\int_{S^2} p(\omega, \omega') \ d \omega' = 1\]根据归一化条件可以发现，相位函数实际上定义了球面上任意方向$\omega’$散射到$\omega$方向的概率分布。
把发射和内散射结合起来就得到了能量增加项：
\[d L_o(\text{p}, \omega) = L_s(\text{p}, \omega) \ dt\] \[L_s(\text{p}, \omega) = L_e(\text{p}, \omega) + \sigma_s(\text{p}, \omega) \int_{S^2} p(\text{p}, \omega, \omega') \ L_i(\text{p}, \omega') \ d \omega'\]其中内散射对应$\sigma_s(\text{p}, \omega) \int_{S^2} p(\text{p}, \omega, \omega’) \ L_i(\text{p}, \omega’) \ d \omega’$。同时不难发现$L_s$与散射方程非常相似，主要区别在于相位函数是直接定义在radiance上的，因此无需余弦项来计算irradiance。
Phase Functions
对于大多数自然界中的介质，相位函数是两个方向$\omega_i$和$\omega_o$夹角的函数，因此也常常记为$p(\cos \theta)$。这样的介质也称为各向同性(isotropic)介质，显然各向同性介质的相位函数满足旋转不变性。同时各向同性介质还满足对称性，即对称$\omega_i$和$\omega_o$函数的值不变，即$p(\cos \theta) = p(\cos -\theta)$。
而对于各向异性(anisotropic)介质，它的相位函数则需要使用包含两个方向的四元函数来表示。这里需要注意的是各向异性介质的相位函数同样需要满足对称性。
在PBRT中使用PhaseFunction类定义了所有相位函数的接口。
class PhaseFunction {
public:
virtual Float p(const Vector3f &wo, const Vector3f &wi) const = 0;
virtual Float Sample_p(const Vector3f &wo, Vector3f *wi,
const Point2f &u) const = 0;
};
其中PhaseFunction::p()方法用来计算给定两个方向相位函数的值。
最为常用的相位函数是Henyey–Greenstein相位函数。通过对真实介质散射行为的测量，拟合得到介质的相位函数为：
\[p_{\text{HG}} (\cos \theta) = \frac{1}{4 \pi} \frac{1 - g^2}{(1 + g^2 + 2g \cos \theta)^{3/2}}\]其中参数$g$位于区间(-1, 1)，用来控制散射光的分布。当$g$取负值时表示介质会出现反向散射(back-scattering)的现象，即大部分散射光线会散射到后方上(蓝色实线)；而$g$取正值时散射光线只会指向前方(橙色虚线)。同时$g$的绝对值越大，散射方向越接近入射方向(或其反方向)。
在PBRT中通过PhaseHG()函数实现了Henyey–Greenstein相位函数的计算：
inline Float PhaseHG(Float cosTheta, Float g) {
Float denom = 1 + g * g + 2 * g * cosTheta;
return Inv4Pi * (1 - g * g) / (denom * std::sqrt(denom));
}
在此基础上，Henyey–Greenstein相位函数由HenyeyGreenstein类实现：
class HenyeyGreenstein : public PhaseFunction {
public:
HenyeyGreenstein(Float g) : g(g) { }
Float p(const Vector3f &wo, const Vector3f &wi) const;
Float Sample_p(const Vector3f &wo, Vector3f *wi, const Point2f &sample) const;
private:
const Float g;
};
Float HenyeyGreenstein::p(const Vector3f &wo, const Vector3f &wi) const {
return PhaseHG(Dot(wo, wi), g);
}
实际上Henyey–Greenstein相位函数的$g$值还有明确的数学意义，它表示$\cos \theta$在$p$分布下的平均值(数学期望)。
\[g = \int_{S^2} p(-\omega, \omega') \ (\omega \cdot \omega') \ d \omega' = 2\pi \int_0^{\pi} p(-\cos \theta) \cos \theta \sin \theta \ d \theta\]显然对于各向同性的相位函数$g = 0$。
对于更加复杂的相位函数我们可以使用一些简单的相位函数通过加权叠加得到：
\[p(\omega, \omega') = \sum_{i=1}^n w_i p_i(\omega \rightarrow \omega')\]其中权重项$w_i$需要满足归一化条件$\sum_{i=1}^n w_i = 1$。
Media
在PBRT中介质使用Medium类作为通用接口，其中Medium::Tr()函数用来计算介质的透射系数。
class Medium {
public:
virtual ~Medium() { }
virtual Spectrum Tr(const Ray &ray, Sampler &sampler) const = 0;
virtual Spectrum Sample(const Ray &ray, Sampler &sampler,
MemoryArena &arena, MediumInteraction *mi) const = 0;
};
同时对于复杂的场景可能存在多种不同的介质，因此在PBRT中使用MediumInterface来表示两种不同介质之间的界面。GeometricPrimitive实例会存储一个MediumInterface指针用来存储两种介质的信息。
struct MediumInterface {
MediumInterface(const Medium *medium)
: inside(medium), outside(medium) { }
MediumInterface(const Medium *inside, const Medium *outside)
: inside(inside), outside(outside) { }
bool IsMediumTransition() const { return inside != outside; }
const Medium *inside, *outside;
};
因此在计算光线与场景交点时需要同时添加交点处的界面信息：
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;
// Initialize SurfaceInteraction::mediumInterface after Shape intersection
if (mediumInterface.IsMediumTransition())
isect->mediumInterface = mediumInterface;
else
isect->mediumInterface = MediumInterface(r.medium);
return true;
}
在很多情况下我们需要处理由参与介质构成几何模型的情况。此时光线与模型相交后不会改变前进路线，而是进入到介质中继续传播。为了处理这种情况，在PBRT中允许表面的材质指针为nullptr，这表示相交表面不会影响光线的传输；同样地，此时相交对象的BSDF指针SurfaceInteraction::bsdf也为nullptr，表示光线不会在表面发生散射。
在此基础上我们可以实现包含参与介质的场景与光线求交函数Scene::IntersectTr()，它会在计算光线和场景交点的基础上考虑参与介质对光线能量的影响：
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);
}
}
Medium Interactions
为了描述体散射我们还需要更新一些交点类Interaction的构造函数：
Interaction(const Point3f &p, const Vector3f &wo, Float time,
const MediumInterface &mediumInterface)
: p(p), time(time), wo(wo), mediumInterface(mediumInterface) { }
Interaction(const Point3f &p, Float time,
const MediumInterface &mediumInterface)
: p(p), time(time), mediumInterface(mediumInterface) { }
bool IsMediumInteraction() const { return !IsSurfaceInteraction(); }
如果交点设置的法线，还可以利用它来判断光线所处的介质：
const Medium *GetMedium(const Vector3f &w) const {
return Dot(w, n) > 0 ? mediumInterface.outside :
mediumInterface.inside;
}
对于光线位于物体内侧的情况则无需判断直接返回介质：
const Medium *GetMedium() const {
Assert(mediumInterface.inside == mediumInterface.outside);
return mediumInterface.inside;
}
把它们组合到一起就得到了考虑参与介质的交点类MediumInteraction：
class MediumInteraction : public Interaction {
public:
MediumInteraction(const Point3f &p, const Vector3f &wo, Float time,
const Medium *medium, const PhaseFunction *phase)
: Interaction(p, wo, time, medium), phase(phase) { }
bool IsValid() const { return phase != nullptr; }
const PhaseFunction *phase;
Homogeneous Medium
对于homogeneous medium，我们需要$\sigma_a$和$\sigma_s$来描述光线的衰减。同时我们使用Henyey–Greenstein相位函数来描述散射的行为，因此还需要一个参数$g$。这样就可以得到HomogeneousMedium类的完整定义：
class HomogeneousMedium : public Medium {
public:
HomogeneousMedium(const Spectrum &sigma_a, const Spectrum &sigma_s, Float g)
: sigma_a(sigma_a), sigma_s(sigma_s), sigma_t(sigma_s + sigma_a), g(g) { }
Spectrum Tr(const Ray &ray, Sampler &sampler) const;
Spectrum Sample(const Ray &ray, Sampler &sampler, MemoryArena &arena, MediumInteraction *mi) const;
private:
const Spectrum sigma_a, sigma_s, sigma_t;
const Float g;
};
由于$\sigma_t$是一个常量，我们只需要将它和光线传播的距离乘起来就得到了光学深度进而计算透射系数。当然在编写程序时还需要注意一些数值计算的精度和溢出问题，因此在PBRT中还对路径的长度进行了截断：
Spectrum HomogeneousMedium::Tr(const Ray &ray, Sampler &sampler) const {
return Exp(-sigma_t * std::min(ray.tMax * ray.d.Length(), MaxFloat));
}
3D Grids
对于inhomogeneous medium，PBRT使用GridDensityMedium类来存储空间中不同位置上介质的密度。
class GridDensityMedium : public Medium {
public:
GridDensityMedium(const Spectrum &sigma_a, const Spectrum &sigma_s,
Float g, int nx, int ny, int nz, const Transform &mediumToWorld,
const Float *d)
: sigma_a(sigma_a), sigma_s(sigma_s), g(g), nx(nx), ny(ny), nz(nz),
WorldToMedium(Inverse(mediumToWorld)),
density(new Float[nx * ny * nz]) {
memcpy((Float *)density.get(), d, sizeof(Float) * nx * ny * nz);
// Precompute values for Monte Carlo sampling of GridDensityMedium
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;
}
Float Density(const Point3f &p) const;
Float D(const Point3i &p) const {
Bounds3i sampleBounds(Point3i(0, 0, 0), Point3i(nx, ny, nz));
if (!InsideExclusive(p, sampleBounds))
return 0;
return density[(p.z * ny + p.y) * nx + p.x];
}
Spectrum Sample(const Ray &ray, Sampler &sampler, MemoryArena &arena,
MediumInteraction *mi) const;
Spectrum Tr(const Ray &ray, Sampler &sampler) const;
private:
const Spectrum sigma_a, sigma_s;
const Float g;
const int nx, ny, nz;
const Transform WorldToMedium;
std::unique_ptr<Float[]> density;
Float sigma_t;
Float invMaxDensity;
};
GridDensityMedium::Density()函数通过三线性插值计算空间中任意点处的密度。需要注意的是这里$p$点使用的是规范化坐标[0, 1]3，因此需要首先转换到物体坐标系中。
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);
}
GridDensityMedium::D()函数用来返回采样点位置处的密度。
Float D(const Point3i &p) const {
Bounds3i sampleBounds(Point3i(0, 0, 0), Point3i(nx, ny, nz));
if (!InsideExclusive(p, sampleBounds))
return 0;
return density[(p.z * ny + p.y) * nx + p.x];
}
The 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, \omega_i) \vert \cos{\theta_i} \vert \ d \omega_i \ dA\]实际上我们可以使用体散射来描述这一过程，因此整个渲染的核心在于如何设置BSSRDF项$S$。
在PBRT中各类BSSRDF的通用接口为BSSRDF类：
class BSSRDF {
public:
BSSRDF(const SurfaceInteraction &po, Float eta)
: po(po), eta(eta) { }
virtual Spectrum S(const SurfaceInteraction &pi, const Vector3f &wi) = 0;
virtual Spectrum Sample_S(const Scene &scene, Float u1, const Point2f &u2,
MemoryArena &arena, SurfaceInteraction *si, Float *pdf) const = 0;
protected:
const SurfaceInteraction &po;
Float eta;
};
在初始化时BSSRDF必须提供出射点po以及材质的折射率eta，而在求值时会调用BSSRDF::S()函数计算pi位置以ei角度入射后材质返回的BSSRDF函数值。
Separable BSSRDFs
可以想象即使是对平面或是球面这种简单的几何形状想要精确地计算BSSRDF都是十分困难的。为了不失一般性，PBRT使用了一种简化的模型并将BSSRDF拆分成3个独立的部分，这样的简化模型通过SeparableBSSRDF类来表示。
其中第一项$(1 - F_r(\cos \theta_o))$对应从出射点通过透射离开材质的比例；第二项$S_p(p_o, p_i)$表示光线从入射点到出射点进行传输的比例；第三项$S_\omega(\omega_i)$则是从入射位置进入材质的能量比例。对于高 albedo的材质，它们往往表现出更强的各向同性而且Fresnel项对最终的方向相关的结果影响更大，因此这一近似方法较为准确；而对于低albedo的材质这样的近似则不够准确。
把近似后的$S(p_o, \omega_o, p_i, \omega_i)$带入次表面散射方程可以得到：
\[\begin{aligned} L_o (p_o, \omega_o) &= \int_A \int_{H^2(n)} S(p_o, \omega_o, p_i, \omega_i) L_i(p, \omega_i) \vert \cos{\theta_i} \vert \ d \omega_i \ dA \\ &= (1 - F_r(\cos \theta_o)) \int_A S_p(p_o, p_i) \int_{H^2(n)} S_\omega(\omega_i) \ L_i(p, \omega_i) \vert \cos{\theta_i} \vert \ d \omega_i \ dA(p_i) \end{aligned}\]我们把$S_\omega(\omega)$定义为规范化后的Fresnel透射：
\[S_\omega(\omega) = \frac{1 - F_r(\cos \theta)}{c \pi}\]其中常数$c$需要满足在半球面上的归一化条件：
\[\int_{H^2} S_\omega(\omega) \cos \theta \ d \omega = 1\]整理后可以得到：
\[\begin{aligned} c &= \int_0^{2 \pi} \int_0^{\frac{\pi}{2}} \frac{1 - F_r(\eta, \cos \theta)}{\pi} \sin \theta \cos \theta \ d \theta \ d \phi \\ &= 1 - 2 \int_0^{\frac{\pi}{2}} F_r(\eta, \cos \theta) \sin \theta \cos \theta d \theta \end{aligned}\]其中的积分项称为Fresnel反射函数的一阶矩。类似地，$i$阶矩可以定义为：
\[\bar{F}_{r, i} (\eta) = \int_0^{\frac{\pi}{2}} F_r(\eta, \cos \theta) \sin \theta \cos^i \theta d \theta\]对于空间项$S_p(p_o, p_i)$，PBRT将其简化为只与两点空间位置有关的函数：
\[S_p(p_o, p_i) \approx S_r (\Vert p_o - p_i \Vert)\]其中$S_r$函数可由继承的子类自行实现。我们把整个流程梳理一下就得到了SeparableBSSRDF类的完整定义：
class SeparableBSSRDF : public BSSRDF {
public:
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) { }
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);
}
Spectrum Sw(const Vector3f &w) const {
Float c = 1 - 2 * FresnelMoment1(1 / eta);
return (1 - FrDielectric(CosTheta(w), 1, eta)) / (c * Pi);
}
Spectrum Sp(const SurfaceInteraction &pi) const {
return Sr(Distance(po.p, pi.p));
}
Spectrum Sample_S(const Scene &scene, Float u1, const Point2f &u2,
MemoryArena &arena, SurfaceInteraction *si, Float *pdf) const;
Spectrum Sample_Sp(const Scene &scene, Float u1, const Point2f &u2,
MemoryArena &arena, SurfaceInteraction *si, Float *pdf) const;
Float Pdf_Sp(const SurfaceInteraction &si) const;
virtual Spectrum Sr(Float d) const = 0;
virtual Float Sample_Sr(int ch, Float u) const = 0;
virtual Float Pdf_Sr(int ch, Float r) const = 0;
private:
const Normal3f ns;
const Vector3f ss, ts;
const Material *material;
const TransportMode mode;
};
Tabulated BSSRDF
和BRDF类似，我们同样可以使用一个表格来存储材质的BSSRDF，这在PBRT中使用TabulatedBSSRDF类来表示。
class TabulatedBSSRDF : public SeparableBSSRDF {
public:
TabulatedBSSRDF(const SurfaceInteraction &po,
const Material *material, TransportMode mode, Float eta,
const Spectrum &sigma_a, const Spectrum &sigma_s,
const BSSRDFTable &table)
: SeparableBSSRDF(po, eta, material, mode), table(table) {
sigma_t = sigma_a + sigma_s;
for (int c = 0; c < Spectrum::nSamples; ++c)
rho[c] = sigma_t[c] != 0 ? (sigma_s[c] / sigma_t[c]) : 0;
}
Spectrum Sr(Float distance) const;
Float Pdf_Sr(int ch, Float distance) const;
Float Sample_Sr(int ch, Float sample) const;
private:
const BSSRDFTable &table;
Spectrum sigma_t, rho;
};
需要注意的是当材质确定时$S_r$实际上只是一个一维函数。更确切地说$S_r$是一个关于材质折射率$\eta$、散射各向异性$g$、反照率(albedo)$\rho$、衰减率$\sigma_t$以及距离$r$的函数。对于这样一个高维的函数直接进行离散化是不现实的，因此需要对这些参数进行简化。简化后的$S_r$使用固定的$\eta$和$g$，同时将物理距离$r$转换为无量纲的光学距离$r_\text{optical} = \sigma_t r$。这样$S_r$就化简为：
\[S_r(\eta, g, \rho, \sigma_t, r) = \sigma_t^2 S_r(\eta, g, \rho, 1, r_\text{optical})\]这样我们就可以使用$\rho$和$r_\text{optical}$组成的表格来表示$S_r$了。在PBRT中使用BSSRDFTable结构体来存储这样的表格。
struct BSSRDFTable {
const int nRhoSamples, nRadiusSamples;
std::unique_ptr<Float[]> rhoSamples, radiusSamples;
std::unique_ptr<Float[]> profile;
std::unique_ptr<Float[]> rhoEff;
std::unique_ptr<Float[]> profileCDF;
BSSRDFTable(int nRhoSamples, int nRadiusSamples);
inline Float EvalProfile(int rhoIndex, int radiusIndex) const {
return profile[rhoIndex * nRadiusSamples + radiusIndex];
}
};
还需要注意的是TabulatedBSSRDF::rho表示的是单次散射的能量损失，而材质的albedo则是经过多次散射后能量损失。为了区分这两个概念我们把后者记为$\rho_\text{eff}$，它存储在BSSRDFTable::rhoEff中：
利用TabulatedBSSRDF实例就可以实现TabulatedBSSRDF::Sr()函数完成BSSRDF了：
Spectrum TabulatedBSSRDF::Sr(Float r) const {
Spectrum Sr(0.f);
for (int ch = 0; ch < Spectrum::nSamples; ++ch) {
// Convert r into unitless optical r_radius
Float rOptical = r * sigma_t[ch];
// Compute spline weights to interpolate BSSRDF 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))
continue;
// Set BSSRDF value Sr[ch] using tensor spline interpolation
Float sr = 0;
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
Float weight = rhoWeights[i] * radiusWeights[j];
if (weight != 0)
sr += weight * table.EvalProfile(rhoOffset + i,
radiusOffset + j);
}
}
// Cancel marginal PDF factor from tabulated BSSRDF profile
if (rOptical != 0)
sr /= 2 * Pi * rOptical;
Sr[ch] = sr;
}
// Transform BSSRDF value into world space units
Sr *= sigma_t * sigma_t;
return Sr.Clamp();
}
Subsurface Scattering Materials
最后我们使用SubsurfaceMaterial材质来封装BSSRDF的计算：
class SubsurfaceMaterial : public Material {
public:
SubsurfaceMaterial(Float scale, const std::shared_ptr<Texture<Spectrum>> &Kr,
const std::shared_ptr<Texture<Spectrum>> &Kt,
const std::shared_ptr<Texture<Spectrum>> &sigma_a,
const std::shared_ptr<Texture<Spectrum>> &sigma_s,
Float g,
Float eta,
const std::shared_ptr<Texture<Float>> &uRoughness,
const std::shared_ptr<Texture<Float>> &vRoughness,
const std::shared_ptr<Texture<Float>> &bumpMap,
bool remapRoughness)
: scale(scale), Kr(Kr), Kt(Kt), sigma_a(sigma_a), sigma_s(sigma_s),
uRoughness(uRoughness), vRoughness(vRoughness), bumpMap(bumpMap), eta(eta),
remapRoughness(remapRoughness), table(100, 64) {
ComputeBeamDiffusionBSSRDF(g, eta, &table);
}
void ComputeScatteringFunctions(SurfaceInteraction *si, MemoryArena &arena,
TransportMode mode, bool allowMultipleLobes) const;
private:
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;
};
计算散射性质时可以调用SubsurfaceMaterial::ComputeScatteringFunctions()函数，它会根据材质的性质为交点添加BSDF以及BSSRDF：
void SubsurfaceMaterial::ComputeScatteringFunctions(
SurfaceInteraction *si, MemoryArena &arena, TransportMode mode,
bool allowMultipleLobes) const {
// Perform bump mapping with bumpMap, if present
if (bumpMap)
Bump(bumpMap, si);
// Initialize BSDF for SubsurfaceMaterial
Spectrum R = Kr->Evaluate(*si).Clamp();
Spectrum T = Kt->Evaluate(*si).Clamp();
Float urough = uRoughness->Evaluate(*si);
Float vrough = vRoughness->Evaluate(*si);
// Initialize bsdf for smooth or rough dielectric
si->bsdf = ARENA_ALLOC(arena, BSDF)(*si, eta);
if (R.IsBlack() && T.IsBlack())
return;
bool isSpecular = urough == 0 && vrough == 0;
if (isSpecular && allowMultipleLobes) {
si->bsdf->Add(ARENA_ALLOC(arena, FresnelSpecular)(R, T, 1.f, eta, mode));
} else {
if (remapRoughness) {
urough = TrowbridgeReitzDistribution::RoughnessToAlpha(urough);
vrough = TrowbridgeReitzDistribution::RoughnessToAlpha(vrough);
}
MicrofacetDistribution *distrib = isSpecular ? nullptr :
ARENA_ALLOC(arena, TrowbridgeReitzDistribution)(urough, vrough);
if (!R.IsBlack()) {
Fresnel *fresnel = ARENA_ALLOC(arena, FresnelDielectric)(1.f, eta);
if (isSpecular)
si->bsdf->Add(ARENA_ALLOC(arena, SpecularReflection)(R, fresnel));
else
si->bsdf->Add(ARENA_ALLOC(arena, MicrofacetReflection)(R, distrib, fresnel));
}
if (!T.IsBlack()) {
if (isSpecular)
si->bsdf->Add(ARENA_ALLOC(arena, SpecularTransmission)(T, 1.f, eta, mode));
else
si->bsdf->Add(ARENA_ALLOC(arena, MicrofacetTransmission)(T, distrib, 1.f, eta, mode));
}
}
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);
}