《Fundamentals of Computer Graphics》第五版 第十二章 图形学数据结构

259 阅读11分钟

有些数据结构,由于回答了诸如曲面、空间、场景结构等基本问题,而在图形学中反复出现。本章主要介绍一些最常见、最有用的基础数据结构:网格结构(mesh structure)、空间数据结构、场景图(scene graph)、分片多维数组(tiled multidimensional array)。

对于网格,除了用于存储和传输静态网格的基本存储方案之外,还介绍了翼边结构(winged-edge structure)和半边结构(half-edge structure),它们常用于管理动态模型。

由于在管理物体及其变换方面非常有用,场景图数据结构在图形学应用中很普遍。新设计的图形 API 都可以很好地支持场景图。

空间数据结构主要介绍了包围体层级结构(bounding volume hierarchy)、层次空间细分(hierarchical space subdivision)、均匀空间细分(uniform space subdivision),以及 BSP 树隐面消除算法。这些方法也适用于几何剔除(geometry culling)、碰撞检测(collision detection)。

最初为了提升分页性能(paging performance)而发展分片多维数组的,是那些需要和磁盘交换图形数据的应用。现在这种结构对内存局部性(memory locality)至关重要。

三角网格(Triangle Meshes)

真实世界中的大多数模型都是由共用顶点的三角形构成,通常称为三角网格(triangular meshes,或 triangle meshes)、不规则三角网络(triangular irregular networks, 简称 TINs),因此高效地处理它们对图形程序的性能至关重要。哪种效率最重要取决于具体应用场景,常见的有:存储空间、传输数据(网络传输、CPU 向图形系统传输)所占用的带宽等;而对于需要操作网格——比如细分(subdivision)、网格编辑(mesh editing)、网格压缩(mesh compression)等——的应用,高效地访问邻接信息至关重要。

三角网格并不是一系列无关的三角形,而是互相共用顶点和边的三角形所构成的网络,因而处理起来可以更加高效。

一个三角网格所需的最基本信息是一组三角形——顶点三元组(triples of vertices),及其顶点坐标。此外,许多程序还需要在顶点、边、面上存储额外数据。最常用的是顶点数据,然后借助顶点数据对三角形线性插值,有时还需要在每条边、每个面上存储数据。

网格拓扑(Mesh Topology)

网格的类表面(surface-like)性可以表述为对网格拓扑(mesh topology)——三角形互相连接的方式——的限制。最简单、最严格的要求是:网格必须为流形(manifold)。粗略地说,二维流形就是一个曲面,上面任一点的邻域都可以用一个平面来近似。

manifold.png

验证一个网格是否为流形等价于检查下面两个条件是否满足:

  • 每条边仅被两个三角形共用
  • 每个顶点只被一个三角形闭环(complete loop of triangles)围绕

上述对流形的限制过于严格,有时也会放宽限制,允许网格是带边流形(manifold with boundary),此时检查条件为:

  • 每条边仅被一个或两个三角形共用
  • 每个顶点只与一个边连通三角形集(edge-connected set of triangles)邻接

manifold-with-boundary.png

许多应用依赖于曲面方向(orientation of surface)。曲面方向根据顶点顺序来定义:按三角形顶点顺序右手螺旋,大拇指方向为正面(front)或外面(outside),反之为背面(back)或里面(inside);另一种等价表述是,从正面观测三角形,顶点沿逆时针方向。如果一个连通网格的所有三角形正面都保持一致,就称该网格为方向一致的(consistently oriented),这等价于每一对相邻三角形方向一致。

有些系统使用顺时针顶点方向定义三角形正面。

orientation.png

方向一致的网格也称为可定向的(orientable)。流形或带边流形可以是不可定向的,但不会影响实际应用。

mobius.png

索引网格存储(Indexed Mesh Storage)

每个三角形可以存储为:

Triangle {vector3 vertexPosition[3]}\begin{aligned} &\mathrm{Triangle}\ \{ \\ &\quad\mathrm{vector3\ vertexPosition[3]} \\ &\} \end{aligned}

separate-shared.png

每个三角形都存储自身顶点数据会导致存储空间的浪费。因此可以让三角形共用顶点数据以减少空间占用,也就是使用共用顶点网格(shared-vertex mesh)。

Triangle {Vertex v[3]}Vertex {vector3 position// or other vertex data}\begin{aligned} &\mathrm{Triangle}\ \{ \\ &\quad\mathrm{Vertex\ v[3]} \\ &\} \\ \\ &\mathrm{Vertex}\ \{ \\ &\quad\mathrm{vector3\ position}\quad \text{// or other vertex data} \\ &\} \end{aligned}

数组 v\mathrm{v} 中存储的是 Vertex\mathrm{Vertex} 对象的引用或指针。在实现中,顶点和三角形通常用数组存储,三角形对顶点的引用通过数组索引表示:

IndexedMesh {int tInd[nt][3]vector3 verts[nv]}\begin{aligned} &\mathrm{IndexedMesh}\ \{ \\ &\quad\begin{aligned} &\mathrm{int\ tInd}[n_{t}][3] \\ &\mathrm{vector3\ verts}[n_{v}] \end{aligned} \\ &\} \end{aligned}

这种存储共用顶点的网格表示称为索引三角网格(indexed triangle mesh)。

假定浮点数、指针、整数所占空间相同,那么上述两种存储方式所需空间分别为:

  • Triangle\mathrm{Triangle},每个三角形存储三个向量,总共占据 9nt9n_{t} 个存储单元。
  • IndexedMesh\mathrm{IndexedMesh},每个顶点存储一个向量,每个三角形存储三个整数,总共占据 3nv+3nt3n_{v}+3n_{t} 个存储单元。

根据经验,大网格中每个顶点平均连接约 6 个三角形,因此 nt2nvn_{t}\approx 2n_{v}。这使得索引三角网格节省了一半的空间,如果有其它顶点属性,节约的空间会更多。

三角带(Triangle Strips)与三角扇(Triangle Fans)

索引网格是最常见的三角网格表示,也常用于网络传输以及应用和图形流水线间数据传输。如果需要更紧凑的存储空间,可以使用三角带(triangle strips)或三角扇(triangle fans)。

triangle-fan.png

在三角扇中,有一个公共顶点被所有三角形共用,其它顶点生成一组类似风扇扇叶的三角形。顶点序列中的第一个顶点表示三角扇中心,后面每一对连续顶点对应一个三角形。

triangle-strip.png

三角带与三角扇类似,但它适用于范围较宽的网格。三角带就是在沿着一个线状条带上下交替地排列顶点。顶点序列中,每三个连续顶点构成一个三角形。为了保持方向一致,相邻两个三角形需要反转顶点顺序。

strips-indices.png

在三角带和三角扇中,n+2n+2 个顶点足以描述 nn 个三角形。相比于索引三角网格,顶点索引的空间占用降至 (n+2)/(3n)(n+2)/(3n),如果三角带足够长,大约为 1/31/3

relative-size.png

可以看到,(n+2)/(3n)(n+2)/(3n) 衰减得很快。因此,即使是无结构的网格,使用贪婪算法把它们聚合成较短的三角带也是有价值的。

网格连通性(Mesh Connectivity)

索引网格、三角带、三角扇都可以紧凑地表示静态网格,但不适合修改网格。网格编辑需要更复杂的数据结构来高效地解决连通性相关的问题,比如:

  • 找出指定三角形的三个邻接三角形
  • 找出共用指定边的两个三角形
  • 找出所有共用指定顶点的面
  • 找出所有共用指定顶点的边

连通性本质上是一种拓扑性质,而网格的连通性完全由邻接关系决定。

由于许多应用中网格规模都非常大,因此高效的表示至关重要。最直接的实现方式是存储三种类型 Vertex\mathrm{Vertex}Edge\mathrm{Edge}Triangle\mathrm{Triangle} 的数据以及所有关系:

Triangle {Vertex v[3]Edge e[3]}Edge {Vertex v[2]Triangle t[2]}Vertex {Triangle t[ ]Edge e[ ]}\begin{aligned} &\mathrm{Triangle}\ \{ \\ &\quad\mathrm{Vertex\ v}[3] \\ &\quad\mathrm{Edge\ e}[3] \\ &\} \\ \\ &\mathrm{Edge}\ \{ \\ &\quad\mathrm{Vertex\ v}[2] \\ &\quad\mathrm{Triangle\ t}[2] \\ &\} \\ \\ &\mathrm{Vertex}\ \{ \\ &\quad\mathrm{Triangle\ t}[\ ] \\ &\quad\mathrm{Edge\ e}[\ ] \\ &\} \end{aligned}

但由于这些信息互相关联,存储冗余过多。而且上述 Vertex\mathrm{Vertex} 中涉及变长数组,一般没有高效的实现方式。因此,最好的做法是定义一个类接口,基于更高效的数据结构来解决上面关于连通性的问题,而不需要显式存储所有关系。

虽然上述 Edge\mathrm{Edge}Triangle\mathrm{Triangle} 中存储的都是定长数组,但多边形网格(polygon mesh)中每个多边形可以有任意多的边和顶点。这使得许多传统数据结构都基于边。对于只有三角形的网格,也可以将连通性存储在面上。

一个好的网格数据结构应该相当紧凑,而且可以高效地进行邻接查询,即查询时间不依赖网格规模。

三角形近邻结构(The Triangle-Neighbor Structure)

在基本的共用顶点网格中,对每个三角形扩充三个到邻接三角形的指针,对每个顶点扩充一个到任一邻接三角形的指针,便可以得到一种基于三角形的紧凑的网格数据结构:

Triangle {Triangle nbr[3]Vertex v[3]}Vertex {// ... per-vertex data ...Triangle t // any adjacent tri}\begin{aligned} &\mathrm{Triangle}\ \{ \\ &\quad\mathrm{Triangle\ nbr}[3] \\ &\quad\mathrm{Vertex\ v}[3] \\ &\} \\ \\ &\mathrm{Vertex}\ \{ \\ &\quad\text{// ... per-vertex data ...} \\ &\quad\mathrm{Triangle}\ t\ \text{// any adjacent tri} \\ &\} \end{aligned}

Triangle.nbr[k]\mathrm{Triangle.nbr}[k] 指向共用顶点 kkk+1k+1 的邻接三角形。这种结构称为三角形近邻结构(triangle-neighbor structure)。

triangle-neighbor-structure.png

利用这种数据结构,除了三角形的顶点和邻接三角形之外,还可以找到顶点的邻接三角形。如果顶点 vv 在三角形 tt 中的索引为 kk,那么三角形 t.nbr[k]t.\mathrm{nbr}[k] 就是绕顶点 vv 顺时针方向的下一个三角形。下面是顶点的邻接三角形遍历算法:

TrianglesOfVertex(v) {t=v.tdo {find i such that (t.v[i]==v)t=t.nbr[i]} while (t != v.t)}\begin{aligned} &\mathrm{TrianglesOfVertex}(v)\ \{ \\ &\quad\begin{aligned} &t = v.t \\ &\mathrm{do}\ \{ \\ &\quad\begin{aligned} &\text{find}\ i\ \text{such that}\ (t.\mathrm{v}[i] == v) \\ &t = t.\mathrm{nbr}[i] \end{aligned} \\ &\}\ \mathrm{while}\ (t\ \text{!=}\ v.t) \end{aligned} \\ &\} \end{aligned}

上述算法需要搜索中心顶点在三角形中的索引,这引入了额外分支。为减少分支,可做如下改进,不再存储到邻接三角形的指针,而是存储到邻接三角形中相应邻边的指针:

Triangle {Edge nbr[3]Vertex v[3]}Edge { // the i-th edge of triangle tTriangle tint i // in {0, 1, 2}}Vertex {// ... per-vertex data ...Edge e // any edge leaving vertex, whose// triangle is counterclockwise}\begin{aligned} &\mathrm{Triangle}\ \{ \\ &\quad\mathrm{Edge\ nbr}[3] \\ &\quad\mathrm{Vertex\ v}[3] \\ &\} \\ \\ &\mathrm{Edge}\ \{\ \text{// the $i$-th edge of triangle $t$} \\ &\quad\mathrm{Triangle}\ t \\ &\quad\mathrm{int}\ i\ \text{// in \{0, 1, 2\}} \\ &\} \\ \\ &\mathrm{Vertex}\ \{ \\ &\quad\text{// ... per-vertex data ...} \\ &\quad\mathrm{Edge}\ e\ \text{// any edge leaving vertex, whose} \\ &\qquad\qquad \text{// triangle is counterclockwise} \\ &\} \end{aligned}

Edge\mathrm{Edge} 表达的实际含义是三角形 tt 的第 ii 条边。

在实践中,可以从三角形索引 tt 中借出两个比特来存储边的索引 ii,以此来实现 Edge\mathrm{Edge},这让存储空间保持不变。借助 Edge\mathrm{Edge},可以快速地找到上一个或下一个邻接三角形,比如,对于任意三角形 tt 的索引为 jj 的边,下式恒成立:

t.nbr[j].t.nbr[t.nbr[j].i].t==tt.\mathrm{nbr}[j].t.\mathrm{nbr}[t.\mathrm{nbr}[j].i].t == t

邻接三角形遍历算法也可以改进为:

TrianglesOfVertex(v) {{t, i}=v.edo {{t, i}=t.nbr[i]i=(i+1) mod 3} while (t != v.e.t)}\begin{aligned} &\mathrm{TrianglesOfVertex}(v)\ \{ \\ &\quad\begin{aligned} &\{t,\ i\} = v.e \\ &\mathrm{do}\ \{ \\ &\quad\begin{aligned} &\{t,\ i\} = t.\mathrm{nbr}[i] \\ &i = (i + 1)\ \mathrm{mod}\ 3 \end{aligned} \\ &\}\ \mathrm{while}\ (t\ \text{!=}\ v.e.t) \end{aligned} \\ &\} \end{aligned}

三角形近邻结构相当紧凑,对于只有顶点坐标属性的网格,大约总共占据 4nv+6nt16nv4n_{v}+6n_{t}\approx 16n_{v} 个存储单元,而基础索引网格占用 9nv9n_{v} 个存储单元。

为了让上述算法适用于带边流形,可以为边界三角形的邻居引入合适的哨兵值(sentinel value)——比如 1-1,并保证边界顶点指向逆时针方向最后一个邻接三角形。

翼边结构(The Winged-Edge Structure)

翼边结构(winged-edge structure)是一种广泛使用的将连通性存储在边(edge)上的网格数据结构。在这种结构中,边是一等公民,每条边上需要存储:

  • (head)、(tail)两个顶点的指针
  • (left)、(right)两个面的指针
  • 在左、右面上沿逆时针方向遍历的前一条边和后一条边

此外,每个顶点、每个面也要扩充一个到任一条关联的边的指针。

Edge {Edge lprev, lnext, rprev, rnextVertex head, tailFace left, right}Face {// ... per-face data ...Edge e // any adjacent edge}Vertex {// ... per-vertex data ...Edge e // any incident edge}\begin{aligned} &\mathrm{Edge}\ \{ \\ &\quad\begin{aligned} &\mathrm{Edge\ lprev,\ lnext,\ rprev,\ rnext} \\ &\mathrm{Vertex\ head,\ tail} \\ &\mathrm{Face\ left,\ right} \end{aligned} \\ &\} \\ \\ &\mathrm{Face}\ \{ \\ &\quad\text{// ... per-face data ...} \\ &\quad\mathrm{Edge}\ e\ \text{// any adjacent edge} \\ &\} \\ \\ &\mathrm{Vertex}\ \{ \\ &\quad\text{// ... per-vertex data ...} \\ &\quad\mathrm{Edge}\ e\ \text{// any incident edge} \\ &\} \end{aligned}

winged-edge-structure.png

上图中 rnext 和 rprev 标反了!

翼边数据结构可以快速找出一个面或一个顶点的所有的边,借助这些边可以进一步找到相邻的所有顶点和面:

EdgesOfVertex(v) {e=v.edo {if (e.tail==v)e=e.lprevelsee=e.rprev} while (e != v.e)}EdgesOfFace(f) {e=f.edo {if (e.left==f)e=e.lnextelsee=e.rnext} while (e != f.e)}\begin{aligned} &\mathrm{EdgesOfVertex}(v)\ \{ \\ &\quad\begin{aligned} &e = v.e \\ &\mathrm{do}\ \{ \\ &\quad\begin{aligned} &\mathrm{if}\ (e.\mathrm{tail} == v) \\ &\quad e = e.\mathrm{lprev} \\ &\mathrm{else} \\ &\quad e = e.\mathrm{rprev} \end{aligned} \\ &\}\ \mathrm{while}\ (e\ \text{!=}\ v.e) \end{aligned} \\ &\} \\ \\ &\mathrm{EdgesOfFace}(f)\ \{ \\ &\quad\begin{aligned} &e = f.e \\ &\mathrm{do}\ \{ \\ &\quad\begin{aligned} &\mathrm{if}\ (e.\mathrm{left} == f) \\ &\quad e = e.\mathrm{lnext} \\ &\mathrm{else} \\ &\quad e = e.\mathrm{rnext} \end{aligned} \\ &\}\ \mathrm{while}\ (e\ \text{!=}\ f.e) \end{aligned} \\ &\} \end{aligned}

上述算法和数据结构同样适用于多边形网格(polygon mesh),而不局限于三角网格,这也是基于边的结构的重要优势之一。

和任何数据结构一样,翼边数据结构可以在时间和空间上进行各种权衡。比如,删掉 prev\mathrm{prev} 引用可以节省空间,但却让顺时针绕面遍历或者逆时针绕顶点遍历变得更加困难;想要知道前一条边,只能使用 next\mathrm{next} 引用绕一圈来得到。

半边结构(The Half-Edge Structure)

翼边结构中遍历所有的边时需要检查边的方向,和三角形近邻结构类似,用半边(half-edge)代替边,可以避免这种检查。原先存储在边上的数据分摊到了两条半边上,它们方向相反,且分别从属于两侧的面。每条半边包含了:

  • 到所属面的指针
  • 到头顶点的指针
  • 到所属面上邻边的指针
  • 到另一侧的半边的指针

和翼边结构类似,可以把前一条和后一条邻边指针都记录下来,也可以只记录其中一条。下面给出的是只记录后一条邻接半边的数据结构:

HEdge {HEdge pair, nextVertex vFace f}Face {// ... per-face data ...HEdge h // any h-edge of this face}Vertex {// ... per-vertex data ...HEdge h // any h-edge pointing toward this vertex}\begin{aligned} &\mathrm{HEdge}\ \{ \\ &\quad\begin{aligned} &\mathrm{HEdge\ pair,\ next} \\ &\mathrm{Vertex}\ v \\ &\mathrm{Face}\ f \end{aligned} \\ &\} \\ \\ &\mathrm{Face}\ \{ \\ &\quad\text{// ... per-face data ...} \\ &\quad\mathrm{HEdge}\ h\ \text{// any h-edge of this face} \\ &\} \\ \\ &\mathrm{Vertex}\ \{ \\ &\quad\text{// ... per-vertex data ...} \\ &\quad\mathrm{HEdge}\ h\ \text{// any h-edge pointing toward this vertex} \\ &\} \end{aligned}

half-edge-structure.png

边的遍历算法如下:

EdgesOfVertex(v) {h=v.hdo {h=h.next.pair} while (h != v.h)}EdgesOfFace(f) {h=f.hdo {h=h.next} while (h != f.h)}\begin{aligned} &\mathrm{EdgesOfVertex}(v)\ \{ \\ &\quad\begin{aligned} &h = v.h \\ &\mathrm{do}\ \{ \\ &\quad h = h.\mathrm{next}.\mathrm{pair} \\ &\}\ \mathrm{while}\ (h\ \text{!=}\ v.h) \end{aligned} \\ &\} \\ \\ &\mathrm{EdgesOfFace}(f)\ \{ \\ &\quad\begin{aligned} &h = f.h \\ &\mathrm{do}\ \{ \\ &\quad h = h.\mathrm{next} \\ &\}\ \mathrm{while}\ (h\ \text{!=}\ f.h) \end{aligned} \\ &\} \\ \end{aligned}

由于半边通常成对出现,因此很多实现可以移除 pair\mathrm{pair} 指针。

场景图(Scene Graphs)

一个三角网格可以表示一个物体,而图形应用中另一个常见问题是如何摆放物体,这可以通过变换来实现。大多数场景允许层次化组织,因此可以使用场景图(scene graph)根据这种层次结构来管理这些变换。场景图中每个节点都对应一个物体或一个坐标系,上面还记录了当前物体相对于父节点坐标系的变换矩阵。施加在物体上的变换矩阵就是场景图中沿着从物体到根节点的路径上所有矩阵的乘积。

scene-graph-ferry.png

栅格系统

场景图在栅格系统中可以借助矩阵栈(matrix stack)来高效地实现,许多 API 都支持这种数据结构。矩阵栈的 push\mathrm{push}pop\mathrm{pop} 操作表示在乘积的右手边添加、移除矩阵。借助矩阵栈可以对场景图做深度优先遍历:

function traverse(node)push(Mlocal)draw object using composite matrix from stacktraverse(left child)traverse(right child)pop()\begin{aligned} &\textbf{function}\ \text{traverse(node)} \\ &\quad\begin{aligned} &\text{push}(\mathbf{M}_{\text{local}}) \\ &\text{draw object using composite matrix from stack} \\ &\text{traverse(left child)} \\ &\text{traverse(right child)} \\ &\text{pop()} \end{aligned} \end{aligned}

场景图可能有很多变形,但基本思想是一致的。

射线追踪

射线追踪可以非常自然地应用矩阵变换,而无需更改几何表示。物体实例可以用原始几何和变换矩阵来表示,渲染时再处理显式构造。

射线追踪的优势是可以选择计算交点的空间。射线 a+tb\vec{a}+t\vec{b} 与变换后物体的交点等价于原始空间中逆变换射线(inverse-transformed ray)与未变换物体(untransformed object)的交点。在原始空间中计算交点可能有下面两个优势:

  1. 未变换物体计算交点的程序可能更加简单
  2. 许多变换后的物体都共用同一个未变换物体,因此节约存储空间

intersection-ray-tracing.png

对于 Surface\mathrm{Surface} 及其子类实例,结合法向量变换,可以给出 hit\mathrm{hit} 函数的伪代码,用于计算射线与变换后物体的交点:

instance::hit(Ray a+tb, real t0, real t1, HitRecord rec)Ray ray=M1a+tM1bif (baseObjecthit(ray, t0, t1, rec)) thenrec.n=(M1)Trec.nreturn trueelsereturn false\begin{aligned} &\mathrm{instance::hit}(\mathrm{Ray}\ \vec{a}+t\vec{b},\ \mathrm{real}\ t_{0},\ \mathrm{real}\ t_{1},\ \mathrm{HitRecord\ rec}) \\ &\quad\begin{aligned} &\mathrm{Ray\ ray} = \mathbf{M}^{-1}\vec{a} + t\mathbf{M}^{-1}\vec{b} \\ &\mathbf{if}\ (\mathrm{baseObject}\rightarrow\mathrm{hit}(\mathrm{ray},\ t_{0},\ t_{1},\ \mathrm{rec}))\ \mathbf{then} \\ &\quad\mathrm{rec}.\vec{n} = (\mathbf{M}^{-1})^{T}\mathrm{rec}.\vec{n} \\ &\quad\mathbf{return}\ \mathrm{true} \\ &\mathbf{else} \\ &\quad\mathbf{return}\ \mathrm{false} \end{aligned} \end{aligned}

上述伪代码要求射线方向 b\vec{b} 不能局限于单位向量。

空间数据结构(Spatial Data Structures)

很多图形应用都需要快速定位到某一空间区域中的几何物体,比如射线追踪器中计算交点、交互应用中寻找视口范围内的物体、游戏和物理仿真中的碰撞检测等。这些需求都可以通过各种空间数据结构(spatial data structures)来实现。

将物体层次化组织起来称为物体划分(object partitioning)策略,物体被分成不相交的组,但不同组可能会空间交叠。将空间分割为不相交区域称为空间划分(space partitioning)策略,空间被分成独立的区域,一个物体可能跨越多个分区。空间划分可以是规则的——空间被分成形状一致的部分,也可以是不规则的——空间被自适应地分成不规则的部分。

下面以射线追踪为例来讨论这些结构,它们同样适用于视区剔除(view culling)、碰撞检测(collision detection)等。和大多数搜索问题一样,计算射线与物体的交点也可以用分治(divide and conquer)思想在亚线性(sub-linear)时间复杂度内解决。

包围盒(Bounding Box)

对于大多数加速计算交点的方法,关键操作都是计算射线和包围盒的交点。下面以 2D 情形为例进行说明,3D 情形类似。2D 包围盒围成的区域可以表示为:

[xmin,xmax]×[ymin,ymax][x_{\text{min}}, x_{\text{max}}] \times [y_{\text{min}}, y_{\text{max}}]

假定射线与包围盒边界所在直线的交点参数分别为 txmint_{\text{xmin}}txmaxt_{\text{xmax}}tymint_{\text{ymin}}tymaxt_{\text{ymax}},区间 [txmin,txmax][t_{\text{xmin}},t_{\text{xmax}}] 表示射线上介于直线 x=xminx=x_{\text{min}}x=xmaxx=x_{\text{max}} 之间的点的参数,y 方向同理。易知,射线与包围盒相交等价于区间 [txmin,txmax][t_{\text{xmin}}, t_{\text{xmax}}][tymin,tymax][t_{\text{ymin}}, t_{\text{ymax}}] 有交集

ray-box-intersection.png

判断射线与包围盒是否相交的伪代码如下:

txmin=(xminxe)/xdtxmax=(xmaxxe)/xdtymin=(yminye)/ydtymax=(ymaxye)/ydif (txmin>tymax) or (tymin>txmax) thenreturn falseelsereturn true\begin{aligned} &t_{\text{xmin}} = (x_{\text{min}} - x_{e})/x_{d} \\ &t_{\text{xmax}} = (x_{\text{max}} - x_{e})/x_{d} \\ &t_{\text{ymin}} = (y_{\text{min}} - y_{e})/y_{d} \\ &t_{\text{ymax}} = (y_{\text{max}} - y_{e})/y_{d} \\ &\mathbf{if}\ (t_{\text{xmin}}>t_{\text{ymax}})\ \mathrm{or}\ (t_{\text{ymin}}>t_{\text{xmax}})\ \mathbf{then} \\ &\quad\mathbf{return}\ \mathrm{false} \\ &\mathbf{else} \\ &\quad\mathbf{return}\ \mathrm{true} \end{aligned}

上述伪代码只考虑了射线的方向向量分量均为正的情况,对于非正的情况,以 x 方向为例,借助 IEEE 浮点数规范,考虑到 0-0,做如下处理即可:

a=1/xdif (a0) thentxmin=a(xminxe)txmax=a(xmaxxe)elsetxmin=a(xmaxxe)txmax=a(xminxe)\begin{aligned} &a = 1/x_{d} \\ &\mathbf{if}\ (a\geqslant 0)\ \mathbf{then} \\ &\quad t_{\text{xmin}} = a(x_{\text{min}} - x_{e}) \\ &\quad t_{\text{xmax}} = a(x_{\text{max}} - x_{e}) \\ &\mathbf{else} \\ &\quad t_{\text{xmin}} = a(x_{\text{max}} - x_{e}) \\ &\quad t_{\text{xmax}} = a(x_{\text{min}} - x_{e}) \end{aligned}

这些伪代码其实都是在判断直线与包围盒是否相交,对于射线,应该在计算参数 txmint_{\text{xmin}}tymint_{\text{ymin}} 时,将负值截断到 00;计算 txmaxt_{\text{xmax}}tymaxt_{\text{ymax}} 时,根据 txmint_{\text{xmin}}tymint_{\text{ymin}} 决定是否提前 return\mathrm{return}

层次包围盒(Hierarchical Bounding Box)

假定用二叉树来描述包围体层级结构(bounding volume hierarchy)。其中,每个叶节点包含一个图元,每个非叶节点包含一颗或两颗子树,并带有一个包围盒。包围盒包含了层级在该节点下的所有物体,但不一定包含与它存在空间交叠的所有物体。由于非叶节点下的两颗子树没有几何顺序,因此射线可能同时与这两颗子树相交。

hierarchical-bounding-boxes.png

包围体层级树结构中的节点 bvh-node\text{bvh-node}Surface\text{Surface} 的子类:

class bvh-node subclass of Surfacevirtual bool hit(Ray e+td, real t0, real t1, HitRecord rec)virtual box boundingBox(object a)Surface-pointer leftSurface-pointer rightbox bbox\begin{aligned} &\textbf{class}\ \text{bvh-node subclass of Surface} \\ &\quad\begin{aligned} &\text{virtual bool hit}(\text{Ray}\ \vec{e}+t\vec{d},\ \text{real}\ t_{0},\ \text{real}\ t_{1},\ \text{HitRecord rec}) \\ &\text{virtual box boundingBox}(\text{object}\ a) \\ &\text{Surface-pointer left} \\ &\text{Surface-pointer right} \\ &\text{box bbox} \end{aligned} \end{aligned}

bvh-node\text{bvh-node}hit\text{hit} 方法伪代码如下:

function bool bvh-node::hit(Ray e+td, real t0, real t1, HitRecord rec)if (bbox.hitbox(e+td, t0, t1)) thenHitRecord lrec, rrecleftHit=(leftNULL) and (lefthit(e+td, t0, t1, lrec))rightHit=(rightNULL) and (righthit(e+td, t0, t1, rrec))if (leftHit and rightHit) thenif (lrec.t<rrec.t) thenrec = lrecelserec = rrecreturn trueelse if (leftHit) thenrec = lrecreturn trueelse if (rightHit) thenrec = rrecreturn trueelsereturn falseelsereturn false\begin{aligned} &\textbf{function}\ \text{bool bvh-node::hit}(\text{Ray}\ \vec{e}+t\vec{d},\ \text{real}\ t_{0},\ \text{real}\ t_{1},\ \text{HitRecord rec}) \\ &\quad\begin{aligned} &\textbf{if}\ (\text{bbox.hitbox}(\vec{e}+t\vec{d},\ t_{0},\ t_{1}))\ \textbf{then} \\ &\quad\begin{aligned} &\text{HitRecord lrec, rrec} \\ &\text{leftHit} = (\text{left}\neq\text{NULL})\ \text{and}\ (\text{left}\rightarrow\text{hit}(\vec{e}+t\vec{d},\ t_{0},\ t_{1},\ \text{lrec})) \\ &\text{rightHit} = (\text{right}\neq\text{NULL})\ \text{and}\ (\text{right}\rightarrow\text{hit}(\vec{e}+t\vec{d},\ t_{0},\ t_{1},\ \text{rrec})) \\ &\textbf{if}\ (\text{leftHit and rightHit})\ \textbf{then} \\ &\quad\begin{aligned} &\textbf{if}\ (\text{lrec}.t < \text{rrec}.t)\ \textbf{then} \\ &\quad \text{rec = lrec} \\ &\textbf{else} \\ &\quad \text{rec = rrec} \\ &\textbf{return}\ \text{true} \end{aligned} \\ &\textbf{else if}\ (\text{leftHit})\ \textbf{then} \\ &\quad \text{rec = lrec} \\ &\quad\textbf{return}\ \text{true} \\ &\textbf{else if}\ (\text{rightHit})\ \textbf{then} \\ &\quad \text{rec = rrec} \\ &\quad\textbf{return}\ \text{true} \\ &\textbf{else} \\ &\quad\textbf{return}\ \text{false} \end{aligned} \\ &\textbf{else} \\ &\quad\textbf{return}\ \text{false} \end{aligned} \end{aligned}

如果树构建得足够合适,可以省略掉 left\text{left}NULL\text{NULL} 检查。对于 right\text{right},用 left\text{left} 来代替 NULL\text{NULL},可以省略掉 right\text{right}NULL\text{NULL} 检查;但这导致 left\text{left} 被检查了两遍,这种做法是否值得取决于树构建时的细节。

对于包围盒层级树,除了二叉(binary)和大致平衡(roughly balanced)之外,一般还要保证兄弟节点的包围盒之间尽可能少地交叠。一种构建方法是,在划分一组物体之前先沿坐标轴对物体排序,伪代码如下:

function bvh-node::create(object-array A, int AXIS)N = A.lengthif (N==1) thenleft = A[0]right = NULLbbox = boundingBox(A[0])else if (N==2) thenleft = A[0]right = A[1]bbox = combine(boundingBox(A[0]), boundingBox(A[1]))elsesort A by the object center along AXISleft = new bvh-node(A[0..N/21], (AXIS + 1) mod 3)right = new bvh-node(A[N/2..N1], (AXIS + 1) mod 3)bbox = combine(leftbbox, rightbbox)\begin{aligned} &\textbf{function}\ \text{bvh-node::create(object-array A, int AXIS)} \\ &\quad\begin{aligned} &\text{N = A.length} \\ &\textbf{if}\ (\text{N} == 1)\ \textbf{then} \\ &\quad\begin{aligned} &\text{left = A[0]} \\ &\text{right = NULL} \\ &\text{bbox = boundingBox(A[0])} \end{aligned} \\ &\textbf{else if}\ (\text{N} == 2)\ \textbf{then} \\ &\quad\begin{aligned} &\text{left = A[0]} \\ &\text{right = A[1]} \\ &\text{bbox = combine(boundingBox(A[0]), boundingBox(A[1]))} \end{aligned} \\ &\textbf{else} \\ &\quad\begin{aligned} &\text{sort A by the object center along AXIS} \\ &\text{left = new bvh-node(A[0..N/2$-$1], (AXIS + 1) mod 3)} \\ &\text{right = new bvh-node(A[N/2..N$-$1], (AXIS + 1) mod 3)} \\ &\text{bbox = combine(left$\rightarrow$bbox, right$\rightarrow$bbox)} \end{aligned} \end{aligned} \end{aligned}

其中,x、y、z 轴分别用整数 001122 表示。为了提高树的构建质量,可以在每次递归时选择那个让两颗子树的包围盒体积之和最小的坐标轴。这对于物体排布不规则的场景有显著提升。只做划分不做全排也可以让上述代码效率更高。

另一种可能更好的构建方法是,让两颗子树包含的空间大小——而不是物体个数——尽可能相等,伪代码如下:

function bvh-node::create(object-array A, int AXIS)N = A.lengthif (N==1) thenleft = A[0]right = NULLbbox = boundingBox(A[0])else if (N==2) thenleft = A[0]right = A[1]bbox = combine(boundingBox(A[0]), boundingBox(A[1]))elsefind the midpoint m of the bounding box of A along AXISpartition A into lists with lengths k and (Nk) surrounding mleft = new bvh-node(A[0..k], (AXIS + 1) mod 3)right = new bvh-node(A[k+1..N1], (AXIS + 1) mod 3)bbox = combine(leftbbox, rightbbox)\begin{aligned} &\textbf{function}\ \text{bvh-node::create(object-array A, int AXIS)} \\ &\quad\begin{aligned} &\text{N = A.length} \\ &\textbf{if}\ (\text{N} == 1)\ \textbf{then} \\ &\quad\begin{aligned} &\text{left = A[0]} \\ &\text{right = NULL} \\ &\text{bbox = boundingBox(A[0])} \end{aligned} \\ &\textbf{else if}\ (\text{N} == 2)\ \textbf{then} \\ &\quad\begin{aligned} &\text{left = A[0]} \\ &\text{right = A[1]} \\ &\text{bbox = combine(boundingBox(A[0]), boundingBox(A[1]))} \end{aligned} \\ &\textbf{else} \\ &\quad\begin{aligned} &\text{find the midpoint $m$ of the bounding box of A along AXIS} \\ &\text{partition A into lists with lengths $k$ and (N$-k$) surrounding $m$} \\ &\text{left = new bvh-node(A[0..$k$], (AXIS + 1) mod 3)} \\ &\text{right = new bvh-node(A[$k+1$..N$-$1], (AXIS + 1) mod 3)} \\ &\text{bbox = combine(left$\rightarrow$bbox, right$\rightarrow$bbox)} \end{aligned} \end{aligned} \end{aligned}

虽然这种方法让二叉树不再平衡,但它对空白空间的遍历更加简单,而且构建速度较快。

均匀空间细分(Uniform Spatial Subdivision)

均匀空间细分(uniform spatial subdivision)是将整个场景分成齐轴盒子(axis-aligned box)——也称单元格(cell),它们大小相同,但不必是立方体。

网格本身是 Surface\text{Surface} 的子类,而且应该用三维 Surface\text{Surface} 指针数组来实现。空单元格的指针为 NULL\text{NULL};只包含一个物体的单元格,指针指向该物体;包含多个物体的单元格,指针指向一个列表、另一个网格、或是其它数据结构——如包围体层级结构。

uniform-spatial-subdivision.png

射线沿着网格中的单元格遍历,直到命中一个物体。考虑到网格由一系列等间距平行平面构成,射线可以增量方式遍历单元格。以二维情形为例,首先需要找到射线 e+td\vec{e}+t\vec{d} 命中的第一个单元格的索引 (i,j)(i,j),然后再根据下一组平面交点参数 txnextt_{\text{xnext}}tynextt_{\text{ynext}} 来决定移动 ii 还是移动 jj。检查单元格内射线与物体是否相交时,必须将参数 tt 的范围限制在单元格内部。

三维 Surface\text{Surface} 指针数组可以通过分片来改善局部性(locality)。

traversal-uniform-spatial-subdivision.png

齐轴空间二分(Axis-Aligned Binary Space Partitioning)

空间二分树(binary space partitioning tree,简称 BSP 树)是一种空间层级结构。树中每个节点都包含一个截平面(cutting plane)和左右两颗子树,每颗子树分别包含截平面一侧的所有物体,对于跨平面的物体,两颗子树均存储。BSP 树节点:

class bsp-node subclass of Surfacevirtual bool hit(Ray e+td, real t0, real t1, HitRecord rec)virtual box boundingBox(object a)Surface-pointer leftSurface-pointer rightreal D\begin{aligned} &\textbf{class}\ \text{bsp-node subclass of Surface} \\ &\quad\begin{aligned} &\text{virtual bool hit}(\text{Ray}\ \vec{e}+t\vec{d},\ \text{real}\ t_{0},\ \text{real}\ t_{1},\ \text{HitRecord rec}) \\ &\text{virtual box boundingBox}(\text{object}\ a) \\ &\text{Surface-pointer left} \\ &\text{Surface-pointer right} \\ &\text{real}\ D \end{aligned} \end{aligned}

假定截平面为 x=Dx=D,射线求交的起点为 p=e+t0d\vec{p}=\vec{e}+t_{0}\vec{d},需考虑以下四种情况:

  1. xp<Dx_{p}<Dxd<0x_{d}<0。只检查左半空间。
  2. xp<Dx_{p}<Dxd>0x_{d}>0。先检查左半空间,如果没有命中,再检查右半空间。
  3. xp>Dx_{p}>Dxd>0x_{d}>0。只检查右半空间。
  4. xp>Dx_{p}>Dxd<0x_{d}<0。先检查右半空间,如果没有命中,再检查左半空间。

BSP-cutting-plane.png

射线求交伪代码如下:

function bool bsp-node::hit(Ray e+td, real t0, real t1, HitRecord rec)xp=xe+t0xdif (xp<D) thenif (xd<0) thenreturn (leftNULL) and (lefthit(e+td,t0,t1,rec))t=(Dxe)/xdif (t>t1) thenreturn (leftNULL) and (lefthit(e+td,t0,t1,rec))if (leftNULL) and (lefthit(e+td,t0,t,rec)) thenreturn truereturn (rightNULL) and (righthit(e+td,t,t1,rec))elseanalogous code for cases 3 and 4\begin{aligned} &\textbf{function}\ \text{bool bsp-node::hit}(\text{Ray}\ \vec{e}+t\vec{d},\ \text{real}\ t_{0},\ \text{real}\ t_{1},\ \text{HitRecord rec}) \\ &\quad\begin{aligned} &x_{p} = x_{e} + t_{0}x_{d} \\ &\textbf{if}\ (x_{p} < D)\ \textbf{then} \\ &\quad\begin{aligned} &\textbf{if}\ (x_{d} < 0)\ \textbf{then} \\ &\quad\textbf{return}\ (\text{left}\neq\text{NULL})\ \text{and}\ (\text{left}\rightarrow\text{hit}(\vec{e}+t\vec{d}, t_{0}, t_{1}, \text{rec})) \\ &t = (D - x_{e})/x_{d} \\ &\textbf{if}\ (t > t_{1})\ \textbf{then} \\ &\quad\textbf{return}\ (\text{left}\neq\text{NULL})\ \text{and}\ (\text{left}\rightarrow\text{hit}(\vec{e}+t\vec{d}, t_{0}, t_{1}, \text{rec})) \\ &\textbf{if}\ (\text{left}\neq\text{NULL})\ \text{and}\ (\text{left}\rightarrow\text{hit}(\vec{e}+t\vec{d}, t_{0}, t, \text{rec}))\ \textbf{then} \\ &\quad\textbf{return}\ \text{true} \\ &\textbf{return}\ (\text{right}\neq\text{NULL})\ \text{and}\ (\text{right}\rightarrow\text{hit}(\vec{e}+t\vec{d}, t, t_{1}, \text{rec})) \end{aligned} \\ &\textbf{else} \\ &\quad\text{analogous code for cases 3 and 4} \end{aligned} \end{aligned}

调用上述方法前,需要先计算射线与根节点包围盒的交点以确定 t0t_{0}t1t_{1}。考虑到齐轴截平面的方向有三种,可以给 bsp-node\text{bsp-node} 添加一个整数索引 axis\text{axis} 表示截平面的垂直轴;假定点可以通过索引来访问分量,则上述伪代码中的 x 分量用索引代替即可,如:用 a[axis]a[\text{axis}] 替换 xax_{a}

虽然单个 bsp-node\text{bsp-node} 处理起来比 bvh-node\text{bvh-node} 更快,但是一个物体在 BSP 树中可能存储于不只一颗子树,这让 BSP 树占据更多空间。树的构建质量决定了最终的快慢,BSP 树的构建和 BVH 树类似,可以循环更换截平面的垂直轴并且每次让截平面平分整个区域,也可以更加细致地控制截平面。

BSP 树在隐面消除中的应用

空间数据结构还可用来确定任意视角下物体的遮挡关系。使用非齐轴空间二分树(non–axis-aligned binary space partitioning tree)可以给出一种优雅的隐面消除算法,称为 BSP 树算法(BSP tree algorithm)。BSP 树隐面消除算法的关键点在于,事先创建一个在任意视角下都有用的数据结构。

BSP 树算法是画家算法的一个例子。下面假定场景中所有模型均由三角形构成,但算法也适用于任意平面多边形(planar polygon)构成的模型。

以三角形 T1T_{1} 为根节点构建一颗二叉树,f1f_{1} 是三角形 T1T_{1} 所在平面的隐式方程。负分支(negative branch)包含了顶点均满足 f1(p)<0f_{1}(\vec{p})<0 的所有三角形,正分支(positive branch)包含了顶点均满足 f1(p)>0f_{1}(\vec{p})>0 的所有三角形。

BSP-draw.png

绘制整个场景的伪代码如下:

function draw(bsptree tree, point e)if (tree.empty) thenreturnif (ftree.root(e)<0) thendraw(tree.plus,e)rasterize tree.triangledraw(tree.minus,e)elsedraw(tree.minus,e)rasterize tree.triangledraw(tree.plus,e)\begin{aligned} &\textbf{function}\ \text{draw}(\text{bsptree tree},\ \text{point}\ \vec{e}) \\ &\quad\begin{aligned} &\textbf{if}\ (\text{tree.empty})\ \textbf{then} \\ &\quad\textbf{return} \\ &\textbf{if}\ (f_{\text{tree.root}}(\vec{e}) < 0)\ \textbf{then} \\ &\quad\begin{aligned} &\text{draw}(\text{tree.plus}, \vec{e}) \\ &\text{rasterize tree.triangle} \\ &\text{draw}(\text{tree.minus}, \vec{e}) \end{aligned} \\ &\textbf{else} \\ &\quad\begin{aligned} &\text{draw}(\text{tree.minus}, \vec{e}) \\ &\text{rasterize tree.triangle} \\ &\text{draw}(\text{tree.plus}, \vec{e}) \end{aligned} \end{aligned} \end{aligned}

将递归终止条件提高一个层级可以稍微提高一点代码效率。平面方程选择三点式还是一般式,以及是否需要存储法向量,取决于具体实现上对时间和空间的权衡。

树的构建

上述绘制算法假定了任何一个三角形不会横跨任何其它三角形所在平面。当三角形横跨分割面时,可以将其分割成三个三角形:

T1=(a,b,A)T2=(b,B,A)T3=(A,B,c)\begin{aligned} &T_{1} = (\vec{a}, \vec{b}, \vec{A}) \\ &T_{2} = (\vec{b}, \vec{B}, \vec{A}) \\ &T_{3} = (\vec{A}, \vec{B}, \vec{c}) \end{aligned}

顶点顺序保证了分割前后三角形的法向一致。

cross-dividing-plane.png

考虑到顶点可能离分割面非常近,而且这种情况对于共用顶点的三角形模型很常见,因此需要对数值精度做处理。树的构建算法伪代码如下:

tree-root=node(T1)for i{2,,N} dotree-root.add(Ti)function add(triangle T)fa=f(a)fb=f(b)fc=f(c)if (abs(fa)<ϵ) thenfa=0if (abs(fb)<ϵ) thenfb=0if (abs(fc)<ϵ) thenfc=0if (fa0 and fb0 and fc0) thenif (negative subtree is empty) thennegative-subtree=node(T)elsenegative-subtree.add(T)else if (fa0 and fb0 and fc0) thenif (positive subtree is empty) thenpositive-subtree=node(T)elsepositive-subtree.add(T)elsecut triangle into three triangles and add to each side\begin{aligned} &\text{tree-root} = \text{node}(T_{1}) \\ &\textbf{for}\ i\in\{2,\cdots,N\} \ \textbf{do} \\ &\quad\text{tree-root.add}(T_{i}) \\ &\textbf{function}\ \text{add}(\text{triangle}\ T) \\ &\quad\begin{aligned} &\text{fa} = f(\vec{a}) \\ &\text{fb} = f(\vec{b}) \\ &\text{fc} = f(\vec{c}) \\ &\textbf{if}\ (\text{abs(fa)}<\epsilon)\ \textbf{then} \\ &\quad\text{fa} = 0 \\ &\textbf{if}\ (\text{abs(fb)}<\epsilon)\ \textbf{then} \\ &\quad\text{fb} = 0 \\ &\textbf{if}\ (\text{abs(fc)}<\epsilon)\ \textbf{then} \\ &\quad\text{fc} = 0 \\ &\textbf{if}\ (\text{fa}\leqslant 0\ \text{and}\ \text{fb}\leqslant 0\ \text{and}\ \text{fc}\leqslant 0)\ \textbf{then} \\ &\quad\begin{aligned} &\textbf{if}\ (\text{negative subtree is empty})\ \textbf{then} \\ &\quad\text{negative-subtree} = \text{node}(T) \\ &\textbf{else} \\ &\quad\text{negative-subtree.add}(T) \end{aligned} \\ &\textbf{else if}\ (\text{fa}\geqslant 0\ \text{and}\ \text{fb}\geqslant 0\ \text{and}\ \text{fc}\geqslant 0)\ \textbf{then} \\ &\quad\begin{aligned} &\textbf{if}\ (\text{positive subtree is empty})\ \textbf{then} \\ &\quad\text{positive-subtree} = \text{node}(T) \\ &\textbf{else} \\ &\quad\text{positive-subtree.add}(T) \end{aligned} \\ &\textbf{else} \\ &\quad\text{cut triangle into three triangles and add to each side} \end{aligned} \end{aligned}

切割三角形

树的构建作为预处理过程,一般不需要特别高效,而应尽量保证代码干净紧凑。一个技巧是让所有情况收缩为如下情形:c\vec{c} 在平面的一侧,其余两点在另一侧。为简单考虑,假定子树非空,上面 add\text{add} 方法的最后一种情况的伪代码如下:

if (fafc0) thenswap(fb, fc)swap(b,c)swap(fa, fb)swap(a,b)else if (fbfc0) thenswap(fa, fc)swap(a,c)swap(fa, fb)swap(a,b)compute Acompute BT1=(a,b,A)T2=(b,B,A)T3=(A,B,c)if (fc0) thennegative-subtree.add(T1)negative-subtree.add(T2)positive-subtree.add(T3)elsepositive-subtree.add(T1)positive-subtree.add(T2)negative-subtree.add(T3)\begin{aligned} &\textbf{if}\ (\text{fa}\ast\text{fc}\geqslant 0)\ \textbf{then} \\ &\quad\begin{aligned} &\text{swap(fb, fc)} \\ &\text{swap}(\vec{b}, \vec{c}) \\ &\text{swap(fa, fb)} \\ &\text{swap}(\vec{a}, \vec{b}) \end{aligned} \\ &\textbf{else if}\ (\text{fb}\ast\text{fc}\geqslant 0)\ \textbf{then} \\ &\quad\begin{aligned} &\text{swap(fa, fc)} \\ &\text{swap}(\vec{a}, \vec{c}) \\ &\text{swap(fa, fb)} \\ &\text{swap}(\vec{a}, \vec{b}) \end{aligned} \\ &\text{compute}\ \vec{A} \\ &\text{compute}\ \vec{B} \\ &T_{1} = (\vec{a}, \vec{b}, \vec{A}) \\ &T_{2} = (\vec{b}, \vec{B}, \vec{A}) \\ &T_{3} = (\vec{A}, \vec{B}, \vec{c}) \\ &\textbf{if}\ (\text{fc}\geqslant 0)\ \textbf{then} \\ &\quad\begin{aligned} &\text{negative-subtree.add}(T_{1}) \\ &\text{negative-subtree.add}(T_{2}) \\ &\text{positive-subtree.add}(T_{3}) \end{aligned} \\ &\textbf{else} \\ &\quad\begin{aligned} &\text{positive-subtree.add}(T_{1}) \\ &\text{positive-subtree.add}(T_{2}) \\ &\text{negative-subtree.add}(T_{3}) \end{aligned} \end{aligned}

上述伪代码中可能会切出一个面积接近于零的三角形,将 fa\text{fa}fb\text{fb}fc\text{fc} 中是否恰好有一个为零作为特例来处理,把原始三角形切为两个,可避免出现面积为零的三角形。但是栅格化本身也会处理面积为零的三角形,因此,也可以直接忽略这种情况。

交点坐标可以将直线参数方程 p(t)=a+t(ca)\vec{p}(t)=\vec{a}+t(\vec{c}-\vec{a}) 代入平面隐式方程 np+D=0\vec{n}\cdot\vec{p}+D=0 来得出:

A=a+tA(ca)tA=na+Dn(ca)\begin{gathered} \vec{A} = \vec{a} + t_{A}(\vec{c} - \vec{a}) \\ t_{A} = -\frac{\vec{n}\cdot\vec{a} + D}{\vec{n}\cdot(\vec{c} - \vec{a})} \end{gathered}

BSP 树的遍历效率正比于它的节点数,也就是三角形的个数,包括切割出的三角形。而三角形的个数依赖于建树时三角形的添加顺序。可以随机选取一些三角形顺序来尝试,将效果最好的保留下来。

切割三角形时,把一个三角形切成一个三角形和一个凸四边形会更加高效。对于只包含三角形的模型,这种做法没有什么价值,但它易于支持多边形构成的模型。

分片多维数组(Tiling Multidimensional Arrays)

高效利用存储层级结构(memory hierarchy)是在现代计算机架构中设计算法的关键。通过分片(tiling),或称为分块(bricking),可以确保多维数组中的数据按较好的方式排列。传统的高维数组都会存储为一维数组,这使得列(或行)方向的内存局部性(memory locality)较差。使用 tile 可以让行列方向的内存局部性更加平等。tile 应该和机器上的存储单元——如缓存行(cache line)——大小相仿。由于存在更加粗粒度的存储单元——如分页(pages),层次化分片(hierarchical tiling)也是有用的。

2D-array-contrast.png

二维数组的一级分片

假定 Nx×NyN_{x}\times N_{y} 数组分解为 n×nn\times n 方形 tile,tile 的数目为:

Bx=Nx/nBy=Ny/n\begin{aligned} B_{x} &= N_{x}/n \\ B_{y} &= N_{y}/n \end{aligned}

nn 无法整除 NxN_{x}NyN_{y} 时,需要填充(pad)数组。

tiled-2D-array.png

容易计算,元素 (x,y)(x,y) 所在 tile 的索引为:

bx=x÷nby=y÷n\begin{aligned} b_{x} &= x\div n \\ b_{y} &= y\div n \end{aligned}

元素 (x,y)(x,y) 在 tile 内的索引:

x=x mod ny=y mod n\begin{aligned} x' &= x\ \text{mod}\ n \\ y' &= y\ \text{mod}\ n \end{aligned}

其中,÷\divmod\text{mod} 分别是求商、求余运算。将上述公式整合起来,可以得到分片二维数组中元素 (x,y)(x,y) 的索引公式:

index=n2(Bxby+bx)+yn+x=n2[(Nx/n)(y÷n)+x÷n]+(y mod n)n+x mod n=Fx(x)+Fy(y)\begin{aligned} \text{index} &= n^{2}(B_{x}b_{y} + b_{x}) + y'n + x' \\ &= n^{2}[(N_{x}/n)(y\div n) + x\div n] + (y\ \text{mod}\ n)n + x\ \text{mod}\ n \\ &= F_{x}(x) + F_{y}(y) \\ \end{aligned}

其中,

Fx(x)=n2(x÷n)+x mod nFy(y)=n2(Nx/n)(y÷n)+(y mod n)n\begin{aligned} F_{x}(x) &= n^{2}(x\div n) + x\ \text{mod}\ n \\ F_{y}(y) &= n^{2}(N_{x}/n)(y\div n) + (y\ \text{mod}\ n)n \end{aligned}

Fx(x)F_{x}(x)Fy(y)F_{y}(y) 表格化之后可用于计算 index\text{index},而不必实时计算很多乘除法、模运算。

三维数组的二级分片

将三维数组按宏块(macrobrick)划分,每个宏块由 m×m×mm\times m\times m(brick)构成,每个块由 n×n×nn\times n\times n 个元素构成。可以得到分块三维数组中元素 (x,y,z)(x,y,z) 的索引公式:

index=Fx(x)+Fy(y)+Fz(z)\text{index} = F_{x}(x) + F_{y}(y) + F_{z}(z)

其中,

Fx(x)= [x÷(mn)](mn)3+[(x÷n) mod m]n3+(x mod n)Fy(y)= [y÷(mn)](mn)3Nxmn+[(y÷n) mod m]mn3+(y mod n)nFz(z)= [z÷(mn)](mn)3NxmnNymn+[(z÷n) mod m]m2n3+(z mod n)n2\begin{aligned} F_{x}(x) =\ &[x\div (mn)](mn)^{3} \\ &+ [(x\div n)\ \text{mod}\ m]n^{3} \\ &+ (x\ \text{mod}\ n) \\ F_{y}(y) =\ &[y\div (mn)](mn)^{3}\frac{N_{x}}{mn} \\ &+ [(y\div n)\ \text{mod}\ m]mn^{3} \\ &+ (y\ \text{mod}\ n)n \\ F_{z}(z) =\ &[z\div (mn)](mn)^{3}\frac{N_{x}}{mn}\frac{N_{y}}{mn} \\ &+ [(z\div n)\ \text{mod}\ m]m^{2}n^{3} \\ &+ (z\ \text{mod}\ n)n^{2} \end{aligned}