three-mesh-bvh 源码阅读(4) rayCast通过前面的几篇分析，对于 BVH 节点树的构建过程应该已经很

前言

通过前面的几篇分析，对于 BVH 节点树的构建过程应该已经很清楚了，本篇我们继续探索，看下有了几何体加速结构后如何加速计算，直观地感受下这种把问题规模缩小带来的计算速度提升。

本篇没有涉及到太多新的知识点，更多是技巧的使用，也就是之前在构建过程中反复提到的遍历，通过遍历已构建好的加速结构去实现加速，然后和之前一样需要从 TypedArray 里对数据进行操作。截图为这个库的封面图，在 BVH 核心算法完成的情况下这个效果其实也比较简单, 无非是根据碰撞结果实时绘制出线和点。最后花了一点篇幅介绍了两种碰撞检测的底层算法。

写作的时候使用的版本是最新的 "version": "0.7.3",, 事实上由于核心的东西不会变化，更多的是 API 形式 上的调整，因而系列写作过程中并没有去强调版本，系列的读者可以轻易地识别变化。

解析

入口文件是 /examples/raycast.html , 这个入口文件对应的 js：

作者通过在 ExtensionUtilitiess 写一个 acceleratedRaycast 方法来替换 Mesh 上的 raycast 方法，并和原来保持 API 的兼容来使用加速结构加速碰撞：

convertRaycastIntersect 这个函数的作用是把结果乘以模型的世界矩阵转换成正确的坐标，否则结果就是基于原点的结果。

再来看下 three 的 Mesh 代码的位置, 查看高频使用的 Raycaster 类，都是非常熟悉的 APIs，注册在类方法上：

点击跳转到 three 的 Raycaster.js, 查看最后的 intersectObject 方法:

原来根据参数，Raycaster 会决定是否递归调用 Mesh 类方法 .raycast。

最后经过梳理，把整个调用关系画出来：

这里面有一个 rayCastFirst 方法，字面意思就是首次碰撞，它和 rayCast 的区别，我们下文会分析。

碰撞检测算法

查看 three.js 源码的 Mesh 类，最终你会找到检测碰撞的方法是 Ray 类的 intersectTriangle 方法，至此，你会发现 three.js 自带的碰撞检测低效的原因是因为没有使用任何加速结构，仅仅对点是否在几何体的 Box 和 Sphere 包围盒内做检测（Mesh.raycast），前者只需要判断点是否在 aabb内，后者只需要判断点和圆心距离。这样虽然计算量小，但是显然对于高面数，复杂拓扑结构的几何体十分低效。

所以把加速结构应用到 three 原有的碰撞检测上那么效率应该就会大幅提升。通过前文的分析我们已经清楚地知道作者是通过保持兼容的方式替换了 Mesh 类的 raycast 方法从而应用已生成的加速结构。

接下来我们通过分析 MeshBVH 类上的 rayCast 和 rayCastFirst 方法然后一步步过渡到最基本的射线和三角形求交算法，这个算法被巧妙地设计过，相比于射线和平面方程求交再判断位于三角形内的算法，计算复杂度显著降低，这在许多场景非常重要，比如各种运用到光线追踪的渲染算法，检测算法等，在 AI 的帮助下我得以快速地了解这整个过程，由于篇幅原因，这个算法的详细介绍写在文末。

对比下rayCast 和 rayCastFirst ：

一下子就能看出来两函数的作用和区别了，raycastFirst 是返回射线和 bvh 的 buffer数据（三角形）的首个碰撞结果，而 raycast 是会寻找所有求交的结果。后面的 .face.materialIndex 是当该 Mesh 使用了 arrayMaterials （材质组）的时候把求教的结果（三角面）所使用的相应的材质编号记录下来，为后面的光线追踪程序提供便利。显然，文张开头的那个动画效果只需要求解首个结果，求解所有结果无非是遍历 buffer 记录所有结果，接下来我们把精力放到 raycastFirst 上，继续往下走。

我们定位到 raycastFirstFunc, 这里会判断 indirect, 看过系列文章的读者或许可以想起这个作用，它是为 BVH 数据额外生成了一个 buffer 避免直接在 geometryBuffer 里维持 BVH 结构。点击函数定位到如下位置，这个 generated 和 template 字样很熟悉，之前的我文章也有提过，作用是让 rollup 根据 template 内容条件编译不同的代码出来，作者大量使用这个方法来区分 direct 和 indirect 不同的代码。

raycastFirst.generated.png

这里对于 buffer 数据的遍历如果看过之前文章的朋友应该已经非常熟悉了，这里不做大篇幅重新分析，让我们来看看几个关键字： BufferStack, offset, count 分别代表当前 BVH数据（分组 Mesh 有多颗树），buffer 里的数据起点和计数（三角形）。这个闭包 _raycastFirst 函数对左右节点进行遍历，遇到第一个碰撞结果停止遍历。看下左右节点的遍历：

射线与 AABB 求交：

常规做法

常规做法是对 AABB 的六个面分别做碰撞检测，那么本质上就是解是解线性方程组，需要解 6 次再根据射线方向俩判断 t_min 和 t_max。那么对于这样的计算，我们简单分析一下时间复杂度：

求 t: $O + t\vec{d} = Ax + By + Cz$ ，需要若干次加法和一次浮点除法，复杂度为 $O(n)$
有 6（n） 个面, 复杂度为 $O(n)$

所以常规做法的复杂度是 $O(n^2)$ 。

高效的方法

标红的是射线与 AABB 包围盒的碰撞检测，这个方法我们看这个图：

写过比较多 3D 程序的读者比较熟悉的一个话题了，这里再复习一下。

上图里对于射线与 AABB 包围盒的碰撞检测，二维的规则可以拓展到三维。这个方法是怎么样的来简单分析下：以 XOY 平面为例，x 轴向的两个平面分别是 $x = x0$ 和 $x = x1$ , 这两平面称为 slabs 平面。y 轴同理。

首先，需要求 $t_{minx}, t_{miny}, t_{maxx}, t_{maxy}$ ，由于是 AABB 轴向包围盒，那么很容易可以求得：

\begin{cases} & x_{min} = x_p + t_{xmin} x_d \\ & t_{xmin} = (x_{min} - x_p) / x_d \\ & y_{min} = y_p + t_{ymin} y_d \\ & t_{ymin} = (y_{min} - y_p) / y_d \\ \end{cases}

其中 $(x_{p}, y_{p})$ 是射线原点的坐标。

由于如果与盒子碰撞，会穿过盒子产生两个交点，分别把离射线原点近的标注为 $t_{min}$ , 远端的标注为 $t_{max}$ ，那么可以观察到射线总是先碰到 slab 平面再接触到盒子，所以有：

\begin{cases} t_{min} &= max(t_{ymin}, t_{xmin}) \\ t_{max} &= min(t_{ymax}, t_{xmax}) \\ \end{cases}

最后，判断出界的情况：

\begin{cases} & \text{ if } t_{min} > t_{max} , \space \text{超出范围} \\ & \text{ if } t_{max} < 0, \space \text{反方向出界} \\ \end{cases}

最后算下时间复杂度，每次比较的复杂度只有几次加法，两次比较，最坏的情况是检测 3(n)次, 所以时间复杂度是 $O(n)$ 。

显然第二种方法高效得多，还省去了一次浮点除法。

遍历的顺序

看闭包函数 _rayCast 函数的遍历代码：

	...
	else {
	
			// consider the position of the split plane with respect to the oncoming ray; whichever direction
			// the ray is coming from, look for an intersection among that side of the tree first
			const splitAxis = SPLIT_AXIS( nodeIndex32, uint32Array );
			const xyzAxis = _xyzFields[ splitAxis ];
			const rayDir = ray.direction[ xyzAxis ];
			const leftToRight = rayDir >= 0;
	
			// c1 is the child to check first
			let c1, c2;
			if ( leftToRight ) {
	
				c1 = LEFT_NODE( nodeIndex32 );
				c2 = RIGHT_NODE( nodeIndex32, uint32Array );
	
			} else {
	
				c1 = RIGHT_NODE( nodeIndex32, uint32Array );
				c2 = LEFT_NODE( nodeIndex32 );
	
			}
	
			const c1Intersection = intersectRay( c1, float32Array, ray, _boxIntersection );
			const c1Result = c1Intersection ? _raycastFirst( c1, bvh, side, ray ) : null;
	
			// if we got an intersection in the first node and it's closer than the second node's bounding
			// box, we don't need to consider the second node because it couldn't possibly be a better result
			if ( c1Result ) {
	
				// check if the point is within the second bounds
				// "point" is in the local frame of the bvh
				const point = c1Result.point[ xyzAxis ];
				const isOutside = leftToRight ?
					point <= float32Array[ c2 + splitAxis ] : // min bounding data
					point >= float32Array[ c2 + splitAxis + 3 ]; // max bounding data
	
				if ( isOutside ) {
	
					return c1Result;
	
				}
	
			}
	
			// either there was no intersection in the first node, or there could still be a closer
			// intersection in the second, so check the second node and then take the better of the two
			const c2Intersection = intersectRay( c2, float32Array, ray, _boxIntersection );
			const c2Result = c2Intersection ? _raycastFirst( c2, bvh, side, ray ) : null;
	
			if ( c1Result && c2Result ) {
	
				return c1Result.distance <= c2Result.distance ? c1Result : c2Result;
	
			} else {
	
				return c1Result || c2Result || null;
	
			}
	
		}
	}

这里面有一个 leftToRight 的标志位判断，由于我用 GeoGebra 画了个示意图，蓝色和黄色分别代表某非叶子节点的左节点和右节点，MN 代表了射线方向 rayDir > 0 的情况，这种情况下左节点碰撞优先级高于右节点，所以由射线方向上距离射线原点更近的节点优先遍历，反之亦然。

还有一点非常重要，那就是当 c1Result 也就是优先遍历的这个节点如果判断产生了碰撞那么就不需要遍历另一个节点了，提升了遍历效率。

结合注释，这段代码也就搞明白了。

射线与三角形求交

最后来到判断交点的最后步骤，射线与三角形交点的计算。由于为了保持兼容和非侵入式设计，这个库的调用链还蛮长的，由 intersectClosestTri 函数点进去会一路到达 ThreeRayIntersectUtilities.js 文件，这本是 three.js 底层工具函数，作者导出来做了适配，而最终执行射线与三角形交点的计算是在 Ray.js 这个 three.js 的文件里。

ray.intersectTriangle.png

传统办法

大致分为两步，第一步是求交点，也就是求方程 $O + t\vec{d} = Ax + By + Cz$ ，这个和上面的一样，这里更具体地给出步骤方便读者比较。步骤及其复杂度：

计算平面方程：O(1)

两个向量叉乘：约9次乘法和6次加法
计算D值：3次乘法和1次加法

求解 t 值：O(1)

点乘计算：6次乘法和5次加法
1次除法

计算交点坐标：O(1)

3次乘法和3次加法

判断点是否在三角形内：O(1)

使用边加权法：约27次乘法和18次加法

总体复杂度：O(1)，大约需要45-50次乘法，30-35次加法，1次除法

Möller-Trumbore 算法

推导

用边加权形式：

\mathbf{Q} + t\mathbf{D} = b_1\mathbf{E}_1 + b_2\mathbf{E} \quad ..... \quad(1)

其中：

$\mathbf{Q}$ 是射线起点到三角形一个顶点的向量, $\mathbf{Q} = \mathbf{O} - \mathbf{V_0}$
$\mathbf{D}$ 是射线方向, $\mathbf{O}$ 是射线原点， $t$ 是射线参数
$\mathbf{V_0}$ 是三角形某个顶点
$\mathbf{E}_1$ 和 $\mathbf{E}_2$ 是三角形的两条边向量
$b_1$ 和 $b_2$ 是交点的边加权分量

写成行列式：

\left[\begin{array}{lll} -\mathbf{D} & \mathbf{E_1} & \mathbf{E_2} \end{array}\right] \left[\begin{array}{l} \mathbf{t} \\ \mathbf{b1} \\ \mathbf{b2} \end{array}\right]=\mathbf{Q} \quad ..... \quad (2)

此时三个未知量，三个等式，可以用熟悉的克莱姆法则求解。这里通过变换来对 $(1)$ 等式进行变换分别求出三个未知量，以 $t$ 为例：

首先，等式两边同时叉乘 $(\mathbf{E_1} \times \mathbf{E_2})$ ：

(\mathbf{Q} + t\mathbf{D} ) \times (\mathbf{E_1} \times \mathbf{E_2}) = b_1\mathbf{E_1} + b_2\mathbf{E_2} \times (\mathbf{E_1} \times \mathbf{E_2})

右边等于 $\mathbf{0}$ ，等式化为：

\mathbf{Q} \cdot (\mathbf{E_1} \times \mathbf{E_2}) + t(\mathbf{D} \cdot (\mathbf{E_1} \times \mathbf{E_2})) = \mathbf{0}

令

$\mathbf{N} = \mathbf{E_1} \times \mathbf{E_2}$ 则

t = - \mathbf{Q} \cdot \mathbf{N} / \mathbf{D} \cdot \mathbf{N}

由于 $\mathbf{D} ， \mathbf{N}$ 可能不同号，引入 $sign$ 来表达，这个符号返回 1 或者 -1, 两向量同符号返回 1，不同号则返回 -1, 于是有：

$|\mathbf{D} \cdot \mathbf{N}| \cdot t = -\text{sign}(\mathbf{D} \cdot \mathbf{N}) \cdot \mathbf{Q} \cdot \mathbf{N}$

另外两个等式通过叉乘 $\mathbf{E_1}$ 或者 $\mathbf{E_2}$ 也可以变成以下形式：

$|\mathbf{D} \cdot \mathbf{N}| \cdot b_1 = \text{sign}(\mathbf{D} \cdot \mathbf{N}) \cdot \mathbf{D} \cdot (\mathbf{Q} \times \mathbf{E_2})$

$|\mathbf{D} \cdot \mathbf{N}| \cdot b_2 = \text{sign}(\mathbf{D} \cdot \mathbf{N}) \cdot \mathbf{D} \cdot (\mathbf{E_1} \times \mathbf{Q})$

这也就是 three.js 里这部分源码注释的公式来源：

three.moller-trumbore.png

它的步骤及其复杂度如下：

计算关键向量：O(1)

1次向量叉乘：9次乘法和6次加法
2次向量点乘：6次乘法和4次加法

计算 t, b1, b2：O(1)

3次标量除法
约9次乘法和6次加法

总体复杂度：O(1)，大约需要24次乘法，16次加法，3次除法

计算量这部分内容由 AI 辅助写作，如有错误，欢迎勘误，水平有限，不胜感激。

可以看到通过代数变换的 Möller-Trumbore 显著减少了计算量，避免使用计算量很大的克莱姆法则求解，这在数值稳定性很有帮助（浮点数算得步骤越多精度越不稳定）。

结语

前面的系列文章：

写文章不易，点赞收藏是最好的支持！