面向开发人员的体积渲染：基础（中）_参与介质渲染

本文翻译自Volume Rendering for Developers: Foundations，如有错误，欢迎指正。

三维密度场的体绘制(Volume Rendering of a 3D Density Field)

到目前为止，我们只学习了如何渲染同质体积对象。在空间中散射系数和吸收系数恒定的物体。这很好，但有点无聊，而且并不是事物本质上的样子。如果你观察蒸汽火车中冒出的烟雾，这些烟雾的体积是不均匀的。有些部分比其他部分更不透明。那么我们如何渲染异构体积对象呢？

在现实世界中，我们会说吸收或散射系数在空间中变化。吸收系数越高，体积越不透明。并且散射和吸收可以在空间中彼此独立地变化。不过，我们通常选择一种更实用的方法，即对给定体积对象的散射系数和吸收系数使用恒定值，并使用密度参数(density parameter)来调节体积体在空间中的外观。想象一下密度参数（在编程术语中只是一个实数，如浮点数或双精度数）随空间变化。然后我们可以做这样的事情：


实际上原文的公式是连在一块的，但是我感觉作者意思应该是要分开介绍，由读者们自己取舍了：

这里的 ${\sigma }_s^{\prime }$ 就是由空间变化密度参数调节的散射系数。比尔定律方程中使用的消光系数 ${\sigma }_t$ ，即吸收系数和散射系数之和，也受到空间变化密度参数的调节。而 $density(p)$ 是一种返回空间中p点密度的函数。通过这样做，我们可以生成体积对象的图像，例如上图中的火山烟柱。

现在的问题是我们如何生成这个密度场(density field)？您可以使用两种技术：

1. 程序化(Procedural)：您可以使用 3D 纹理（例如 Perlin 噪声函数）以程序化的方式创建空间变化的密度场。
2. 模拟(Simulation)：您也可以使用流体模拟程序来模拟流体（例如烟雾）的运动，以产生随空间变化的密度场。

在本章中，我们将使用第一种方法。在下一章中，我们将学习如何渲染流体模拟的结果。

使用噪声函数生成密度场

图 1：我们可以使用 Perlin 噪声来程序化生成密度场。噪声函数将一个点作为参数，并返回该点在 [-1,1] 范围内的噪声值。

---

在本课程中，我们将使用 Perlin 噪声函数程序化生成 3D 密度场。如果您不熟悉程序噪声生成的概念，我们建议您阅读以下两课：Value Noise and Procedural Patterns (Part 1) 和 Perlin Noise (Part 2)。

什么是程序化噪声函数（在本例中为 Perlin 噪声函数）？它是一个通过 3D 空间按程序生成噪声模式的函数（在该术语的编程意义上）。我们可以使用这种模式来生成一个密度场，其值随空间变化。噪声函数将一个点作为参数，并返回该点的 3D 噪声纹理值（实数，如浮点数或双精度数）。该值绑定在[-1,1]范围内。密度可以是 0（无体积）或正值，因此我们需要裁剪或重新映射噪声函数的值以获得正的密度值。在以下代码片段中，我们将值从 [-1,1] 重新映射到 [0,1]：

float density = (noise(pSample) + 1) * 0.5;
1

其中 pSample 是当我们穿过体积时样本沿相机光线的位置(position)。

对于噪声函数，我们将使用 Ken Perlin 本人提供的改进 Perlin 噪声的实现（链接）。同样，如果您有兴趣，可以在专门介绍 Perlin Noise (Part 2) 的课程中了解此代码的工作方式和原因。但在本课中，我们将假设您熟悉该函数。如果没有，也不必太担心。您需要关心的只是将 3D 空间中您想要对函数求值的点的位置传递给函数，并且它会为该点返回 [-1,1] 范围内的值。以下是供参考的代码（查看源代码部分中提供的文件以获取完整的实现）：

int p[512]; // permutation table (see source code)
 
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 ))));
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

最后，我们将在光线步进的程序中使用它：

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; // fully transmissive to start with
    for (size_t n = 0; n < numSteps; ++n) {
        float t = tMin + stepSize * (n + 0.5);
        Point p = ray.orig + ray.dir * t;
        // **change**：density is no longer a constant value. It varies through space.
        float density = (1 + noise(p)) / 2;
        float sampleAtt = exp(-density * sigma_t * stepSize);
        // transmission is attenuated by sample opacity
        transmission *= samplAtt;
        ...
    }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

图 2：密度沿射线变化。

图 3：异质体积(heterogeneous volume)的结果，其密度场由噪声函数生成。

---

正如您所看到的，与我们在前一章中用来渲染同质体积(homogeneous volume)对象的程序相比，这是一个相当简单的更改。我们将密度变量的声明移动到光线行进循环内，该变量不再是一个常量值，而是一个空间变化的参数。图 2 直观地显示了所发生的情况。当我们沿着射线行进时，我们对密度场进行采样，其中再次使用噪声函数来生成该密度场。对于沿射线的每个样本，我们使用样本位置作为函数输入参数来评估噪声函数，并将结果用作该点的密度值。图 3 显示了应用于体积球体的结果。

为了确保您了解这里发生的情况，我们使用相同的密度噪声函数绘制了一段距离上的噪声函数和相同距离上的透射值(transmission value)。结果如下图所示。我们还绘制了如果该体积具有恒定密度（绿色）与随空间变化的密度时比尔兰伯特定律(Beer’s Lambert law)的样子。如您所见，绿色曲线（恒定密度）非常平滑，而红色曲线（非均匀）则不然。请注意，当噪声函数返回接近 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);
}
1
2
3
4
5
6
7
8
9
10
11
12
13

但有一个问题。该代码用于计算体积透射值(transmission value)，但我们用来计算内散射贡献的代码（还记得 Li 项吗？）却不会按原样计算。现在，我们将解释为什么以及需要进行哪些更改才能使其与异构参与介质一起工作。

异质参与介质(heterogenous participating medium)的内散射

图 4：光线穿过异质参与介质。

---

希望通过查看图 4，您能了解问题所在。对于均质体积物体，当涉及到光线时，我们所要做的就是找到在光线方向上从样本点到物体边界的距离。然后使用该距离（我们称之为 Dl）和体积消光系数（ sigma_t = sigma_a + sigma_s ）应用比尔定律来找出在穿过到采样点的体积。简单的：

Color lighRayContrib = exp(-Dl * sigma_t * density) * lightColor;
1

异质体积的问题是，这不再有效，因为密度也沿着光线变化，如图 4 所示。请注意：我们解决的问题与我们在第 2 章中使用前向光线解决的问题不同。到目前为止，我们使用光线步进的原因是为了估计沿相机光线的内散射项。没有其他的。

然而，光线步进技术在这里将再次有用：用以估计内散射项以及估计相机和光线在穿过异构参与介质时的透射(transmission)。尽管技术相同（在这种特殊情况下是前向光线行进，这是随机采样方法的一种形式），但不是同样的问题（估计内散射与估计光线的透射）。我们需要将光线分成一系列片段，并估计每个段的透射，假设在片段的长度上（在由步长定义的小体积元素内），体积元素的密度是均匀的，然后将总透射率值乘以沿射线移动时样本的透射率。在伪代码中，可以按如下方式实现：

// compute light ray transmission in heterogeneous medium
float transmission = 1;
float stepSize = Dl / numSteps;
for (n = 0; n < numSteps; ++n) {
    float t = stepSize * (n + 0.5);
    float sampleAtt = exp(-density(evalDensity(t) * stepSize * sigma_t);
    transmission *= samplAtt;
}
1
2
3
4
5
6
7
8

将此代码与我们在前几章中使用的代码进行比较，以计算当我们沿着光线从 t0 行进到 t1 时相机光线透射值的值。

float sigma_t = sigma_a + sigma_s;
float density = 0.1; // density is constant. Used to scale sigma_t
float transparency = 1; // initialize transparency to 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; // sample position (middle of the step) 
 
// **heightlight**   
	// compute sample transparency using Beer's law
    float sample_transparency = exp(-step_size * sigma_t * density); 
 
    // attenuate global transparency by sample transparency
    transparency *= sample_transparency; 
// **end heightlight**
    // In-scattering. 
    if (hitObject->intersect(sample_pos, light_dir, isect_vol) && isect_vol.inside) { 
        ...
        result += ... 
    } 
 
    // finally attenuate the result by sample transparency
    result *= sample_transparency; 
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

这两个代码片段正在做同样的事情。我们可以利用一个数学小技巧，直到现在我们才谈到它（因为现在正是这样做的好时机）。这里是：

如果你看一下代码，你会发现透射值本质上是一系列指数相互相乘。如果展开光线步进循环（片段 2），您会得到如下内容：

// dx = stepSize, and noise(x) is in the range [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)) * 
                     ...;
1
2
3
4
5
6
7
8
9

如果我们使用我们刚刚学到的指数的数学属性重写这段代码，我们会得到：

float tau = noise(t0) + noise(t1) + noise(t2) + ...;
float transmission = exp(-tau * sigma_t * dx);
1
2

换句话说，当我们沿着光线步进时（就像我们对相机光线所做的那样），我们需要做的就是累积沿着光线的每个样本的密度值，然后使用这个总和来计算光线衰减/透射单次调用指数函数的值（这确实节省了一些时间）。下图说明了这个概念。

详细信息
从技术上讲，我们可以对相机光线透射值执行相同的操作，但请注意下面提供的代码，我们在沿着相机光线步进时使用透射值来衰减 Li 项。当我们在体积中前进时，我们需要光线透射的中间值，这就是为什么我们不只是将密度汇总成一个变量并在最后计算最终的光线透射值，就像我们计算光线的透射率一样。

“tau”这个名字并不是错误选择的。您经常会在文献中看到它被用来表示称为光学深度的量。该数量通常使用两个希腊字母：tau ( $\tau$ ) 或 rho ( $\rho$ )。我们不会在本章中给出光学深度的正式定义，因为此时这可能会令人困惑。但我们在Volume Rendering: Summary, Equations / Theory一章中进行了介绍。

就是这样！现在您已拥有渲染异构体积对象的准确图像所需的一切。最后一张图显示了给定噪声曲线的传输曲线。

实际（和功能）实施

让我们对我们的程序进行必要的调整，以演示这在实践中如何运作。请记住，该算法现在的工作原理如下：

1. 沿着摄像机光线进行光线步进。在沿相机射线的每个样本处，估计样本位置处的密度，以计算样本透射率并估计内散射贡献。
2. 沿着光线进行光线步进：之前我们只沿着相机光线进行光线行进，现在我们也需要沿着光线方向进行光线步进。对于均匀体积物体，我们只需要沿着相机光线行进即可估计内散射项。而对于异质，我们需要沿着相机光线这样做来估计内散射项（这不会改变），但现在还要评估沿着光线的密度项。我们需要沿着光线进行光线步进，以评估沿着这些光线的密度函数。

换句话说，我们现在需要沿着相机和光线行进，并沿着这些光线评估每个样本的密度函数。目前是 Perlin 噪声函数。这是大量的操作，正如您将看到的，将我们的球体渲染为异质介质将比其同质介质花费更长的时间。

// [comment]
// This function is now called by the integrate function to evaluate the density of the 
// heterogeneous volume sphere at sample position p. It returns the value of the Perlin noise
// function at that 3D position remapped to the range [0,1]
// [/comment]
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; // absorption coefficient
    float sigma_s = 0.5; // scattering coefficient
    float sigma_t = sigma_a + sigma_s; // extinction coefficient
    float g = 0; // henyey-greenstein asymetry factor
    uint8_t d = 2; // russian roulette "probability"

    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; // initialize transmission to 1 (fully transparent)
    vec3 result{ 0 }; // initialize volumetric sphere color to 0

    // The main ray-marching loop (forward, march from t0 to t1)
    for (int n = 0; n < ns; ++n) {
        // Jittering the sample position
        float t = isect.t0 + stride * (n + distribution(generator));
        vec3 sample_pos = ray_orig + t * ray_dir;

        // [comment]
        // Evaluate the density at the sample location (space varying density)
        // [/comment] float eval_density(const vec3& p)
        <span style="color: red; font-weight: bold; background-color: rgba(255,0,0,0.1); border: 1px none rgba(255,0,0,0.3);">float density = eval_density(sample_pos);</span>
        float sample_attenuation = exp(-step_size * density * sigma_t);
        transparency *= sample_attenuation;

        // In-scattering.
        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;
            // [comment]
            // Ray-march along the light ray. Store the density values in the tau variable.
            // [/comment] float eval_density(const vec3& p)
            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;
                <span style="color: red; font-weight: bold; background-color: rgba(255,0,0,0.1); border: 1px none rgba(255,0,0,0.3);">tau += eval_density(light_sample_pos);</span>
            }
            float light_ray_att = exp(-tau * stide_light * sigma_t);
            result += light_color *       // light color
                      light_ray_att *     // light ray transmission value
                      phaseHG(-ray_orig, light_dir, g) * // phase function
                      sigma_s *           // scattering coefficient
                      transparency *      // ray current transmission value
                      stride *            // dx in our Riemann sum
                      density;            // volume density at the sample location
        }

        // Russian roulette
        if (transparency < 1e-3) {
            if (distribution(generator) > 1.f / d)
                break;
            else
                transparency *= d;
        }
    }

    // combine background color and volumetric sphere color
    return background_color * transparency + result;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85

在此实现中，用于相机和光线的步长是相同的。这没有必要。为了加快速度，您可以使用更大的步长来估计光线传输值。另外，现在我们使用程序纹理，我们可能会遇到一些过滤问题。如果对噪声函数进行采样的频率太低，您可能会错过程序纹理中的一些细节，并且最终会出现锯齿问题。同样，这是一个过滤问题，我们现在不会深入研究，但请注意步长、噪声频率和图像分辨率以某种方式相互关联（从采样的角度来看）。

您可能会感到惊讶，程序输出（右）看起来并不更像云。然而，正如您在噪声模式图像（左）中看到的那样，构成噪声模式的“块”默认情况下相当平滑，这就是体积也具有平滑外观的原因。制作云状程序噪声的艺术需要以不同的方式调整噪声函数的结果，以获得视觉上更有趣的结果（我们将在下面进一步展示）。

您可以尝试以下几个简单的变化：从噪声函数中删除负值（左）并获取噪声函数的绝对值（右）。下面进一步解释衰减参数。

无论多么不起眼，这正是电影《接触》开场镜头（1997 年）中用于创建壮观体积效果的技术。这些图像是使用 Pixar 渲染器渲染并由 Sony Picture Imageworks 创建的。这个序列在当时是一项巨大的技术任务。如前所述，他们使用了与本课程中提供的技术相同的技术。不同之处在于，他们使用一些几何体（而不是基本球体）来定义体积物体的形状，并使用一些分形图案来赋予星云云状纹理。我们将在本课的最后一章中讨论前一种技术（生产中的体积渲染）。至于云状纹理，现在让我们看看我们能做什么…

使用密度函数（创建更有趣的外观和动画）

编写光线步进器是一回事。创建类似云的程序纹理是另一回事。前者是一门科学，而后者更是一门艺术：使用一系列数学工具来塑造程序纹理，并花费大量时间调整参数，直到最终得到令你满意的东西。我们的目标不是创建令人信服的令人惊叹的图像，而是为您提供工具或砖块，首先帮助您理解事物是如何工作的，其次，您最终可以将自己重新组合到更复杂的系统中，以制作一些“令人惊叹”的图像，如果您想要。但我们会将其留给您…如果您创造了一些很酷的东西，请与我们分享！

现在，您可以使用一些技巧来创建更有趣的云状“球体”。这只是其中一些工具的简短选择。

Smoothstep

smoothstep 函数是您熟悉的函数，因为我们已经将它用于多种用途，包括噪声函数。它在两个值之间创建“平滑”过渡。这是该函数的一种可能的实现：

float smoothstep(float lo, float hi, float x)
{
    float t = std::clamp((x - lo) / (hi - lo), 0.f, 1.f);	//t=(x-lo)/(hi-lo)返回0-1之间的数
    return t * t * (3.0 - (2.0 * t));	//返回3t²-2t³
}
1
2
3
4
5

我们可以使用此函数在球体边界附近创建衰减。为此，我们将对 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); // smooth transition from 0 to 1 as distance goes from 0.1 to 1
    return (1 - falloff);
}
1
2
3
4
5
6
7

同样，此技术对于在体积到达球体边界之前淡出体积非常有用。我们使用 smoothstep 函数来控制衰减开始（并最终结束，但对于衰减效果，应保持为 1）距离边缘或侧面多远。

fBm

fBm（在很久以前也称为等离子体plasma或“混沌chaos”纹理）代表分形布朗运动(Fractal Brownian Motion)。在计算机图形学领域，它具有与数学不同的含义。我们将考虑这对我们图形工程师和艺术家意味着什么。它是一种由程序化噪声层的总和组成的分形图案，其频率和幅度因层而异。您可以在“程序生成”部分的课程中找到有关此模式的一些信息。

典型的 fBm 程序纹理的代码可以构造如下：

float eval_density(const vec3& p, ...)
    vec3 vp = p - sphere_center;
    ...
    // build an fBm fractal pattern
    float frequency = 1;
    vp *= frequency; // scale the initial point value if necessary
    size_t numOctaves = 5; // number of layers
    float lacunarity = 2.f; // gap between successive frequencies
    float H = 0.4; // fractal increment parameter  分形增量参数
    float value = 0; // result of the fBm (use this for our density)
    for (size_t i = 0; i < numOctaves; ++i) {
        value += noise(vp) * powf(lacunarity, -H * i);
        vp *= lacunarity;
    }

    // clip negative values
    return std::max(0.f, value) * (1 - falloff);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

构建 fBm 模式需要两行代码（在各层上循环）。可以从此基本版本构建许多变体，例如取噪声函数绝对值（一种称为湍流(turbulence)的模式）等。如果您想了解有关此主题的更多信息（并了解如何正确过滤此模式，请再次检查程序生成部分）。

Bias

偏差会移动中心值 (0.5)，同时保持 0 和 1 相同。(Bias shifts where the center (0.5) value will be while leaving 0 and 1 the same.)

float eval_density(...) 
{
    float bias = 0.2;
    float exponent = (bias - 1.0) / (-bias - 1.0)
    // assuming exponent > 0
    return powf(noisePattern, exponent);
}
1
2
3
4
5
6
7

Rotate the noise pattern using the primitive local frame

另一种技术在很久以前（由于计算能力的提高，在流体模拟成为主流之前）被广泛使用，对体积基本体（球体、立方体等）内的噪声模式进行动画处理。我们传递给密度函数的样本位置是在世界空间中定义的，但我们可以在附加到球体的参考系统中定义该点位置。换句话说，我们将样本点从世界空间变换到对象空间。要在球体参考系统内定义样本点，我们需要做的是：

vec3 sample_object_space = sample_pos - center;
1

更通用的解决方案包括使用矩阵转换球状基本体。然后，我们将使用对象到世界矩阵的逆矩阵将样本点从世界空间转换为球体局部参考系统。然而，为了简单起见，我们在示例程序中没有使用矩阵。但您可以利用 Scratchapixel 上提供的信息轻松地自行完成此操作。

该技术很有用，原因如下：如果我们移动球体，球体的局部参考系中定义的点的坐标将不会移动。这可以用来确保噪声图案与球体保持一致，无论其变换如何。类似于您在手指之间滚动的玻璃弹珠的纹理。在这个例子中，我们选择来说明这个想法（尽管它与我们刚刚解释的有点不同，但概念和结果是相同的）：我们将围绕原始球体 y 轴旋转物体空间中的点（一个更好的解决方案是使用对象到世界矩阵对球体进行评价(rate the sphere)，然后使用世界到对象矩阵将样本点从世界空间转换到对象空间，但我们懒得在这里这样做，所以我们决定旋转在对象空间中的点）。结果，我们可以看到分形图案也在旋转，有点像我们看到玻璃弹珠旋转一样。使用此方法可以应用更复杂的变形（线性或非线性）。

float eval_density(const vec3& p, const vec3& center, const float& radius)
{ 
    // transform the point from world to object space
    vec3 vp = p - center;
    vec3 vp_xform;

    // rotate our sample point in object space (frame is a global variable going from 1 to 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);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

请注意我们现在如何清楚地看到分形图案的 3D 结构。

And many more…

“程序化”技巧的列表可以一直列下去：位移（displacement ，我们使用噪声函数置换体积的边缘）、不同类型的噪声（波浪billow、时空space-time，这是一种动画类型的噪声等）。我们最终将在本课程的未来修订中扩展此列表。

到目前为止我们学到了什么？

请参阅最后一章，了解本课所学习技术的总结。

不过，您现在可以做的是使用光线行进来渲染单散射异质体积对象。这已经是一个很大的成就了。现在应该清楚的另一件事是，通过光线步进，我们将体积对象分解为由步长定义的更小的体积元素。这项技术之所以有效，是因为我们假设这些微小的体积元素或样本足够小，本身是均匀的。简而言之，您将异构对象分解为小“砖块”，您可以将其视为“同质”。尽管砖块本身可能彼此不同。有点像用乐高积木搭建一个物体。

源代码

本课程的源代码部分提供了实现此技术的全功能程序的源代码。您可以使用 eval_density() 函数来创建不同的外观。

基于三维体素网格(3D Voxel Grids)的体绘制

在本章结束时，您将能够生成以下图像序列：

如前一章所述，您可以使用两种技术创建密度场：程序化或使用流体模拟软件。在本章中，我们将讨论后者。请注意，在本课程中，我们不会学习流体模拟的工作原理。我们将学习如何渲染流体模拟生成的数据。我们承诺，将来我们会学习流体模拟。

步骤 1：使用 3D 网格存储密度值

存在许多用于模拟流体（烟雾、水等）的技术，但通常，在过程中的某个时刻，模拟结果存储在 3D 网格中。为了简单起见，在本课中，我们将假设这些网格在所有三个维度上具有相同的分辨率（例如，x、y 和 z 坐标分别为 32x32x32），并且该分辨率是 2 的幂（8、16、32、 64、128 等）。这只是为了本课的目的。实际上，情况不一定如此。另外，我们在本课中不会改变网格。因此，我们假设网格是一个轴对齐的盒子：我们可以将网格视为轴对齐的边界框（AABB），这将简化我们的射线盒相交测试。使解决方案通用（又名支持未轴对齐的网格，可以通过使用网格的世界到对象变换矩阵将相机光线变换到网格对象空间来轻松完成，如本课程中所述）。

图 1：8x8x8 网格，存储密度值。

---

网格很适合模拟流体的运动，因为网格的体素（这些是构成网格的小体积元素，我们也可以称为单元）设置有一些初始密度（想象它们充满了烟雾），并且随着时间的推移，该密度分布在相邻单元之间。它们从一个体素移动到另一个体素的方式由纳维-斯托克斯方程(Navier-Stokes equation)决定。但同样，这个话题留到另一课。在本课程中，您需要了解的是，流体模拟的结果存储在由存储密度值（0 或大于 0 的值）的体素组成的 3D 网格中。网格中的每个体素都存储一个唯一的密度值（Each voxel in the grid stores one unique density value，密度“填充”体素的体积）。 3D 网格过于流畅的模拟，就像位图之于图像一样。密度场被量子化。在任何编程语言中定义存储密度值的 3D 网格都相当简单，因为我们将在本章末尾研究它们的微妙之处。

除了网格分辨率（grid resolution，它在任何维度（如 32x32x32）中包含的体素数量）之外，我们还需要定义世界空间中网格的大小(size of the grid)。这是场景中物体的大小。我们可以通过不同的方式来做到这一点，但为了方便起见，我们将通过定义网格在世界空间中的最小和最大边界或范围来实现这一点，如下所示：(-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) 不行。以这种方式定义网格世界空间大小很方便，因为最小和最大范围可以直接插入光线盒相交例程（我们在 A Minimal Ray-Tracer: Rendering Simple Shapes (Sphere, Cube, Disk, Plane, etc.) 中学习到的）。因此，我们已经知道如何计算与立方体相交的任何射线的 t0 和 t1 值（类似于我们为球状基本体所做的方式）。

您可以修改前面章节中的代码来渲染立方体而不是球体。现在我们还可以使用矩阵从更有趣的角度渲染场景。你应该得到这样的东西：

struct Grid {
    float *density;
    size_t dimension = 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) {
        // ray march from t0 to 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;
    // allocate memory to read the data from the cache file
    size_t numVoxels = grid.resolution * grid.resolution * grid.resolution;
    // feel free to use a unique_ptr if you want to (the sample program does)
    grid.density = new float[numVoxels];
    std::ifstream cacheFile;
    cacheFile("./cache.0100.bin", std::ios::binary);
    // read the density values in memory
    cacheFile.read((char*)grid.density, sizeof(float) * numVoxels);

    Point origin(0);
    for (each pixel in the frame) {
        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;
    }

    // free memory
    delete [] grid.density;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

现在不用太担心缓存文件的内容。我们稍后会研究这个问题。但从上面的代码中您应该看到，我们所做的就是创建一个连续的内存块来存储与网格中的体素一样多的浮点值。然后将值从磁盘移至内存（第 36 行）。没有什么花哨。

这是您应该得到的图像。现在我们知道我们将使用 3D 网格来存储密度场值，该网格具有分辨率和大小，并且我们知道如何计算光线与框（也称为网格）相交的点，让我们看看当我们沿着射线行进时如何读取网格密度值。

步骤2：光线行进穿过网格

正如本课程中所解释的，使用网格作为加速结构非常复杂。我们必须从射线相交的网格中找到每个单元，以查看每个相交单元包含的几何图形（三角形）是否也与射线相交。虽然不是超级复杂，但这也不是一个超级琐碎的任务（特别是如果你想让它更快）。所以你可能会担心我们必须再次经历这个过程。好消息：一点也不。对于光线行进，这是一项微不足道的任务，因为我们关心的是沿着光线的样本点，如下图所示。

由于我们知道这些点的位置，我们剩下的就是找出这些点重叠的网格中体素的位置（这是下一项的主题）并读取存储在这些体素中的密度值。不涉及网格遍历。超级简单。现在，让我们看看如何从采样点转到体素坐标并使用这些坐标从内存中读取值。

步骤 3：沿着光线步进时读取密度值

我们需要对迄今为止使用的代码进行更改的是 evalDensity 函数。我们将读取该位置处网格的内容，而不是评估空间中任何给定点处由程序噪声函数定义的密度场。这是一个非常小的变化。其他一切都一样。通常我们知道样本点总是包含在网格的边界内（因为它应该位于射线与盒子感兴趣的点之间的某个位置）。因此，这应该不是问题，但为了使代码更加健壮，您可能仍然希望根据网格边界检查样本点坐标，以防止崩溃。

现在，从网格读取密度值可以有不同的名称。我们通常称之为查找(lookup)。在 OpenVDB（一个旨在高效存储和读取体积数据的库）中，这称为探测(probing)（尽管是一个相当不常见的名称）。好的，那么我们如何读取数据呢？

我们需要做的是找到样本点和它所属网格中的体素之间的对应关系。原始样本点是在世界空间中定义的。所以我们需要将该样本点转换为网格离散坐标。为此，我们需要：

• 将点从世界空间变换到对象空间。在我们的情况下，我们需要做的就是从样本位置移除网格最小范围(没有负数的坐标)。现在，样本点是相对于网格的“左下”角定义的。
• 然后我们需要将对象空间中的点除以网格大小（grid size，在我们的示例中为 10）。该点的坐标现已标准化。我们可以说标准化坐标（点的坐标包含在 [0,1] 范围内）。
• 最后，我们需要将标准化坐标乘以网格分辨率。我们称这个空间为体素空间(voxel space)。此时，坐标仍然定义为浮点数。为了计算网格坐标并从内存中读取这些坐标的值，我们需要将这些浮点数字四舍五入到最接近的整数。为了防止崩溃，我们需要确保体素空间中样本的坐标不大于网格分辨率减1（如果任何标准化坐标恰好为 1，则可能会发生这种情况）。我们现在可以读取该位置的密度值。

如果您查看上面的 (2D) 示例，您可以看到该点在体素空间中的坐标为 (3.36, 3.28)。从那里我们知道，找到内存中体素的索引(index of the voxel)只需将 y 坐标的整数部分乘以 8（网格分辨率）加上 3（x 坐标的整数部分）即可。

这是代码：

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;		//转换到体素空间

	//对体素空间的xyz取整，为了方便后期确定在体素空间的的索引
    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));

    // nearest neighbor 使用最近邻方法来采样网格密度
    return grid->density[(zi * grid->resolution + yi) * grid->resolution + xi]; 
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14

在这里（就像 Perlin 噪声函数一样，它也需要将点从世界空间转换为格子空间(lattice space)），我们使用floor函数，因为如果归一化点的坐标低于 0（我们将看到为什么会发生这种情况）等一下），我们希望函数返回 0 而不是 -1（简单的 static_cast 最终会这样做）。

我们将此方法称为最近邻搜索 (nearest-neighbor search) 插值，因为您本质上所做的是检索样本点落入的体素的坐标（因此它不是插值）。不过，还存在其他方法，它们都是从计算该体素的坐标开始的。这些其他方法涉及对我们接下来将研究的相邻体素的值进行插值。但在我们研究该主题之前，我们现在已经拥有渲染流体模拟的第一张图像所需的一切（源代码部分中提供了各种流体缓存，供您与示例程序一起下载）。让我们看看它是什么样子的。

如您所知，我们的目标是提供不依赖于外部库的示例。因此，我们不使用 OpenVDB（或 Field3D，OpenVDB 的祖先）等行业格式来存储网格数据。相反，我们使用类似于以下内容的最基本方式将流体模拟的数据转储到二进制文件：

// our binary cache format
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));
        }
    }
}
1
2
3
4
5
6
7
8
9

请注意，这和OpenVDB 的工作原理相同，但是，它们支持压缩以及多分辨率和稀疏体积(sparse volumes)等功能（我们将在本课末尾讨论）。我们的方法的优点是简单明了，这正是我们在教学时所寻求的。因此，在上面的代码中需要注意的是，我们沿着 z 轴从后到前存储体素切片。

• 请记住（对于右手坐标系，这里应该是指在世界空间中）网格的最小范围位于网格的背面，左下角（下图中的洋红色体素）。
• 在体素空间中，网格中的第一个体素的坐标为 (0,0,0)，在我们的示例中，最后一个体素（网格前面，沿正 z 轴，位于右上角）的坐标为 ( 8,8,8)（下图中黄色体素，本例中网格分辨率为8）。

因此，如果我们在一个简单的浮点数组中映射体素，就像我们将数据转储到文件时（以及当我们从文件中将它们读回内存时）一样，该数组中的第一个体素为索引 0，最后一个体素为索引 0。 voxel 的索引为 511（它是一个 0-based array，因此 888 减 1）。要从体素空间坐标计算体素的索引，我们需要做的就是：

size_t index = (zi * gridResolution + yi) * gridResolution + xi;
1

这是获得结果的最基本方法。然而，我们可以使用三线性插值来提高图像质量。让我们看看这是如何工作的。

步骤 4：通过三线性插值(Trilinear Interpolation)改进结果

基本思想如下：体素代表单个密度值，但是样本点可以位于任何给定体素体积内的任何位置。因此，假设该位置的密度值应该以某种方式混合样本点所在的体素的密度以及与其相邻的体素的密度值，如图所示，这是有意义的：

对于三线性插值，我们简单地混合 8 个体素的值。我们如何混合它们？简单地通过计算样本点到这 8 个体素中每一个体素的距离，并使用这些距离来权衡这 8 个体素中每一个对最终结果的贡献。

该过程是线性插值，但是是 3D 的。我们用来插值 2 个值的方程是 $a\ast (1-weight)+b\ast weight$ ，其中 $weight$ 从 0 到 1 变化。在 2D 中，此过程称为双线性插值，需要 4 个像素。在 3D 中，我们称之为三线性插值，它需要 8 个体素。请注意，用于在 a 和 b 之间进行插值的方程称为插值。在线性插值的情况下，插值是一阶多项式。

我们假设给定体素中存储的密度是采样点恰好位于该体素中间时应采用的值。由于三线性插值方案的性质，为了获得该结果，我们需要将体素空间中的样本位置移动 -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;
    //**highlight**
    Vector pLattice(pVoxel.x - 0.5, pVoxel.y - 0.5, pVoxel.z - 0.5);
    //**end highlight**
    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;

    // trilinear interpolation
    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(xi + i, yi + j, zi + k);
            }
        }
    }

    return value;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

现在让我们看看这是如何（以及为什么）起作用的。

上图中说明了这两种技术。在左上角，您有一个最近邻插值方法的示例。这里没什么特别的。我们有 4 个体素（在 2D 情况下，在 3D 情况下有 8 个），样本点（蓝点）落在其中之一。因此，采样点呈现存储在该体素中的密度值。

三线性插值使事情变得有点复杂。您可能会同意，如果我们的采样点位于四个体素的中间，我们可能应该返回这四个体素中存储的密度值的平均值。也就是说，在此特定示例中，体素密度值的总和除以 4。

result = (0.9 + 0.14 + 0.08 + 0.63) / 4
1

现在，如果该点正好位于四个体素的交叉处（这是一种特殊情况），则可以正常工作。所以我们需要找到一个通用的解决方案。

首先（稍后我们将解释为什么需要这样做），我们需要将采样点偏移半个体素空间：（-0.5，-0.5，-0.5）。如果我们保留恰好落在四个体素中间的样本点的示例，则应用了偏移的样本点现在位于左下体素的正中间。我们接下来要做的就是计算从这一点到每个体素边界的距离。在我们的插图（2D 情况）中，这由距离 dx0（从样本点 x 坐标到 x0（体素左下 x 坐标）的距离）、dy0（从样本点 y 坐标到 y0（体素）的距离来表示。左下 y 坐标）、dx1（样本点 x 坐标到体素右上 x 坐标 x1 的距离）和 dy1（样本点 y 坐标到体素右上 y 坐标 y1 的距离）协调）。这些在技术上称为 曼哈顿距离(Manhattan distances) 。

曼哈顿距离：沿直角轴测量的两点之间的距离。在 p1 位于 (x1, y1) 且 p2 位于 (x2, y2) 的平面中，它是 |x1 - x2| + |y1 - y2|。

实际上，正如我们将看到的，我们只需要 dx0、dy0 和 dz0。这是三线性插值的总体思路（如果您想要更完整的解释，请查看专门介绍插值的课程）。

我们有一个 2x2x2 体素块。我们称它们为 v000、v100（向右移动 1 个体素）、v010（向上移动 1 个体素）、v110（向右 1 个、向上 1 个）、v001（向前移动 1 个体素）、v101（现在你明白了)、v011 和 v111。这个想法是使用 x0 作为权重对 v000-v100、v010-v110、v001-v1001 和 v011-v111 进行线性插值。这给了我们 4 个值，然后我们使用 y0 作为权重进行线性插值（再次成对）。我们留下 2 个值，然后我们最终使用 z0 作为权重进行线性插值。请注意，无论您使用 x0、y0 还是 z0 开始前 4 次线性插值，只要每次使用不同的维度对连续结果进行线性插值，都没有任何区别。请记住，我们的线性插值是 $a\ast (1-\omega )+b\ast \omega$ ，其中 $\omega$ 是权重。在代码中，我们得到这样的结果：

float result = 
    (1 - z0) * (  // blue
                (1 - y0) * (v000 * (1 - x0) + v100 * x0) + // green and red
                     y0  * (v010 * (1 - x0) + v110 * x0)   // green and red
                ) + 
          z0 * ( // blue
                (1 - y0) * (v001 * (1 - x0) + v101 * x0) + // green and red
                     y0  * (v011 * (1 - x0) + v111 * x0)); // green and red
1
2
3
4
5
6
7
8

希望该布局可以帮助您了解插值的三个嵌套级别。具有 4（红色）+ 2（绿色）+ 1（蓝色）插值。现在，如果我们分别用 wx0、wy0、wz0、wx1、wy1 和 wz1 替换 (1-x0)、(1-y0)、(1-z0)、x0、y0 和 z0，则展开并重新排列前面的代码片段（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
2
3
4
5
6
7
8
9

这是代码片段 1 中的代码（尽管在代码片段 1 中，权重是在嵌套循环中动态计算的）。如果您使用上图中的 2D 示例进行数学计算，其中 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
1
2
3

这就是我们正在寻找的结果。很棒的三线性插值（在 2D 情况下为双线性）有效。您现在还了解了为什么我们需要在体素空间中将采样点偏移 -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
1
2
3
4
5
6
7
8

换句话说，该样本返回的值将是存储在该体素中的值（如果我们假设偏移之前的样本点位于左下体素的中间，则在我们的示例中为 0.9），这正是您所期望的。

这就是三线性插值的全部内容。现在让我们看看结果与最近邻方法相比有何不同。

正如您所看到的，使用三线性插值（右）渲染的缓存图像更柔和。尽管这种改进是有代价的，因为使用此方法时渲染时间会显着增加（与最近邻方法相比）。三线性插值大约是最近邻查找所需时间的 5 倍（假设没有优化）。

现在你可以（再次）阅读其他人的代码：OpenVDB

让我们看看 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)); 
    }; 
 
    // Trilinear interpolation:
    // The eight surrounding lattice values are used to construct the result. \n
    // 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 
BoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord, 
                   typename TreeT::ValueType& result) 
{ 
    ... 
    const Vec3i inIdx = local_util::floorVec3(inCoord); 
    const Vec3R uvw = inCoord - inIdx; 
 
    // Retrieve the values of the eight voxels surrounding the
    // fractional source coordinates.
    ValueT data[2][2][2]; 
 
    const bool hasActiveValues = BoxSampler::probeValues(data, inTree, Coord(inIdx)); 
 
    result = BoxSampler::trilinearInterpolation(data, uvw); 
 
    ... 
} 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

sample() 方法首先获取存储在 8 个体素中的值。这发生在方法 probeValues 中。然后，如上所述（在 trilinearInterpolation 方法中）对体素进行插值。函数 _interpolate 只不过是代码中更常称为 lerp 的线性插值。

进一步了解体积缓存(volume caches)

我们只是列出一系列我们认为值得一提的主题，以供参考。为了使本课程保持相当简短，我们现在不会讨论太多细节。这些主题中的大多数将在本课程的未来修订中进行讨论，或者在单独的课程中进行研究。

其他插值方案：cubic and filter kernels

如前所述，线性插值是一阶多项式。可以使用高阶多项式来获得更平滑的结果。例如，您可以使用三次插值。在 2D 中，我们需要采样 4x4 像素来执行双三次插值。在 3D 中，我们需要 $2^4=64$ 体素。正如您可以想象的那样，结果看起来会更平滑，但渲染时间会增加。

您还可以使用filter kernels，类似于图像滤波中使用的滤波核（例如三角形、米切尔或高斯滤波器内核）。与图像一样，滤波器的范围越大，处理时间就越长。

这两种技术将在本课程的未来修订版中得到全面解释。

存储更多数据和密度

生产中使用的卷缓存格式通常旨在让您在缓存中存储所需的任何通道以及密度，例如温度（用于渲染火焰）、速度（用于渲染 3D 运动模糊）等。

滤波和砖图(Filtering and brick maps)

如果您出于生产目的而渲染体积缓存，您可能应该关心过滤和 LOD。渲染体积缓存可能遇到的问题与渲染纹理相同。纹理具有固定的分辨率（像素数，如 512、1024 等），如果您从超远距离渲染应用了给定纹理的对象，那么您会假设您的对象在图像中只有几个像素宽，当对象移动或相机移动时，每次渲染新帧时，都会渲染纹理中的不同像素或纹素。因此，纹理会在帧与帧之间出现抖动（至少可以这么说，但更多时候它看起来只是一种噪音）。为了尽量减少这个问题，我们制作了较小版本的纹理，与原始纹理数据一起存储；我们将这些原始纹理的低分辨率版本称为 mipmap。当对象较远（图像中较小）时使用较低分辨率的纹理可以消除噪声（技术上称为混叠）。我们还没有触及 mipmap 主题，但我们很快就会触及。

同时，如果您已经熟悉 mipmap 的概念、它们是什么、它们的用途以及如何创建它们，您应该知道我们在过滤 2D 纹理时遇到的问题也是我们在处理 2D 纹理时遇到的问题3D 缓存（或 3D 纹理）。因此，我们可以将mipmap方法应用于3D缓存。为此，您需要创建原始网格数据的缩小版本，并根据图像中缓存的大小使用正确的级别（缩小版本）。虽然该过程与 mipmap 非常相似，并且您也可以将这些缩小版本称为 mipmap，但在 3D 中，我们将它们称为砖块贴图( brick maps )。它们看起来像《我的世界》中的物体。

详细信息
brick maps这个术语是由皮克斯 Renderman 团队创造的，我们相信它，但我们不太确定它的起源…如果有人知道，请写信给我们。但一般来说，让缓存存储其数据的多种分辨率的想法并不是特别新或不常见。所有生产渲染器和格式（例如 OpenVDB）可能都支持它。

构建这些 brick maps 与构建纹理的 mipmap 级别非常相似，但是，我们取 8 个体素的平均值（而不是取 4 个像素或纹素的平均值）。这就是为什么处理分辨率正好是 2 的幂的网格很方便，因为在这种情况下缩小级别变得微不足道。

下面是一些代码，展示了如何从我们原始的流体模拟缓存创建这些级别（级别 0 的基本分辨率为 128）。这只是一个简单的例子，直到我们有机会正确解决这个主题：

size_t baseResolution = 128;
size_t numLevels = log2(baseResolution); /* float to size_t implicit cast */
std::unique_ptr<Grid []> gridLod = std::make_unique<Grid []>(numLevels - 2); // ignore res. 2 and 4

// load level 0 data
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));
			}
		}
	}
}
...
// to render level 3 of the cache for example (resolution: 16)
trace(ray, L, transmittance, rc, gridLod[3]);
...
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

我们现在不会解释如何选择正确的级别，但我们会在以后的课程修订中解释。与此同时，至少您知道问题所在以及解决方案是什么。

这是在不同级别渲染的缓存图像（级别 0 是原始缓存分辨率，在我们的示例中为 128）：

运动模糊(Motion blur)

3D 运动模糊怎么样？为了渲染这种效果，我们需要存储具有流体运动的体素。这通常表示为称为运动矢量的方向。该矢量的方向指示流体在网格中移动的平均方向；它的大小是流体移动的速度。使用此信息，可以在渲染时模拟 3D 运动模糊。

该主题将在以后的课程中研究。

平流(Advection)

虽然平流与体渲染没有直接关系（它更多地与流体模拟有关），但我们将其添加到此处作为提醒（以便我们将来可以编写有关它的课程（如果可以的话））。可以在渲染时完成，将细节添加到现有的流体模拟中。

Sparse volumes: 它们是什么？

很多时候，网格中只有一小部分体素存储大于 0 的密度值。这就留下了很多“空”体素；储存它们可能会造成巨大的浪费。同样，您可能有一组 8 个体素存储相同的密度值，例如 0.138。这也是一种空间浪费。Sparse volumes 旨在解决这个问题。

一般思路是这样的：我们创建一个体素，其大小相当于一个 2x2x2 体素块。因此，这个较大的体素是原始体素大小的两倍。下一步：如果 2x2x2 体素块中的所有体素都具有相同的密度值（可以是 0 或任何大于 0 的值，例如 0.139），然后我们删除 2x2x2 体素组并将单个值存储在较大的体素中（例如 0.139 以遵循我们的示例），否则，较大的体素指向 2x2x2 体素的原始块。

此过程是递归的。我们从最高分辨率的网格（最高级别）开始，然后自下而上构建较低的级别。指向更多体素块的体素是内部节点( internal nodes )，而存储值的体素是叶节点( leaf nodes )。下图显示了稀疏 2D 体积的示例。

当一个块没有合并时，我们存储 1（较大的体素）+ 8 个体素（2x2x2 体素块），但是当一个块被合并时，我们只存储 1 个体素。在大多数流体模拟中，大部分原始体素要么为空，要么共享相似的密度值。例如，云(clouds)通常就是这种情况。云的核心通常是非常均匀的。密度在云的边缘基本上不同。因此，使用稀疏(sparse)表示形式对这些卷进行编码可能是减小磁盘和内存中缓存大小的好方法。

请注意，使用 2x2x2 体素块不是强制性的。包括 OpenVDB 在内的许多系统为层次结构的前几个上层选择八叉树结构，然后将以下级别存储到 32x32x32 体素块中（例如）。以这种方式组织数据可能会提高缓存一致性和数据访问（与使用较小的块大小相比）。

最后，稀疏体积(sparse volumes)在渲染中也很有用。可以跳过空的大体素。具有均匀密度的大体素可以呈现为均匀体积。加在一起，这些可以节省大量时间。

稀疏体积（与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) {
            ...
        }
    }
}
1
2
3
4
5
6
7

而不是那样：

for (size_t z = 0; z < resolution; ++z) {
    for (size_t y = 0; y < resolution; ++y) {
        for (size_t x = 0; x < resolution; ++x) {
            ...
        }
    }
}
1
2
3
4
5
6
7

这是因为数据在内存中的布局方式。通过首先迭代 x 和最后 z，我们按顺序迭代存储在内存中的值。通过首先迭代 z 和最后 x，我们将跳转到内存的不同部分。这可能会产生大量缓存未命中，从而降低性能。当然，这取决于格式。我们将数据存储在缓存中，首先遍历 x，然后遍历 y，然后遍历 z，但其他格式可能使用不同的约定。因此，请注意代码的这一特定方面，因为这可能是寻找优化的地方。

烘培光照

您可以在网格的体素中烘焙光照，并在渲染时从体素读取此数据，以跳过在摄像机光线的每个光线行进步骤中计算 Li 项的步骤。这将加快渲染速度。当然，您仍然需要在预通道（烘焙通道）中为网格中的每个体素渲染此 Li 项，这仍然需要一段时间。此外，每当照明或模拟发生变化时，您都需要重新开始烘焙。我们提到这种技术主要是为了参考和历史原因。

深度阴影贴图(Deep shadow maps)

阴影贴图即将成为过去的产物，尽管出于完整性和历史原因，我们认为提及深阴影贴图会很棒，它在某种程度上与将光照烘焙到网格体素中的想法有相似之处。这里的想法是计算场景中每个光源的阴影贴图。但是，深阴影贴图不是存储任何给定像素从光线到场景中最近对象的深度（这是阴影贴图的用途），而是将体积的密度存储为距离的函数(deep shadow maps store the density of the volume instead as a function of distance.)。换句话说，深阴影贴图中的每个像素都存储了一条曲线，该曲线表示通过体积的透射(transmission)如何随着我们移动而变化。这项技术是由皮克斯的Tom Lokovic和Eric Veach于2000年开发的（如果你想了解更多，你可以在这里找到他们的论文）。该技术存在几种变体。

源代码（外部链接GitHub）

像往常一样，您可以在本课程的源代码部分找到本章的源代码。我们还提供了一个大约 100 个缓存文件的序列（下载 cachefiles.zip 文件并将其解压缩到本地驱动器上以及程序源代码），您可以使用这些文件来渲染本章开头所示的动画。

本章源代码地址