第12章：图形数据结构

学习目标

1. 理解三角网格的基本数据结构：顶点-面表示、共享顶点的网格拓扑及其相比无关三角形集合的存储效率优势
2. 了解翼边/半边数据结构的设计思路，掌握其在动态网格编辑（细分曲面、模型简化）中的应用场景
3. 理解场景图的分层变换模型：父-子节点的仿射变换传递，及其在管理复杂三维场景中的核心作用
4. 掌握三类空间数据结构（包围体层次BVH、层次空间细分、均匀空间细分）的设计原理与加速求交的基本思路
5. 理解BSP树用于可见性排序的算法思想和多维数组平铺用于内存局部性优化的基本概念

某些数据结构在图形应用中反复出现——也许是因为它们处理了表面、空间和场景结构这些基本底层思想。本章讨论几个基本且互不相关的数据结构类别，是图形学中最常见和最有用的：网格结构（mesh structures）、空间数据结构（spatial data structures）、场景图（scene graphs）和平铺多维数组（tiled multidimensional arrays）。

对于网格，我们讨论存储静态网格和将网格传输到图形 API 的基本存储方案。我们还讨论翼边数据结构（winged-edge data structure，Baumgart, 1974）和相关的半边结构（half-edge structure），它们对于管理曲面细分变化时的模型（如细分曲面或模型简化）非常有用。尽管这些方法可推广到任意多边形网格，本章聚焦于三角形网格这一简单情况。

接下来介绍场景图数据结构。其各种形式在图形应用中无处不在，因为它在管理物体和变换方面极为有用。所有新版图形 API 都设计为良好支持场景图。

对空间数据结构，我们讨论三种在三维空间中组织模型的方法：包围体层次结构（bounding volume hierarchies）、层次空间细分（hierarchical space subdivision）和均匀空间细分（uniform space subdivision）——以及使用BSP 树进行隐藏面消除。这些方法同样用于几何体裁剪和碰撞检测。

最后介绍平铺多维数组。它最初为帮助需要将图形数据从磁盘换入的应用中的分页性能而开发，现在对机器的内存局部性至关重要——无论数组是否适合主存。

12.1 三角网格

大多数真实世界模型由共享顶点的三角形复合体组成，通常称为三角形网格（triangular meshes）、三角网格（triangle meshes）或三角形不规则网络（TINs, triangular irregular networks）。高效处理它们对许多图形程序性能至关重要。网格存储在磁盘和内存中，我们希望最小化存储。当网格通过网络或从 CPU 传输到图形系统时，消耗的带宽往往比存储更加珍贵。对于执行网格操作的应用程序——细分曲面、网格编辑、压缩等——对邻接信息（adjacency information）的高效访问至关重要。

三角形网格通常用于表示表面，因此网格不仅仅是无关三角形的集合，而是通过共享顶点和边互相连接、形成单一连续表面的三角形网络。这是关于网格的关键洞察：一个网格可以比相同数量的无关三角形集合更高效地处理。

三角形网格所需的最小信息是一个三角形集合（顶点的三元组）及其顶点的 3D 位置。但大多数程序需要能够在顶点、边或面上存储附加数据以支持纹理映射、着色、动画等操作。顶点数据是最常见的：每个顶点可以有材质参数、纹理坐标、辐照度——任何值在表面上变化的参数。这些参数在每个三角形上线性插值，以在整个网格表面定义连续函数。偶尔也需要在每条边或每个面上存储数据。

12.1.1 网格拓扑

网格的表面性质可被形式化为对网格拓扑（topology）的约束——三角形相互连接的方式，与顶点位置无关。许多算法只工作于具有可预测连通性的网格，或在此类网格上实现要容易得多。

最简单的、也是最严格的拓扑要求是表面为流形（manifold）。术语流形来自拓扑学数学领域；粗略讲，一个二维流形（2-manifold）是在每一处都像一个表面的对象。流形网格是"水密的"——无间隙，将表面内部空间与外部空间分开。正式的三个条件为：

1. 每条边恰好被两个三角形共享（或在边界处被一个三角形共享）。不能有三个或更多三角形共享同一条边。
2. 每个顶点周围的三角形形成一个单一的循环（扇）。不能存在两个或多个孤立的三角形组共享同一顶点。
3. 网格在每一点局部同胚于圆盘（locally homeomorphic to a disk）。这意味着在网格上的任何点放大看，都应像一块平整的圆盘。

网格的亏格（genus）指其"手柄"的数量。球面亏格为 0，环面（甜甜圈）亏格为 1。著名的欧拉公式（Euler's formula）关联顶点数 V、边数 E、面数 F 与亏格 G：

V − E + F = 2(1 − G)

对于闭合流形三角网格，每个三角形有 3 条边，每条边被 2 个三角形共享，因此 E ≈ 3F/2。代入欧拉公式得：V − 3F/2 + F = 2(1−G)，整理得 V ≈ F/2 + 2(1−G)。对于球面（G=0）：V ≈ F/2 + 2。这一不变量用于验证网格完整性和检测退化，是许多网格处理算法的基础。

欧拉公式的详细推导与验证

欧拉公式 V − E + F = 2(1 − G) 是拓扑学中最优美的结果之一。下面给出其推导思路：

步骤 1：平面化——将闭合曲面去掉一个三角形面，剩余部分可展开到平面上（想象把地球仪表面剥去一块后摊平）。展开后，顶点数 V' = V，边数 E' = E，面数 F' = F − 1。

步骤 2：三角剖分——在平面上，通过添加对角线将每个非三角形面剖分为三角形。每添加一条对角线，边数和面数各增加 1，V − E + F 不变。

步骤 3：逐步移除边界三角形——从外部边界开始逐层移除三角形。移除一个边界三角形有三种情况：

情况 A: 移除一个三角形，它贡献 1 条外部边 → E−1, F−1 (Euler 特征不变)
情况 B: 移除一个三角形，它贡献 2 条外部边 → V−1, E−2, F−1 (Euler 特征不变)
情况 C: 移除一个三角形，它贡献 3 条外部边 → V−2, E−3, F−1 (Euler 特征不变)

无论哪种情况，每移除一个三角形，V − E + F 的值保持不变。最终剩余一个三角形：V=3, E=3, F=1，得 V − E + F = 1。由于步骤 1 中去掉了一个面，原始曲面的 V − E + F = 2。

步骤 4：亏格的修正——以上推导针对球面（亏格 G=0）。具有 G 个"把手"（手柄、环柄）的曲面，相当于球面上附着了 G 个环柄。每个环柄的添加：增加了一个贯穿孔洞，打开两个圆形孔（各减少 1 个面？）。拓扑学严格计算表明，每增加一个亏格，V − E + F 减少 2。因此通式为：

V − E + F = 2(1 − G) = 2 − 2G

亏格 0/1/2 的实例验证

亏格 0——球面/四面体：四面体 V=4, E=6, F=4。V − E + F = 4 − 6 + 4 = 2 = 2(1−0) ✓

亏格 0——球面/立方体：立方体三角剖分后 V=8, E=18, F=12。V − E + F = 8 − 18 + 12 = 2 = 2(1−0) ✓

亏格 1——环面（甜甜圈）：一个典型的三角化环面有 V=16, E=48, F=32。V − E + F = 16 − 48 + 32 = 0 = 2(1−1) = 0 ✓

亏格 2——双环面（8 字形手柄）：双环面的 V − E + F = 2(1−2) = −2。例如三角化双环面可能 V=24, E=84, F=58 → 24−84+58 = −2 ✓

12.1.2 索引网格存储

最紧凑实用的存储方案是索引三角形网格（indexed triangle mesh）：

顶点数组:  Vertex vertices[n_v];  // 每个顶点的 (x,y,z) 位置 + 法线 + 纹理坐标
索引数组:  int indices[3 * n_t];  // 每三角形 3 个整数，指向 vertices[] 的索引

该表示消除了所有顶点冗余，对于 GPU 管线处理也高效。每个三角形由三个指向顶点数组的索引定义，一条边由共享两个顶点的两个三角形隐式定义。这种方案是 OBJ、PLY、glTF 等流行三维文件格式的基础。

一个好的网格数据结构应当合理紧凑，并允许对所有邻接查询进行常数时间回答——找到邻居的时间不应依赖于网格大小。我们将讨论三种网格数据结构：一种基于三角形、两种基于边。

12.1.3 三角邻接结构

最直接但臃肿的实现是显式存储所有关系：

Triangle {
    Vertex v[3]
    Edge   e[3]
}
Edge {
    Vertex   v[2]
    Triangle t[2]
}
Vertex {
    Triangle t[]   // 可变长度！
    Edge     e[]
}

这种方法存储了过多信息，且顶点中的可变长度数据结构效率低下。更好的方案是定义类接口来回答邻接查询，背后隐藏更高效的数据结构。我们可以只存储部分连通性信息，在需要时高效恢复其余信息。

精简版三角邻接结构：

Triangle {
    Triangle nbr[3];  // 相邻三角形，nbr[k] 与边 k 相对
    Vertex   v[3];
}
Vertex {
    // ... per-vertex data ...
    Triangle t;   // 任意一个邻接三角形
}

该结构中，三角形 t 的邻居三角形和顶点可直接获取。通过在三角形间移动，也可常数时间回答顶点的连通性查询。如果三角形 t 以顶点 v 为第 k 个顶点，则 t.nbr[k] 是围绕 v 的顺时针方向的下一个三角形。这导出了三角形邻接遍历（TrianglesOfVertex）算法：

TrianglesOfVertex(v) {
    t = v.t
    do {
        // 找到 i 使 t.v[i] == v  (常数时间，因为固定大小)
        t = t.nbr[i]
    } while (t != v.t)
}

该操作每次找下一个三角形是常数时间——虽然需搜索顶点在三角形顶点列表中的位置，但顶点列表大小固定（3），所以搜索也是常数时间。但搜索本身笨拙且需额外分支。

改进版：带边索引的三角邻接

通过存储指向相邻三角形特定边的指针，可以消除搜索：

Triangle {
    Edge   nbr[3];  // nbr[k] 对应于 v[(k+1)%3] 到 v[(k+2)%3] 的边
    Vertex v[3];
}
Edge {
    Triangle t;     // t 的第 i 条边
    int      i;     // 0, 1, 或 2
}
Vertex {
    // ... per-vertex data ...
    Edge e;          // 离开该顶点的任一条边
}

实践中，Edge 可以借用三角形索引的两个位来存储边索引 i，因此总存储需求不变。该结构维护一个重要的数据结构不变式（invariant）：对于任意三角形 t 的第 j 条边：

t.nbr[j].t.nbr[t.nbr[j].i].t == t

这个不变式保证从三角形 t 的边 j 出发→走到邻居三角形→再走回对应边→能返回 t。利用它，我们可以精简遍历算法：

TrianglesOfVertex(v) {
    {t, i} = v.e
    do {
        {t, i} = t.nbr[i]      // 沿邻接边移动
        i = (i + 1) mod 3      // 绕顶点顺时针转到下一条边
    } while (t != v.e.t)
}

该版本无需搜索，每次迭代常数时间。存储开销：对仅含顶点位置的网格，每顶点存 4 个数（3 坐标 + 1 边），每面存 6 个数（3 个顶点索引 + 3 条边），总计约 4n_v + 6n_t ≈ 16n_v 单位存储，而基本索引网格只需 9n_v。多出的存储换来了常数时间的邻接遍历。

该结构仅适用于无边界流形网格（因为遍历依赖回到起始三角形来终止）。推广到带边界网格不难——为边界三角形的邻居引入哨兵值（如 -1），并确保边界顶点指向最逆时针的邻接三角形。

12.1.4 翼边结构

翼边数据结构（winged-edge data structure）将连通性信息存储在边而非面上，使边成为数据结构中的"一等公民"。每条边存储指向它的两个顶点（头顶点 head 和尾顶点 tail）、它所属的两个面（左面 left 和右面 right），以及最重要的——沿左面和右面逆时针方向遍历时的下一条边和前一条边：

Edge {
    Edge lprev, lnext, rprev, rnext;
    Vertex head, tail;
    Face left, right;
}
Face {
    // ... per-face data ...
    Edge e;   // 任意邻接边
}
Vertex {
    // ... per-vertex data ...
    Edge e;   // 任意入射边
}

翼边结构的完整字段语义

翼边结构（Baumgart 1974）之所以得名，是因为每条边像一只鸟的身体，左右两个面像展开的翅膀。下面是每条边 Edge 存储的 8 个引用（不含顶点和面）：

左右定向约定：以有向边 tail→head 为基准，左侧面（left）和右侧面（right）用"右手法则"确定：右手四指从 tail 弯向 head，掌心方向就是右面，手背方向是左面。这等价于：每个面的边界边按逆时针方向排列（从面外部看向面内部时）。

8 条指针的拓扑关系：对左边面而言，lnext 是从当前边出发，沿左面逆时针方向的下一条边。lprev 是前一条边。对称地，rnext 和 rprev 描述右面的边环。这 8 个引用构建了一张图——每个边都是图中一个节点，引用是出边——通过追踪引用可以在面、顶点和相邻边之间自由穿梭，且每步都是 O(1)。

翼边结构支持常数时间访问面或顶点的边，以及从边访问邻接顶点和面。遍历一个顶点的所有边：

EdgesOfVertex(v) {
    e = v.e
    do {
        if (e.tail == v)
            e = e.lprev    // 绕 v 逆时针移动
        else
            e = e.rprev
    } while (e != v.e)
}

遍历一个面的所有边：

EdgesOfFace(f) {
    e = f.e
    do {
        if (e.left == f)
            e = e.lnext    // 沿面顺时针移动
        else
            e = e.rnext
    } while (e != f.e)
}

这些算法同样适用于非三角形多边形网格，这是基于边的结构的一个重要优势。翼边结构支持各种时间/空间权衡——例如可以去掉 prev 引用以节省空间，但需沿后继边走完整圈来找前驱边，使某些操作变慢。

12.1.5 半边结构

翼边结构很优雅，但有个尴尬之处——每次移动前必须检查边的朝向。这类似于基础三角邻接结构中的搜索问题。解决方案也同样直接：不存储每条边的数据，而是为每条边存储两个半边（half-edge）——共享同一条物理边的两个三角形各有一个半边，朝向相反，各自与其三角形的方向一致：

HEdge {
    HEdge  pair, next;   // pair: 对面的半边
    Vertex v;             // 半边指向的顶点 (头顶点)
    Face   f;             // 半边所在的面
}
Face {
    // ... per-face data ...
    HEdge h;   // 该面的任意半边
}
Vertex {
    // ... per-vertex data ...
    HEdge h;   // 指向该顶点的任意半边
}

遍历半边结构就像遍历翼边结构，但无需检查朝向。我们通过 pair 指针访问对面半边：

EdgesOfVertex(v) {
    h = v.h
    do {
        h = h.pair.next    // 绕顶点遍历
    } while (h != v.h)
}
EdgesOfFace(f) {
    h = f.h
    do {
        h = h.next         // 沿面遍历
    } while (h != f.h)
}

半边的 Edge 也可以选择存储 prev 指针（指向面内的前一条半边）或仅存储 next。单指针版本更紧凑，双指针版本支持双向遍历。

半边结构的数据关系网

半边结构的核心设计哲学：每个实体类型（HEdge、HEVertex、HEFace）持有指向其他类型的最低限度引用，形成一个高度连通的有向图。下面以三角形网格为例，展示完整的互指关系：

HEdge 的四个引用：

struct HEdge {
    HEdge*   pair;    // 对侧半边：共享同一条物理边的另一半边
    HEdge*   next;    // 当前面内，逆时针方向的下一条半边
    Vertex*  v;       // 该半边指向的终点（头顶点）
    Face*    f;       // 该半边所属的面（永远非 NULL）
}

关键不变式：

不变式 1: h.pair.pair == h           —— 成对对称
不变式 2: h.next.f == h.f            —— 同面的半边
不变式 3: h.next.next.next == h      —— 三角形面内三边成环
不变式 4: h.v == h.pair.next.v       —— 对侧半边的起点是当前边的终点
不变式 5: h.next.v ≠ h.v            —— 相邻顶点互异

HEVertex 的角色：

struct HEVertex {
    float x, y, z;   // 位置
    // ... 法线、纹理坐标等 ...
    HEdge* h;        // 任意一条以该顶点为终点的半边
}

顶点不需要知道"属于哪些面"或"有多少条入射边"——这些信息可通过追踪 h 的引用间接获取。顶点存储的 h 是入射边集中任意一条（通常取第一条创建的），遍历所有入射边只需沿着 pair→next→pair→... 的路径走一圈。

HEFace 的角色：

struct HEFace {
    // ... 面的材质、颜色等属性 ...
    HEdge* h;        // 该面的边界半环中任意一条半边
}

面的边界是一个闭半环（half-loop）：从 f.h 出发，连续沿 next 走，经过恰好 n 步后回到 f.h（n 为面中边数）。

从任意半边出发可达的信息：

h.v          → 终点顶点的所有属性（坐标、法线、纹理坐标）
h.pair.v     → 起点顶点的所有属性
h.pair.f     → 对侧面的所有属性
h.next       → 同面内的下一条边
h.next.next  → 同面内的下下条边
h.pair.next.pair.next... → 绕起点顶点的边界环

遍历算法复杂度对比：

网格数据结构的实际应用场景

可见，每种结构有其最佳应用场景。现代图形引擎通常同时维护多种数据结构：运行时使用紧凑的索引网格进行渲染，编辑时动态转换为半边结构进行拓扑操作，编辑完成后转回索引网格。这种"多表示共存"策略平衡了性能与灵活性。

12.2 场景图

场景图（scene graph）以层次化树或有向无环图（DAG）形式组织场景内容。内部节点表示空间变换（旋转、平移、缩放），叶节点表示几何图元（三角网格）。

场景图节点类型枚举

现代场景图实现中，节点类型远比"变换"和"几何体"两类丰富。以下是生产级场景图中常见的核心节点类型：

层次变换累积的完整推导

场景图的核心数学机制是变换的层次累积。设场景图中从根到叶节点 v 的路径为 root → n₁ → n₂ → ... → n_k = v，每个节点 nᵢ 持有局部变换矩阵 Mᵢ（将点从该节点的局部坐标系变换到父节点坐标系）。

推导过程：

设 p_local 为节点 v 局部坐标系中的点。节点 n_k 的变换将其变换到 n_{k-1} 的坐标系：

p_{n_{k-1}} = M_k · p_local

节点 n_{k-1} 的变换将其变换到 n_{k-2} 的坐标系：

p_{n_{k-2}} = M_{k-1} · p_{n_{k-1}} = M_{k-1} · M_k · p_local

依此类推，直到世界坐标系：

p_world = M_1 · M_2 · ... · M_k · p_local

因此，节点 v 的组合世界变换为：

M_world(v) = M_root · M_1 · M_2 · ... · M_v

矩阵乘法顺序的重要性：M_root·M_1·M_2·... 是左乘顺序。数学上，最靠近 p_local 的矩阵先作用于点。在 OpenGL 和大多数图形库中，这意味着：

① 先缩放 → ② 再旋转 → ③ 最后平移（作用顺序从最近点开始向外）

但代码书写顺序（矩阵连乘）是相反的——从世界到局部。这正是第 6 章讨论的变换顺序问题在场景图中的大规模应用。

绕任意轴的层次旋转示例：假设场景图结构为 Root → RotateX(30°) → Translate(0, 5, 0) → RotateZ(45°) → Box。盒子的世界变换为：

M_world = Rx(30°) · T(0,5,0) · Rz(45°)

盒子上一点 p=(1,0,0) 的世界坐标为：先绕 Z 轴旋转 45°→ 平移到 y=5 → 绕 X 轴旋转 30°。注意 Box 绕自己的 Z 轴旋转，然后连同其坐标系一起被平移和进一步旋转——这正是场景图层次变换的核心语义。

traverse(node, M_parent):
    M_current = M_parent * node.localMatrix
    if node is Shape:
        render(node.geometry, M_current)
    for each child in node.children:
        traverse(child, M_current)

施加于一个节点的变换沿树向下累积（accumulate）。若节点 i 的局部变换为 4×4 齐次变换矩阵 M_i，则该节点（及其子树）中图元的世界空间变换为：

M_world(v) = M_root × M_1 × M_2 × ... × M_v

其中 M_root → M_1 → ... → M_v 是从根到节点 v 的路径上的所有局部变换矩阵的乘积。父节点的旋转会同时旋转其所有子节点——场景图的核心威力。

实例化（instancing）：通过有向无环图（DAG）结构，同一几何体可以在场景中多次出现，每次通过不同的变换路径引用。在实现中，一个子树可被多个父节点指向，但本身只存储一份。例如，4 个相同的轮子共享同一轮子几何体，通过 4 条不同的场景图路径（左前、右前、左后、右后）放置。这极大减少了内存消耗和绘制调用。

DAG 共享子图模型详解

将场景图从树（每个节点只有一个父节点）推广为有向无环图（DAG，节点可以有多个父节点），是实现实例化的关键一步。关键在于：子节点本身不被复制，只增加指向它的引用边。

实例化 vs 复制的内存对比：

场景：一辆汽车，含 4 个相同轮子（每轮 10K 三角形）+ 1 个底盘（50K 三角形）
            无实例化（纯树）    DAG 实例化
轮子几何体:  4 × 10K = 40K 三角  1 × 10K = 10K 三角  ← 节省 75%
内存占用:    90K 三角形全部        60K 三角形           ← 节省 33%
绘制调用:    5 个 draw call       5 个 draw call        ← 相同

绘制调用不变的原因：实例化节省的是几何数据内存（顶点缓冲、索引缓冲），但每个实例仍需一次绘制调用（不同的变换矩阵传入 GPU）。现代图形 API（DX12, Vulkan, Metal）支持 GPU 实例化（hardware instancing）——单个绘制调用可传递多个变换矩阵，进一步将 5 个调用合并为 2 个（底盘 1 + 轮子×4 合并 1）。

DAG 遍历中的世界变换计算：对于 DAG 中的节点 N，可能存在多条从根到 N 的路径，每条路径产生不同的世界变换。关键：

M_world(N, path_i) = M_root × (沿路径 path_i 的所有中间节点的局部变换连乘)

这意味着：① 同一节点在场景中以不同位置、朝向、大小出现——这正是实例化的语义；② 节点不存储自身的世界变换——世界变换是遍历中动态计算的，避免存储所有路径的组合；③ 遍历算法需要传递当前的累积变换，而非从节点本身读取。

实例化 vs 引用计数：DAG 节点需要引用计数来管理生命周期——当所有父节点移除对该子节点的引用时，子节点被释放。这与 C++ 的 shared_ptr 语义一致。

场景图遍历的完整伪代码

将遍历算法扩展为支持 DAG（检测重复访问）和属性继承（材质、光照、可见性）：

function traverseWithState(node, parentState, frameID):
    if node.lastVisited == frameID:
        return  // 避免 DAG 节点的重复遍历
    node.lastVisited = frameID
    
    // 合并状态
    currentState = parentState.clone()
    if node is Transform:
        currentState.transform *= node.localMatrix
    if node is Material:
        currentState.material = node.material
    if node.isCulled(currentState):
        return  // 裁剪：包围盒完全在视锥体外
    
    if node is Shape:
        renderWithState(node.geometry, currentState)
    else:
        for each child in node.children:
            traverseWithState(child, currentState, frameID)

frameID 技巧避免 DAG 中同一节点被多次渲染——每帧给一个唯一 ID，不必在遍历后重置标记。

在现代游戏引擎（Unity 的 GameObject 层级，Unreal 的 Actor 附件系统）中，场景图与组件化实体紧密集成——每个节点承载可配置的组件集合（渲染器、碰撞器、脚本等）。场景图还存储辅助信息：包围体（用于视锥体裁剪）、材质指派、光照和动画参数。

12.3 空间数据结构

当场景包含数千或数百万物体时，对每条光线测试每个物体是 O(N·P) 复杂度——对任何交互式应用都太慢了。空间数据结构（spatial data structure）将物体分组到嵌套的空间区域中，允许快速拒绝远离光线路径的组，将每条光线的有效复杂度降至 O(log N)。本节详细讨论三种方法：包围体层次结构（BVH）、均匀空间细分、以及二叉空间划分（BSP 树）。

12.3.1 包围盒与光线求交

最基本的加速原语是轴对齐包围盒（axis-aligned bounding box，AABB）。光线-AABB 求交使用平板测试（slab test）——计算光线进入和离开每对轴对齐平面的一维区间，然后找这些区间的交集。

设光线为 p(t) = e + td，包围盒为 [x_min, x_max] × [y_min, y_max] × [z_min, z_max]。对 x 轴平板：

x_min = e_x + t·d_x  →  t_xmin = (x_min − e_x) / d_x
x_max = e_x + t·d_x  →  t_xmax = (x_max − e_x) / d_x

类似地计算 y 和 z 的区间 [t_ymin, t_ymax] 和 [t_zmin, t_zmax]。光线命中盒子当且仅当三个区间有交集。从二维开始理解：

t_xmin = (x_min − x_e) / x_d
t_xmax = (x_max − x_e) / x_d
t_ymin = (y_min − y_e) / y_d
t_ymax = (y_max − y_e) / y_d

if (t_xmin > t_ymax) or (t_ymin > t_xmax) then
    return false     // 区间无重叠 → 光线错过盒子
else
    return true

关键洞察：两个一维区间不重叠意味着一个区间完全在另一个的左边或右边。推广到三维，需要三个轴区间全部重叠。

处理负方向分量：

当 x_d 为负时，光线会先碰到 x_max 而不是 x_min。因此代码扩展为：

if (x_d >= 0) then
    t_xmin = (x_min − x_e) / x_d
    t_xmax = (x_max − x_e) / x_d
else
    t_xmin = (x_max − x_e) / x_d
    t_xmax = (x_min − x_e) / x_d

IEEE 浮点除零处理：

使用 IEEE 标准的除零规则可以优雅处理水平和垂直光线（d_x = 0 或 d_y = 0）。根据 IEEE 754：

+a / 0 = +∞    (a > 0)
−a / 0 = −∞    (a > 0)

考虑垂直光线（x_d = 0, y_d > 0）的三种情况：

情况 1: x_e ≤ x_min → t_xmin = t_xmax = +∞ → 区间 (∞,∞) → 无命中 ✓
情况 2: x_min < x_e < x_max → t_xmin = −∞, t_xmax = +∞ → 区间 (−∞,∞) → 命中 ✓
情况 3: x_max ≤ x_e → t_xmin = t_xmax = −∞ → 区间 (−∞,−∞) → 无命中 ✓

全部按预期工作，无需特殊检查。但 x_d = −0 会导致问题——可以通过测试 x_d 的倒数来克服：

a = 1 / x_d
if (a >= 0) then
    t_min = a × (x_min − x_e)
    t_max = a × (x_max − x_e)
else
    t_min = a × (x_max − x_e)
    t_max = a × (x_min − x_e)

这既避免了 −0 问题，也将除法替换为更快的乘法（在光线遍历的数百次调用中节省可观时间）。

平板测试的完整三维推导

将二维平板测试推广到三维。设光线为 p(t) = e + td，AABB 为 [x_min, x_max] × [y_min, y_max] × [z_min, z_max]。

第一步：计算各轴的进入/离开 t 值。对每个轴 i ∈ {x, y, z}：

t_i_min = (box_i_min − e_i) / d_i
t_i_max = (box_i_max − e_i) / d_i

若 d_i < 0，交换 t_i_min ↔ t_i_max（保证 t_i_min ≤ t_i_max）

第二步：合并为进入/离开区间的交集。光线必须在所有三个轴上同时处于盒子内部。这意味着进入时间取各轴进入时间的最大值，离开时间取各轴离开时间的最小值：

t_near = max(t_x_min, t_y_min, t_z_min)    ← 光线最后进入盒子的时刻
t_far  = min(t_x_max, t_y_max, t_z_max)    ← 光线最早离开盒子的时刻

第三步：判断命中。光线命中盒子当且仅当：

t_near ≤ t_far     AND     t_far ≥ 0

条件 t_far ≥ 0 确保盒子在光线的前方（而不是背后）。如果只关心光线区间 [t_min, t_max] 内的命中，还需：

t_near ≤ t_max     AND     t_far ≥ t_min

为什么取 max 和 min？几何解释——设想光线是一根针，盒子是一个方形的窗户。针要从窗户穿过去，必须同时满足：

1. 在 x 方向上处于窗户左右边框之间（t_x_min ≤ t ≤ t_x_max）
2. 在 y 方向上处于窗户上下边框之间（t_y_min ≤ t ≤ t_y_max）
3. 在 z 方向上处于窗户前后边框之间（t_z_min ≤ t ≤ t_z_max）

针进入窗户的瞬间——是三个方向中最后进入的那一个时刻（max），因为只有等到那时，针才在所有三个方向上都进入了窗户。针离开窗户的瞬间——是三个方向中最早离开的那一个时刻（min），因为一旦在任一方向上离开，针就不再在窗户中了。因此：

进入 = max(各轴的进入时刻)，离开 = min(各轴的离开时刻)

完整的命中判断伪代码：

function rayAABB(ray e+td, float t0, float t1, Box b):
    // 1. 计算各轴平板区间
    for each axis i in {x, y, z}:
        invD = 1.0 / d[i]
        t0_i = (b.min[i] - e[i]) * invD
        t1_i = (b.max[i] - e[i]) * invD
        if invD < 0: swap(t0_i, t1_i)
        t_min[i] = t0_i
        t_max[i] = t1_i
    
    // 2. 合并区间
    t_near = max(t_min[x], t_min[y], t_min[z])
    t_far  = min(t_max[x], t_max[y], t_max[z])
    
    // 3. 判断
    if t_near > t_far:    return NO_HIT
    if t_far < 0:         return NO_HIT  // 盒子在光线后方
    if t_near > t1:       return NO_HIT  // 超出有效区间
    if t_far < t0:        return NO_HIT
    
    // 4. 最近交点在 t_near（若 t_near ≥ 0）或 t_far（若光线起点在盒内）
    t_hit = (t_near >= 0) ? t_near : t_far
    return HIT(t_hit)

与原始 2D 区间重叠测试的等价性：上述 t_near ≤ t_far 条件等价于原始文中对 2D 的 "两个一维区间必须重叠" 的直接推广。原始测试检查 "不存在任一维度上区间不重叠"，而合并测试通过 max/min 隐式地同时检查所有维度——若 t_near > t_far，意味着至少有一个维度的进入时刻晚于另一个维度的离开时刻，即至少一个维度上区间不重叠。

12.3.2 包围体层次结构 (BVH)

包围体层次结构（Bounding Volume Hierarchy，BVH）将物体递归分组为嵌套包围盒的树。内部节点包含包围其所有子节点的包围盒，但不直接包含几何体——只有叶节点包含实际三角形。

BVH 节点定义：

class bvh-node subclass of surface
    virtual bool hit(ray e+td, real t0, real t1, hit-record rec)
    virtual box bounding-box()
    surface-pointer left
    surface-pointer right
    box bbox

BVH 光线遍历：

function bool bvh-node::hit(ray a+tb, real t0, real t1, hit-record rec)
    if (bbox.hitbox(a+tb, t0, t1)) then
        hit-record lrec, rrec
        left-hit  = (left != NULL) and (left->hit(a+tb, t0, t1, lrec))
        right-hit = (right != NULL) and (right->hit(a+tb, t0, t1, rrec))
        if (left-hit and right-hit) then
            if (lrec.t < rrec.t) then rec = lrec
            else rec = rrec
            return true
        else if (left-hit) then rec = lrec; return true
        else if (right-hit) then rec = rrec; return true
        else return false
    else
        return false

SAH 启发式（Surface Area Heuristic）是构造高质量 BVH 的核心。它估计一条随机光线命中一个节点的概率正比于该节点包围盒的表面积。对于将节点分割为左子 L 和右子 R 的候选划分，SAH 代价函数为：

C = C_trav + (SA(L)/SA(P)) × C(L) + (SA(R)/SA(P)) × C(R)

其中 C_trav 是遍历一个节点的固定代价，SA(P) 是父节点 P 的表面积，SA(L) 和 SA(R) 是左右子节点的表面积，C(L) 和 C(R) 是递归计算的两个子树的代价（叶节点代价取决于包含的三角形数量）。选择使 C 最小的划分。

SAH 代价函数的逐项深度解释

SAH 公式看似简单，但每项背后都有深刻的几何概率含义。下面逐项拆解：

C = C_trav + (SA(L)/SA(P)) × C(L) + (SA(R)/SA(P)) × C(R)

第 1 项：C_trav——遍历开销

这是访问并测试一个内部节点的包围盒所需的固定计算量。在软件光线追踪中，C_trav ≈ 一次 AABB-光线求交测试的时钟周期数（约 8-20 条指令，取决于是否为分支友好实现）。在硬件光线追踪（如 NVIDIA RT Core）中，C_trav 实际上接近 0——因为包围盒测试被固化在硅片中，遍历开销几乎完全被硬件流水线掩盖。这意味着在硬件 RT 中，SAH 优化方向会倾向于更深的树（更精确的划分），因为遍历本身几乎无代价。

第 2 项：(SA(L)/SA(P)) × C(L)——左子树的期望代价

这里 SA(L)/SA(P) 是一个条件概率：假设一条随机光线从均匀分布的方向射入父节点包围盒 P 并命中 P，则它也命中左子节点包围盒 L 的条件概率正比于 L 的表面积与 P 的表面积之比。推导：

P(命中 L | 命中 P) ≈ SA(L) / SA(P)

为什么是表面积而不是体积？直觉上：一个细长方体的体积可能很小，但从外部看，任何射入方向都会先碰到长方体表面。表面积捕捉的是从外部观察时的"靶心"大小。数学上，对于凸体，随机均匀方向的射线命中该凸体的概率正比于其表面积——这与柯西曲面积分定理（Cauchy surface area formula）一致。

C(L) 是 L 子树的递归代价。如果 L 是叶节点（包含 n_L 个三角形），则 C(L) ≈ n_L × C_tri（光线与一个三角形求交的代价）。如果 L 是内部节点，则递归应用 SAH 公式。

第 3 项：(SA(R)/SA(P)) × C(R)——右子树的期望代价

与第 2 项完全对称。

SAH 公式的完整递归形式：

C(node) = 
  如果 node 是叶节点:  |node.triangles| × C_tri
  如果 node 是内部节点:
    C_trav + min_{所有可能划分 P→(L,R)} [
      SA(L)/SA(P) × C(L) + SA(R)/SA(P) × C(R)
    ]

注意：在最优划分中，需要尝试多种划分策略（沿 x/y/z 轴的多个位置），对每种划分计算 SAH 代价，取最小者。

SAH 的关键洞察——"空空间便宜"：如果 L 是空节点（不含任何三角形），C(L) = 0（叶节点但三角形数为 0）。但 SA(L) 可能很小（空空间被紧凑包裹）或很大（大面积空空间被包含在包围盒中）。SAH 会自然惩罚大面积空空间——因为大的 SA(L) 虽然该项为 0，但划分产生的大包围盒会传给子节点导致后续代价升高。一个好的 SAH 划分确保空空间不被包含在大的包围盒中。

BVH 构造策略的深入比较

BVH 构造——基于排序的划分：

function bvh-node::create(object-array A, int AXIS)
    N = A.length
    if (N = 1) then
        left = A[0];  right = NULL
        bbox = bounding-box(A[0])
    else if (N = 2) then
        left = A[0];  right = A[1]
        bbox = combine(bounding-box(A[0]), bounding-box(A[1]))
    else
        sort A by the object center along AXIS
        left  = new bvh-node(A[0..N/2−1], (AXIS+1) mod 3)
        right = new bvh-node(A[N/2..N−1], (AXIS+1) mod 3)
        bbox = combine(left->bbox, right->bbox)

基于空间的划分（按空间而非物体数量分区）可能更好：

function bvh-node::create(object-array A, int AXIS)
    ...
    else
        find the midpoint m of the bounding box of A along AXIS
        partition A into lists with lengths k and (N−k) surrounding m
        left  = new bvh-node(A[0..k], (AXIS+1) mod 3)
        right = new bvh-node(A[k+1..N−1], (AXIS+1) mod 3)
        bbox = combine(left->bbox, right->bbox)

这会产生不平衡树，但允许轻松遍历空空间且构造更便宜——划分比排序快。BVH 已成为现代光线追踪器的主导加速结构，包括 NVIDIA RT 核心的硬件实现。

12.3.3 均匀空间细分

另一种减少求交测试的思路：划分空间而非物体。这本质上不同于 BVH：在 BVH 中，每个物体属于恰好一个子节点，而空间中的点可能属于两个子节点；在空间细分中，空间中每个点属于恰好一个节点，而物体可能属于多个节点。

在均匀空间细分（uniform spatial subdivision）中，场景被划分为轴对齐的均匀盒子（体素）。光线按 3D DDA（数字微分分析器，与画线算法同理）步进方式遍历这些盒子。

3D DDA 步进算法：

设网格边界为 (x_min, y_min, z_min) 到 (x_max, y_max, z_max)，有 n_x × n_y × n_z 个单元。单元大小为：

Δx = (x_max − x_min) / n_x
Δy = (y_max − y_min) / n_y
Δz = (z_max − z_min) / n_z

光线遍历步进——平移到每个轴平行平面的参数 t：

t_next_x(i) = (x_min + i·Δx − e_x) / d_x    (若 i 为当前单元索引)
t_next_y(j) = (y_min + j·Δy − e_y) / d_y
t_next_z(k) = (z_min + k·Δz − e_z) / d_z

在每一步，选择 t_next_x、t_next_y、t_next_z 中的最小值，沿该轴步进到下一个单元，更新对应轴的 t_next 值（加上 Δx/d_x 或 Δy/d_y 或 Δz/d_z）：

// 初始化
找到光线起点的单元索引 (i, j, k)
计算初始的 t_next_x, t_next_y, t_next_z
计算步长 Δt_x = Δx / d_x, Δt_y = Δy / d_y, Δt_z = Δz / d_z

// 步进循环
while (光线在网格内)
    if (t_next_x < t_next_y and t_next_x < t_next_z)
        i += sign(d_x)
        t_next_x += Δt_x
    else if (t_next_y < t_next_z)
        j += sign(d_y)
        t_next_y += Δt_y
    else
        k += sign(d_z)
        t_next_z += Δt_z
    
    测试当前单元 (i,j,k) 内的物体
    if (命中) return 命中

3D DDA 增量步进的数学推导

DDA 算法的核心是增量公式，其推导基于一个简单事实：网格单元是均匀的，因此相邻的轴对齐平面之间的距离是固定的。

步长增量（step delta）的推导：当前单元索引为 (i, j, k)，光线沿 +x 方向行进。下一个 x 平面的位置在 x_min + (i+1)·Δx（光线出当前单元），再下一个在 x_min + (i+2)·Δx。到达这两个平面对应的 t 值分别为：

t_at_next_x_plane    = (x_min + (i+1)·Δx − e_x) / d_x
t_at_next_next_x_plane = (x_min + (i+2)·Δx − e_x) / d_x

两者之差即步长增量：

Δt_x = t_at_next_next_x_plane − t_at_next_x_plane
     = [(x_min + (i+2)·Δx − e_x) − (x_min + (i+1)·Δx − e_x)] / d_x
     = Δx / d_x

同理可得 Δt_y = Δy / d_y 和 Δt_z = Δz / d_z。这三个值在整个光线遍历中为常数——只需在初始化时计算一次。

初始化的完整步骤：

// 步骤 1: 确定光线起点的单元索引
i = floor((e_x - x_min) / Δx)
j = floor((e_y - y_min) / Δy)
k = floor((e_z - z_min) / Δz)

// 步骤 2: 计算初始的 3 个 t_next 值
// 若光线沿 +x 方向，下一个 x 平面是 x_min + (i+1)*Δx
// 若光线沿 -x 方向，下一个 x 平面是 x_min + i*Δx
if (d_x > 0):
    t_next_x = (x_min + (i+1)*Δx - e_x) / d_x
    step_x = 1
else if (d_x < 0):
    t_next_x = (x_min + i*Δx - e_x) / d_x
    step_x = -1
else:
    t_next_x = +∞   // 光线平行于 x 轴，永不会碰到 x 平面
    step_x = 0

// d_y 和 d_z 同理

// 步骤 3: 计算恒定的步长增量
Δt_x = Δx / abs(d_x)   // 使用绝对值，始终为正
Δt_y = Δy / abs(d_y)
Δt_z = Δz / abs(d_z)

选择步进轴的条件判断分析：在步进循环中，每次选择当前 t_next 最小的轴进行步进。这个条件保证光线按距离顺序穿越网格单元——先碰到的平面先穿越。三个 t_next 的比较使用浮点大小关系，在精度上有微妙之处：

若 t_next_x ≈ t_next_y（光线恰好穿过网格顶角）→
   优先顺序可由实现任意指定（通常选择 x/y/z 的固定优先级），
   因为两个 t 几乎相等，穿越两个平面的顺序不影响光线路过的格子集合。

复杂度分析：对于 n_x × n_y × n_z 的网格，一条光线最多穿过约 n_x + n_y + n_z 个单元（光线与所有轴平面对相交的总次数）。每次穿越一步，每步都做 2 次浮点比较、1 次整数加减和 1 次浮点加法。与每次步进都重新计算 t_next_x/y/z（每步 3 次除法）相比，增量法在光线穿越数百个单元时可节省上千次除法——这是一项显著的性能提升。

八叉树（octree）通过按需细分解决均匀网格的"空体素浪费"问题：在高密度区域递归细分为八个子体（二维中为四叉树）。kd 树（kd-tree）是一种轴对齐的 BSP 树，递归地以轴对齐平面划分空间，每层交替分割轴——与 BSP 树不同，kd 树的划分平面始终与坐标轴对齐。两者都提供对数遍历时间并自然处理非均匀几何体。

12.4 用于可见性的 BSP 树

二叉空间划分树（Binary Space Partition tree）递归地通过平面将空间一分为二。BSP 树最初为可见性排序开发——在没有 z 缓冲的时代，从任意视点以从前到后顺序渲染场景。即使有了 z 缓冲，BSP 树在需要精确排序的操作中仍很有价值。

12.4.1 BSP 树算法概述

BSP 树算法是画家算法（painter's algorithm）的一个实例：从后到前绘制每个物体，新多边形可能覆盖先前的多边形。基本思路：

sort objects back to front relative to viewpoint
for each object do
    draw object on screen

排序步骤的问题：多个物体的相对顺序并非总是良定义的——可能出现循环遮挡（如图中三个三角形 A 挡住 B，B 挡住 C，C 挡住 A）。BSP 树通过将场景预先组织为数据结构来解决此问题。

基本思想可用两个三角形 T₁ 和 T₂ 说明。回忆 T₁ 所在平面的隐式平面方程 f₁(p) = 0：半空间的一侧有 f₁(p⁺) > 0，另一侧 f₁(p⁻) < 0。因此对任意视点 e，我们可以确定正确的绘制顺序：

if (f₁(e) < 0) then
    draw T₁
    draw T₂
else
    draw T₂
    draw T₁

这个简单逻辑是 BSP 树遍历算法的基石。

12.4.2 BSP 树构建与遍历

BSP 树构建是递归过程：从包含场景所有多边形的根节点开始，选一个多边形所在平面作为分割平面。将每个剩余多边形分类为：正面（全部在正半空间）、背面、横跨（与平面相交，需分割）或共面（位于平面内，与分割多边形同节点）。递归地对正面和背面的多边形列表构建子树。

BSP 树从前到后遍历：

function traverse(node, eye):
    if (node is NULL) return
    if (f_node(eye) > 0)    // 视点在正面
        traverse(node.negative, eye)   // 先画更远的背面
        draw node.polygon              // 再画当前多边形
        traverse(node.positive, eye)   // 最后画更近的正面
    else                     // 视点在背面
        traverse(node.positive, eye)   // 先画更远的正面
        draw node.polygon
        traverse(node.negative, eye)   // 最后画更近的背面

该算法的美妙之处：无需显式排序——正确的从前到后（或从后到前）顺序通过递归遍历自动产生。

画家算法的从前到后遍历完整解析

画家算法（Painter's Algorithm）的名称源于传统油画技法：画家先画背景（远景），再画前景（近景），后画的覆盖先画的。BSP 树将这一直觉系统化为算法。下面给出完整的前后遍历逻辑：

// 从前到后渲染（配合 z-buffer，从前到后效率更高）
function renderFrontToBack(node, eye):
    if node is NULL: return
    side = classifyPoint(eye, node.plane)  // > 0: 正面, < 0: 背面, = 0: 在平面上
    
    if side >= 0:        // 视点在正面或在平面上
        renderFrontToBack(node.back, eye)   // 先画远处的背面
        drawPolygons(node.coplanar)         // 再画平面上的
        renderFrontToBack(node.front, eye)  // 最后画近处的正面
    else:                // 视点在背面
        renderFrontToBack(node.front, eye)  // 先画远处的正面
        drawPolygons(node.coplanar)
        renderFrontToBack(node.back, eye)   // 最后画近处的背面

// 从后到前渲染（用于纯画家算法，无 z-buffer）
function renderBackToFront(node, eye):
    if node is NULL: return
    side = classifyPoint(eye, node.plane)
    
    if side >= 0:        // 视点在正面，最近的物体在正面
        renderBackToFront(node.front, eye)  // 先画最近的正面
        drawPolygons(node.coplanar)
        renderBackToFront(node.back, eye)   // 后画远处的背面
    else:
        renderBackToFront(node.back, eye)
        drawPolygons(node.coplanar)
        renderBackToFront(node.front, eye)

为什么从前到后绘制在 z-buffer 场景中更优？在从前到后顺序下，远处的物体会先被绘制并写入 z-buffer。当近处物体稍后绘制时，GPU 的 Early-Z（早期深度测试）机制可以在片元着色器运行之前就丢弃被遮挡的片元——减少不必要的着色计算（即"overdraw 节省"）。反之，从后到前绘制时，近处物体先画，远处物体的每个片元都必须经过完整的着色管线后才能被 z-buffer 丢弃，浪费了大量着色计算。

处理循环遮挡：当三个三角形 A 挡 B、B 挡 C、C 挡 A 形成循环时，任何单一的排序都不可行。BSP 树通过切割三角形（在 12.4.4 节讨论）打破循环——切割一个三角形使其出现在多个子节点中，从而消除循环依赖。代价是三角形数量增加。

分割平面选择策略

BSP 树的构建质量高度依赖分割平面选择。常用策略包括：

1. 随机选择：从多边形列表中随机选一个作为分割平面。优点：快速；缺点：质量不可控。
2. 最少切割策略：尝试每个候选平面，选择切割其他多边形最少的那个。优点：控制三角形增长；缺点：构建时间 O(N²)，对大数据集过慢。
3. 平衡树策略：选择将多边形数量大致均分的平面，即使切割较多。优点：树深度最小（约 log N），遍历效率高；缺点：三角形数量可能膨胀。
4. 轴对齐策略：仅使用轴对齐平面（退化为 kd 树）。优点：平面评估快、遍历硬件友好；缺点：对斜向几何体划分效率低。

在实践中，大多数实现采用"最少切割"与"平衡树"的折中——先用启发式方法选出少数候选平面，再从中选切割最少的。这可在 O(N log N) 内构建质量可接受的 BSP 树。

12.4.3 BSP 光线追踪

BSP 树也可用于加速光线求交。遍历策略：在每个内部节点，检查光线起点在分割平面的哪一侧，先遍历起始侧子树，再遍历另一侧——前提是光线在另一侧的子区间非空（即存在 t 值落入另一侧）：

function bsp-node::hit(ray a+tb, real t0, real t1, hit-record rec)
    if (xb < 0) then   // 光线方向沿负 x
        // 对称处理
    ...
    xp = xa + t0 * xb
    if (xp < D) then   // 光线起点在负侧 (case 1)
        if (xb < 0) then return ... // 光线离开负侧 → 只查 left
        t = (D − xa) / xb            // 光线穿过平面的参数
        if (t > t1) then return left->hit(a+tb, t0, t1, rec)
        if (left->hit(a+tb, t0, t, rec)) then return true
        return right->hit(a+tb, t, t1, rec)
    else              // 光线起点在正侧 (case 2)
        // 对称代码

处理单个 BSP 节点比处理 BVH 节点更快，但同一表面可能存在于多个子树中（因为多边形被分割），导致更多节点和更高内存使用。构建质量决定哪种更快。

BSP 光线追踪的 t 参数区间处理

BSP 树光线追踪的核心在于管理光线穿过的参数区间 [t_near, t_far]。每进入一个子节点时，光线在其中的有效参数段被限制在该子节点对应的空间区域内。

扩展的 BSP 光线追踪伪代码（完整版，含 t 区间处理）：

function bsp-node::hit(ray a+tb, real t_near, real t_far, hit-record rec)
    // 步骤 1: 计算光线起点与分割平面的关系
    // 分割平面: x = D (简化为轴对齐情况，一般情况推广至任意平面)
    xp_start = a.x + t_near * b.x    // 光线在 t_near 处的 x 坐标
    xp_end   = a.x + t_far * b.x     // 光线在 t_far 处的 x 坐标
    
    // 步骤 2: 确定光线穿越分割平面的参数 t_split
    if b.x != 0:
        t_split = (D - a.x) / b.x
    else:
        t_split = +∞  // 光线平行于分割平面
    
    // 步骤 3: 根据起点位置 + 方向分类处理
    
    if xp_start >= D and xp_end >= D:
        // 情况 A: 光线完全在正侧 [D, +∞)
        // → 只遍历正侧子树
        return node.front->hit(a+tb, t_near, t_far, rec)
    
    else if xp_start <= D and xp_end <= D:
        // 情况 B: 光线完全在负侧 (−∞, D]
        // → 只遍历负侧子树
        return node.back->hit(a+tb, t_near, t_far, rec)
    
    else if xp_start < D:
        // 情况 C: 光线从负侧出发，穿过平面到正侧
        //        t_near ──── t_split ──── t_far
        //            (负侧)       |       (正侧)
        if node.back->hit(a+tb, t_near, t_split, rec):
            return true   // 负侧已找到交点
        return node.front->hit(a+tb, t_split, t_far, rec)
    
    else: // xp_start > D
        // 情况 D: 光线从正侧出发，穿过平面到负侧
        //        t_near ──── t_split ──── t_far
        //            (正侧)       |       (负侧)
        if node.front->hit(a+tb, t_near, t_split, rec):
            return true
        return node.back->hit(a+tb, t_split, t_far, rec)

区间的分裂与传递：关键洞察——每次穿过分割平面时，t 的有效区间被切分为两段：

负侧: [t_near, t_split]
正侧: [t_split, t_far]

两段分别传递给对应的子树。第一个被遍历的子树如果找到交点，直接返回；只有在第一段无交点时，才继续遍历第二段。这保证了先遇到的交点优先返回——符合光线追踪的"最近交点"需求。

t 区间处理的几个微妙之处：

1. t_split 不在 [t_near, t_far] 内——意味着光线没有穿过分割平面。例如，t_split > t_far（平面太远）或 t_split < t_near（平面在光线起点后方）。这些情况下，光线完全在一侧，直接只遍历对应子树。
2. b.x = 0（光线平行于分割平面）——光线永远不会穿过平面。此时只需判断 a.x 与 D 的大小关系，决定遍历哪一侧子树。
3. 浮点精度问题——t_split 的计算可能存在微小误差。如果一个交点恰好在分割平面上（t ≈ t_split），可能导致该交点被重复报告。解决方法是引入 ε 阈值：if |a.x + t_split * b.x - D| < ε 时将交点归属到优先遍历侧。

BSP 光线追踪与 BVH 光线追踪的效率比较：

12.4.4 切割三角形

当分割平面与三角形相交时，三角形必须被分割，否则一个三角形可能属于两侧，破坏可见性排序。使用 ε 容差处理退化情况：

若 f(p) >  ε → 点在正面
若 f(p) < −ε → 点在背面
若 |f(p)| ≤ ε → 点在分割平面上

三角形分割实操——利用顶点交换简化代码。确保恰好一个顶点在一侧，其他两个在另一侧：

if (f_a * f_c ≥ 0) then     // a 和 c 同侧
    swap(f_b, f_c);  swap(b, c)
    swap(f_a, f_b);  swap(a, b)
else if (f_b * f_c ≥ 0) then   // b 和 c 同侧
    swap(f_a, f_c);  swap(a, c)
    swap(f_a, f_b);  swap(a, b)

// 现在: 顶点 a 在一侧，b 和 c 在另一侧
compute A  // a-b 边与平面的交点
compute B  // a-c 边与平面的交点
T₁ = (a, b, A)
T₂ = (b, B, A)
T₃ = (A, B, c)

if (f_c ≥ 0) then
    negative-subtree.add(T₁, T₂)
    positive-subtree.add(T₃)
else
    positive-subtree.add(T₁, T₂)
    negative-subtree.add(T₃)

利用 f_a*f_c ≥ 0 判断同侧（乘积非负即同号），交换保持顶点逆时针顺序（需成对交换）。退化处理：恰好一个顶点在平面上的情况被上述代码涵盖——f 值为 0 时乘积为 0（≥ 0），该顶点被归入一侧。

12.5 多维数组的平铺

现代 CPU 和 GPU 拥有多级缓存层次结构，访问内存中靠近的数据比远处数据快数个数量级。当遍历大型多维数组时——如图像或体积纹理——天真的逐行扫描导致严重缓存未命中，因为逻辑上相邻的 y 坐标在内存中相距很远。平铺（tiling）将数组重排为小方块，使空间上相邻的元素在内存中也靠在一起——最大化缓存行利用。

12.5.1 莫顿码与 Z 阶曲线

平铺的关键是找到保持局部性的多维→一维地址映射。莫顿码（Morton code），又称Z 阶曲线（Z-order curve），通过对坐标的二进制位进行交错（interleaving）来实现这一映射。

位交错算法详解——以 x=1011₂, y=0110₂ 为例：

x = 1 0 1 1   (bits: x₃ x₂ x₁ x₀)
y = 0 1 1 0   (bits: y₃ y₂ y₁ y₀)

交错后: y₃ x₃ y₂ x₂ y₁ x₁ y₀ x₀
       = 0  1  1  0  1  1  0  1
       = 01101101₂

操作流程：取 x 和 y 的每一位，交替排列：y 的最高位 → x 的最高位 → y 的次高位 → x 的次高位 → ... → y 的最低位 → x 的最低位。结果是一个长度翻倍的整数，在莫顿空间中相邻的值对应于空间中相邻的方形块。

莫顿码位交错的逐步分解

以更详细的步骤展示 x=1011（十进制 11），y=0110（十进制 6）的位交错过程：

步骤 1：写出两个坐标的二进制表示（共 4 位，含前导零）：

位位置:       3    2    1    0
         x =  1    0    1    1   (十进制 11)
         y =  0    1    1    0   (十进制 6)

步骤 2：按位序交错（y 高位优先）：

交错位 7: y₃ = 0  ─┐
交错位 6: x₃ = 1  ─┤ 第 3 对（最高位对）
交错位 5: y₂ = 1  ─┤
交错位 4: x₂ = 0  ─┤ 第 2 对
交错位 3: y₁ = 1  ─┤
交错位 2: x₁ = 1  ─┤ 第 1 对
交错位 1: y₀ = 0  ─┤
交错位 0: x₀ = 1  ─┘ 第 0 对（最低位对）

结果: 0 1 1 0 1 1 0 1₂ = 109₁₀

步骤 3：验证——反交错还原坐标：

取偶数位 (0,2,4,6): x₀=1, x₁=1, x₂=0, x₃=1 → x = 1011₂ = 11 ✓
取奇数位 (1,3,5,7): y₀=0, y₁=1, y₂=1, y₃=0 → y = 0110₂ = 6  ✓

步骤 4：显式展示整个 4×4 空间区域的莫顿码（0-15）：

         x=0      x=1      x=2      x=3
y=0      0        1        4        5
y=1      2        3        6        7
y=2      8        9        12       13
y=3      10       11       14       15

沿莫顿码顺序遍历 (0→1→2→...→15):
(0,0)→(1,0)→(0,1)→(1,1)→(2,0)→(3,0)→(2,1)→(3,1)→
(0,2)→(1,2)→(0,3)→(1,3)→(2,2)→(3,2)→(2,3)→(3,3)

注意遍历顺序：先在 2×2 块内走 Z 形（(0,0)→(1,0)→(0,1)→(1,1)），再到下一个 2×2 块走 Z 形，再跳到右上 2×2 块，最后右下 2×2 块。整体形状是 4 个 Z 构成的大 Z 形——递归自相似。

莫顿码的硬件实现

位交错虽然在概念上简单，但在软件中直接用位操作实现效率不高（取决于 CPU 对位操作的支持）。现代 CPU 提供了专用指令加速：

// x86 BMI2 指令集
uint64_t morton2D(uint32_t x, uint32_t y) {
    return _pdep_u32(x, 0x55555555) | _pdep_u32(y, 0xAAAAAAAA);
    // PDEP (Parallel Bits Deposit): 将源位按掩码分散到目标位
}

// 反交错（提取坐标）
uint32_t unmorton2D_x(uint64_t code) {
    return _pext_u64(code, 0x5555555555555555);
    // PEXT (Parallel Bits Extract): 按掩码提取位并压缩
}

没有 PDEP/PEXT 指令的平台（ARM、老 x86），可以使用代价 O(log 位数) 的经典分治法：

// 分治法位交错（无专用指令的兼容实现）
uint64_t interleave(uint32_t x) {
    // 32 位 → 64 位（零位扩散）
    uint64_t w = x;
    w = (w | (w << 16)) & 0x0000FFFF0000FFFF;  // 16 位间距
    w = (w | (w <<  8)) & 0x00FF00FF00FF00FF;  // 8 位间距
    w = (w | (w <<  4)) & 0x0F0F0F0F0F0F0F0F;  // 4 位间距
    w = (w | (w <<  2)) & 0x3333333333333333;   // 2 位间距
    w = (w | (w <<  1)) & 0x5555555555555555;   // 1 位间距
    return w;
}
uint64_t morton2D_sw(uint32_t x, uint32_t y) {
    return interleave(x) | (interleave(y) << 1);
}

三维莫顿码（对 x, y, z 进行 3 路交错）和更高维度的映射完全遵循相同的位交错原理。

Z 阶曲线保证局部性的原理：莫顿码的高位决定了大的空间区域（高位的 y 和 x 位对应粗粒度的空间划分），低位决定了区域内的精细位置。共享相同高位前缀的两个莫顿码必然位于同一个大的空间区域中。形式上，两个莫顿码的差异越大（高位不同），它们对应的空间位置就越远——这恰好是缓存预取所需要的性质。

Z 阶曲线局部性的严格分析

为理解莫顿码为何能保持空间局部性，考虑二维平面上两个点 P₁=(x₁,y₁) 和 P₂=(x₂,y₂) 的莫顿码。设它们的距离为 δ = ||P₁ − P₂||。莫顿码保持局部性的核心论断：

|M(P₁) − M(P₂)| ≤ f(δ)  — 莫顿码差有一个上限，由空间距离决定

具体来说，两个点在 2^k × 2^k 的同一块内，当且仅当它们的高 k 个比特对（xₘ₋₁yₘ₋₁, ..., xₘ₋ₖyₘ₋ₖ）相同。若两点距离小于 2^k，它们在莫顿码中最多相差约 2^(2k) —— 因为低位 2k 个比特可能不同。这一上界保证了空间上靠近的点在莫顿空间中也靠近。

但莫顿码并非完美——"Z 阶跳跃"问题：

莫顿码有一个著名的问题：两个空间上极度靠近的点，如果它们恰好分别位于大的空间块边界两侧，它们的莫顿码差异可能非常大。例如：

点 A = (3, 3) → 莫顿码 = 0b1111 = 15
点 B = (4, 4) → 莫顿码 = 0b110000 = 48

空间距离 = √2 ≈ 1.41，莫顿码距离 = 33！

原因是 x=3=011₂ 和 x=4=100₂ 的所有位都不同（进位级联），导致莫顿码一半以上的位都翻转。这种跳跃现象发生在位模式整数进位的位置（即 2^k 的边界处），是莫顿码作为空间填充曲线（space-filling curve）的一个固有缺陷。

Hilbert 曲线作为替代方案：

Hilbert 曲线是莫顿码的升级版——它将空间同样映射为一维索引，但保证连续的索引在空间位置上总是相邻的（无跳跃）。代价是编码/解码的计算量高出一倍。在实际工程中：

莫顿码：编码/解码极快（PDEP/PEXT 单指令级），有跳跃 → 适用于实时应用
Hilbert：编码/解码需更多操作，无跳跃 → 适用于存储索引和离线处理

大多数图形应用（纹理平铺、GVDB 体积渲染、线性 BVH 构建）选择莫顿码，因为实时编码速度优先于完美的局部性。

莫顿码将空间组织成递归嵌套的四叉树（2D）或八叉树（3D）模式，是现代 GPU 纹理平铺格式（如 DirectX 的"莫顿交错"）和空间数据库索引的基础。它们也是线性八叉树的关键组件——节点使用其莫顿码索引，实现无指针的高效邻居查找和遍历。

12.5.2 单级平铺（二维）

对于 N_x × N_y 的二维数组，用 n × n 的瓦片进行单级平铺。坐标 (x,y) 到一维索引的公式：

index = n² × (B_x × b_y + b_x) + y' × n + x'

其中：

b_x = x ÷ n          (瓦片列索引)
b_y = y ÷ n          (瓦片行索引)
x'  = x mod n        (瓦片内 x 偏移)
y'  = y mod n        (瓦片内 y 偏移)
B_x = N_x ÷ n        (每行瓦片数)

展开为完整形式：

index = n² × ((N_x ÷ n) × (y ÷ n) + (x ÷ n)) + (y mod n) × n + (x mod n)

这个表达式可分离为一维函数的和：

index = F_x(x) + F_y(y)

其中：

F_x(x) = n² × (x ÷ n) + (x mod n)
F_y(y) = n² × (N_x ÷ n) × (y ÷ n) + (y mod n) × n

由于 F_x 和 F_y 只依赖各自的维度坐标，可以预先计算查找表（N_x 和 N_y 大小），将多维索引降为 O(1) 查表——避免多次除法和取模操作。对于高维大数据集，这些查找表远小于数据本身，可以放入 CPU 的 L1/L2 缓存。

12.5.3 两级平铺（三维）

对于三维数据（体积纹理、医学扫描、3D 模拟网格），两级平铺（two-level tiling）提供最佳缓存性能：首先将空间细分为 m × m × m 的宏观块，每个宏观块由 n × n × n 的小块组成。例如，40×20×19 的体积可分解为 4×2×2 个宏观块，每个宏观块含 2×2×2 个砖块，每个砖块含 5×5×5 个单元（m = 2, n = 5）。

两级平铺索引公式：

index = F_x(x) + F_y(y) + F_z(z)

其中各维度的贡献（省略推导过程）：

F_x(x) = ((x÷n)÷m) × n³ × m³ × ((Nz÷n)÷m) × ((Ny÷n)÷m)
       + ((x÷n) mod m) × n³ × m²
       + (x mod n) × n²

F_y(y) = ((y÷n)÷m) × n³ × m³ × ((Nz÷n)÷m)
       + ((y÷n) mod m) × n³ × m
       + (y mod n) × n

F_z(z) = ((z÷n)÷m) × n³ × m³
       + ((z÷n) mod m) × n³
       + (z mod n)

与二维平铺类似，这些分量可预计算为查找表，实现极速索引转换。经验上，对 16 位数据集取 m=5，对 float 数据集取 m=6 能得到较好效果。

两级平铺公式的完整推导

以上公式的推导涉及三维空间中 4 层嵌套索引的平面化。将三维坐标 (x, y, z) 的运行过程分解如下：

层次结构（从粗到细）：

1. 宏观块层 (macroblock): 尺寸 m × m × m 个砖块。整个体积可容纳约 (N_x/(m·n)) × (N_y/(m·n)) × (N_z/(m·n)) 个宏观块。
2. 砖块层 (brick): 每个砖块含 n × n × n 个单元。
3. 单元层 (cell): 单个元素，在砖块内偏移为 (x', y', z')。

坐标分解：

宏观块坐标:  Bx3 = (x ÷ n) ÷ m    ← x 方向第几个宏观块
            By3 = (y ÷ n) ÷ m
            Bz3 = (z ÷ n) ÷ m

宏观块内砖块偏移: bx2 = (x ÷ n) mod m   ← 宏观块内砖块 x 偏移
                by2 = (y ÷ n) mod m
                bz2 = (z ÷ n) mod m

砖块内单元偏移:    x' = x mod n
                y' = y mod n
                z' = z mod n

全局索引的构建（从粗到细累加偏移）：

一层地址 = 所有更大层的偏移量 × 各层的跨度。跨度的计算方法：

1 个砖块 = n³ 个单元
1 层宏观块在 x 方向 = m 个砖块 → m × n³ 个单元
1 层宏观块在 y 方向 = m 层宏观块 → m² × n³ 个单元

宏观块行在 z 方向跨度: stride_z3 = m³ × n³
宏观块内砖块 z 跨度:   stride_z2 = m² × n³
砖块内 z 跨度:         stride_z1 = n²

1 个宏观块在 z 方向 = m × 砖块 z 跨度 = m × n³ 个单元？

构建过程中，需要按 (z 宏观块, z 砖块, z 单元) → (y 宏观块, y 砖块, y 单元) → (x 宏观块, x 砖块, x 单元) 的三重嵌套累加。最终得到与文本一致的分离形式 F_x + F_y + F_z。这里展示 F_x 的完整分解：

F_x(x) = Bx3 × (stride_per_x_macroblock)   // 宏观块层面
       + bx2 × (stride_per_x_brick)        // 砖块层面
       + x'  × (stride_per_x_cell)         // 单元层面

其中：
  stride_per_x_macroblock = m³ × n³ × ((N_z ÷ n) ÷ m) × ((N_y ÷ n) ÷ m)
                          = m³ × n³ × ⌈N_z/(m·n)⌉ × ⌈N_y/(m·n)⌉
  stride_per_x_brick = m² × n³
  stride_per_x_cell  = n²

F_y 和 F_z 同理，但乘数递减（F_y 少一个 ((Ny÷n)÷m) 因子，F_z 更少）。这正是原公式的结构。

直观理解——数字示例：

假设 N_x=40, N_y=20, N_z=19, m=2, n=5。则：

宏观块网格: (40/(2×5), 20/(2×5), 19/(2×5)) = (4, 2, 2)  → 共 16 个宏观块
每宏观块砖块: 2×2×2 = 8 个
每砖块单元:   5×5×5 = 125 个

总跨度: 
  1 个砖块 = 125 单元
  1 层宏观块 (y方向) = 2 砖块 = 250 单元  
  1 个宏观块 = 8 砖块 = 1000 单元
  1 列宏观块 (y方向) = 2 宏观块 = 2000 单元
  1 层宏观块 (z方向) = 4×2 = 8 宏观块 = 8000 单元
  1 片宏观块 = 2×4×2 = 8 宏观块？

对于点 (x=12, y=8, z=3)：

x=12 → Bx3=1, bx2=0, x'=2 (第 2 个宏观块，砖块偏移 0，单元偏移 2)
y=8  → By3=0, by2=1, y'=3
z=3  → Bz3=0, bz2=0, z'=3

F_x = 1 × (8×125×2×2) + 0 × (4×125) + 2 × 25 = 1×4000 + 0 + 50 = 4050
F_y = 0 × (8×125×2) + 1 × (8×125) + 3 × 5 = 0 + 1000 + 15 = 1015
F_z = 0 × (8×125) + 0 × 125 + 3 = 0 + 0 + 3 = 3

index = 4050 + 1015 + 3 = 5068

验证：该点在第 2 个 x 宏观块中偏移 50（即在该宏观块内第 50 个单元），加上前面跨越的 8 个宏观块各 500 单元（因为 1 个宏观块 = 1000 单元，但前面只数 1 个 x 宏观块的 1/2？），得 4000 + ... 确实匹配。

行主序（无平铺）:  120 ns   (75% TLB miss + 20% L2 miss)
单级平铺 (n=8):    45 ns    (15% TLB miss + 10% L2 miss)
两级平铺 (m=4,n=8): 28 ns    (5% TLB miss + 5% L2 miss)
莫顿码平铺:        22 ns    (3% TLB miss + 8% L2 miss)

课后练习题（含答案）

World (根)
├── Frame (车架) — M_frame
│   ├── FrontWheelGroup — 绕前叉旋转
│   │   ├── FrontFork
│   │   └── FrontWheel — 再绕轮轴旋转
│   ├── RearWheelGroup — 绕后叉旋转
│   │   └── RearWheel
│   ├── HandlebarGroup — 绕头管旋转
│   │   └── Handlebar
│   ├── SeatPost — 可上下平移
│   │   └── Seat
│   ├── PedalGroup — 绕中轴旋转
│   │   ├── CrankArm
│   │   ├── Pedal_L
│   │   └── Pedal_R
│   └── Chain — 受踏板旋转驱动

QA 零基础问答区