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

这一章描述了GJK和EPA算法。首先描述了市面上大部分资料说明的通俗的GJK和EPA算法。然后给出了GJK论文的解析。

GJK算法

Minkowski和与差

首先定义“Minkowski和”与“Minkowski差”：

设有任意两个几何体（对凹凸性不作要求）$P_1$, $P_2$。

$$ P_1 \oplus P_2 = \{p_1 + p_2 | p_1 \in P_1, p_2 \in P_2 \} $$

即两个几何体中所有点和的集合。

而Minkowski差则相反：

$$ P_1 \ominus P_2 = \{p_1 - p_2 | p_1 \in P_1, p_2 \in P_2 \} $$

即两个几何体中所有点差的集合。

这两个定义在碰撞检测中十分有用。Minkowski和会在后面扫略体算法中广泛使用。而Minkowski差则在GJK中被使用。

Minkowski差有三个性质：

如果两个物体相交，其Minkowski差一定包含原点

显然，总可以在两个物体上分别找到两个点，他们在相交区域内且坐标相同，那么他们的Minkowski差就是原点。

如果两个物体不相交，其距离为Minkowski差到原点的距离

因为距离的定义是$||\vec{a} - \vec{b}||$而Minkowski差为$\vec{a} - \vec{b}$。那么两几何体的距离显然是$\min(||\vec{a} - \vec{b}||) = \min(||(\vec{a}-\vec{b}) - (0, 0, 0)||)$为最后的Minkowski差中到原点的距离。

凸体的Minkowski和/差也是凸体

$$ \begin{aligned} & \forall x_1,x_2 \in Polygon \Rightarrow \\ & x = \lambda x_1 + (1 - \lambda) x_2 \in Polygon, \lambda \in [0, 1] \end{aligned} $$$$ \begin{aligned} & \forall a_1, a_2 \in Polygon_1, \forall b_1, b_2 \in Polygon_2 \\ & a_1 - b_1 = p_1 \in Polygon_1 \oplus Polygon_2 \\ & a_2 - b_2 = p_2 \in Polygon_1 \oplus Polygon_2 \\ & 且 \\ & \lambda p_1 + (1-\lambda) p_2 = \lambda(a_1 + b_1) + (1 - \lambda)(a_2 + b_2) = [\lambda a_1 + (1-\lambda)a_2] + [\lambda b_1 + (1-\lambda)b_2] = a^* + b^* \in Polygon_1 \oplus Polygon_2 \end{aligned} $$

正如上面所写，$Polygon_1 \oplus Polygon_2$内的点满足凸体定义，那么$Polygon_1 \oplus Polygon_2$就是凸体。

Minkowski差也同理可证。

有了这些性质，我们就可以有构造两个凸体的Minkowski和/差的算法（这个算法的动画过程可以看“参考3”）。Minkowski和的构造是：

按一定顺序（比如顺时针），在凸体的中心朝其顶点的方向（两个凸体都要），找$Polygon_1$和$Polygon_2$的Support点，将其做Minkowski和。最后得到的点集全部在$Polygon_1 \oplus Polygon_2$所成凸体的顶点。并且这些顶点也是按照指定顺序（这里是顺时针）排列。这样我们按顺序将这些点连起来就找到最后的凸体了。

这算法很容易理解：Support点是两个凸体某个方向最远的点，所以他们的和一定在这个方向最大，从而构成Minkowski和凸体中的顶点。而只找凸体中心朝顶点的方向是为了让算法在有限步骤内完成。

那么Minkowski差只是反过来：

按一定顺序，在凸体的中心朝其顶点的方向（设为$\vec{d}$），找$Polygon_1$中在$\vec{d}$上的Support点以及和$Polygon_2$在$-\vec{d}$的Support点，将其做Minkowski差。最后得到的点集全部在$Polygon_1 \ominus Polygon_2$所成凸体的顶点。并且这些顶点也是按照指定顺序排列。这样我们按顺序将这些点连起来就找到最后的凸体了。

小注意点：严格来说，这里找Support点是隐含了用上一节说的Support点算法寻找。因为从数学意义上来说，Support点不一定是凸体端点（比如垂直于边的方向找到的点是在边上），从而两个点相加不构成最终凸体的顶点。

但我们也不会真的用这个算法求出最终的凸体，然后算原点到凸体的距离。GJK只用了这里面的思想。

GJK算法

先说明2D GJK算法，3D的只是其扩展：

随意选取一个方向，找到此方向上Support点的Minkowski差$P_1$
然后选取朝向原点的方向$\vec{P_1O}$，找到这个方向上Support点的Minkowski差$P_2$
以$P_1$和$P_2$构成的线段的朝向原点方向的法线为新方向，找到此方向上Support点的Minkowski差$P_3$
此时，$P_1,P_2,P_3$构成一个2D Simplex（也就是三角形）。判断此Simplex是否包含原点。如果包含，则两物体相交
如果不包含原点，找到三角形中离原点最近的边（假设为$P_1,P_2$），以此边朝向原点方向的垂线作为新方向，找到此方向上的Support点的Minkowski差$P_4$
$P_1,P_2,P_4$会构成新的Simplex（这时$P_3$被丢弃）。使用此Simplex回到步骤4
如果构造新Simplex时，新找的Support点已经在原本的Simplex中了，此时说明原点总是在Minkowski差外，那么两物体不相交

可以看到。GJK算法的核心是“贪心地”构造Simplex，直到无法构造。

3D GJK算法一样，先如法炮制找到三个点构成一个三角形。后面的步骤不再找线到原点的垂线，而是用三角面朝向原点那侧的法线作为新方向，从而构造3D的Simplex（四面体）。然后判断原点是否在四面体内。

GJK中较为迷惑的是为什么总要取朝向原点的方向，因为相交要求Minkowski差包含原点，所以他不停地往原点的方向构造Simplex。算是一种贪心算法（其实是迭代下降法）

GJK的优点：

大部分时候不需要遍历两个物体的所有顶点
算法仅要求能够找到两个物体在某方向上的Support点。不需要物体顶点，面法线等乱七八糟信息。严格意义上来说GJK的输入参数是两个Support点获取函数。

GJK的退化

通过GJK的步骤可以看到，他是先构建点，再构建线段，再是三角面再是四面体。那如果其中任意一个构建中，新找到的Support点无法支持构建下一步的Simplex，这要怎么办呢？或者在还没构建出完整Simplex时，新点就已经在旧的Simplex上又要怎么处理呢。这种就是GJK的退化。

退化的处理在Gino的GJK论文中有说。我也有写解析文章。

EPA算法

GJK能够告诉我们物体是否相交。但对于计算MTD来说，我们还想要知道两个物体的挤出深度（相交深度）和挤出方向。

某些引擎中，会使用Simplex到原点的距离以及对应垂线方向作为挤出深度和挤出方向。这是一种近似做法。因为最后的Simplex不一定在Minkowski差所产生的多边形的边上。理论上，只有边上离原点的距离才是最正确的：

GJK退化成线段

这里Simplex是$\Delta_{BGE}$。分离向量是$\vec{OH}$，距离为$\sqrt{2}$。但$GE$边并不在Minkowski差的几何体边上。这意味着这条边并不是正确的。因为相交深度一定是可以将两个几何体刚好分开的长度，而此长度会导致物体的边和边刚好分开。在Minkowski差中就是差的几何体上的边到原点的方向和大小才对。

而EPA算法正是解决此问题的关键。他通过不停扩充Simplex为新的多边形（多面体）来逐步逼近到Minkowski差的几何体边。

其输入时GJK返回的Simplex。以及两个几何体寻找Support点的函数。

首先找到Simplex中离原点最近的边（记为$E$）
找到$E$的外法线（从原点指向边的垂直方向，如图中${OH}$），在这个方向寻找新Support点的Minkowski差（记为$P$）
将$P$加入Simplex扩充Simplex为多面体
重复步骤1~3，直到没有新Support点加入。此时，距离原点最近的边的法线和长度就是挤出深度和相交深度

EPA流程

注意：图中画的扩展Simplex是多边形。实际上是使用多个三角形表示此多边形。

这里比较棘手的地方是，当增加新点时，如何增加对应的新面。

各家物理引擎的做法都差不多，使用面缝合(Suture)方法。假设当找到某个面$S$法线方向的Support点$P$:

删除$S$
找到$S$相邻的所有面$C_i$，然后判断$P$在$C_i$的外侧还是内侧
1. 如果在内侧，用$P$和$S$面留下的边连起来构建三个新的面，扩展结束。
2. 如果在外侧，说明$S_i$也需要被删除，那么做一样的事情：删除$S_i$，判断$P$是否在所有和$S_i$相邻的面的内或外。这样递归地构建直到无面可删。此时将剩下的边和$P$点连接构成新面

那么这里看一下经典的EPA扩展Simplex存储结构：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
struct EPAPolyhedron {
    using uint_type = SmallestUInt<EPA_MAX_FACE_COUNT>;

    struct Face {
        Vector3 m_normal{};
        real m_dist{};
        std::array<uint_type, 3> m_vertices;
        std::array<uint_type, 3> m_adj_facets;  // neighbor facet indices
        std::array<uint_type, 3> m_adj_edges;   // which edge in each neighbor
        bool m_obsolete = false;
    };

    std::array<Vector3, static_cast<size_t>((EPA_MAX_FACE_COUNT + 4) * 0.5)>
        m_pts;
    std::array<Face, EPA_MAX_FACE_COUNT> m_faces;
    uint_type m_vertex_count = 0;
    uint_type m_face_count = 0;
};

这里Face::m_obsolete用于标记面是否被删除。我们在删除面的时候不会移动面在EPAPolyhedron::m_faces里面的位置。只是简单标记为删除。新面也是直接加到m_faces后面。这样可能会更快一些。

而这里面没有任何的内存分配，用EPA_MAX_FACE_COUNT指定了最大的面上限（这里是128（抄的PhysX，Bullet是64，Jolt是动态分配内存的不定大小））。

这个面上限是个经验数字。达到上限之后就直接返回离EPAPolyhedron最近面的方向作为挤出方向就行（近似结果）。所以大部分物理引擎都不支持顶点数特别多的Convex（PhysX是256个顶点256个面）。

GJK和EPA论文解析

GJK和EPA的论文非常值得一看。里面记录的方法并不是现代广为流传的重心坐标法（但数学上是等价的）。

我写了一些论文解析文章，请按顺序阅读：