从0开始制作游戏物理引擎（七）

本章说明了射线检测Raycast算法。

大部分Raycast都是使用向量方程的数值解直接得出，这种的公式大全^[1]都有记载不再赘述。

本文主要叙述：

射线与凸多面体求交的Slab方法(Cyrus-Beck, Liang-Barsky和Kay-Kajiya)
Gino编写的对任意凸体的Raycat算法的论文^[2]解释。此算法也用来做任意凸体之间的Sweep算法。

射线与三角形求交

使用Möller–Trumbore算法^[3]，源于《Fast Minimum Storage Ray-Triangle Intersection》：

一般的做法是：

求射线与三角形所在平面的交点$p$
看$p$在不在三角形内（可使用重心坐标法）

而Möller–Trumbore算法等于是将两步并做一步，并且论文中给了实现，此实现将所有的变量延后到需要计算的时候才算，并且给出了剔除三角形背面相交和双面三角形相交的。论文说对于TriangleMesh，可以减少25%~50%的内存（取决于顶点共享形式）。此算法可以同时得到交点和交点在三角面上的UV。

对于射线$O + tD$以及三角形$V_1,V_2,V_3$来说，其公共点为（三角形用重心坐标表示）：

$$ O + tD = (1 - u - v)V_0+ uV_1 + vV_2 $$

用矩阵表示为：

$$ \begin{bmatrix} -D & V_1 - V_0 & V_2 - V_0 \end{bmatrix} \begin{bmatrix} t \\ u \\ v \end{bmatrix} = O - V_0 $$

然后令$E_1 = V_1 - V_0, E_2 = V_2 - V_0, T = O - V_0$。根据Cramer法则^[4]可得：

$$ \begin{bmatrix} t \\ u \\ v \\ \end{bmatrix} = \frac{1}{\begin{vmatrix}-D & E_1 & E_2\end{vmatrix}} \begin{bmatrix} \begin{vmatrix} T & E_1 & E_2 \end{vmatrix}\\ \begin{vmatrix} -D & T & E_2 \end{vmatrix}\\ \begin{vmatrix} -D & E_1 & T \end{vmatrix} \end{bmatrix} $$

又由线性代数中:

$$ \begin{vmatrix} A & B & C \end{vmatrix} = -(A\times C)\cdot B = -(C \times B)\cdot A $$

可以将上式变化为：

$$ \begin{bmatrix} t \\ u \\ v \\ \end{bmatrix} = \frac{1}{(D\times E_2)\cdot E_1} \begin{bmatrix} (T\times E_1)\cdot E_2 \\ (D\times E_2)\cdot T \\ (T\times E_1)\cdot D \\ \end{bmatrix} = \frac{1}{P \cdot E_1} \begin{bmatrix} Q\cdot E_2 \\ P\cdot T \\ Q\cdot D \\ \end{bmatrix} $$

这里$P = (D\times E_2), Q = T \times E_1$。

之后就是使用这些变量进行计算。论文中给出了详细代码，十分通俗易懂，这里就不再细说了。

注：剔除背面的算法中，默认规定在右手系下，顶点顺时针（CW）绕序的一侧是背面。

射线与AABB求交

其实是图形学里面的方法，用于做线段裁切的。但在物理引擎求线段/射线/直线和凸体交点也很好用。

本文一共有如下算法：

由Cyrus-Beck在Generalized Two- and Three-Dimensional Clipping中编写的通用凸体裁剪方法
Kay-Kajiya于1986年在SIGGRAPH上发布的Ray Tracing Complex Scenes中提出的基于Cyrus-Beck算法的优化版本^[5]
而Kay-Kayjiya算法由Amy Williams在An Efficient and Robust Ray–Box Intersection Algorithm^[6]中提出一种更高效更鲁棒的实现，Github上也有代码^[7]

下文中使用的定义：

定义射线为$\vec{s} + \vec{d}t$。其中$\vec{s}$为起点，$\vec{d}$为方向
半空间$H_i(\vec{n}, \vec{p})$：$\vec{p}$为平面上一点，$\vec{n}$为法向量

Cyrus-Beck算法

将凸体的各个面视为半空间（朝向内侧视为正半空间），于是凸体变为这些半空间$H_i(\vec{n}, \vec{p})$的交集。

首先看射线射入半空间是“射入”还是“射出”。只需要看$\vec{d}$与$\vec{n}$是否同向。如果同向就是射入，反向就是射出。

论文的思想是：通过裁切射线，只保留射入和射出这一段，对所有半空间做完后，如果射线没有全被裁切完，那显然就和凸体相交。

首先计算射线和半空间平面的所有交点，以射入和射出进行分类。但由于交点都在射线上，那么只需求得$t$即可，于是有

$$ t_{in_i},t_{out_i} $$

那么显然，我最后裁切的线段一定满足$t \in [\max \{t_{in_i}\}, \min \{t_{out_i}\}]$。

那么只要

$$ \max \{t_{in_i}\} \lt \min \{t_{out_i}\} $$

就说明射线和凸体相交了。

而交点则是$\max \{t_{in_i}\}$和$\min \{t_{out_i}\}$这两个$t$所在点。

真正实现中有两个提前退出算法的小技巧：

只要发现任意$t_{in_i} \ge t_{out_i}$就直接退出，不需要真的算出全部$t$
当射线起点在平面外侧，且方向远离平面，说明射线永远不可能射入平面，那么肯定不相交，直接返回
- 先判断点是否在平面半空间$H(\vec{n}, \vec{p})$外：$(\vec{d} - \vec{p}) \cdot \vec{n} \le 0$
- 再看方向是否远离：$\vec{d} \cdot \vec{n} \le 0$

Kay-Kajiya算法和Williams的无分支版本

是Cyrus-Beck算法在3D空间中对盒体的特化。现在广泛用于图形学光线追踪中计算射线和BVH的交点（物理这边也会用到BVH）。

由于提出Slab概念，又叫Slab算法。

算法中只考虑射线与AABB盒子相交。与OBB相交可以通过仿射变换变成与AABB相交问题。

Slab就是：盒体的一对对面所夹的空间。所以盒体由三个slab的交集构成，每个轴各有两个垂直于此轴的Slab。我们称沿着负轴更的Slab为$Slab^-$，而沿着正轴更远的那个为$Slab^+$。

其算法理论就是Cyrus-Beck的理论，而其在Cyrus-Beck算法上的优化是工程上的。假设AABB盒子的中心是$\vec{c}$，半长为$\vec{hl}$。我们不需要再通过计算判断射线是射入面还是射出面：对于一对Slab，只需看$\vec{c} - \vec{s}$的$x,y,z$分量的符号，如果是正的，表示是从$Slab^-$射入$Slab^+$射出。负的就是反过来。

Kay-Kajiya算法的性能上的突破在于一条射线和许多AABB判交的场景下。

接下来看一下Williams的实现⁷。注意：Williams论文中并没有给出SIMD写法，他只是提了一嘴说可以利用SIMD。但PhysX，JoltPhysics中都有SIMD版本。SIMD版本主要用于一条射线和多个AABB（一般是4个因为SSE指令只有4-lane）求交，一般用于BVH空间划分算法而非Raycast中（这个后面要用到再看吧，裸SIMD指令看着也费劲）。如果只和一个AABB求交一般还是用标量版本。

首先是Ray的实现，里面缓存了一些信息：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
class Ray {
  public:
    Ray() { }
    Ray(Vector3 o, Vector3 d) {
      origin = o;
      direction = d;
      // 注意：这里没判断d的各个分量是不是0，这是故意的！
      // 这是Williams实现的核心设计：如果分量为0，结果是+/-inf（IEEE754标准规定）
      inv_direction = Vector3(1/d.x(), 1/d.y(), 1/d.z());
      sign[0] = (inv_direction.x() < 0);
      sign[1] = (inv_direction.y() < 0);
      sign[2] = (inv_direction.z() < 0);
    }
    Ray(const Ray &r) {
      origin = r.origin;
      direction = r.direction;
      inv_direction = r.inv_direction;
      sign[0] = r.sign[0]; sign[1] = r.sign[1]; sign[2] = r.sign[2];
    }

    Vector3 origin;
    Vector3 direction;
    Vector3 inv_direction;  // direction中各分量的倒数
    int sign[3];	// direction中各分量的符号
};

然后是判交算法：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
// Box只有一个成员变量Vector3 parameters[2]，其中parameters[0]存储min点，parameters[1]存储max点

// 这里参数t0和t1是射线自身的有效参数区间。通过这两个t可以处理所有直线类型的求交：
//  * 直线：t0 = -FLT_MAX, t1 = FLT_MAX
//  * 射线：t0 = 0, t1 = FLT_MAX
//  * 线段: t0 = 0, t1 = 1（或射线长度）
bool Box::intersect(const Ray &r, float t0, float t1) const {
  float tmin, tmax, tymin, tymax, tzmin, tzmax;

  // 求射线到x轴Slab的t
  // 这里一行解决了三种情况：
  // 1. r.origin在左Slab左边或者右边，这里都能算得最小值而不用判断射线方向
  // 2. 如果r.direction平行于面，那inv_direction.x()是+inf，那么：
  //    1. origin.x < min.x：tmin = +inf, tmax = +inf
  //    2. origin.x > max.x: tmin = -inf, tmax = -inf
  //    2. origin.x在slab内: tmin = -inf, tmax = +inf
  // 3. 如果r.direction.x()是inf，而乘法左边是0，那么tmin = tmax = NaN
  tmin = (parameters[r.sign[0]].x() - r.origin.x()) * r.inv_direction.x();
  tmax = (parameters[1-r.sign[0]].x() - r.origin.x()) * r.inv_direction.x();
    
  // 求射线到y轴Slab的t
  tymin = (parameters[r.sign[1]].y() - r.origin.y()) * r.inv_direction.y();
  tymax = (parameters[1-r.sign[1]].y() - r.origin.y()) * r.inv_direction.y();
    
  // 提前退出，无论值是inf, NaN还是正常值都能正确处理，可以自己试一试
  if ( (tmin > tymax) || (tymin > tmax) ) 
    return false;
  if (tymin > tmin)
    tmin = tymin;
  if (tymax < tmax)
    tmax = tymax;
    
  // 再求z轴Slab的t
  tzmin = (parameters[r.sign[2]].z() - r.origin.z()) * r.inv_direction.z();
  tzmax = (parameters[1-r.sign[2]].z() - r.origin.z()) * r.inv_direction.z();
    
  // 继续判断
  if ( (tmin > tzmax) || (tzmin > tmax) ) 
    return false;
  if (tzmin > tmin)
    tmin = tzmin;
  if (tzmax < tmax)
    tmax = tzmax;
  return ( (tmin < t1) && (tmax > t0) );
}

这就是为什么说Williams的算法是鲁棒的，他可以处理任意的数值。

Gino的论文解释

Gino在论文²中提出了一种基于GJK的，针对任意凸体的射线检测算法。但此算法只能得到$t$最小的那个交点及对应法线。

为了和论文统一下面还是使用论文的符号。

定义射线$R = \{ \vec{s} + \lambda \vec{r} : \lambda \ge 0\}$。其中$\vec{s}$是起点，$\vec{r}$是方向。

用$C$表示任意凸体。注意凸体上某点的法线不一定是唯一的（比如凸多面体顶点处法线）。只有某个点$\vec{p}$附近光滑（或者用分析的说法是可微）的此点的法线才是唯一的。

通过法线定义可以得到如下方程：

$$ \vec{n} \cdot (\vec{x} - \vec{p}) \le 0, \forall \vec{x} \in C $$

那么反过来就是所有不在$C$中的点$\vec{x}$可以定义为：

$$ \begin{aligned} & \vec{n} \cdot (\vec{x} - \vec{p}) \le 0 & &(1) \end{aligned} $$

那么将$\vec{x}$换成射线上任意点$\vec{x} = \vec{s} + \lambda \vec{r}$，就可以得到既在射线上又不在$C$上的点的条件：

$$ \begin{aligned} \lambda \vec{n} \cdot \vec{r} \gt \vec{n} \cdot (\vec{p} - \vec{s}) & & (2) \end{aligned} $$

显然，射线与$C$的交点也满足这个方程，或者说的更明确一些：

如果$\vec{n} \cdot \vec{r} \gt 0$，那么所有满足$\lambda \gt \frac{\vec{n}\cdot(\vec{p} - \vec{s})}{\vec{n}\cdot \vec{r}}$的点就都不是交点。此时$\frac{\vec{n}\cdot(\vec{p} - \vec{s})}{\vec{n}\cdot \vec{r}}$是$\lambda$的上界
如果$\vec{n} \cdot \vec{r} \lt 0$，那么所有满足$\lambda \lt \frac{\vec{n}\cdot(\vec{p} - \vec{s})}{\vec{n}\cdot \vec{r}}$的点就都不是交点。此时$\frac{\vec{n}\cdot(\vec{p} - \vec{s})}{\vec{n}\cdot \vec{r}}$是$\lambda$的下界
如果$\vec{n} \cdot \vec{r} = 0$且$\vec{n} \cdot (\vec{p} - \vec{s}) \lt 0$，那么这条射线就不和$C$相交

然后Gino又使用他的经典逼近思路：通过不停地裁切$\lambda$可能取值的区间，将交点不可能存在的区间从中移除，直到找到交点。逼近的方法如下：

假设现在的$\lambda$下界是$\lambda_i$，那么有点$\vec{x_i} = \vec{s} + \lambda_i \vec{r}$。此时$\lambda \in [0, \lambda_i)$的所有点就被排除了。然后令$\vec{c_i} \in C$是距离$\vec{x_i}$的最近点。那么显然：

如果$\vec{c_i} = \vec{x_i}$，那么最近点就是$\vec{c_i}$（算法结束条件）
如果$\vec{c_i}$只是$C$边界上一点，那么$\vec{n_i} = \vec{x_i} - \vec{c_i}$是$\vec{c_i}$这一点的法线。那么满足$(2)$式：
$$ \begin{aligned} & \lambda \vec{n_i} \cdot \vec{r} \gt \vec{n_i} \cdot (\vec{c_i} - \vec{s}) \\ & \lambda \vec{n_i} \cdot \vec{r} \gt \lambda_i \vec{n_i} \cdot \vec{r_i} - ||\vec{n_i}^2|| \\ \end{aligned} $$
那么：
1. 如果$\vec{n_i} \cdot \vec{r} = 0$，那么射线就不和$C$相交（因为不等式右边只有$-||\vec{n}_i||^2 \le 0$恒成立）。
2. 如果$\vec{n}\cdot\vec{r} \ne 0$。令$\lambda_{i+1} = \lambda_i - \frac{||\vec{n}_i||^2}{\vec{n_i}\cdot \vec{r}}$。此时：
  1. 若$\vec{n} \cdot \vec{r} \gt 0$，显然$\lambda_{i+1} \lt \lambda_i$，那么根据假设$\lambda \in [0, \lambda_i)$中的点都不是交点，所以$\lambda_{i+1}$早已被剔除，此时射线和$C$就没有交点。
  2. 若$\vec{n} \cdot \vec{r} \lt 0$，显然$\lambda_{i+1} \gt \lambda_i$，那么$\lambda_{i+1}$等于是更大的$\lambda$下界，这样就增大了下界。
  那么就可以通过不停迭代$\lambda_i$来不停增大下界，直到$\vec{x_i} = \vec{c_i}$找到交点或判断根本不相交。

于是算法流程就是：

$$ \begin{aligned} & \lambda = 0 \\ & \vec{x} = \vec{s} \\ & \vec{n} = \vec{0} \\ & \vec{c} = \text{"the point of C closest to }\vec{x}\text{"} \\ & \text{while not "}\vec{x} \text{ is not close enough to }\vec{c}\text{" do} \\ & \text{begin} \\ & \quad \vec{n} = \vec{x} - \vec{c} \\ & \quad \text{if }\vec{n} \cdot \vec{r} \ge 0 \text{ then return false} \\ & \quad \text{else} \\ & \quad \text{begin} \\ & \qquad \lambda = \lambda - \frac{||\vec{n}||^2}{\vec{n}\cdot \vec{r}} \\ & \qquad \vec{x} = \vec{s} + \lambda \vec{r} \\ & \qquad \vec{c} = \text{"the point of C closest to }\vec{x}\text{"} \\ & \quad \text{end} \\ & \text{end} \\ \end{aligned} $$

原论文随后证明了$\lambda_i \lt \lambda_{i+1} \le \lambda_{hit}$是全局收敛的充要条件。这里就不写了。

$\vec{x_i}$足够接近$\vec{c_i}$的条件可以使用他两的距离去判断，而每轮迭代中找到$\vec{c}$就可以使用GJK。

但每轮迭代都要从头算GJK，这是一个很慢的事。Gino随后将GJK与这整个算法合并得到一个更高效的形式：

首先，$\vec{n_i} = \vec{x_i} - \vec{c_i}$可以改成用Minkowski差和支撑映射的计算方式：$\vec{n_i} = v(\{\vec{x_i}\} - C)$。

然后，可以看到原算法并不能给出点$\vec{c_i}$处的法线（他只是靠法线逐渐逼近）。但我们可以通过每一步迭代出的$\vec{v}_k$逐渐逼近$\vec{n_i}$。此处$\vec{v_k}$是Simplex上距离原点最近点，那么通过此点找到的支撑点记为$\vec{p_k} = s_C(\vec{v_k})$。显然，$\vec{v_k}$是$\vec{p_k}$的法线，因为其满足$(1)$式中对法线的定义$\vec{v_k} \cdot (x - \vec{p_k}) \le 0$。那么将其带入$(2)$式就可以得到被排除的点的新条件：

$$ \lambda \vec{v_k}\cdot \vec{r} \gt \vec{v_k} \cdot (\vec{p_k} - \vec{s}) $$

令$\vec{x_i} = \vec{s} + \lambda_i \vec{r}$是目前最大下界，那么将$\vec{s} = \vec{x_i} - \lambda_i \vec{r}$带入上式可得：

$$ \begin{aligned} & \lambda \vec{v_k}\cdot \vec{r} \gt \lambda_i \vec{v_k}\cdot \vec{r} -\vec{v_k}\cdot \vec{w_k} \\ & 其中\vec{w_k} = s_{\{x_i\} - C}(-\vec{v_k}) = \vec{x_i} - \vec{p_k} \end{aligned} $$

那么如法炮制对$\vec{v_k}\cdot \vec{r}$进行分析：

$\vec{v_k}\cdot \vec{r} = 0$且$\vec{v_k}\cdot \vec{w_k} \gt 0$，那么射线无交点
$\vec{v_k}\cdot \vec{r} \gt 0$那么$\lambda_i - \frac{\vec{v_k}\cdot \vec{w_k}}{\vec{v_k}\cdot \vec{r}}$同样位于不相交区间。那么此时射线无交点
$\vec{v_k}\cdot \vec{r} \lt 0$那么$\lambda_i - \frac{\vec{v_k}\cdot \vec{w_k}}{\vec{v_k}\cdot \vec{r}}$就是新下界。

那么现在算法就变成（这里要对照GinoGJK的论文过程看一下可能才能看明白怎么融合的）：

$$ \begin{aligned} & \lambda = 0 \\ & \vec{x} = \vec{s} \\ & \vec{n} = \vec{0} \\ & \vec{v} = \vec{x} - \vec{p}, \text{any } \vec{p} \in C \\ & P = \emptyset \\ & \text{while } ||\vec{v}||^2 \gt \epsilon^2 \text{ do} \\ & \text{begin} \\ & \quad \vec{p} = s_C(\vec{v}) \\ & \quad \vec{w} = \vec{x} - \vec{p} \\ & \quad \text{if } \vec{v}\cdot \vec{w} \gt 0 \text{ then} \\ & \quad \text{begin} \\ & \qquad \text{if } \vec{v} \cdot \vec{r} \ge 0 \text{ then return false} \\ & \qquad \text{else} \\ & \qquad \text{begin} \\ & \qquad \quad \lambda = \lambda - \frac{\vec{v}\cdot \vec{w}}{\vec{v} \cdot \vec{r}} \\ & \qquad \quad \vec{x} = \vec{s} + \lambda \vec{r} \\ & \qquad \quad \vec{n} = \vec{v} \\ & \qquad \text{end} \\ & \quad \text{end} \\ & \quad Y = P \cup \{\vec{p}\} \\ & \quad \vec{v} = v(\operatorname{co}(\{\vec{x}\} - Y)) \\ & \quad P = \text{"smallest }X \subseteq Y \text{ such that } \vec{v} \in \operatorname{co}(\{\vec{x}\} - X) \\ & \text{end} \\ & \text{return true} \\ \end{aligned} $$

论文最后提到还能做如下优化：

原始GJK中，同一支撑点在多次迭代中有可能重复出现（ill condition），而本算法不能再依赖于这个检查。Gino在32位单精度实现中发现不需要再进行这个检查。
在Johnson子距离算法中，行列式$\Delta_i(Y)$可能由于数值精度问题出错。此时就直接退出循环，将当前算出的$\vec{x}$作为结果返回（类似GinoGJK中的处理方法）
终止容差$\epsilon$（第6行的那个）可以用相对容差$||v_k|^2 \le \epsilon_{tol} \max \{||\vec{x} - \vec{p}||^2: \vec{p} \in P_k\}$，其中$\epsilon_{tol}$是std::numeric_limits<float>::epsilon() * 10。
算法可以提前退出：在Raycast中射线都是有长度的，一旦$\lambda$超过长度$\lambda_{max}$就直接退出算法，表示无交点
帧间相关性：判断不相交时返回的$\vec{v}$是射线和$C$的分离轴，若物体连续移动，则可以作为下帧的热启动（这在连续碰撞检测中有用）。注意此时$\vec{v_0}$不一定在$\vec{x_0} - C$内，所以得跳过第一轮的$while$判断。

此算法也可以作为两个物体$C_1, C_2$的Sweep的算法。显然，此问题可以转化为一原点为起点，方向为$\vec{d}$的射线与$C_1 - C_2$物体的射线检测。而绝大部分物体的Sweep算法都是此算法（除了有数值解的球和球，球和胶囊）。

在线演示

下面是目前文章提到的各种查询方法演示：