体素地形：物理
在上一篇文章中，我们讨论了 Roblox 中体素地形的体素数据定义与存储的具体细节。在此基础上，许多其他系统会从存储中读取和写入数据，并以不同的方式进行解释。每个系统（渲染、网络、物理）的实现都是完全独立的，并不依赖于存储或其他系统做出的任何决策，因此我们可以独立研究它们。
虽然从逻辑上讲，接下来研究网格生成器（mesher）似乎更为合理（这是我们对那个能够将一盒体素数据转换为带有材质属性的三角形数据——即地形表面——的组件的称呼），但由于该组件同时被物理和渲染系统所使用，其算法相当复杂且包含不少“魔法”，因此我们将留待以后再探讨，今天先来研究物理系统。
初始原型
物理支持对地形至关重要，因为 Roblox 上的大量内容都依赖于稳健的物理行为，无论是在游戏中还是在编辑器中。特别是对于地形而言，实现物理支持意味着需要实现地形与我们使用的所有其他形状之间的碰撞检测，并高效支持光线投射。我们非常关注性能和内存消耗，因为我们假设某些世界将高度依赖地形及其物理特性。
尽管我们的物理引擎是自研的，但在广相/窄相阶段（特别是针对复杂凸体和凸分解，依赖于Bullet的GJK实现及其他算法）仍使用了Bullet Physics的部分组件，因此从基于Bullet的解决方案开始原型开发是合乎逻辑的。
与体素存储类似，我们将整个世界划分为多个区块，并将每个区块表示为一个碰撞对象；这种划分至关重要，因为每个区块都将成为物理数据更新的基本单位。由于地形可能随时发生变化，区块大小需要在更新成本（如果区块过大，每次更新该区块内的体素时，就必须付出更新整个区块的高昂代价。 目前我们认为增量更新过于复杂难以实现，因为体素变化可能导致生成的网格拓扑结构发生改变），以及区块开销（如果一个区块包含 2^3 个体素，那么我们需要花费大量时间和内存来管理这些区块）之间的平衡。 为了在这些因素之间取得平衡，我们最终确定了 8^3 个区块（提醒一下，一个角色的高度略高于 1 个体素，这应该能让你对规模有个概念）。
虽然我们可以尝试直接利用底层体素表示进行碰撞检测，但我们使用的网格化算法较为复杂，且在不运行算法的情况下难以准确预测表面位置。由于我们的体素体积较大，为了消除视觉伪影，我们希望渲染和物理模拟的表示能高度一致，因此决定采用多边形表示进行碰撞检测。
因此，原型程序会处理每个块，对其运行网格生成器以生成三角网格，然后使用 btBvhTriangleMeshShape 创建一个 Bullet 碰撞对象。生成的对象会被插入到通用广相结构中，与世界中的其他对象一同存在。每当某个对象（例如一个球）与其中一个块发生相交时，我们会生成一个接触对象，然后运行 Bullet 的算法来确定接触点。
虽然这个原型让我们迈出了第一步，但也凸显了几个关键的改进方向：
- 广相数据过于粗略：每当任何物体与 8^3 体素块发生相交时，我们都需要对该接触点运行窄相算法；这导致了大量冗余的接触点，这些接触点实际上并未产生任何实际的接触点，尤其是在洞穴场景中，悬浮在空中的物体可能与相对较大的块边界框发生重叠，却从未触碰到任何几何体。
- 窄相位数据体积过大：对于每个块，我们最终需要存储三角网格（其紧凑性不如体素数据，因为需要存储顶点位置/索引），以及用于加速碰撞检测的 Bullet BVH 结构，后者同样相当庞大。
- 窄阶段数据的生成速度太慢：虽然我们的网格生成算法经过了深度优化，但 Bullet 对生成的网格进行处理的速度相对较慢（比生成网格慢达 3 倍）。
这表明我们需要调整广相和窄相的实现方案。让我们看看最终采取了哪些措施。
位宽阶段：构建
由于广相阶段的核心问题在于精度，我们决定让广相阶段了解每个物体的结构，并基于实际体素数据进行早期剔除。每当物体移动时，我们不再与每个重叠的块建立接触，而是先检查该块中的体素数据，以确定物体的AABB是否与任何体素数据相交。
尽管我们的体素数据相当紧凑，但我们仍希望将广相位数据的内存占用降至最低。此外，由于我们的网格生成算法特性，单个体素生成的几何体并不局限于该体素内部，可能会溢出到相邻体素中，因此我们实际上还需要检查相邻体素。 为了使这一切高效运行，我们决定为每个区块存储一个位掩码（其中每个体素对应一位），该掩码将告知我们每个体素是否存在可能发生碰撞的几何体。
能够不依赖网格生成算法来生成该掩码至关重要（因为对整个地形运行网格生成算法耗时过长，且根据我们代码的结构，必须为整个世界提供广相数据）， 因此我们采用近似方法：假设一个体素为实心，理论上它可能在其任何相邻体素（包括对角线相邻体素，共计 3^3=27 个体素）内部生成几何体，并将掩码中这些体素的位置全部设为 1；这一过程称为膨胀。 这同样要求我们考察该区块的邻近体素，因此我们的输入是一个 10^3 体素的盒子，输出是一个 8^3 的位掩码。
最后，我们有两种接触类型：固体和水（我们利用基元与水之间的接触来计算浮力），因此生成两个位掩码。该过程大致如下（请注意，为便于说明，本文中的图片假设使用的是4³个区块而非8³个）：
为了加快这个过程，我们使用位运算进行扩展。首先生成 10^3 位掩码，并将其存储在一个由 10^2 个 16 位整数组成的数组中，然后像这样对每个整数进行水平扩展：
data[y][z] = (data[y][z] << <span style="color: #ff0000;">1</span>) | data[y][z] | (data[y][z] >> <span style="color: #ff0000;">1</span>);
接着，我们分两轮沿另外两个轴进行扩展，具体如下：
temp[y][z] = data[y][z - <span style="color: #ff0000;">1</span>] | data[y][z] | (data[y][z + <span style="color: #ff0000;">1</span>];
最终，我们得到 10^2 个 10 位扩展位掩码。 接着，我们提取 8^2 个 8 位位掩码（对应块中的体素，并舍弃冗余的边界体素），并将它们存储为广相数据。在此过程中，我们还会过滤掉冗余的块。如果 10^3 体积内的所有体素都充满空气，则该处不存在几何体； 此外，如果 10^3 体积内的所有体素都充满固体材料或水，那么我们的网格生成算法将永远不会生成任何多边形，但我们仍需将该块保留在广相中（例如，以便判断物体是否完全浸没在水中），因此我们会以特殊方式对其进行标记。
总体而言，这带来了极低的存储开销。最坏情况是每体素2位（1位用于固体掩码，1位用于水掩码），但许多区块会被丢弃，因为它们要么是空的，要么是满的。而在这种情况下，我们只需存储区块结构，而不需要存储掩码。
位广相：重叠测试
现在我们已经生成了广域阶段数据（该操作在关卡加载时完成，且当块内任何体素发生变化时，该块的广域阶段数据会重新生成），我们可以探讨如何使用它。
每当物体移动时，我们会获取物体在旧位置和新位置的 AABB，查询广相数据以检测重叠，并进行接触管理。如果物体在旧位置与某个区块存在固体接触，但在新位置与该区块不存在固体接触，则可以移除该接触。 我们为每个物体-区块对最多追踪 1 个固体接触和 1 个水体接触，因此一个非常大的物体最终可能与地形产生多个接触点（而且每个接触点都可能因我们稍后将讨论的窄相处理而生成多个接触点）。
查询分为两步。首先，我们获取对象的 AABB，将其扩展 1 个体素以覆盖因溢出至相邻体素而与几何体发生接触的情况，将其投影到区块空间（通过将坐标除以 8 个体素），并将最小/最大值转换为整数（使用 floor/ceil 函数），从而得到区块范围。 接着，我们遍历每个区块并检查位数据，以判断对象的AABB是否触及任何标记为固体/水体的位。虽然后者可以通过查询对象AABB内的每个体素来实现，但我们利用位数据进行了如下优化。
请注意，每个区块包含 8^3 位；我们将该区块的每个 Y 切片（Y 轴向上）存储在一个 64 位整数中，其中每组连续的 8 位决定 X 行中的数据。因此，每个掩码都是一个简单的数组，我们分别存储一个用于表示固体和一个用于表示水的掩码：
<span style="color: #007788;">uint64_t</span> solid[<span style="color: #ff0000;">8</span>];
<span style="color: #007788;">uint64_t</span> water[<span style="color: #ff0000;">8</span>];
现在，我们获取物体的AABB，将其投影到区块空间并确定体素空间中的边界。然后，我们观察物体在与区块相交时的XZ方向范围，并生成一个64位掩码，其中物体在XZ平面上与体素相交的位置被设置为1：
然后，我们遍历对象覆盖的所有 Y 切片，并对每个切片进行简单的位测试：
<span style="color: #006699;">if</span> (chunk.solid[y] & mask)
touchesSolid = <span style="color: #336666;">true</span>;
这使我们能够通过一条简单指令（在 32 位架构上，此测试约需 3 条指令）检查一个数据块的整个 XZ 切片（多达 64 个体素！）是否存在重叠，这使得对相当大的对象进行精确查询变得非常高效。 为了优化小对象的开销，我们确保以最快速度生成 XZ 范围的掩码。虽然可以通过遍历 XZ 范围并设置位来实现，但我们注意到，该掩码实际上是代表一个垂直条带和一个水平条带的两个掩码的交集；对于每个方向，我们维护一个查找表，并将查找结果组合以获得最终的掩码：
<span style="color: #007788;">uint64_t <span style="color: #000000;">mask = masksVer[cmin.x][cmax.x] & masksHor[cmin.z][cmax.z];
</span></span>
窄阶段：分析与计划
既然广相数据已经整理妥当，接下来我们来看看窄相。从宏观层面来看，基于Bullet的窄相工作原理如下：
- 每个块存储一个顶点数组（每个顶点包含 3 个浮点数用于位置和 1 个字节用于材质）、一个索引数组（每个三角形有 3 个 16 位索引）以及一个 BVH（即 AABB 树）。
- 在构建树时，三角形列表会被反复拆分为节点，直到叶节点仅包含一个三角形
- 进行碰撞检测时，通过简单的AABB查询从树中提取一组三角形；每个三角形会与目标基元进行碰撞检测，方法是使用针对该组合的专用算法（如三角形-球体碰撞）或通用的GJK/EPA算法。这些碰撞产生的点会被输入到一个单纯形结构中，该结构最多保存4个接触点，并试图最大化接触面积。
- 在执行光线投射时，会运行一个简单的光线-AABB树查询；每个匹配的叶节点都会与光线进行相交运算。保留最近的点并将其作为结果返回。
我们决定保留碰撞检测算法和整体结构，但将所有其他组件替换为更适合我们用例的版本。为了降低内存开销，我们开始采用懒加载方式生成碰撞对象。我们仅在建立接触点或对块进行射线投射时才生成三角形/树数据，并缓存生成的对象，以确保内存开销保持在可控范围内。
这显著改善了窄阶段的内存消耗，但生成成本依然高得令人望而却步。 Bullet 的 BVH 三角网格树有两个版本：非量化版和量化版。两者均维护 AABB 树，但存储方式不同（一个使用 32 位浮点数，另一个使用 16 位整数）。由于每个三角形都有一个叶节点，且树是二叉树，因此所需的节点数量大约是三角形数量的两倍。
原始网格的基础内存开销约为每个顶点 13 字节（3 个浮点坐标和 1 字节材质数据），以及每个三角形 6 字节（用于 3 个 16 位索引），由于平均三角形数量是顶点数量的两倍，因此每个三角形总计约 12.5 字节。
未量化的 Bullet BVH 每个节点占用 64 字节，因此每个三角形约需 128 字节（节点结构本身仅需 44 字节，但被填充至 64 字节），这大约是三角形数据内存开销的 10 倍。此外，生成该树所需的时间约为生成网格（使用我们那套将体素数据转换为三角形数据的实现方案）的 3 倍。
量化版 Bullet BVH 每个节点占用 16 字节，折合每个三角形约 32 字节（约为三角形数据内存开销的 2.5 倍）。其生成速度比非量化版慢，与生成网格相比耗时约 5 倍。
鉴于这两种方案在内存占用和数据生成时间方面均表现不佳（由于采用延迟生成机制，构建时间至关重要），我们决定用自定义的 kD 树替代原有树结构，该树采用每个节点包含两个平面的设计。
松散 kD 树：构建
kD树是一种二叉树，每个节点沿一个轴将空间分割为两部分。通常kD树每个节点仅有一个分割平面，但这会导致处理那些未完全包含在某个子节点内的三角形变得复杂（必须用分割平面将其截断），因此我们采用了松散kD树。 每个节点沿同一轴线拥有两个分割平面，其中所有左子节点都包含在由其中一个平面定义的子空间内，所有右子节点都包含在由另一个平面定义的子空间内。这些平面紧密贴合子节点的内容，从而产生两种可能的平面配置：
虽然总体而言 AABB 树在空间定位方面比 kD 树稍好一些，但 kD 树的构建速度更快，而且正如我们稍后将讨论的那样，在递归遍历过程中可以恢复大部分信息。 其主要优势在于每个节点只需存储 2 个值，而非完整的 AABB。我们还决定在每个叶节点中存储多个三角形——通常情况下，存储 2 个三角形而非 1 个并不会显著影响查询质量——最终形成了如下结构：
<span style="color: #006bb8;">union</span> KDNode {
<span style="color: #006bb8;"> struct</span> {
<span style="color: #007788;"> float</span> splits[<span style="color: #ff0000;">2</span>];
<span style="color: #007788;"> unsigned int</span> axis: <span style="color: #ff0000;">2</span>; <span style="color: #999999;">// 0=X, 1=Y, 2=Z</span>
<span style="color: #007788;"> unsigned int</span> childIndex: <span style="color: #ff0000;">30</span>; <span style="color: #999999;">// children are at childIndex+0,1</span>
} branch;
<span style="color: #006bb8;">struct</span> {
<span style="color: #007788;"> unsigned int</span> triangles[<span style="color: #ff0000;">2</span>]; <span style="color: #999999;">// up to two triangles per leaf</span>
<span style="color: #007788;"> unsigned int</span> axis: <span style="color: #ff0000;">2</span>;<span style="color: #999999;"> // must be 3; same offset as branch.axis </span>
<span style="color: #007788;"> unsigned int</span> triangleCount: <span style="color: #ff0000;">30</span>;
} leaf;
};
树结构相对简单；所有节点都存储在一个大型数组中，以便我们可以通过索引来引用节点。轴索引的第4个值用作标记以区分叶节点，并且只存储一个子节点索引，因为每个分支节点总是同时拥有两个子节点。 我们在每个叶节点中最多存储 2 个三角形；虽然空间足以容纳 4 个，但存储 4 个三角形而非 2 个会使碰撞处理稍慢，因此我们选择了 2 个。
请注意，每个节点占用12字节。由于每个叶节点存储2个三角形，因此总节点数大约等于三角形总数，所以内存开销最终为12字节/三角形。此处显然还有进一步优化的空间。 我们可以通过使用16位整数压缩顶点位置来降低顶点数据的内存开销，并以类似方式优化kD树节点，从而将每个kD节点降至约6字节，最终导致每个三角形的总内存开销为3.5b顶点 + 6b索引 + 6b树 = 15.5b。但当前约24.5b/三角形的结果已足够满足发布需求。
树的构建过程相对标准。我们从一个三角形数组开始，对其进行递归细分，在此过程中生成分支节点（当每个节点包含的三角形数达到 2 个或更少时停止递归）。 每次划分时，我们通过选取当前 AABB 的最长轴来确定分割轴，将分割点置于该轴上所有三角形中点的平均位置，利用中点作为判定条件筛选该平面左侧/右侧的三角形，随后利用所有 3 个三角形顶点重新计算左右两侧的分割平面。 如果最终的分布在三角形数量上过于偏斜（目前定义为任一节点中的三角形占比低于25%），我们会重新划分，并将（按中点排序的列表中）一半的三角形放入一个子树，另一半放入另一个子树。这在空间连贯性和树深度之间保持了平衡，限制了树的深度，并确保了过程能够终止。
最终的构建过程比Bullet的方案更为简洁，实现效率也更高，因此比Bullet最快的树构建算法快约3-4倍。这使得网格和树数据的大小以及生成时间大致相当，达到了良好的平衡（或者说，这意味着若要显著提升整体流程的性能，必须对两者都进行大幅优化 :D）。
松散 kD 树：查询
如前所述，我们仅需两种查询：AABB查询（需收集给定AABB内的所有三角形，用于窄阶段）和射线投射查询（需收集与射线相交的所有三角形，或仅收集交点最近的那个）。 这两种查询均采用无栈遍历实现。由于树的深度有上限，因此很容易预先计算树的深度，并为给定的遍历预分配临时空间。无栈遍历未必能为我们节省大量空间或时间，但它有助于理解性能分析结果，因为分支和叶节点中所有遍历产生的开销都集中在同一个函数中，使得处理起来相对容易一些。
kD树不像AABB树那样完全掌握每个节点的范围信息，但我们可以在遍历过程中恢复这些信息。对于AABB查询，当遇到分支节点时，我们只沿AABB位于右半空间的分支向下遍历。由于采用了分层遍历，这最终只会遍历那些体积与AABB重叠的节点：
<span style="color: #006699;">if</span> (node.branch.splits[<span style="color: #ff0000;">1</span>] <= aabbMax[axis])
buffer[offset++] = childIndex + <span style="color: #ff0000;">1</span>; <span style="color: #999999;">// push right child</span>
<span style="color: #006699;">if</span> (node.branch.splits[<span style="color: #ff0000;">0</span>] >= aabbMin[axis])
buffer[offset++] = childIndex + <span style="color: #ff0000;">0</span>; <span style="color: #999999;">// push left child
</span>如果查询的 AABB 位于整个 kD 树边界之外，这种方法可能效率较低，因此我们为每个 kD 树存储一个 AABB 以实现早期拒绝。
对于射线查询，我们选择采用区段树遍历，其中区段由射线起点/方向以及参数 t 的两个边界 tmin 和 tmax 定义，这些边界包含每个节点所定义子空间内的所有点。当遇到分支时，我们需要将分割平面与射线进行相交运算（由于平面与轴对齐，此操作简单且快速），并调整后续遍历的 t 边界：
<span style="color: #007788;">float</span> sa = raySource[axis];
<span style="color: #007788;">float</span> da = rayDir[axis];
<span style="color: #007788;">float</span> t0 = (node.branch.splits[i0] - sa) / da;
<span style="color: #007788;">float</span> t1 = (node.branch.splits[i1] - sa) / da;
<span style="color: #006699;">if</span> (t1 <= rn.tmax)
buffer[offset++] = { childIndex + i1, max(t1, rn.tmin), rn.tmax };
<span style="color: #006699;">if</span> (t0 >= rn.tmin)
buffer[offset++] = { childIndex + i0, rn.tmin, min(t0, rn.tmax) };
与 AABB 查询类似，如果射线未与 kD 树的完整边界相交，遍历效率可能较低，因此我们通过将射线与存储的 AABB 进行相交运算来计算初始区间边界。
此外，为了加速仅需第一个点的射线查询，我们将遍历顺序调整为优先访问定义了沿射线方向更早出现子空间的分支（这会影响上文代码片段中的 i0/i1 索引）。 如果发现的交点位于光线方向上早于给定线段的最小点 t，则可以终止遍历。遗憾的是，由于浮点精度问题，我们需要在每个分支上稍微扩展线段。
通过这种方式，我们实现了针对 AABB 和光线投射的高效树查询。AABB 查询性能与 Bullet 实现相当，但光线投射查询更快；这既是因为我们对每个分支的工作量较少（仅需让光线与两个平面相交），也是因为若找到合适的交点，我们可以提前终止遍历。
未来工作
虽然我们构建的算法运行效果相当不错，但无疑还有进一步改进的空间。 窄阶段的内存消耗尚有优化空间；此外，我们目前在 kD 树中存储的是三角形。Christer Ericson 建议，若将四边形作为第一类基本图形存储，光线投射效率可提升至 2 倍，因为 Möller-Trumbore 算法能以极小的额外计算量处理四边形（我们的地形由四边形构建，其中部分为平面，因此该方案可行）。
我们尚未探索的一个领域是窄阶段中仍沿用自Bullet的部分。或许可以采用更优的算法来处理三角形凸碰撞或缩减接触点流形；此外，我们目前使用了一些权宜之计来处理内部边碰撞，而实际上可以生成关于每条边是外部还是内部的数据，并在生成碰撞点/法线时加以利用。
最后，我们在原型开发阶段曾探索过但未投入正式生产的功能是“不相连的地形区域”。每当修改体素数据时，我们的原型会执行基础连通性分析，并将所有不相连的体素区域转换为自由移动的对象。要对此进行实时碰撞计算，最可能的方法是用凸包近似其形状，尽管可能存在其他方案，我们日后或许会对此进行探索。
