SIGGRAPH 2023 | 《地平线 西之绝境》DLC「燃烧海岸」中的体积云系统NUBIS 3

导语：最新的SIGGRAPH 2023高级实时渲染专场已经部分Release，本篇为大家介绍Guerrilla Games的首席环境艺术家兼技术美术工程师Andrew Schneider的Nubis 3系统。简单概括，Nubis是一套完整的体积云解决方案，包括引擎内渲染，资源数据生成，编辑器和天气驱动等模块。Nubis 3的实际开发大约花了不到一年的时间，但系统所涵盖的设计理念与实现思路却经过了将近10年的沉淀【2014~2023】。对于云建模大体经历了“高度梯度插值”（SIG 15）, “垂直剖面模型”(SIG 22), “3D网格距离场”(SIG 22~23)三个版本，启发了非常多In-House Engine。

今年的课程PPT一共194页，喜欢直接阅读第一手资料可参看：

https://advances.realtimerendering.com/s2023/index.html

或者选择阅读本文的总结，这里先给一下PPT的阅读建议：

• 1~49页是SIG 22的快速回顾，如果已经读过或者看过我写的文章可以直接略过

• 51~58页是方案开发背景，便于大家理解系统重构时面临的情况与考虑，Nubis 3的Voxel Cloud是为《地平线：西之绝境》DLC“燃烧海岸”准备的方案【暗示了小规模试验性应用的前提】，因为场景设计与玩法对云提出了更严格的考验，原来的方案面对新的要求比较吃力

• 59~109页讲的时体素云的密度场表征

• 110~133是对光照的设计解读 ，也就是适应体素渲染的照明积分方程建模，作者对这部分着墨并不多。可能是方案仍未十分完善的缘故

• 134~160页介绍的是步进控制，特别由于引入3D距离场，专门介绍了3D距离场的解编码

• 161~189页是该方案的用例介绍，主要是设计上的思维启发。特别的，花大篇幅介绍了新的雷暴云建模与游戏玩法的配合

• 190~194页扩展总结

可根据自己的需要跳到对应部分阅读。

接下来我们依旧从分布与造型方案开始分析。

一、体积密度场构建

1.1 体素云建模

首先是云建模，也就是云造型的数据来源。这部分方案有很多，但最好的手段是离线生成VDB数据。但注意，在这里只是使用支持VDB的DCC工具生成原始造型，但引擎里所谓的体素化是指对场景空间和Lighting空间体素化，并不是直接使用VDB的稀疏体素。生成造型的最终目的是烘焙Dimentional Profile【可以由Mesh SDF得到，注意这里的Profile是3D的，而在这之前Dimentional Profile是2.5D】，在CPU端Dimentional Profile是一份独立的NVDF【Nubis Volume Data Field】字段，其在GPU端对应的则是编码后的体积纹理。

VDB的云生成部分则是影视工作流的主场了，作者的思路是在Houdini中使用自己编写的流体模拟器来“生长”云，这个小型求解器被称为“Aero”，是Andrew 2014年之前还在Blue Sky Studio时开发的。

Aero相比较于复杂的燃烧式求解器要精简很多，因为它设计之初只用于模拟云和烟雾在空气中的升腾作用

解决了生成模拟，还剩下编辑工具。针对这部分需要又开发了“Atlas Tool”【也是基于Houdini做的】，Atlas提供多个单体云的合成【可以是体素、点云、Mesh和之前的2.5D作家系统数据】，也提供对生成云的二次编辑比如挤压，侵蚀调整和打洞。在云上打洞可以将云体本身做为一个独立场景，进而带来玩法上的启发【后面会介绍】。

网格转VDB也是一种重要的构建单体云的思路，《穿靴子的猫》中就使用了这种流程，视效部门关注完整的VDB云，而其余部门只需要使用占位用Mesh Proxy。Atlas Tool提供Mesh转VDB，而VDB又可以生成网格距离场，这样美术工作者可以选择性的将任何场景元素转化为体积云：

这里额外扩展一下，直接Mesh生成网格距离场也可以，通常使用的算法是“遍历体素化”和“3D跳洪算法”。

1.2 分布与大型

云在场景中的分布规模目前还没有很好的解决，空间体素化通常意味着规模不能太大。数据的管理是个很大的问题，特别是基于光线步进的方案中这种循环内频繁的数据读取与指令计算的切换对内存并不友好。所以能够发现大多数使用体素的体积渲染方案都是局部体素化，无论是视锥体素还是世界空间的局部体素。对于DLC来说，场景规模在4km * 4km，高度范围为500m。考虑体素网格还需要有足够的精度来捕捉玩家飞过云层的细节【游戏中飞行座骑的翼展越2m】。体素网格被控制在2048 * 2048 * 256，这样勉强够用：

看到这里你肯定会想，为什么不用稀疏体素呢？答案是时间不够了，当前内存占用如下：

为了降低性能，体素内存储的并不是最终的密度，而是Dimentional Profile数据。

最终密度还是采用3D噪声侵蚀Profile，Profile的另一个重要作用是做为一种“概率密度场”表征照明，后文会提到这部分。除了Dimentional Profile，还有Detail Type和Density Type两个重要NVDF字段，它们三个共同组成了体积云的Modeling Data。简单说明下作用：

• Dimentional Profile：构建大型，提供梯度信息

• Detail Type：描述Billow和Wispy两种细节形态在云体上的分布，后面的1.3小节会介绍这两种细节特征

• Density Type：提供密度调制

作为Modeling Data的NVDF‘s对应GPU端的数据是“VoxelCloudModelingDataTexture”，这是一个RGB三通道体积纹理，R-mDimentionalProfile, G-mDetailType, B-mDensityScale，使用BC6压缩，分辨率为512 * 512 * 64，每纹素占用一个字节，总体大概16.777Mib:

这样体素步进过程中采样Dimentional Profile，返回非零点意味着接触云，则使用3D噪声侵蚀并计算照明。

1.3 细节侵蚀

细节侵蚀部分需要考虑云的两种升腾作用，一种是向上升腾过程中遇到冷空气时会受到一种反向的挤压力，这形成了云上部波涛汹涌的结构，我们称之为Billow；另外一种是云体向低密度空间散发形成的一缕缕的絮状结构，我们称之为Wispy。云的升腾过程会同时受到两种效应的影响，因此云体会同时存在Billow与Wisp两种特征，并且存在一定方向性。为了在新的单体云上模拟这种 效用，重新设计了侵蚀用细节噪声。它是一个4通道3D纹理，分辨率128 * 128 * 128，RGBA分别存储R：Low Freq "Curl-Alligator", G:High Freq "Curl-Alligator", B:Low Freq "Alligator", A: High Freq "Alligator"。

这里解释下为这么要这样设计，我们已经知道云体表面波涛汹涌的细节可以使用Inv-Worley做出来，但一份Inv-Worley噪声只能侵蚀出“扁球形”的单一细节特征，通常我们会使用三个频率倍增的Inv-Worley噪声按强度递减规律构建分型噪声去侵蚀。但为了尽可能减少体积采样，我们希望将细节侵蚀用的噪声全部塞进一张Texture。于是这里使用了一种被称为“Alligator”的噪声，它的信号特征足够丰富而不需要再特意构建分型。Alligator的生成算法可参考Houdini，这里给出代码：

/*

 * Copyright (c) 2023

 *      Side Effects Software Inc.  All rights reserved.

 *

 * Redistribution and use of Houdini Development Kit samples in source and

 * binary forms, with or without modification, are permitted provided that the

 * following conditions are met:

 * 1. Redistributions of source code must retain the above copyright notice,

 *    this list of conditions and the following disclaimer.

 * 2. The name of Side Effects Software may not be used to endorse or

 *    promote products derived from this software without specific prior

 *    written permission.

 *

 * THIS SOFTWARE IS PROVIDED BY SIDE EFFECTS SOFTWARE `AS IS' AND ANY EXPRESS

 * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES

 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN

 * NO EVENT SHALL SIDE EFFECTS SOFTWARE BE LIABLE FOR ANY DIRECT, INDIRECT,

 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT

 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,

 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF

 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING

 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,

 * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

 *

 *----------------------------------------------------------------------------

 */

/// Alligator Noise is provided by Side Effects Software Inc. and is licensed

/// under a Creative Commons Attribution-ShareAlike 4.0 International License.

///

/// {

///  "dct:Title"  : "Alligator Noise",

///  "dct:Source" : "http://www.sidefx.com/docs/hdk15.0/alligator_2alligator_8_c-example.html"

///  "license"    : "http://creativecommons.org/licenses/by-sa/4.0/",

///  "cc:attributionName" : "Side Effects Software Inc",

/// }

///

///

/// This code is intended to be reference implementation of Houdini's Alligator

/// Noise algorithm.  It is not an optimal implementation by any means.

///

/// @note When using the Houdini's HDK, it's much easier to simply call @code

/// #include

/// static UT_Noise     alligatorNoise(0, UT_Noise::ALLIGATOR);

/// static double

/// alligator(const UT_Vector3 &pos)

/// {

///    return alligatorNoise.turbulence(pos, 0);

/// }

/// @endcode

#include

#include

#include

#include

#include

#include

#include

namespace HDK_Sample {

static double

frandom(std::size_t hash)

{

#if 1

    // For VEX, we need something that produces the same results in a

    // thread-safe fashion (thus random() and even random_r() are not

    // possibilities).

    uint        seed = hash & (0xffffffff);

    return SYSfastRandom(seed);

#else

    srandom(hash);

    return double(random()) * (1.0/RAND_MAX);

#endif

}

static double

hash(int ix, int iy, int iz)

{

    std::size_t hash = 0;

    hboost::hash_combine(hash, ix);

    hboost::hash_combine(hash, iy);

    hboost::hash_combine(hash, iz);

    return frandom(hash);

}

// These hash functions permute indices differently to make different values

static double hash_1(int ix, int iy, int iz) { return hash(ix, iy, iz); }

static double hash_2(int ix, int iy, int iz) { return hash(iy, iz, ix); }

static double hash_3(int ix, int iy, int iz) { return hash(iz, ix, iy); }

static double hash_4(int ix, int iy, int iz) { return hash(ix, iz, iy); }

static double

rbf(double d)

{

    // Radial basis function

    auto smooth = [](double d) { return d*d*(3 - 2*d); };

    return d < 1 ? smooth(1-d) : 0;

}

static void

getCenter(int ix, int iy, int iz, const int ipos[3], double center[3])

{

    // Compute the center point for the given noise cell

    // hash_* is a function which takes the seeds and returns a random double

    // between 0 and 1.

    center[0] = hash_1(ix+ipos[0], iy+ipos[1], iz+ipos[2]) + ix;

    center[1] = hash_2(ix+ipos[0], iy+ipos[1], iz+ipos[2]) + iy;

    center[2] = hash_3(ix+ipos[0], iy+ipos[1], iz+ipos[2]) + iz;

}

double

noiseValue(int ix, int iy, int iz, const int ipos[3])

{

    return hash_4(ix+ipos[0], iy+ipos[1], iz+ipos[2]);

}

static double

distance(const double *p0, const double *p1)

{

    // Distance from one point to another

    double      sum = 0;

    for (int i = 0; i < 3; ++i)

        sum += (p0[i]-p1[i]) * (p0[i]-p1[i]);

    return SYSsqrt(sum);

}

static double

alligator(const double pos[3])

{

    int                         ipos[3]; // Integer coordinates

    double                      fpos[3]; // Fractional coordinates

    double                      vpos[3]; // Sparse point

    int                         idx = 0;

    std::priority_queuenvals;

    for (int i = 0; i < 3; ++i)

    {

        ipos[i] = floor(pos[i]);

        fpos[i] = pos[i] - ipos[i];

    }

    for (int ix = -1; ix <= 1; ++ix)

    {

        for (int iy = -1; iy <= 1; ++iy)

        {

            for (int iz = -1; iz <= 1; ++iz, ++idx)

            {

                getCenter(ix, iy, iz, ipos, vpos);

                double  d = distance(fpos, vpos);

                if (d < 1)

                {

                    // Scale value by noise associated with the point.

                    double      v = noiseValue(ix, iy, iz, ipos) * rbf(d);

                    nvals.push(v);

                }

            }

        }

    }

    if (!nvals.size())

        return 0;

    double max = nvals.top();

    nvals.pop();

    if (nvals.size())

        max -= nvals.top();

    return max;

}

/// VEX callback to implement: float alligator(vector pos);

static void

alligator_Evaluate(int, void *argv[], void *)

{

    auto result = (float *)argv[0];

    auto fpos = ((const UT_Vector3 *)argv[1])->data();

    double      pos[3];

    std::copy(fpos, fpos+3, pos);

    *result = alligator(pos);

}

}

using namespace HDK_Sample;

void

newVEXOp(void *)

{

    new VEX_VexOp("alligator@&FV",

                    alligator_Evaluate,

                    VEX_ALL_CONTEXT,

                    nullptr,

                    nullptr);

}

[ 上下滑动查看代码 ]

另一个观察规律是在云体边缘往往需要较低频的信号，而密度高的地方需要更为高频的信号来塑造细节，因此需要高低频两张Alligator噪声去侵蚀不同密度处的形体。现在也能做到了，由于Mesh SDF自带从外到内的梯度信息，所以很容易就能控制中心较为致密的区域使用High Freq Alligator，而外围较低松散区域使用Low Freq Alligator。

RG通道解决了Billow特征，BA同理来解决Wispy效应。Perin-InvWorley经过验证只适合塑造大型体但不适合用来构建这种一缕缕的信号分布，这里使用的Curl Noise去扭曲Inv-Alligator Noise。其实思路类似，Curl Noise包含了一缕缕的信号特征，Alligator包含了良好的丛聚状分布，两者结合一下就能用来描述Wispy。同样要区分高低频。这里的一点优化是根据距离选择不同MipLevel来采样Detail Noise，距离越远，Level越高。根据Andrew的说法有效控制Mip层级能节省大概15%的性能。

搞清楚侵蚀用的数据后，我们来看下侵蚀算法，它被描述为“Up-prez"。通过调用GetUprezzedVoxelCloudDensity函数实现，代码如下：

float GetUprezzedVoxelCloudDensity(CloudRenderingRaymarchInfo inRaymarchInfo, float3 inSamplePosition, float inDimensionalProfile, float inType, float inDensityScale, float inMipLevel, bool inHFDetails)

{

    // Step1-Apply wind offset ：注意只有水平偏移

    inSamplePosition -= float3(cCloudWindOffset.x, cCloudWindOffset.y, 0.0) * voxel_cloud_animation_speed;

    // Step2-Sample noise ：注意MipLevel层级随距离而增加，根据Andrew的说法有效控制Mip层级能节省大概15%的性能

    float mipmap_level = GetVoxelCloudMipLevel(inRaymarchInfo, inMipLevel);

    // 4通道RGBA体积纹理，分辨率128 * 128 * 128，未压缩尺寸为 2Bytes/Texel

    // R：Low Freq "Curl-Alligator", G:High Freq "Curl-Alligator", B:Low Freq "Alligator", A: High Freq "Alligator"

    float4 noise = Cloud3DNoiseTextureC.SampleLOD(Cloud3DNoiseSamplerC, inSamplePosition * 0.01, mipmap_level);

    // Step3-Define Detail Erosion ：一个观察规律是高度密度区域需要较高频细节，低密度区需要较低频细节

    // 由于DimensionalProfile【准确得说是mProfile】已经包含了云从表面到内部的密度梯度变化【密度逐渐增大】

    // 定义wispy噪声

    float wispy_noise = lerp(noise.r, noise.g, inDimensionalProfile);

    // 定义billowy噪声

    float billowy_type_gradient = pow(inDimensionalProfile, 0.25);

    float billowy_noise = lerp(noise.b * 0.3, noise.a * 0.3, billowy_type_gradient);

    // 定义High Noise composite - blend to wispy as the density scale decreases.

    float noise_composite = lerp(wispy_noise, billowy_noise, inType);

    // 是否启用超高频细节

    if(inHFDetails)

    {

        // 使用采样获得的高频噪声调制出更高频的噪声，以弥补从云体上表面飞过摄像机非常靠近云时细节信息的缺乏

        float hhf_wisps = 1.0 - pow(abs(abs(noise.g * 2.0 - 1.0) * 2.0 - 1.0), 4.0);

        float hhf_billows = pow(abs(abs(noise.a * 2.0 - 1.0) * 2.0 - 1.0), 2.0);

        float hhf_noise = saturate(lerp(hhf_wisps, hhf_billows, inType));

        float hhf_noise_distance_range_blender = ValueRemap(inRaymarchInfo.mDistance, 50.0, 150.0, 0.9, 1.0); // 50，150的单位是m

        noise_composite = lerp(hhf_noise, noise_composite, hhf_noise_distance_range_blender);

    }

    // 组合噪声

    float uprezzed_density = ValueErosion(inDimensionalProfile, noise_composite);

    float powered_density_scale = pow(saturate(inDensityScale), 4.0);

    uprezzed_density *= powered_density_scale;

    // 这一行是一个技巧性做法，这样的power操作可以使得低密度区进一步锐化呈现更清晰的结果

    uprezzed_density = pow(uprezzed_density, lerp(0.3, 0.6, max(EPSILON, powered_density_scale)));

    if(inHFDetails)

    {

        float distance_range_blender = GetFractionFromValue(inRaymarchInfo.mDistance, 50.0, 150.0);

        uprezzed_density = pow(uprezzed_density, lerp(0.5, 1.0, distance_range_blender)) * lerp(0.666, 1.0, distance_range_blender);

    }

    // 返回带有软边缘的最终噪声

    return uprezzed_density;

}

[ 上下滑动查看代码 ]

其中GetVoxelCloudMipLevel用来评估细节噪声纹理的采样MipLevel，代码中的Cloud3DNoiseTextureC就是上述介绍的4通道细节噪声，也能够看到依据DimentionalProfile插值高低频Billow, Wispy噪声的过程。这里有额外一点设计需要说明，代码中的HFDetail分支是为了解决超近视距下细节不够的缺陷设计的，但为了避免采样次数这里并没有使用更高频的噪声，而是通过“数值折叠”的方式合成频率更高的噪声：

然后类似原始细节噪声那样使用mDetailType【Modeling Data的三个NVDF字段之一】进行插值。这里的近视距被认定为50m之类，50~150m使用high Detail Noise与high high Detail Noise的插值，150m以外使用采样直接得到high Noise依据mDetailType进行插值。得到细节噪声后，侵蚀Profile大型，同SIG 22一致也是使用ValueErosion实现的。返回侵蚀后的细节后再乘上mDensityScale得到最终细节。

1.4 全局距离场加速

加速Raymarching计算的核心有两个，一是尽可能减少0密度点采样，而是使用可变步幅步进。Nubis 3中同样使用了2种手段，关于自适应步幅我们1.5再议，先看第一点如何解决。

现在想象体素网格内步进的情景，由于暂未引入稀疏结构加速，体素域内的步进类似于世界空间步进。Raymarching的过程是步进采样，然后判断是否返回非零密度值，如果非零意味着有可能接触到云体，那么进一步采样细节噪声并计算光照。光照计算涉及二次子步进的过程，这方面也有很多优化手段我们后面再说明。虽然基于Profile的思路已经尽可能减少3D采样次数了，但为了判断是否进入云内，Profile本身的采样还是无法避免的。虽然一些特定的步进过程可以省去这一步，比如SIG 22中的地形云使用数值计算粗略判断大型范围等，但对于大多数情况Profile采样牺牲是必须的。问题就在于如何找到云这一步骤上，假如能够一步到位的定位到云表面位置，那么就可以最大程度减少Profile采样。于是距离场技术再次被请了出来，通过在编辑过程将多个Mesh SDF合成Global SDF , GPU仅通过一次采样就能找到距离云表面最近的点：

一旦艺术家完成云景的编辑，就可以生成一个全局的3D SDF。它被做为一个单独的NVDF【Field Data NVDF】, 这主要考虑到距离场对精度的要求需要重新设计压缩与解编码算法。Field Data NVDF将-256m~4096m的距离范围重新映射到0~1，距离云体越远值越大，云表面为0。它在GPU对应的状态是EncodedScalarTexture的3D Texture，分辨率为512 * 512 * 64。并且使用自定义的BC1压缩，而前面介绍的Modeling NVDF统一使用BC6压缩。区分开来是考虑距离场对精度的要求较高，使用Hugh Malan开发的压缩算法，提供16位精度的方式压缩SDF。经过测试，与未压缩的纹理相比可节省10~30%。压缩算法如下：

RGB565Color gConvertScalarToDXT1Endpoint_RoundDown(float inScalar, rcVec3 inUnpackDot)

{

    // gDXT1_5BitChannelValues and gDXT1_6BitChannelValues are lookup tables holding the exact value of each of the 32 R/B levels, and 64 G levels.

    int index_r = sFindLowerBoundIndex(gDXT1_5BitChannelValues, inScalar / inUnpackDot.mX);

    float val_r = gDXT1_5BitChannelValues[index_r] * inUnpackDot.mX;

    float residual_after_r = inScalar - val_r;

    float query_g = residual_after_r / inUnpackDot.mY;

    int index_g = sFindLowerBoundIndex(gDXT1_6BitChannelValues, query_g);

    float val_g = gDXT1_6BitChannelValues[index_g] * inUnpackDot.mY;

    float residual_after_g = residual_after_r - val_g;

    float query_b = residual_after_g / inUnpackDot.mZ;

    int index_b = sFindLowerBoundIndex(gDXT1_5BitChannelValues, query_b);

    return gCreate565Color(index_r, index_g, index_b);

}

[ 上下滑动查看代码 ]

基本思想是将 16 位值表示为 565 RGB 颜色：即红色通道具有最高有效的 5 位，G 具有接下来的 6 位，B 具有最低有效的 5 位。在运行时，着色器将对纹理进行双线性采样，并使用 float3(1.0, 1.0/32, 1.0/(32*64)) 的点积将 RGB 重新转换为标量：

着色器采样编码后的距离场功能被封装在GetVoxelCloudDistance函数中，实现如下：

float GetVoxelCloudDistance(float3 inSampePosition)

{

    // Decompress BC1

    // 通过减少频繁访问的体素数据的内存瓶颈，与未压缩的纹理相比可以节省 10-30%

    float3 sampled_color = EncodedScalarTexture.Sample(BilinearSampler, uv).rgb;

    float result = dot(sampled_color, float3(1.0, 0.03529415, 0.00069204));

}

[ 上下滑动查看代码 ]

这样一来，加上用于造型的3个Modeling NVDF字段，一共需要4个NVDFs，分别以BC6、Custom BC1方式压缩传给GPU。

1.5 自适应Raymarching

加速Raymarching的另一个思路，可能是更为主流的思路是设计可变步幅或者叫自适应步幅。

与以往不同，Nubis3中的Raymarching使用球面步进来判断与表面的交点，对于云内使用步长逐渐增大的自适应步幅，并且对于250m以内使用Jitter扰动方向来减少“切片式伪影”。整套Raymarching实现过程可参考如下代码：

// Raymarching结构

void RaymarchVoxelClouds(CloudRenderingGlobalInfo inGlobalInfo, float2 inRaySegment, inout CloudRenderingPixelData ioPixelData)

{

    // 初始化 raymarching 信息结构并设置到射线段开头的起始距离

    CloudRenderingRaymarchInfo raymarch_info;

    raymarch_info.mDistance = inRaySegment[0];

    // 主步进循环

    while(ioPixelData.mTransmittance > view_ray_transmittance_limit && raymarch_info.mDistance < min(inRaySegment[1], inGlobalInfo.mMaxDistance))

    {

        // 定位到摄像机方向上的云表面位置

        float3 sample_pos = inGlobalInfo.mCameraPosition.xyz + inGlobalInfo.mViewDirection * raymarch_info.mDistance;

        // 云体内使用自适应步幅（最小不小于1m）

        float adaptive_step_size = max( 1.0, max(sqrt(raymarch_info.mDistance), EPSILON) * 0.08);

        // 采样距离场，找到最近体素云表面距离

        raymarch_info.mCloudDistance = GetVoxelCloudDistance(sample_pos);

        // 取有符号距离和自适应步长中的最大值，这可以保证云体内时使用自适应步长来做为步进步幅

        raymarch_info.mStepSize = max(raymarch_info.mCloudDistance, adaptive_step_size);

        // 基于步长的抖动来减少摄像机方向上切片式伪影

        float jitter_distance = (raymarch_info.mDistance < 250.0 ? inGlobalInfo.mJitteredHash : inGlobalInfo.mStaticHash) * raymarch_info.mStepSize;

        sample_pos += inGlobalInfo.mViewDirection * jitter_distance;

        // 进入云内

        if (raymarch_info.mCloudDistance < 0.0)

        {

            // 获取采样数据

            VoxelCloudDensitySamples voxel_cloud_sample_data = GetVoxelCloudDensitySamples(raymarch_info, sample_pos, inMipLevel, inHFDetails);

            // 积分方程计算云照明

            // 前两个LightingPos采用传统向着光源的阴影射线步进，从第三个开始采样体素化光场中预计算的数据【256 * 256 * 32，分辨率是密度体素网格的1/8】

            IntegrateCloudSampleData(voxel_cloud_sample_data, inGlobalInfo, raymarch_info, ioPixelData, false);

        }

        // 移动至下一步

        raymarch_info.mDistance += raymarch_info.mStepSize;

    }

}

[ 上下滑动查看代码 ]

其中GetVoxelCloudDistance是上文提到了获取全局距离场的函数，GetVoxelCloudDensitySamples是密度模型的实现封装，包括大型构建与细节侵蚀，IntegrateCloudSampleData是针对体素渲染的散射积分方程。在本文第四部分可以看到完整的实现代码。但在这之前，先让我们补齐最后一环：照明实现。

二、体素云照明

体素云照明部分要考虑两个问题，一是如何构建散射积分方程（是单散射还是多散射？），而是如何优化LightingMarching过程，这可以优化非常多的性能。LightingMarching的一个重大优化是只针对前2步进行传统射线步进，后面则采样预存在Lighting Voxel Grid中的数据：

由此引入第二个体素场景，第一个是用来构建密度场的体素网格分辨率如前文提到的被设计为2048 * 2048 * 256，而用于照明的体素网格分辨率为密度场网格的1/8，只有256 * 256 * 32。

引入Lighting Voxel 的另一个好处是便于渲染长距离的云间阴影【在这之前使用了10个LightSample】，省去8个样本的误差当然是有的，但实践证明在可接受范围内，对比如下：

在散射实现方面，与SIG 22类似，使用DimentionalProfile做为概率密度场模拟内散射，使用Beer定律模拟吸收与外散射，同样基于DimentionalProfile做为概率分布拟合环境光：

// 使用轮廓剖面做为散射概率场

ms_volume = dimensional_profile;

// 获取云景的距离场

cloud_distance = GetVoxelCloudDistance(inSamplePosition);

// 吸收与外散射评估，inSunLightSummedDensitySamples是预先计算的DensityToSun

ms_volume *= exp(-inSunLightSummedDensitySamples * Remap(sun_dot, 0.0, 0.9, 0.25, ValueRemap(cloud_distance, -128.0, 0.0, 0.05, 0.25)));

// 同2.5D云一样，也使用维度剖面来近似环境光

float ambient_scattering = pow(1.0 - dimensional_profile, 0.5);

// 朝向天光的环境光

float ambient_scattering = pow(1.0 - dimensional_profile, 0.5) * exp(-summed_ambient_density);

[ 上下滑动查看代码 ]

同主光源累积类似，环境光也提供类似的预存手段：

有了主光与环境光，还要考虑局部光源【主要是点光】的支持，因为新的Model同样要满足雷暴模拟的需要。局部光源同样使用概率场近似的思路，使用球形定义势能变化，然后将其乘以密度来伪造非均匀的衰减：

// 第二光源支持：主要是做为局部光源的点光源

float potential_energy = pow( 1.0 - (d1 / radius), 12.0);

float pseudo_attenuation = (1.0 - saturate(density * 5.0));

float glow_energy = potential_energy * pseudo_attenuation;

通过引入Lighting Voxel 解耦Raymarching的手段，将帧占用降低了40%。

三、考虑玩法的雷暴云

这一节介绍新的雷暴云实现，SIG 22中Andrew花了大量的篇幅介绍如何基于Vertical Profile构建连接两个云层的Supercell，以及如何进行照明设计。现在，3D Profile至少在造型上更容易实现，而且新的流程能够在云体上挖洞。这意味着Single Cloud本身就能做为一个场景存在，这种新的体验可以引出很多有趣的玩法。比如世界观设计师Elijah Houck设想了这样一个情景：风暴机械鸟做为大型的顶级掠食者而存在，剧情设定中机械鸟会吸引玩家进入并探索它的巢穴——一大片风暴云，而玩家可以骑在坐骑的背上穿梭于云内外，当玩家杀死Boss后雷暴云会消失，天空放晴。这意味着要设计云的外型，内部结构，内外照明（包括闪电）等等，还要考虑如何实现云的消散过程。

Aero流体模拟器和Atlas Tool的存在能够让设计师准确控制云的造型以满足雷暴云的典型特征，并在云体上编辑洞穴结构以引导玩家的前行路线。

每当闪电的能量蕴蓄爆发时，就能照亮内部的结构，整套灯光列表系统的刷新与控制来源于SIG 22 中的Supercell。同样的，在云内部也能看到外部的光照变化。

四、特性支持与代码实现

4.1 特性支持度

Voxel Cloud Model，Vertical Profile与Envelope Model的特性支持与扩展性对比如下：

基于3D Mesh SDF和Voxel的模型相较于垂直剖面与包络模型更灵活，可扩展性也更好。更为重要的是这种思路是实时渲染视角下传统思路向全Voxel方案【实时VDB】过渡的中间版本。

4.2 代码参考

/////////////////////////////////////////////////  SIGGRAPH 2023 Nubis 3 ProtoType Code //////////////////////////////////////////////////////////////

// SDF压缩算法

RGB565Color gConvertScalarToDXT1Endpoint_RoundDown(float inScalar, rcVec3 inUnpackDot)

{

    // gDXT1_5BitChannelValues and gDXT1_6BitChannelValues are lookup tables holding the exact value of each of the 32 R/B levels, and 64 G levels.

    int index_r = sFindLowerBoundIndex(gDXT1_5BitChannelValues, inScalar / inUnpackDot.mX);

    float val_r = gDXT1_5BitChannelValues[index_r] * inUnpackDot.mX;

    float residual_after_r = inScalar - val_r;

    float query_g = residual_after_r / inUnpackDot.mY;

    int index_g = sFindLowerBoundIndex(gDXT1_6BitChannelValues, query_g);

    float val_g = gDXT1_6BitChannelValues[index_g] * inUnpackDot.mY;

    float residual_after_g = residual_after_r - val_g;

    float query_b = residual_after_g / inUnpackDot.mZ;

    int index_b = sFindLowerBoundIndex(gDXT1_5BitChannelValues, query_b);

    return gCreate565Color(index_r, index_g, index_b);

}

// Modeling NVDF's

// BC6 压缩，RGB, 分辨率为512 * 512 * 64】, 1 Byte / Texel，16.777Mb

// R: Dimentional Profile

// G: Detail Type

// B: Density Scale

Texture3D VoxelCloudModelingDataTexture;

SamplerState VoxelDataSampler;

// Field Data NVDF

// DXT1/BC1压缩,分辨率为512 * 512 * 64, 距离范围从-256m~4096m重新钳制到0~1

Texture3D EncodedScalarTexture;

SamplerState BilinearSampler;

// Detail Noise

// 4通道RGBA体积纹理，分辨率128 * 128 * 128，未压缩尺寸为 2Bytes/Texel

// R：Low Freq "Curl-Alligator", G:High Freq "Curl-Alligator", B:Low Freq "Alligator", A: High Freq "Alligator"

Texture3D Cloud3DNoiseTextureC;

SamplerState Cloud3DNoiseSamplerC;

// 云密度建模需要的 Struct

struct VoxelCloudModelingData

{

    float mDimensionalProfile;

    float mDetailType;

    float mDensityScale;

}

// 云密度采样使用的 Struct

struct VoxelCloudDensitySamples

{

    float mProfile; // 大型密度场，也叫轮廓样本

    float mFull;    // 细节密度场，也叫完整样本

}

// 云的射线步进结构

struct CloudRenderingRaymarchInfo

{

    float mDistance;       // 光线行进距离

    float mCloudDistance;  // 到体素云表面的距离

    float mStepSize;       // 步进步幅

}

#define ValueRemapFuncionDef(DATA_TYPE) \

    DATA_TYPE ValueRemap(DATA_TYPE inValue, DATA_TYPE inOldMin, DATA_TYPE inOldMax, DATA_TYPE inMin, DATA_TYPE inMax) \

    { \

        DATA_TYPE old_min_max_range = (inOldMax - inOldMin); \

        DATA_TYPE clamped_normalized = saturate((inValue - inOldMin) / old_min_max_range); \

        return inMin + (clamped_normalized*(inMax - inMin)); \

    }

ValueRemapFuncionDef(float)

ValueRemapFuncionDef(float2)

ValueRemapFuncionDef(float3)

ValueRemapFuncionDef(float4)

float ValueErosion(float inValue, float inOldMin)

{

    // derrived from Set-Range, this function uses the oldMin to erode or inflate the input value. - inValues inflate while + inValues erode

    float old_min_max_range = (1.0 - inOldMin);

    float clamped_normalized = saturate((inValue - inOldMin) / old_min_max_range);

    return (clamped_normalized);

}

float GetFractionFromValue(float inValue, float inMin, float inMax)

{

    return saturate((inValue - inMin) / (inMax - inMin));

}

float GetVoxelCloudDistance(float3 inSampePosition)

{

    // Decompress BC1

    // 通过减少频繁访问的体素数据的内存瓶颈，与未压缩的纹理相比可以节省 10-30%

    float3 sampled_color = EncodedScalarTexture.Sample(BilinearSampler, uv).rgb;

    float result = dot(sampled_color, float3(1.0, 0.03529415, 0.00069204));

}

VoxelCloudModelingData GetVoxelCloudModelingData(float3 inSampePosition, float inMipLevel)

{

    VoxelCloudModelingData modeling_data;

    float3 Modeling_NVDF = VoxelCloudModelingDataTexture.SampleLOD(VoxelDataSampler, inSampePosition, inMipLevel).rgb;

    modeling_data.mDimensionalProfile = Modeling_NVDF.r;

    modeling_data.mDetailType = Modeling_NVDF.g;

    modeling_data.mDensityScale = Modeling_NVDF.b;

    return modeling_data;

}

float GetVoxelCloudMipLevel(CloudRenderingRaymarchInfo inRaymarchInfo, float inMipLevel)

{

    // Apply Distance based Mip Offset

    float mipmap_level = cUseVoxelFineDetailMipMaps ? log2(1.0 + abs(inRaymarchInfo.mDistance * cVoxelFineDetailMipMapDistanceScale)) + inMipLevel : inMipLevel;

    return mipmap_level;

}

float GetUprezzedVoxelCloudDensity(CloudRenderingRaymarchInfo inRaymarchInfo, float3 inSamplePosition, float inDimensionalProfile, float inType, float inDensityScale, float inMipLevel, bool inHFDetails)

{

    // Step1-Apply wind offset ：注意只有水平偏移

    inSamplePosition -= float3(cCloudWindOffset.x, cCloudWindOffset.y, 0.0) * voxel_cloud_animation_speed;

    // Step2-Sample noise ：注意MipLevel层级随距离而增加，根据Andrew的说法有效控制Mip层级能节省大概15%的性能

    float mipmap_level = GetVoxelCloudMipLevel(inRaymarchInfo, inMipLevel);

    // 4通道RGBA体积纹理，分辨率128 * 128 * 128，未压缩尺寸为 2Bytes/Texel

    // R：Low Freq "Curl-Alligator", G:High Freq "Curl-Alligator", B:Low Freq "Alligator", A: High Freq "Alligator"

    float4 noise = Cloud3DNoiseTextureC.SampleLOD(Cloud3DNoiseSamplerC, inSamplePosition * 0.01, mipmap_level);

    // Step3-Define Detail Erosion ：一个观察规律是高度密度区域需要较高频细节，低密度区需要较低频细节

    // 由于DimensionalProfile【准确得说是mProfile】已经包含了云从表面到内部的密度梯度变化【密度逐渐增大】

    // 定义wispy噪声

    float wispy_noise = lerp(noise.r, noise.g, inDimensionalProfile);

    // 定义billowy噪声

    float billowy_type_gradient = pow(inDimensionalProfile, 0.25);

    float billowy_noise = lerp(noise.b * 0.3, noise.a * 0.3, billowy_type_gradient);

    // 定义High Noise composite - blend to wispy as the density scale decreases.

    float noise_composite = lerp(wispy_noise, billowy_noise, inType);

    // 是否启用超高频细节

    if(inHFDetails)

    {

        // 使用采样获得的高频噪声调制出更高频的噪声，以弥补从云体上表面飞过摄像机非常靠近云时细节信息的缺乏

        float hhf_wisps = 1.0 - pow(abs(abs(noise.g * 2.0 - 1.0) * 2.0 - 1.0), 4.0);

        float hhf_billows = pow(abs(abs(noise.a * 2.0 - 1.0) * 2.0 - 1.0), 2.0);

        float hhf_noise = saturate(lerp(hhf_wisps, hhf_billows, inType));

        float hhf_noise_distance_range_blender = ValueRemap(inRaymarchInfo.mDistance, 50.0, 150.0, 0.9, 1.0); // 50，150的单位是m

        noise_composite = lerp(hhf_noise, noise_composite, hhf_noise_distance_range_blender);

    }

    // 组合噪声

    float uprezzed_density = ValueErosion(inDimensionalProfile, noise_composite);

    float powered_density_scale = pow(saturate(inDensityScale), 4.0);

    uprezzed_density *= powered_density_scale;

    // 这一行是一个技巧性做法，这样的power操作可以使得低密度区进一步锐化呈现更清晰的结果

    uprezzed_density = pow(uprezzed_density, lerp(0.3, 0.6, max(EPSILON, powered_density_scale)));

    if(inHFDetails)

    {

        float distance_range_blender = GetFractionFromValue(inRaymarchInfo.mDistance, 50.0, 150.0);

        uprezzed_density = pow(uprezzed_density, lerp(0.5, 1.0, distance_range_blender)) * lerp(0.666, 1.0, distance_range_blender);

    }

    // 返回带有软边缘的最终噪声

    return uprezzed_density;

}

VoxelCloudDensitySamples GetVoxelCloudDensitySamples(CloudRenderingRaymarchInfo inRaymarchInfo, float3 inSamplePosition, float inMipLevel, bool inHFDetails)

{

    // 初始化体素云密度模型结构体

    VoxelCloudModelingData modeling_data;

    // 初始化密度采样结构体

    VoxelCloudDensitySamples density_samples;

    // 获取包围盒内的采样位置，调用函数获取模型数据

    float3 sample_coord = (inSamplePosition - cVoxelCloudBoundsMin) / (cVoxelCloudBoundsMax - cVoxelCloudBoundsMin);

    modeling_data = GetVoxelCloudModelingData(sample_coord, inMipLevel);

    // Assign data

    float dimensional_profile = modeling_data[0];

    float type = modeling_data[1];

    float density_scale = modeling_data[2];

    // 如果剖面距离场数据大于0,则返回细节细节密度场

    // be visible in the reflections anyways.

    if (dimensional_profile > 0.0)

    {

        // 基于SDF的空间剖面 * 调制密度 = 最终大型密度

        density_samples.mProfile = dimensional_profile * density_scale;

        // 调用up-rez函数计算云的最终密度

        if(!cRenderDome)

            density_samples.mFull = GetUprezzedVoxelCloudDensity(inRaymarchInfo, inSamplePosition, dimensional_profile, type, density_scale, inMipLevel, inHFDetails)

                                    * ValueRemap(inRaymarchInfo.mDistance, 10.0, 120.0, 0.25, 1.0);

        else

            density_samples.mFull = density_samples.mProfile;

    }

    else        

    {

        density_samples.mFull = density_samples.mProfile = 0.0;

    }

    // Return result

    return density_samples;

}

// Raymarching结构

void RaymarchVoxelClouds(CloudRenderingGlobalInfo inGlobalInfo, float2 inRaySegment, inout CloudRenderingPixelData ioPixelData)

{

    // 初始化 raymarching 信息结构并设置到射线段开头的起始距离

    CloudRenderingRaymarchInfo raymarch_info;

    raymarch_info.mDistance = inRaySegment[0];

    // 主步进循环

    while(ioPixelData.mTransmittance > view_ray_transmittance_limit && raymarch_info.mDistance < min(inRaySegment[1], inGlobalInfo.mMaxDistance))

    {

        // 定位到摄像机方向上的云表面位置

        float3 sample_pos = inGlobalInfo.mCameraPosition.xyz + inGlobalInfo.mViewDirection * raymarch_info.mDistance;

        // 云体内使用自适应步幅（最小不小于1m）

        float adaptive_step_size = max( 1.0, max(sqrt(raymarch_info.mDistance), EPSILON) * 0.08);

        // 采样距离场，找到最近体素云表面距离

        raymarch_info.mCloudDistance = GetVoxelCloudDistance(sample_pos);

        // 取有符号距离和自适应步长中的最大值，这可以保证云体内时使用自适应步长来做为步进步幅

        raymarch_info.mStepSize = max(raymarch_info.mCloudDistance, adaptive_step_size);

        // 基于步长的抖动来减少摄像机方向上切片式伪影

        float jitter_distance = (raymarch_info.mDistance < 250.0 ? inGlobalInfo.mJitteredHash : inGlobalInfo.mStaticHash) * raymarch_info.mStepSize;

        sample_pos += inGlobalInfo.mViewDirection * jitter_distance;

        // 进入云内

        if (raymarch_info.mCloudDistance < 0.0)

        {

            // 获取采样数据

            VoxelCloudDensitySamples voxel_cloud_sample_data = GetVoxelCloudDensitySamples(raymarch_info, sample_pos, inMipLevel, inHFDetails);

            // 积分方程计算云照明

            // 前两个LightingPos采用传统向着光源的阴影射线步进，从第三个开始采样体素化光场中预计算的数据【256 * 256 * 32，分辨率是密度体素网格的1/8】

            IntegrateCloudSampleData(voxel_cloud_sample_data, inGlobalInfo, raymarch_info, ioPixelData, false);

        }

        // 移动至下一步

        raymarch_info.mDistance += raymarch_info.mStepSize;

    }

}

/////////////////////////////////////////////////   IntegrateCloudSampleData内的部分原型代码 //////////////////////////////////////////////////////////////

// 使用轮廓剖面做为散射概率场

ms_volume = dimensional_profile;

// 获取云景的距离场

cloud_distance = GetVoxelCloudDistance(inSamplePosition);

// 吸收与外散射评估，inSunLightSummedDensitySamples是预先计算的DensityToSun

ms_volume *= exp(-inSunLightSummedDensitySamples * Remap(sun_dot, 0.0, 0.9, 0.25, ValueRemap(cloud_distance, -128.0, 0.0, 0.05, 0.25)));

// 同2.5D云一样，也使用维度剖面来近似环境光

float ambient_scattering = pow(1.0 - dimensional_profile, 0.5);

// 朝向天光的环境光

float ambient_scattering = pow(1.0 - dimensional_profile, 0.5) * exp(-summed_ambient_density);

// 第二光源支持：主要是做为局部光源的点光源

float potential_energy = pow( 1.0 - (d1 / radius), 12.0);

float pseudo_attenuation = (1.0 - saturate(density * 5.0));

float glow_energy = potential_energy * pseudo_attenuation;

[ 上下滑动查看代码 ]

  互动有奖！