Ray-AABB问题：判断线段是否相交于轴对齐边界框(Axially Aligned Bounding Box, AABB)

摘要

Ray-AABB问题：判断线段是否相交于轴对齐边界框(Axially Aligned Bounding Box, AABB)
本文介绍了slab算法的实现，从一个简单实现开始，逐步优化slab算法，算法加速和并进行边界情况处理。

本文转载并翻译自：

FAST, BRANCHLESS RAY/BOUNDING BOX INTERSECTIONS
FAST, BRANCHLESS RAY/BOUNDING BOX INTERSECTIONS, PART 2: NANS

第一部分：简单实现

轴对齐包围盒（AABB）通常用于限制光线跟踪中的有限对象。射线/ AABB交点通常比精确的射线/对象交点计算得更快，并允许构造边界体积层次结构（BVH），减少了每条射线需要考虑的对象数量。（这将在以后的文章中详细介绍BVH。）这意味着光线跟踪器会花费大量时间来计算光线/ AABB交点，因此应高度优化此代码。

进行ray / AABB相交的最快方法是slab方法。想法是将盒子视为三对平行平面内的空间。射线被每对平行平面夹住，并且如果射线的任何部分保留下来，它将与盒子相交。

v1.0 简单实现版本

此算法的简单实现可能看起来像这样（为简洁起见，在两个维度中）：

bool intersection(box b, ray r) {double tmin = -INFINITY, tmax = INFINITY;if (ray.n.x != 0.0) {double tx1 = (b.min.x - r.x0.x)/r.n.x;double tx2 = (b.max.x - r.x0.x)/r.n.x;tmin = max(tmin, min(tx1, tx2));tmax = min(tmax, max(tx1, tx2));}if (ray.n.y != 0.0) {double ty1 = (b.min.y - r.x0.y)/r.n.y;double ty2 = (b.max.y - r.x0.y)/r.n.y;tmin = max(tmin, min(ty1, ty2));tmax = min(tmax, max(ty1, ty2));}return tmax >= tmin;
}

但是，这些部门要花很多时间。由于在进行射线追踪时，同一束射线是针对许多ABB进行测试的，因此预先计算射线方向分量的逆数是有意义的。如果我们可以依靠IEEE 754浮点属性，则这也可以隐式处理方向分量为零的边缘情况-例如，如果射线在射线的范围内，则tx1和tx2值将是相反符号的无穷大。平板，因此tmin和tmax保持不变。如果光线在平板之外，则tx1和tx2将是具有相同符号的无限大，从而使tmin == + inf或tmax == -inf，从而导致测试失败。

v1.1 版本

最终的实现如下所示：

bool intersection(box b, ray r) {double tx1 = (b.min.x - r.x0.x)*r.n_inv.x;double tx2 = (b.max.x - r.x0.x)*r.n_inv.x;double tmin = min(tx1, tx2);double tmax = max(tx1, tx2);double ty1 = (b.min.y - r.x0.y)*r.n_inv.y;double ty2 = (b.max.y - r.x0.y)*r.n_inv.y;tmin = max(tmin, min(ty1, ty2));tmax = min(tmax, max(ty1, ty2));return tmax >= tmin;
}

由于现代浮点指令集可以在没有分支的情况下计算最小值和最大值，因此可以进行无分支或除法的ray / AABB交集测试。

在我写的光线追踪器 Dimension 中对此的实现可以在这里看到。

第二部分：考虑线段与平板重合的情况

在第1部分中，我概述了一种用于计算光线和与轴对齐的边界框之间的交点的算法。依靠IEEE 754浮点属性消除分支的想法可以追溯到[1]中的Brian Smits，该实现被Amy Williams充实了。等。在[2]中。

v2.0

为了快速回顾一下，该想法是替换朴素的平板方法：

bool intersection(box b, ray r) {double tmin = -INFINITY, tmax = INFINITY;for (int i = 0; i < 3; ++i) {if (ray.dir[i] != 0.0) {double t1 = (b.min[i] - r.origin[i])/r.dir[i];double t2 = (b.max[i] - r.origin[i])/r.dir[i];tmin = max(tmin, min(t1, t2));tmax = min(tmax, max(t1, t2));} else if (ray.origin[i] <= b.min[i] || ray.origin[i] >= b.max[i]) {return false;}}return tmax > tmin && tmax > 0.0;
}

v2.1

v2.0的等价写法，但是比v2.0要更快:

bool intersection(box b, ray r) {double tmin = -INFINITY, tmax = INFINITY;for (int i = 0; i < 3; ++i) {double t1 = (b.min[i] - r.origin[i])*r.dir_inv[i];double t2 = (b.max[i] - r.origin[i])*r.dir_inv[i];tmin = max(tmin, min(t1, t2));tmax = min(tmax, max(t1, t2));}return tmax > max(tmin, 0.0);
}

这两种算法真的等效吗？我们已经依靠IEEE 754浮点行为消除了 ray.diri≠0ray.dir_i ≠ 0ray.diri=0 检查。当 ray.diri=±0ray.dir_i = ±0ray.diri=±0时，ray.dir_invi=±∞ray.dir\_inv_i =±∞ray.dir_invi=±∞ 。如果射线原点的iii坐标在框内，则意味着 b.mini<r.origini<b.maxib.min_i <r.origin_i <b.max_ib.mini<r.origini<b.maxi ，我们将得到 t1=−t2=±∞t1 = −t2 =±∞t1=−t2=±∞ 。由于对于所有 nnn 的max（n，−∞）=min（n，+∞）=nmax（n，-∞）= min（n，+∞）= nmax（n，−∞）=min（n，+∞）=n，因此tmin和tmax将保持不变。

另一方面，如果 $i 坐标在框外（坐标在框外（坐标在框外（r.origin_i <b.min_i$ 或 r.origini>b.maxir.origin_i> b.max_ir.origini>b.maxi），则 t1=t2=±∞t1 = t2 =±∞t1=t2=±∞，因此 tmin=+∞或tmax=−∞tmin = +∞或 tmax =-∞tmin=+∞或tmax=−∞ 。这些值之一将贯穿算法的其余部分，从而导致我们返回false。

不幸的是，上述分析有一个警告：如果射线正好位于平板上（r.origini=b.mini或r.origini=b.maxir.origin_i = b.min_i或r.origin_i = b.max_ir.origini=b.mini或r.origini=b.maxi），我们将拥有（说）
t1=(b.mini−r.origini)⋅r.dirinvi=0⋅∞=NaN\begin{aligned} t1 =& (b.min_i−r.origin_i)⋅r.dirinv_i \\ =& 0⋅∞ \\ =& NaN \end{aligned} t1===(b.mini−r.origini)⋅r.dirinvi0⋅∞NaN
这表现得比无穷大得多。正确处理此边缘（字面意义！）的情况取决于min（）和max（）的确切行为。

关于min()和max()

min（）和max（）的最常见实现可能是

#define min(x, y) ((x) < (y) ? (x) : (y))
#define max(x, y) ((x) > (y) ? (x) : (y))

这种形式是如此普遍，以至于被称为SSE / SSE2指令集中的最小/最大指令的行为。使用这些指令是从算法中获得良好性能的关键。话虽如此，这种形式具有涉及NaN的某些奇怪行为。由于所有与NaN的比较都是错误的，
min(x,NaN)=max(x,NaN)=NaNmin(NaN,x)=max(NaN,x)=x\begin{aligned} min(x,NaN) =& max(x,NaN) =& NaN \\ min(NaN,x) =& max(NaN,x) =& x \end{aligned} min(x,NaN)=min(NaN,x)=max(x,NaN)=max(NaN,x)=NaNx
这些操作既不传播也不抑制NaN。相反，当其中一个参数为NaN时，总是返回第二个参数。（关于有符号零，也有类似的奇数行为，但这并不影响此算法。）

相反，IEEE 754指定的最小/最大操作（称为“ minNum”和“ maxNum”）抑制NaN，如果可能的话总是返回一个数字。这也是C99的fmin（）和fmax（）函数的行为。另一方面，Java的Math.min（）和Math.max（）函数传播NaN，与大多数其他对浮点值的二进制运算保持一致。 [3]和[4]对野外的各种最小/最大实现进行了更多讨论。

存在问题

min（）和max（）的IEEE和Java版本提供一致的行为：恰好位于板上的所有光线均被视为不与盒子相交。很容易理解为什么要使用Java版本，因为NaN最终会污染所有计算，并使我们返回false。对于IEEE版本，min(t1,t2)=max(t1,t2)=±∞min(t1,t2)=max(t1,t2)=±∞min(t1,t2)=max(t1,t2)=±∞，与射线完全在盒子外面时相同。

（由于这是一个极端的情况，您可能想知道为什么我们不选择对边界上的光线返回true而不是false。事实证明，使用高效代码来实现此行为要困难得多。）

使用对SSE友好的最小/最大实现，行为是不一致的。平板上的某些光线会相交，即使它们完全位于盒子的另一个维度之外：

在上面的图像中，相机位于立方体的顶面的平面中，并且使用上述算法计算了交点。由于不正确的NaN处理，顶面超出了侧面。

v2.2 解决方案

当最多一个参数是NaN时，我们可以使用
minNum(x,y)=min(x,min(y,∞))maxNum(x,y)=max(x,max(y,−∞))\begin{aligned} minNum(x,y) =& min(x,min(y,∞)) \\ maxNum(x,y) =& max(x,max(y,−∞)) \end{aligned} minNum(x,y)=maxNum(x,y)=min(x,min(y,∞))max(x,max(y,−∞))
Thierry Berger-Perrin在[5]中采用了类似的策略，可以有效地进行计算

tmin = max(tmin, min(min(t1, t2), INFINITY));
tmax = min(tmax, max(max(t1, t2), -INFINITY));

在循环中。这样做也很好。由于CPU处理浮点特殊情况（无穷大，NaN，次正规）的速度较慢，因此改成下面这样速度要比上面快30％左右。

tmin = max(tmin, min(min(t1, t2), tmax));
tmax = min(tmax, max(max(t1, t2), tmin));

出于同样的原因，最好像这样展开循环，以避免处理不必要的更多无穷大：

bool intersection(box b, ray r) {double t1 = (b.min[0] - r.origin[0])*r.dir_inv[0];double t2 = (b.max[0] - r.origin[0])*r.dir_inv[0];double tmin = min(t1, t2);double tmax = max(t1, t2);for (int i = 1; i < 3; ++i) {t1 = (b.min[i] - r.origin[i])*r.dir_inv[i];t2 = (b.max[i] - r.origin[i])*r.dir_inv[i];tmin = max(tmin, min(min(t1, t2), tmax));tmax = min(tmax, max(max(t1, t2), tmin));}return tmax > max(tmin, 0.0);
}

很难理解为什么这个版本是正确的：x坐标的任何NaN都会传播到末端，而其他坐标的NaN将导致tmin≥tmax。在两种情况下，都返回false。

使用-O3的GCC 4.9.2，该实现每秒处理超过9300万条射线，这意味着即使没有向量化，运行时间也大约是30个时钟周期！

v2.3 更好的解决方案

可悲的是，这仍然比没有显式NaN处理的版本慢15％。而且由于遍历有界体积层次结构时通常使用此算法，因此最糟糕的事情是，在退化情况下遍历过多的节点。对于许多应用程序而言，这是非常值得的，并且如果正确实现射线/物体相交功能，则在实践中绝不应导致任何视觉伪影。为了完整起见，这是一个快速实施（1.08亿射线/秒），它不会尝试一致地处理NaN：

bool intersection(box b, ray r) {double t1 = (b.min[0] - r.origin[0])*r.dir_inv[0];double t2 = (b.max[0] - r.origin[0])*r.dir_inv[0];double tmin = min(t1, t2);double tmax = max(t1, t2);for (int i = 1; i < 3; ++i) {t1 = (b.min[i] - r.origin[i])*r.dir_inv[i];t2 = (b.max[i] - r.origin[i])*r.dir_inv[i];tmin = max(tmin, min(t1, t2));tmax = min(tmax, max(t1, t2));}return tmax > max(tmin, 0.0);
}

在[6]中给出了我用来测试各种cross（）实现的程序。在有关该主题的下一篇文章中，我将讨论低级实现的细节，包括矢量化，以从该算法中获得最大的性能。

[1]: Brian Smits: Efficiency Issues for Ray Tracing. Journal of Graphics Tools (1998).
[2]: Amy Williams. et al.: An Efficient and Robust Ray-Box Intersection Algorithm. Journal of Graphics Tools (2005).
[3]: https://groups.google.com/forum/#!topic/llvm-dev/-SKl0nOJW_w
[4]: https://ghc.haskell.org/trac/ghc/ticket/9251
[5]: http://www.flipcode.com/archives/SSE_RayBox_Intersection_Test.shtml
[6]: https://gist.github.com/tavianator/132d081ed4d410c755fd