62 KiB
二 体积渲染实践
一、3D密度场的体积渲染
到目前为止，我们只学习了如何渲染均匀体积物体——即散射系数和吸收系数在空间中保持恒定的物体。这虽然可行，但略显单调，且与自然界中的实际情况并不相符。例如，观察蒸汽火车排出的烟柱，这类体积都是非均匀的，部分区域的不透明度高于其他区域。那么，我们该如何渲染非均匀体积物体呢？
在现实世界中，我们会说吸收系数或散射系数随空间变化：吸收系数越高，体积越不透明；且散射系数和吸收系数的空间变化可以相互独立。不过，我们通常会选择一种更实用的方法：为特定体积物体设定恒定的散射系数和吸收系数，然后使用 密度参数 来调节体积在空间中的外观。可以将密度参数（在编程中通常是 float 或 double 类型的实数）视为 随空间变化 的量，此时我们可以进行如下处理：
\sigma_s' = \sigma_s * density(p)\\ \sigma_t = (\sigma_s + \sigma_a) * density(p)
其中，\sigma_s' 是经空间变化的密度参数调制后的散射系数；比尔定律方程中使用的消光系数 $\sigma_t$（吸收系数与散射系数之和）也会被空间变化的密度参数调制；density(p) 是一个函数，返回空间中点 p 处的密度值。通过这种方式，我们可以生成诸如上图中火山烟柱之类的体积物体图像。
现在的问题是：如何生成这种密度场？主要有两种技术：
-
程序化生成：可以使用 3D 纹理（如佩林噪声函数）来程序化创建空间变化的密度场。
-
模拟生成：也可以使用流体模拟程序（模拟烟雾等流体的运动）来生成空间变化的密度场。
本章将采用第一种方法；下一章将学习如何渲染流体模拟的结果。
使用噪声函数生成密度场
图 1：我们可以使用佩林噪声程序化生成密度场。该噪声函数以点为参数，返回该点在 [-1,1] 范围内的噪声值。
本章将使用柏林噪声函数（Perlin noise function）程序化生成 3D 密度场。
什么是程序化噪声函数（此处特指佩林噪声函数）？从编程角度来说，它是一种在 3D 空间中程序化生成噪声图案的函数。我们可以利用这种图案生成密度随空间变化的密度场。该噪声函数接收一个空间点作为参数，返回该点处 3D 噪声纹理的值（通常是 float 或 double 类型的实数），其取值范围被限制在 [-1,1] 之间。由于密度值只能是 0（无体积）或正数，因此需要对噪声函数的输出值进行裁剪或重映射，以获得正的密度值。在以下代码片段中，我们将噪声值从 [-1,1] 重映射到 [0,1]：
float density = (noise(pSample) + 1) * 0.5;
其中，pSample 是光线步进穿过体积时，相机光线上采样点的位置。
关于噪声函数，我们将使用肯·佩林（Ken Perlin）本人提供的改进版佩林噪声实现。
在本章中，我们假设你已熟悉该函数；如果不熟悉也无需过度担心，你只需知道：向该函数传入需要评估的 3D 空间点位置，它就会返回该点处 [-1,1] 范围内的噪声值。以下是参考代码（完整实现可在源代码部分提供的文件中获取）：
int p[512]; // 置换表（见源代码）
double fade(double t) { return t * t * t * (t * (t * 6 - 15) + 10); }
double lerp(double t, double a, double b) { return a + t * (b - a); }
double grad(int hash, double x, double y, double z){
int h = hash & 15;
double u = h<8 ? x : y,
v = h<4 ? y : h==12||h==14 ? x : z;
return ((h&1) == 0 ? u : -u) + ((h&2) == 0 ? v : -v);
}
double noise(double x, double y, double z){
int X = (int)floor(x) & 255,
Y = (int)floor(y) & 255,
Z = (int)floor(z) & 255;
x -= floor(x);
y -= floor(y);
z -= floor(z);
double u = fade(x),
v = fade(y),
w = fade(z);
int A = p[X ]+Y, AA = p[A]+Z, AB = p[A+1]+Z,
B = p[X+1]+Y, BA = p[B]+Z, BB = p[B+1]+Z;
return lerp(w, lerp(v, lerp(u, grad(p[AA ], x , y , z ),
grad(p[BA ], x-1, y , z )),
lerp(u, grad(p[AB ], x , y-1, z ),
grad(p[BB ], x-1, y-1, z ))),
lerp(v, lerp(u, grad(p[AA+1], x , y , z-1 ),
grad(p[BA+1], x-1, y , z-1 )),
lerp(u, grad(p[AB+1], x , y-1, z-1 ),
grad(p[BB+1], x-1, y-1, z-1 ))));
}
最后，以下是该函数在光线步进程序中的使用方式：
Color integrate(const Ray& ray, ...){
float sigma_a = 0.1;
float sigma_s = 0.1;
float sigma_t = sigma_a + sigma_s;
...
float transmission = 1; // 初始为完全透明
for (size_t n = 0; n < numSteps; ++n) {
float t = tMin + stepSize * (n + 0.5);
Point p = ray.orig + ray.dir * t;
// 密度不再是恒定值，而是随空间变化
float density = (1 + noise(p)) / 2;
float sampleAtt = exp(-density * sigma_t * stepSize);
// 全局透射率被采样点不透明度衰减
transmission *= samplAtt;
...
}
}
图 2：密度沿光线方向变化。
图 3：密度场由噪声函数生成的非均匀体积渲染结果。
与上一章中渲染均匀体积物体的程序相比，你可以看到此处的修改相当简单：我们将密度变量的声明移到了光线步进循环内部，使其不再是恒定值，而是随空间变化的参数。图 2 直观展示了这一过程：我们沿光线步进时对密度场进行采样（密度场仍由噪声函数生成），对于光线上的每个采样点，都以该采样点的位置作为输入参数评估噪声函数，并将结果作为该点的密度值。图 3 展示了将该方法应用于体积球体的渲染结果。
为了确保你理解这一过程，我们绘制了噪声函数在某段距离上的取值，以及基于该噪声函数计算的透射率值；同时，我们还绘制了恒定密度体积（绿色曲线）与非均匀密度体积（红色曲线）的比尔-朗伯定律曲线对比。如图所示，绿色曲线（恒定密度）完全平滑，而红色曲线（非均匀密度）则不然。需要注意的是：当噪声函数返回值接近 0 时，透射率基本保持恒定；当噪声函数返回值较高时，透射率急剧下降（体积密度更大）。这些都是预期结果，但直观展示有助于理解。
float stepSize = 1. / 51.2;
float sigma_t = 0.9;
float t = 0;
float Thomogeneous = 1;
float Theterogeneous = 1;
for (int x = 0; x < 512; x++, t += stepSize) {
float noiseVal = powf((1 + noise(t, 0.625, 0)) / 2.f, 2);
float samplAttHeterogeneous = exp(-noiseVal * stepSize * sigma_t);
Theterogeneous *= samplAttHeterogeneous;
float sampleAttHomogeneous = exp(-0.5 * stepSize * sigma_t);
Thomogeneous *= sampleAttHomogeneous;
fprintf(stderr, "%f %f %f\n", t, Theterogeneous, Thomogeneous);
}
但存在一个问题：这段代码可以正确计算体积的透射率值，但我们之前用于计算内散射贡献（还记得 Li 项吗？）的代码却无法直接使用。接下来，我们将解释原因，以及需要进行哪些修改才能使其适用于非均匀参与介质。
非均匀参与介质的内散射
图 4：光线穿过非均匀参与介质的示意图。
通过观察图 4，你应该能理解问题所在。对于均匀体积物体，计算光线贡献时，只需找到采样点到物体边界沿光线方向的距离（记为 Dl），然后结合该距离和体积消光系数（sigma_t = sigma_a + sigma_s）应用比尔定律，即可得到光线穿过体积到达采样点后的剩余光量（透射光量）。这很简单：
Color lighRayContrib = exp(-Dl * sigma_t * density) * lightColor;
非均匀体积的问题在于，这种方法不再适用——因为密度沿光线方向也会变化（如图 4 所示）。需要注意的是：我们此处要解决的问题，与第 2 章中使用正向光线步进解决的问题并不相同。到目前为止，我们使用光线步进的目的是估计相机光线上的内散射项；而此处，光线步进技术将再次发挥作用，用于估计内散射项以及相机光线和光线在非均匀参与介质中的 透射率。尽管问题不同（估计内散射项 vs 估计光线透射率），但所使用的技术相同（此处特指正向光线步进，一种随机采样方法）。我们需要 将光线划分为一系列小线段，假设每个线段长度（由步长定义的小体积元）内的密度是均匀的（微元思想），然后沿光线推进时，将全局透射率值与每个采样点的透射率相乘。伪代码实现如下：
// 计算非均匀介质中的光线透射率
float transmission = 1;
float stepSize = Dl / numSteps;
for (n = 0; n < numSteps; ++n) {
float t = stepSize * (n + 0.5);
float sampleAtt = exp(-evalDensity(t) * stepSize * sigma_t);
transmission *= samplAtt;
}
将这段代码与上一章中沿光线从 t0 到 t1 推进时计算相机光线透射率的代码进行对比：
float sigma_t = sigma_a + sigma_s;
float density = 0.1; // 密度恒定，用于缩放sigma_t
float transparency = 1; // 初始化透明度为1
for (int n = 0; n < ns; ++n) {
float t = isect.t1 - step_size * (n + 0.5);
vec3 sample_pos= ray_orig + t * ray_dir; // 采样点位置（步长中点）
// 使用比尔定律计算采样点透明度
float sample_transparency = exp(-step_size * sigma_t * density);
// 用采样点透明度衰减全局透明度
transparency *= sample_transparency;
// 内散射计算
if (hitObject->intersect(sample_pos, light_dir, isect_vol) && isect_vol.inside) {
...
result += ...
}
// 最终用采样点透明度衰减结果
result *= sample_transparency;
}
这两段代码的核心逻辑相同。我们可以利用一个数学技巧（现在正是使用它的好时机）：
e^a * e^b = e^{a + b}
观察代码可知，透射率值本质上是一系列指数函数的乘积。如果展开光线步进循环（代码片段 2），会得到如下形式：
// dx = 步长，noise(x)的取值范围为[0,1]
float t0 = dx * (0.5); // n = 0
float t1 = dx * (1 + 0.5); // n = 1
float t2 = dx * (2 + 0.5); // n = 2
...
float transmission = exp(-dx * sigma_t * noise(t0)) *
exp(-dx * sigma_t * noise(t1)) *
exp(-dx * sigma_t * noise(t2)) *
...;
利用我们刚刚学到的指数函数数学性质，可将代码重写为：
float tau = noise(t0) + noise(t1) + noise(t2) + ...;
float transmission = exp(-tau * sigma_t * dx);
换句话说，沿光线推进时（与沿相机光线推进完全相同），我们只需累加光线上每个采样点的密度值，然后利用该总和通过一次指数函数调用计算光线的衰减/透射率值（这确实能节省一些时间）。下图直观展示了这一概念。
从技术上讲，我们也可以对相机光线的透射率进行同样的优化，但请注意，在以下代码中，我们在沿相机光线推进时使用透射率值来衰减 Li 项。由于需要体积推进过程中的透射率中间值，因此我们不能像计算光线透射率那样，将密度值累加后在最后统一计算最终的光线透射率值。
float transmission = 1; // set the camera ray transmission value (full transmission)
vec3 result = 0; // the camera ray radiance (light energy traveling from the volume to the eye)
for (n = 0; n < numSteps; ++n) {
float t = t0 + stepSize * (n + 0.5);
vec3 samplePos = ray.orig + t * ray.dir;
float sampleDensity = evalDensity(samplePos);
// we need this intermediate result to attenuate the Li term (see below)
transmission *= exp(-sampleDensity * sigma_t * stepSize);
// inscattering (Li(x))
if (density > 0 && hit_object->intersect(...) {
float tau = 0;
for (nl = 0; nl < numStepsLight; ++nl) {
float tLight = stepSize * (nl + 0.5);
vec samplePosLight = samplePos + tLight * lightDirection;
tau += evalDensity(samplePosLight);
}
// calculate light ray transmission value at the very end
float lightRayAtt = exp(-tau * sigma_t * stepSize);
result += lightColor * lightRayAtt * sigma_s * ... * transmission;
}
}
“tau” 这个名称并非随意选择——在学术文献中，它常被用来表示“光学深度”（optical depth）这一物理量。有两个希腊字母常被用于表示该物理量：tau（$\tau$）或 rho（$\rho$）。本章不会给出光学深度的正式定义（以免造成混淆），相关定义将在《体积渲染：总结、方程/理论》章节中详细说明。
就是这样！现在你已经掌握了渲染非均匀体积物体准确图像所需的全部知识。最后这幅图展示了特定噪声分布对应的透射率曲线。
实用（且可运行）的实现
让我们对程序进行必要的调整，以演示该方法的实际应用。请记住，算法现在的工作流程如下：
-
沿相机光线步进：在相机光线上的每个采样点，估计该点的密度以计算采样点透射率，并估计内散射贡献。
-
沿光线步进：此前我们只需沿相机光线步进，而现在还需要沿光线步进。对于均匀体积物体，我们直接沿相机光线步进即可估计内散射项；对于非均匀体积物体，沿相机光线步进的目的仍为估计内散射项（这一点没有变化），但还需要评估沿相机光线的密度项；同时，我们还需要沿光线步进，以评估这些光线上的密度函数。
换句话说，我们现在需要同时沿相机光线和光线步进，并在这些光线上的每个采样点评估密度函数（当前为柏林噪声函数）。这涉及大量运算，你会发现：将球体渲染为非均匀介质的耗时，远高于其均匀介质版本。
// [注释]
// 该函数由integrate函数调用，用于评估非均匀体积球体在采样点p处的密度
// 返回该3D位置处佩林噪声函数的值（已重映射到[0,1]范围）
// [/注释]
float eval_density(const vec3& p){
float freq = 1;
return (1 + noise(p.x * freq, p.y * freq, p.z * freq)) * 0.5;
}
vec3 integrate(
const vec3& ray_orig,
const vec3& ray_dir,
const std::vector<std::unique_ptr<Sphere>>& spheres){
...
const float step_size = 0.1;
float sigma_a = 0.5; // 吸收系数
float sigma_s = 0.5; // 散射系数
float sigma_t = sigma_a + sigma_s; // 消光系数
float g = 0; // 亨耶-格林斯坦非对称因子
uint8_t d = 2; // 俄罗斯轮盘赌“概率”参数
int ns = std::ceil((isect.t1 - isect.t0) / step_size);
float stride = (isect.t1 - isect.t0) / ns;
vec3 light_dir{ -0.315798, 0.719361, 0.618702 };
vec3 light_color{ 20, 20, 20 };
float transparency = 1; // 初始化透射率为1（完全透明）
vec3 result{ 0 }; // 初始化体积球体颜色为0
// 主光线步进循环（正向，从t0推进到t1）
for (int n = 0; n < ns; ++n) {
// 采样点位置抖动
float t = isect.t0 + stride * (n + distribution(generator));
vec3 sample_pos = ray_orig + t * ray_dir;
// [注释]
// 评估采样点位置的密度（空间变化的密度）
// [/注释]
float density = eval_density(sample_pos);
float sample_attenuation = exp(-step_size * density * sigma_t);
transparency *= sample_attenuation;
// 内散射计算
IsectData isect_light_ray;
if (density > 0 &&
hit_sphere->intersect(sample_pos, light_dir, isect_light_ray) &&
isect_light_ray.inside) {
size_t num_steps_light = std::ceil(isect_light_ray.t1 / step_size);
float stide_light = isect_light_ray.t1 / num_steps_light;
float tau = 0;
// [注释]
// 沿光线步进：将密度值累加到tau变量中
// [/注释]
for (size_t nl = 0; nl < num_steps_light; ++nl) {
float t_light = stide_light * (nl + 0.5);
vec3 light_sample_pos = sample_pos + light_dir * t_light;
tau += eval_density(light_sample_pos);
}
float light_ray_att = exp(-tau * stide_light [... 655 chars omitted ...]
在该实现中，相机光线和光线使用的步长相同——但这并非必须。为了提高速度，你可以使用更大的步长来估计光线的透射率值。此外，由于我们使用了程序化纹理，可能会遇到滤波问题：如果采样噪声函数的频率过低，可能会丢失程序化纹理的细节，最终导致走样问题。同样，这是一个滤波问题，我们此处不深入探讨，但需要注意：步长、噪声频率和图像分辨率在采样层面是相互关联的。
你可能会感到惊讶：程序的输出结果（右图）看起来并不太像云层。然而，从噪声图案图像（左图）可以看出，默认情况下噪声图案的“块状结构”相当平滑，因此体积的外观也较为平滑。要生成类似云层的程序化噪声，需要通过多种方式调整噪声函数的结果，以获得视觉上更有趣的效果（如下文所示）。
以下是一些你可以立即尝试的简单变体：去除噪声函数的负值（左图），以及取噪声函数的绝对值（右图）。“衰减参数”（falloff parameter）的说明见下文。
尽管效果不够惊艳，但正是这项技术被用于制作 1997 年电影《超时空接触》（Contact）开场片段中的壮观体积特效。这些图像由索尼图形图像运作公司（Sony Picture Imageworks）使用皮克斯（Pixar）的渲染器渲染生成。在当时，这个片段是一项巨大的技术挑战。如前所述，他们使用的技术与本章介绍的方法相同——不同之处在于，他们使用了一些几何图形（而非基础球体）来定义体积物体的形状，并使用分形图案赋予星云类似云层的纹理。我们将在本章的最后一节（生产级体积渲染）中触及前一种技术；至于类似云层的纹理，让我们看看可以如何实现……
调整密度函数（创建更有趣的外观和动画）
编写光线步进器是一回事，创建类似云层的程序化纹理则是另一回事：前者是科学，后者更偏向艺术——利用一系列数学工具塑造程序化纹理，并花费大量时间调整参数，直到获得满意的结果。我们的目标并非创建令人惊叹的逼真图像，而是为你提供工具和基础：首先帮助你理解原理，其次让你能够将这些工具重组为更复杂的系统（如果你愿意的话），从而制作出“惊艳”的图像。不过，这部分工作就交给你了……如果你创作了有趣的作品，记得和我们分享！
以下是一些可用于创建更类似云层的“球体”的技巧——这只是众多工具中的一小部分。
平滑步进函数（Smoothstep）
平滑步进函数是你已经熟悉的函数，我们已在多个场景中使用过它（包括噪声函数）。它能在两个值之间创建“平滑”的过渡。以下是该函数的一种可能实现：
float smoothstep(float lo, float hi, float x){
float t = std::clamp((x - lo) / (hi - lo), 0.f, 1.f);
return t * t * (3.0 - (2.0 * t));
}
我们可以使用该函数在球体边界附近创建衰减效果。为此，需要修改 eval_density 函数，向其传入球体中心和半径作为参数。通过这种方式，我们可以计算采样点到球体中心的归一化距离，并利用该值调整球体的密度，如下所示：
float eval_density(const vec3& sample_pos, const vec3& sphere_center, const float& sphere_radius){
vec3 vp = sample_pos - sphere_center;
float dist = std::min(1.f, vp.length() / sphere_radius);
float falloff = smoothstep(0.8, 1, dist); // 当距离从0.8到1时，平滑过渡从0到1
return (1 - falloff);
}
这项技术适用于在体积到达球体边界前使其渐隐。我们使用平滑步进函数控制衰减开始的位置（距离边缘的距离）——对于衰减效果，过渡的终点应设为 1。
分形布朗运动（fBm）
fBm（全称 Fractal Brownian Motion，在早期也被称为“等离子体纹理”或“混沌纹理”）是一种分形图案，由多个程序化噪声层叠加而成，不同层的频率和振幅各不相同。你可以在“程序化生成”章节的课程中找到关于该图案的更多信息。
典型的 fBm 程序化纹理代码可按如下方式构建：
float eval_density(const vec3& p, ...){
vec3 vp = p - sphere_center;
...
// 构建fBm分形图案
float frequency = 1;
vp *= frequency; // 必要时缩放初始点坐标
size_t numOctaves = 5; // 层数
float lacunarity = 2.f; // 连续频率之间的间隔
float H = 0.4; // 分形增量参数
float value = 0; // fBm的结果（用于密度计算）
for (size_t i = 0; i < numOctaves; ++i) {
value += noise(vp) * powf(lacunarity, -H * i);
vp *= lacunarity;
}
// 裁剪负值
return std::max(0.f, value) * (1 - falloff);
}
构建 fBm 图案只需两行核心代码（循环叠加层）。基于这个基础版本，可以衍生出许多变体——例如取噪声函数的绝对值（一种称为“湍流”的图案）等。如果你想了解更多关于该主题的内容（以及如何正确滤波该图案），请查看“程序化生成”章节。
偏移（Bias）
偏移用于改变中点（0.5）的位置，同时保持 0 和 1 的值不变。
float eval_density(...) {
float bias = 0.2;
float exponent = (bias - 1.0) / (-bias - 1.0);
// 假设指数大于0
return powf(noisePattern, exponent);
}
利用图元局部坐标系旋转噪声图案
在早期（随着计算能力的提升，流体模拟成为主流之前），另一种常用技术是在体积图元（球体、立方体等）内部动画化噪声图案。我们传入密度函数的采样点位置是在世界空间中定义的，但我们也可以在依附于球体的坐标系中定义该点位置——换句话说，将采样点从世界空间转换到物体空间。要在球体坐标系中定义采样点，只需执行以下操作：
vec3 sample_object_space = sample_pos - center;
更通用的解决方案是使用矩阵变换球体图元：通过物体到世界矩阵的逆矩阵，将采样点从世界空间转换到球体的局部坐标系。不过，为了简化示例程序，我们并未使用矩阵——但你可以利用 Scratchapixel 上提供的信息轻松实现这一点。
这项技术的实用价值在于：当我们移动球体时，在球体局部坐标系中定义的点坐标不会改变。这可以确保噪声图案始终“附着”在球体上，无论球体如何变换——类似你在手指间滚动的玻璃弹珠上的纹理。在我们选择的示例中（尽管与我们刚才解释的略有不同，但概念和结果一致）：我们将在物体空间中围绕球体的 y 轴旋转采样点（更好的解决方案是使用物体到世界矩阵旋转球体，然后通过世界到物体矩阵将采样点从世界空间转换到物体空间，但我们在此处为了偷懒，直接在物体空间中旋转点）。这样，我们可以看到分形图案也随之旋转，有点像观察旋转的玻璃弹珠。通过这种方法，还可以应用更复杂的（线性或非线性）变形。
float eval_density(const vec3& p, const vec3& center, const float& radius){
// 将点从世界空间转换到物体空间
vec3 vp = p - center;
vec3 vp_xform;
// 在物体空间中旋转采样点（frame是全局变量，范围从1到120）
float theta = (frame - 1) / 120.f * 2 * M_PI;
vp_xform.x = cos(theta) * vp.x + sin(theta) * vp.z;
vp_xform.y = vp.y;
vp_xform.z = -sin(theta) * vp.x + cos(theta) * vp.z;
float dist = std::min(1.f, vp.length() / radius);
float falloff = smoothstep(0.8, 1, dist);
float freq = 0.5;
size_t octaves = 5;
float lacunarity = 2;
float H = 0.4;
vp_xform *= freq;
float fbmResult = 0;
float offset = 0.75;
for (size_t k = 0; k < octaves; k++) {
fbmResult += noise(vp_xform.x , vp_xform.y, vp_xform.z) * pow(lacunarity, -H * k);
vp_xform *= lacunarity;
}
return std::max(0.f, fbmResult) * (1 - falloff);
}
注意：现在我们可以清晰地看到分形图案的 3D 结构。
还有更多技巧……
“程序化”技巧的列表可以不断延伸：位移（使用噪声函数位移体积的边缘）、不同类型的噪声（波状噪声、时空噪声——一种可动画的噪声等）。我们将在本课程的未来修订版中扩展这个列表。
我们学到了什么？
关于本课程所学技术的总结，见最后一章。
到目前为止，你已经可以使用光线步进渲染单散射非均匀体积物体——这已经是一项重大成就。此外，你现在应该清楚：光线步进技术的核心是将体积物体分解为步长定义的小体积元，该技术之所以有效，是因为我们假设这些微小的体积元（或采样点）足够小，可以视为均匀的。简而言之，你将非均匀物体分解为多个“小砖块”，每个砖块都可视为“均匀的”——尽管不同砖块的属性可能不同，有点像用乐高积木搭建物体。
二、基于 3D 体素网格的体渲染详解
学完本章后，你将能够生成以下序列图像：
如前一章所述，创建密度场有两种技术：过程式生成或使用流体模拟软件。本章我们将聚焦后者。请注意，本课程不会讲解流体模拟的工作原理，而是学习如何渲染流体模拟生成的数据——流体模拟的相关内容我们会在后续课程中介绍，敬请期待。
步骤 1：使用 3D 网格存储密度值
模拟流体（烟雾、水等）的技术有很多，但通常在模拟过程中的某个阶段，结果会被存储在 3D 网格中。为简化教学，本章假设这些网格在三个维度上的分辨率相同（例如 x、y、z 坐标均为 32×32×32），且分辨率为 2 的幂（8、16、32、64、128 等）——这仅为教学方便，实际应用中无需严格遵循。此外，本章暂不涉及网格变换，因此假设网格是轴对齐的立方体：我们可将网格视为轴对齐包围盒（AABB），这会简化光线-立方体相交测试。若要支持非轴对齐网格，只需按照本课程讲解的方法，通过网格的世界空间到对象空间的变换矩阵，将相机光线转换到网格的对象空间即可。
图 1：存储密度值的 8×8×8 网格
网格非常适合模拟流体运动：网格的体素（构成网格的小体积元素，也可称为“单元格”）会被赋予初始密度（可理解为填充烟雾），随着时间推移，密度会在相邻体素间传播。密度在体素间的移动遵循 纳维-斯托克斯方程（Navier-Stokes equation）——这一主题同样留到后续课程讲解。本章你只需了解：流体模拟的结果存储在由体素组成的 3D 网格中，每个体素存储一个密度值（0 或大于 0 的数值），代表“填充”该体素体积的密度。3D 网格之于流体模拟，就如同位图之于图像，密度场被量化处理。 在任何编程语言中，定义存储密度值的 3D 网格都相对简单，但本章末尾会探讨一些需要注意的细节。
除了网格分辨率（每个维度的体素数量，如 32×32×32），我们还需定义网格在世界空间中的尺寸（即场景中对象的大小）。实现方式有多种，为方便起见，我们通过定义网格在世界空间中的最小边界和最大边界（或范围）来确定，例如 (-2,-2,-2) 和 (2,2,2)。需注意，由于目前我们的网格是立方体，边界值不必相对于世界原点对称，但三个维度的世界空间尺寸必须相同：例如 (0.2,0.2,0.2) 与 (10.2,10.2,10.2) 是合法的，而 (-1.2,2.2,3.2) 与 (8.2,4.2,7.2) 则不合法。这种定义方式的优势在于，最小和最大边界可直接代入光线-立方体相交检测程序（我们在《极简光线追踪器：渲染简单图形（球体、立方体、圆盘、平面等）》课程中已学习过）。因此，我们已经掌握了计算任何与立方体相交的光线的 t0 和 t1 值的方法（与球体基元的计算方式类似）。
你可以修改前几章的代码，将渲染球体改为渲染立方体，同时使用矩阵从更有趣的视角渲染场景。最终应得到类似以下的代码：
struct Grid {
float *density;
size_t resolution = 128; // 分辨率
Point bounds[2] = { (-30,-30,-30), (30, 30, 30) }; // 世界空间边界
};
// 光线-立方体相交检测
bool rayBoxIntersect(const Ray& ray, const Point bounds[2], float& t0, float& t1){
...
return true;
}
// 积分计算（光线步进核心）
void integrate(const Ray& ray, const Grid& grid, Color& radiance, Color& transmission){
Vector lightDir(-0.315798, 0.719361, 0.618702); // 光源方向
Color lightColor(20); // 光源颜色
float t0, t1;
if (rayBoxIntersect(ray, grid.bounds, t0, t1)) { // 光线与网格包围盒相交
// 从t0到t1执行光线步进
radiance = ...; // 辐射亮度
transmission = ...; // 透射率
}
}
// 相机到世界空间的变换矩阵
Matrix cameraToWorld{ 0.844328, 0, -0.535827, 0, -0.170907, 0.947768, -0.269306, 0, 0.50784, 0.318959, 0.800227, 0, 83.292171, 45.137326, 126.430772, 1 };
int main(){
Grid grid;
// 分配内存以读取缓存文件中的数据
size_t numVoxels = grid.resolution * grid.resolution * grid.resolution; // 体素总数
grid.density = new float[numVoxels]; // 也可使用unique_ptr管理内存（示例程序采用此方式）
std::ifstream cacheFile;
cacheFile.open("./cache.0100.bin", std::ios::binary); // 打开二进制缓存文件
// 读取密度值到内存
cacheFile.read((char*)grid.density, sizeof(float) * numVoxels);
Point origin(0); // 光线原点
for (帧中的每个像素) {
Vector rayDir = ...; // 计算光线方向
Ray ray;
ray.orig = origin * cameraToWorld; // 转换原点到世界空间
ray.dir = rayDir * cameraToWorld; // 转换方向到世界空间
Color radiance = 0, transmission = 1;
integrate(ray, grid, radiance, transmission); // 执行光线步进积分
pixelColor = radiance; // 像素颜色 = 辐射亮度
pixelOpacity = 1 - transmission; // 像素不透明度 = 1 - 透射率
}
// 释放内存
delete[] grid.density;
}
目前无需过多关注缓存文件的内容（后续会详细说明）。从上述代码可以看出，核心操作是：创建连续的内存块存储与体素数量相等的浮点型密度值，然后将磁盘文件中的数据加载到内存中（第 36 行），逻辑十分简洁。
预期效果
这是你应该得到的图像。现在我们已经明确：将使用 3D 网格存储密度场数据，网格具有固定分辨率和尺寸，且已掌握光线与网格包围盒的相交点计算方法。接下来，我们将学习如何在光线步进过程中读取网格的密度值。
步骤 2：在网格中执行光线步进
如《光线追踪加速结构》课程所述，使用网格作为加速结构较为复杂：需要找到光线相交的所有网格单元格，再判断每个单元格中的几何体（三角形）是否与光线相交。虽然不算特别复杂，但也绝非易事（尤其是追求高效时）。不过好消息是，体渲染中的光线步进完全不需要这一过程！
光线步进的核心是关注光线上的采样点（如下图所示）。由于我们知道这些采样点的位置，只需确定它们所在的网格体素，然后读取该体素的密度值即可——无需遍历网格，极其简单。接下来，我们将学习如何从采样点映射到体素坐标，并通过坐标从内存中读取密度值。
步骤 3：光线步进过程中读取密度值
与前几章的代码相比，我们只需修改 evalDensity 函数：不再通过过程式噪声函数计算空间中任意点的密度，而是从网格中读取对应位置的密度值——这是一个极小的改动，其余逻辑完全不变。
通常情况下，采样点必然位于网格边界内（因为采样范围在光线与立方体的相交点之间），因此无需额外判断。但为了代码健壮性，建议添加采样点坐标与网格边界的校验，避免内存访问越界。
从网格中读取密度值的操作通常称为“查找（lookup）”。在 OpenVDB（一款高效存储和读取体数据的库）中，这一操作被称为“探测（probing）”（较为少见）。具体实现步骤如下：
采样点到体素坐标的转换流程
采样点最初定义在世界空间中，需将其转换为网格的离散坐标，步骤如下：
-
世界空间到对象空间转换：将采样点坐标减去网格的最小边界值，得到相对于网格“左下角”的坐标。
-
归一化：将对象空间坐标除以网格尺寸（示例中为 10），得到归一化坐标（范围在 [0,1] 之间）。
-
体素空间映射：将归一化坐标乘以网格分辨率，得到体素空间坐标（仍为浮点型）。为了读取内存中的密度值，需将浮点型坐标四舍五入为整数（即体素索引）。需注意：体素空间坐标的最大值不能超过“分辨率-1”（避免归一化坐标恰好为 1 时索引越界）。完成以上步骤后，即可读取对应位置的密度值。
2D 示例说明
如上图的 2D 示例所示，采样点的体素空间坐标为 (3.36, 3.28)。此时，该体素在内存中的索引计算方式为：y 坐标的整数部分（3）乘以网格分辨率（8），再加上 x 坐标的整数部分（3），即 3×8+3=27。
代码实现（最近邻插值）
float evalDensity(const Grid* grid, const Point& p){
Vector gridSize = grid->bounds[1] - grid->bounds[0]; // 计算网格尺寸
Vector pLocal = (p - grid->bounds[0]) / gridSize; // 世界空间 → 对象空间（归一化前）
Vector pVoxel = pLocal * grid->resolution; // 归一化 → 体素空间
// 取浮点坐标的整数部分（使用floor避免负索引）
int xi = static_cast<int>(std::floor(pVoxel.x));
int yi = static_cast<int>(std::floor(pVoxel.y));
int zi = static_cast<int>(std::floor(pVoxel.z));
// 最近邻插值：直接读取对应体素的密度值
return grid->density[(zi * grid->resolution + yi) * grid->resolution + xi];
}
关键说明
与 Perlin 噪声中将点转换到晶格空间的逻辑类似，此处使用 std::floor 函数而非直接强制类型转换，是为了避免归一化坐标略小于 0 时（后续会说明原因），索引变为 -1 导致内存访问错误—— std::floor 会确保负坐标的索引为 0，从而返回 0 密度值。
这种方法称为“最近邻插值（nearest-neighbor interpolation）”：本质是直接获取采样点所在体素的密度值，并未进行插值计算。后续我们会学习其他插值方法，但所有方法的第一步都是计算采样点所在的体素坐标。
现在，你已经具备了渲染流体模拟结果的全部条件（源代码区提供了多个流体缓存文件，可与示例程序配合使用）。以下是最近邻插值的渲染效果：
缓存文件格式说明
为了避免依赖外部库，我们没有使用 OpenVDB（或其前身 Field3D）等工业级格式存储网格数据，而是将流体模拟结果以最简洁的二进制格式存储，代码如下：
// 二进制缓存文件的写入逻辑
for (size_t z = 0; z < grid.res; ++z) {
for (size_t y = 0; y < grid.res; ++y) {
for (size_t x = 0; x < grid.res; ++x) {
float density = grid.density[x][y][z];
cacheFile.write((char*)&density, sizeof(float));
}
}
}
OpenVDB 的核心原理与此类似，但支持压缩、多分辨率和稀疏体（本章末尾会简要介绍）。我们的方案优势在于简洁直观，非常适合教学场景。
需特别注意：上述代码中，体素沿 z 轴“从后到前”存储（切片方式）。在右手坐标系中，网格的最小边界位于“后方左下角”（下图中的品红色体素），体素空间中第一个体素的坐标为 (0,0,0)；在 8×8×8 的网格示例中，最后一个体素（z 轴正方向前方右上角）的坐标为 (7,7,7)（下图中的黄色体素），索引为 8×8×8-1=511（0-based 索引）。
体素索引计算
若要通过体素空间坐标（xi, yi, zi）计算其在内存数组中的索引，公式如下：
size_t index = (zi * gridResolution + yi) * gridResolution + xi;
这是最基础的实现方式，但图像质量可通过“三线性插值”进一步提升。
步骤 4：通过三线性插值优化结果
其核心思路如下：一个体素代表单个密度值，但采样点可能位于任意给定体素体积内的任意位置。因此，合理的假设是该位置的密度值应在某种程度上混合采样点所在体素的密度，以及与之直接相邻的体素的密度值（如下图所示）。
对于三线性插值，我们只需混合 8 个体素的值。具体该如何混合？方法很简单：计算采样点到这 8 个体素中每一个体素的距离，并利用这些距离来加权每个体素对最终结果的贡献度。
这个过程本质上是三维空间中的线性插值。我们用于插值两个值的公式为 $a * (1 - weight) + b * weight$，其中
weight的取值范围为 0 到 1。在二维空间中，该过程称为双线性插值，需要 4 个像素；在三维空间中，我们称之为三线性插值，需要 8 个体素。请注意，用于在 a 和 b 之间插值的公式被称为插值函数（interpolant）。在线性插值的情况下，插值函数是一阶多项式。
我们假设存储在某个体素中的密度值，是当采样点恰好位于该体素正中心时应取的值。由于三线性插值方案的特性，要得到这一结果，我们需要将体素空间中的采样位置偏移 -0.5。我们先看代码，再解释其工作原理。
float evalDensity(const Grid* grid, const Point& p){
Vector gridSize = grid->bounds[1] - grid->bounds[0];
Vector pLocal = (p - grid->bounds[0]) / gridSize;
Vector pVoxel = pLocal * grid->baseResolution;
Vector pLattice(pVoxel.x - 0.5, pVoxel.y - 0.5, pVoxel.z - 0.5); // 红色高亮行
int xi = static_cast<int>(std::floor(pLattice.x));
int yi = static_cast<int>(std::floor(pLattice.y));
int zi = static_cast<int>(std::floor(pLattice.z));
float weight[3];
float value = 0;
// 三线性插值
for (int i = 0; i < 2; ++i) {
weight[0] = 1 - std::abs(pLattice.x - (xi + i));
for (int j = 0; j < 2; ++j) {
weight[1] = 1 - std::abs(pLattice.y - (yi + j));
for (int k = 0; k < 2; ++k) {
weight[2] = 1 - std::abs(pLattice.z - (zi + k));
value += weight[0] * weight[1] * weight[2] * grid->operator()(xi + i, yi + j, zi + k);
}
}
}
return value;
}
接下来我们分析这段代码的工作方式及背后的原因。
在上图中，展示了两种插值技术。左上角是最近邻插值法的示例，原理非常简单：我们有 4 个体素（二维场景下，三维则为 8 个），采样点（蓝色圆点）落在其中一个体素内，因此采样点直接采用该体素中存储的密度值。
三线性插值则稍复杂一些。你可能会认同：如果采样点恰好位于四个体素的交叉中心，那么我们应该返回这四个体素中存储的密度值的平均值。也就是说，在这个特定示例中，体素密度值的总和需除以 4：
result = (0.9 + 0.14 + 0.08 + 0.63) / 4
但这仅适用于采样点恰好位于四个体素交叉点的特殊情况，因此我们需要一个通用的解决方案。
首先（稍后解释为何需要这么做），我们需要将采样点在体素空间中偏移半个体素的距离：(-0.5, -0.5, -0.5)。以采样点恰好位于四个体素交叉中心的示例来说，应用偏移后，采样点会落在左下角体素的正中心。接下来，我们计算该点到每个体素边界的距离。在我们的二维示例图中，这一距离表现为 dx0（采样点 x 坐标到左下角体素 x 坐标 x0 的距离）、dy0（采样点 y 坐标到左下角体素 y 坐标 y0 的距离）、dx1（采样点 x 坐标到右上角体素 x 坐标 x1 的距离）和 dy1（采样点 y 坐标到右上角体素 y 坐标 y1 的距离）。这类距离在技术上被称为曼哈顿距离。
曼哈顿距离：沿直角坐标轴测量的两点间距离。在平面中，若 p1 坐标为 (x1, y1)，p2 坐标为 (x2, y2)，则曼哈顿距离为 $|x1 - x2| + |y1 - y2|$。
实际应用中，我们只需用到 dx0、dy0 和 dz0 即可（如下文所示）。以下再次梳理三线性插值的核心思路（若想了解更完整的解释，可参考专门讲解插值的章节）。
我们有一个 2×2×2 的体素块，将其命名为 v000、v100（向右移动 1 个体素）、v010（向上移动 1 个体素）、v110（向右 1 个、向上 1 个）、v001（向前移动 1 个体素）、v101（以此类推）、v011 和 v111。插值的思路是：先以 x0 为权重，对 v000-v100、v010-v110、v001-v101、v011-v111 进行线性插值，得到 4 个值；再以 y0 为权重，将这 4 个值两两线性插值，得到 2 个值；最后以 z0 为权重，对这 2 个值进行线性插值，得到最终结果。请注意，无论先以 x0、y0 还是 z0 开始前 4 次线性插值，只要后续每次插值都使用不同维度的权重，结果都不会受影响。记住，我们的线性插值函数为（原文公式缺失），其中（参数）为权重。用代码表示如下：
float result =
(1 - z0) * ( // 蓝色部分
(1 - y0) * (v000 * (1 - x0) + v100 * x0) + // 绿色和红色部分
y0 * (v010 * (1 - x0) + v110 * x0) // 绿色和红色部分
) +
z0 * ( // 蓝色部分
(1 - y0) * (v001 * (1 - x0) + v101 * x0) + // 绿色和红色部分
y0 * (v011 * (1 - x0) + v111 * x0)); // 绿色和红色部分
希望这样的代码排版能帮助你理解三层嵌套的插值结构：包含 4 次（红色）+ 2 次（绿色）+ 1 次（蓝色）插值操作。若我们将 (1-x0)、(1-y0)、(1-z0)、x0、y0、z0 分别替换为 wx0、wy0、wz0、wx1、wy1、wz1，展开并重新整理上述代码片段（片段 2）中的项，可得到：
float result =
v000 * wx0 * wy0 * wz0 +
v100 * wx1 * wy0 * wz0 +
v010 * wx0 * wy1 * wz0 +
v110 * wx1 * wy1 * wz0 +
v001 * wx0 * wy0 * wz1 +
v101 * wx1 * wy0 * wz1 +
v011 * wx0 * wy1 * wz1 +
v111 * wx1 * wy1 * wz1;
这正是代码片段 1 中的逻辑（尽管片段 1 中权重是在嵌套循环中实时计算的）。若以二维示例图中 dx0 = dy0 = 0.5 的情况代入计算，可得：
result =
0.9 * (1-dx0) * (1-dy0) + 0.14 * dx0 * (1-dy0) + 0.08 * (1-dx0) * dy0 + 0.63 * dx0 * dy0 =
0.9 * 0.25 + 0.14 * 0.25 + 0.08 * 0.25 + 0.63 * 0.25
这正是我们期望的结果。由此可见，三线性插值（二维场景下为双线性插值）是有效的。你现在也能理解为何需要将采样点在体素空间中偏移 -0.5——这是让数学计算成立的必要条件。为了验证这一点，不妨考虑采样点恰好位于体素正中心的情况（偏移前）：应用偏移后，该点会落在体素的左下角（二维场景下），此时 dx0 = dy0（=dz0）=0，而其他所有距离均为 1。代入计算可得：
float result = grid->density(voxelX, voxelY, voxelZ) * (1-0) * (1-0) * (1-0) +
0 + // 0 * (1-0) * (1-0)
0 + // (1-0) * 0 * (1-0)
0 + // 0 * 0 * (1-0)
0 + // (1-0) * (1-0) * 0
0 + // 0 * (1-0) * 0
0 + // (1-0) * 0 * 0
0; // 0 * 0 * 0
换句话说，该采样点返回的值就是该体素中存储的值（若偏移前采样点位于左下角体素正中心，在我们的示例中即为 0.9），这完全符合预期。
以上就是三线性插值的全部原理。接下来我们对比其与最近邻插值法的结果差异。
可以看到，右侧使用三线性插值渲染的缓存图像更平滑。但这种优化是有代价的：使用该方法时，渲染时间会显著增加（相较于最近邻插值法）。在未做优化的情况下，三线性插值的耗时约为最近邻查找的 5 倍。
现在你可以阅读其他开发者的代码了——以 OpenVDB 为例：
template<class ValueT, class TreeT, size_t N> inline bool BoxSampler::probeValues(ValueT (&data)[N][N][N], const TreeT& inTree, Coord ijk) {
bool hasActiveValues = false;
hasActiveValues |= inTree.probeValue(ijk, data[0][0][0]); //i, j, k
ijk[2] += 1;
hasActiveValues |= inTree.probeValue(ijk, data[0][0][1]); //i, j, k + 1
ijk[1] += 1;
hasActiveValues |= inTree.probeValue(ijk, data[0][1][1]); //i, j+1, k + 1
ijk[2] -= 1;
hasActiveValues |= inTree.probeValue(ijk, data[0][1][0]); //i, j+1, k
ijk[0] += 1;
ijk[1] -= 1;
hasActiveValues |= inTree.probeValue(ijk, data[1][0][0]); //i+1, j, k
ijk[2] += 1;
hasActiveValues |= inTree.probeValue(ijk, data[1][0][1]); //i+1, j, k + 1
ijk[1] += 1;
hasActiveValues |= inTree.probeValue(ijk, data[1][1][1]); //i+1, j+1, k + 1
ijk[2] -= 1;
hasActiveValues |= inTree.probeValue(ijk, data[1][1][0]); //i+1, j+1, k
return hasActiveValues;
}
template<class ValueT, size_t N> inline ValueT BoxSampler::trilinearInterpolation(ValueT (&data)[N][N][N], const Vec3R& uvw) {
auto _interpolate = [](const ValueT& a, const ValueT& b, double weight)
{
const auto temp = (b - a) * weight;
return static_cast<ValueT>(a + ValueT(temp));
};
// 三线性插值：
// 使用周围8个晶格点的值构造最终结果。
// result(x,y,z) =
// v000 (1-x)(1-y)(1-z) + v001 (1-x)(1-y)z + v010 (1-x)y(1-z) + v011 (1-x)yz
// + v100 x(1-y)(1-z) + v101 x(1-y)z + v110 xy(1-z) + v111 xyz
return _interpolate(
_interpolate(
_interpolate(data[0][0][0], data[0][0][1], uvw[2]),
_interpolate(data[0][1][0], data[0][1][1], uvw[2]),
uvw[1]),
_interpolate(
_interpolate(data[1][0][0], data[1][0][1], uvw[2]),
_interpolate(data[1][1][0], data[1][1][1], uvw[2]),
uvw[1]),
uvw[0]);
}
template<class TreeT> inline bool [... 428 chars omitted ...]
sample() 方法首先通过 probeValues 方法获取 8 个体素中存储的值，随后在 trilinearInterpolation 方法中按上述逻辑进行插值。_interpolate 函数本质上就是线性插值（在代码中更常被称为 lerp）。
体数据缓存的进阶内容
我们在此列出一系列值得参考的主题。为控制本章篇幅，暂不展开详述。这些主题的大部分内容将在本章后续修订版本中补充，或在独立章节中讲解。
其他插值方案：三次插值与滤波核
如前所述，线性插值函数是一阶多项式。我们也可使用更高阶的多项式以获得更平滑的结果，例如三次插值。在二维场景下，双三次插值需要采样 4×4 个像素；在三维场景下，则需要 2^4 = 64 个体素。可想而知，结果会更平滑，但渲染时间也会相应增加。
你也可以使用滤波核（类似于图像滤波中使用的三角核、Mitchell 核或高斯滤波核）。与图像处理类似，滤波核的范围越大，处理耗时越长。
这两种技术将在本章后续修订版本中完整讲解。
密度之外的附加数据存储
生产环境中使用的体数据缓存格式通常支持在缓存中存储密度之外的任意通道数据，例如 温度（用于渲染火焰）、速度（用于渲染三维运动模糊）等。
滤波与砖块映射（Brick Map）
若你为生产环境渲染体数据缓存，那么滤波和 LOD（细节层次）是需要重点关注的问题。体数据缓存渲染面临的问题与纹理渲染类似：纹理有固定分辨率（如 512、1024 像素等），若从极远的距离渲染贴有该纹理的物体，物体在图像中仅占几个像素；当物体或相机移动时，每帧渲染都会用到纹理中不同的像素/纹素，导致纹理在帧与帧之间出现抖动（轻则如此，重则直接呈现噪点效果）。为缓解这一问题，我们会生成原始纹理的低分辨率版本并与原纹理数据一起存储，这些低分辨率版本被称为多级纹理（mipmap）。当物体距离较远（在图像中尺寸较小）时，使用低分辨率纹理可消除这种噪点（技术上称为走样/aliasing）。我们暂未深入讲解 mipmap，但很快会展开。
如果你已了解 mipmap 的概念、用途及生成方法，那么需要知道：二维纹理的滤波问题同样存在于三维缓存（或三维纹理）中。因此，我们可将 mipmap 方法适配到三维缓存中——为此需要生成原始网格数据的降采样版本，并根据缓存在图像中的尺寸选择合适的层级（降采样版本）。尽管这一过程与 mipmap 非常相似，且你也可将这些降采样版本称为 mipmap，但在三维场景中，我们通常称之为砖块映射（brick map），其形态类似《我的世界》（Minecraft）中的方块。
我们认为 “砖块映射（brick map）” 这一术语是由皮克斯 RenderMan 团队创造的，但我们对其起源并不十分确定…… 如果有人了解相关情况，欢迎来信告知我们。但总体而言，让缓存存储其数据的多分辨率版本这一理念并非特别新颖或罕见，几乎所有生产级渲染器和格式（如 OpenVDB）都可能支持这一理念。
砖块映射的生成过程与纹理 mipmap 层级的生成过程非常相似，但我们需要对 8 个体素取平均值（而非对 4 个像素/纹素取平均）。因此，分辨率为 2 的整数次幂的网格会更便于处理，因为这种情况下的降采样操作会非常简单。
以下代码展示了如何从原始流体模拟缓存生成这些层级（示例中基础分辨率为 128）。这只是一个简易示例，后续我们会专门讲解该主题：
size_t baseResolution = 128;
size_t numLevels = log2(baseResolution); /* 浮点型到size_t的隐式转换 */
std::unique_ptr<Grid []> gridLod = std::make_unique<Grid []>(numLevels - 2); // 忽略分辨率为2和4的层级
// 加载0级数据
gridLod[0].baseResolution = baseResolution;
std::ifstream ifs;
char filename[256];
sprintf_s(filename, "./grid.%d.bin", frame);
ifs.open(filename, std::ios::binary);
gridLod[0].densityData = std::make_unique<float[]>(baseResolution * baseResolution * baseResolution);
ifs.read((char*)gridLod[0].densityData.get(), sizeof(float) * baseResolution * baseResolution * baseResolution);
ifs.close();
for (size_t n = 1; n < numLevels - 2; ++n) {
baseResolution /= 2;
gridLod[n].baseResolution = baseResolution;
gridLod[n].densityData = std::make_unique<float[]>(baseResolution * baseResolution * baseResolution);
for (size_t x = 0; x < baseResolution; ++x) {
for (size_t y = 0; y < baseResolution; ++y) {
for (size_t z = 0; z < baseResolution; ++z) {
gridLod[n](x, y, z) =
0.125 * (gridLod[n - 1](x * 2, y * 2, z * 2) +
gridLod[n - 1](x * 2 + 1, y * 2, z * 2) +
gridLod[n - 1](x * 2, y * 2 + 1, z * 2) +
gridLod[n - 1](x * 2 + 1, y * 2 + 1, z * 2) +
gridLod[n - 1](x * 2, y * 2, z * 2 + 1) +
gridLod[n - 1](x * 2 + 1, y * 2, z * 2 + 1) +
gridLod[n - 1](x * 2, y * 2 + 1, z * 2 + 1) +
gridLod[n - 1](x * 2 + 1, y * 2 + 1, z * 2 + 1));
}
}
}
}
...
// 例如渲染缓存的3级数据（分辨率：16）
trace(ray, L, transmittance, rc, gridLod[3]);
...
我们暂不讲解如何选择合适的层级，但会在本章后续修订版本中补充。至少你现在已了解该问题及其可能的解决方案。
以下是缓存不同层级的渲染效果（0 级为原始缓存分辨率，示例中为 128）：
三维运动模糊
三维运动模糊该如何实现？要渲染这一效果，我们需要为体素存储流体的运动信息——通常用称为 运动矢量 的方向来表示。该矢量的方向代表流体在网格中移动的平均方向，其模长代表流体的移动速度。利用这些信息，可在渲染阶段模拟三维运动模糊。
该主题将在后续章节中讲解。
平流（Advection）
尽管平流与体渲染无直接关联（更多属于流体模拟范畴），但我们在此提及作为备忘（以便后续如有可能，专门撰写章节讲解）。平流可在渲染阶段执行，为现有流体模拟结果增加细节。
稀疏体数据（Sparse Volumes）：是什么？
多数情况下，网格中仅有一小部分体素存储的密度值大于 0，这导致大量“空”体素的存在——存储这些体素会造成显著的空间浪费。此外，可能存在 8 个体素存储相同密度值（例如 0.138）的情况，这同样是空间浪费。稀疏体数据的设计正是为了解决这一问题。
其核心思路如下：我们创建一个尺寸为 2×2×2 体素块大小的大体积素。因此，这个大体积素的尺寸是原始体素的两倍。接下来：若 2×2×2 体素块内的所有体素都具有相同的密度值（可以是 0 或任意大于 0的值，例如 0.139），则删除该 2×2×2 体素块，并将单一值存储在大体积素中（例如示例中的 0.139）；否则，该大体积素指向原始的 2×2×2 体素块。
这一过程是递归的：从最高分辨率网格（最高层级）开始，自底向上构建较低层级。指向更多体素块的体素称为内部节点（internal node），存储具体值的体素称为叶节点（leaf node）。下图展示了二维稀疏体数据的示例：
当体素块未被合并时，我们需要存储 1 个（大体积素）+8 个（2×2×2 体素块）体素；而当体素块被合并时，仅需存储 1 个体素。在大多数流体模拟中，原始体素的很大一部分要么为空，要么具有相似的密度值（例如云层场景中，云的核心区域通常高度均匀，密度变化主要出现在云的边缘）。因此，采用
稀疏表示方式编码这些体数据，可大幅减少缓存在磁盘和内存中的占用空间。
请注意，使用 2×2×2 的体素块并非强制要求。包括 OpenVDB 在内的许多系统，会在层级结构的前几层采用 八叉树（octree） 结构，而后将后续层级存储在 32×32×32 的体素块中（示例）。这种数据组织方式相较于更小的块尺寸，可能提升缓存一致性和数据访问效率。
最后，稀疏体数据在渲染中也十分有用：可跳过空的大体积素，将密度均匀的大体积素作为均质体渲染。这些优化组合可显著节省渲染时间。
稀疏体数据（与 LOD 和砖块映射有相似之处）是一个重要且内容丰富的主题，值得单独撰写章节讲解。若你对该主题感兴趣，也可查阅 Gigavoxels（体素块八叉树）相关资料。
核外渲染（Out-of-Core Rendering）
我们在稀疏体数据后提及核外渲染是有原因的：稀疏体数据和体素块八叉树的设计初衷（除上述优化外），正是为了能够渲染无法全部加载到内存中的超大体积数据缓存。
当无法将整个文件加载到内存中时（你不希望受限于模拟数据的大小），可通过缓存区机制仅加载所需的体素块（例如相机可见的体素块）。当然，这要求体数据缓存被组织为可按需实时加载的体素块集合。该主题将在后续章节中全面讲解。
注意网格数据的遍历方式
若你需要遍历体素，建议采用以下方式：
for (size_t x = 0; x < resolution; ++x) {
for (size_t y = 0; y < resolution; ++y) {
for (size_t z = 0; z < resolution; ++z) {
...
}
}
}
而非：
for (size_t z = 0; z < resolution; ++z) {
for (size_t y = 0; y < resolution; ++y) {
for (size_t x = 0; x < resolution; ++x) {
...
}
}
}
这是由数据在内存中的存储方式决定的：先遍历 x、最后遍历 z 的方式，可按内存中存储的顺序访问值；而先遍历 z、最后遍历 x 的方式，会跳转到内存的不同位置，极易导致大量缓存未命中（cache miss），从而降低性能。当然，这也取决于数据格式：我们的缓存数据是按先 x、再 y、最后 z 的顺序存储的，但其他格式可能采用不同的约定。因此，需留意代码的这一细节——这可能是性能优化的关键点。
光照烘焙
你可将光照烘焙到网格体素中，在渲染阶段直接从体素中读取这些数据，从而跳过在相机光线的每一步光线步进（raymarching）中计算 Li 项的过程，提升渲染速度。当然，你仍需在预处理阶段（烘焙阶段）为网格中的每个体素计算 Li 项，这一过程仍会耗时；此外，每当光照或模拟结果发生变化时，都需要重新烘焙。我们提及该技术主要是为了参考和历史溯源。
深度阴影图（Deep Shadow Maps）
阴影图正逐渐成为历史，但为了内容完整和历史溯源，我们认为有必要提及深度阴影图——其思路与将光照烘焙到体素中的方式有相似之处。该技术的核心是为场景中的每个光源计算一张阴影图，但与传统阴影图存储从光源到场景中最近物体的深度不同，深度阴影图存储的是体密度随距离变化的函数。换句话说，深度阴影图中的每个像素存储一条曲线，代表穿过体数据时的透射率变化。该技术由皮克斯（Pixar）的 Tom Lokovic 和 Eric Veach 于 2000 年提出（若想深入了解，可查阅其论文），目前已有多种衍生版本。
源代码
与往常一样，你可在本章的源代码部分找到对应代码。我们还提供了约 100 个缓存文件序列（下载 cachefiles.zip 并解压到程序源代码所在目录），可用于渲染本章开头展示的动画。