物理引擎学习03-GJK碰撞检测算法基础

GJK是由Gilbert，Johnson，Keerthi 三位前辈发明的，用来计算两个凸多面体之间的碰撞检测，以及最近距离。GJK算法可以在O(M+N)的时间复杂度内，检测出碰撞，算法在每次迭代的过程中，都会优先选择靠近原点的方向，因此收敛速度会很快。算法的证明过程比较复杂，但是原理还是比较容易理解的。

本文作者游蓝海。原创不易，未经许可，禁止任何形式的转载。本文的写作目的，主要是对GJK算法的理解和应用。对算法本身感兴趣的朋友，可以阅读源论文的文献。本系列GJK算法文章共三篇，本篇是第一篇:

• GJK碰撞检测算法基础
• GJK计算多边形之间的最近距离
• GJK和EPA计算穿透向量

1. 基本原理

1.1 直观理解

通俗的讲，GJK算法就是沿着某个方向，从两个多边形上取相距最远的两个点计算差值；然后将方向取反，再次取两个点计算出差值。如果两个多边形发生碰撞，则这两个差值必然有一个大于0，一个小于0。如果不相交，则两个差值是同为正或同为负。

如下图所示，沿着$HI$方向，得到点$G=A-D_1$；沿着$IH$方向，得到点$F=D-A_1$，两个点分别位于原点两侧。

如下图所示，如果两个多边形不相交，则两个点都在原点左侧

当然这种说法只是方便直观的理解，真正的算法要更严谨一些。

1.2 GJK算法原理

GJK算法的结论是：如果两个多边形相交，那么这两个多边形构成的闵可夫斯基差集(Minkowski Difference)，必然会包含原点。就像1.1节所示那样，差集的点，会分布在原点两侧。只不过这里的差集是一个多边形。

1.3 闵可夫斯基差集(Minkowski Difference)

用多边形A的所有点，减去多边形B中所有的点得到的一个点集合。
A – B = { a – b ∣ a ∈ A , b ∈ B } A – B = \{a – b|a∈A, b∈B\} A–B={a–b∣a∈A,b∈B}
闵可夫斯基差集的意义在于，得到两个多边形顶点间的坐标分布关系，如果两个多边形相交，那么差集中点会分布在原点四周，也就是说差集会包含原点。

差集有一些特殊的性质，差集构成的多边形的形状与两个多边形之间的距离没有直接关系。两个多边形距离越大，则差集的中心位置离原点越远；反之，离原点越近。如果相交，则差集多边形会包含原点。

1.4 单形体(Simplex)

计算闵可夫斯基差集是一个非常麻烦的过程，所幸计算碰撞并不需要得到完整的闵可夫斯基差集多边形，我们仅需要计算出一个能够包含原点的差集多边形即可。对于2D空间，需要得到一个三角形；3D空间需要一个四面体。为了方便表示，我们把这样的差集多边形叫做单形体(Simplex)。

1.5 Support函数

为了方便表示，我们把单形体中的点，称作support点；把得到support点的方法称作support函数。

support函数就像1.1节所述的那样，沿着某个方向，从两个多边形上找出距离最远的两个点，然后计算出差值。

2. GJK算法

虽然GJK算法原理理解起来比较困难，但是实现代码却比较简单。基本上手练习一遍，就可以初步掌握GJK算法。如果接着把GJK计算多边形间的最近距离，和计算穿透向量都掌握之后，就算是彻底掌握了GJK算法。

2.1 GJK算法伪代码

bool GJK(Shape shapeA, Shape shapeB)
{
	// 得到初始的方向
    Vector2 direction = findFirstDirection();
    // 得到首个support点
    simplex.add(support(direction));
    // 得到第二个方向
    direction = -direction;
    while(true)
    {
    	Vector2 p = support(direction);
        // 沿着dir的方向，已经找不到能够跨越原点的support点了。
        if (Vector2.Dot(p, direction) < 0)
            return false;
		
        simplex.add(p);
        // 单形体包含原点了
        if (simplex.contains(Vector2(0, 0)))
            return true;

        direction = findNextDirection();
    }
}

这里比较重要的是迭代如何终止，以及下一次迭代的方向选择，其他概念都比较好理解。下面用文字来解释一下算法核心步骤：

1. 随机选取一个初始方向，用support函数得到第一个support点；
2. 将初始方向取反，作为下一次的迭代方向；
3. 迭代循环开始:
4. 用support函数得到一个新的suppport点；
5. 如果新的support点，在迭代方向上的投影小于0，说明在这个方向上，已经无法找到一个能够跨越原点的support点了。也就是说，无法组成一个能够包含原点的单形体了。则两个多边形不相交，检测到此结束；
6. 如果support点达到3个，用这3点组成三角形，如果包含原点，说明发生了碰撞，检测到此结束；
7. 否则，仅保留离原点最近的support边上的两个support点；
8. 此时，将仅剩的两个support点构成一条直线，计算直线的垂线。并选垂线取朝向原点方向，作为下一次的迭代方向；
9. 跳转到步骤3。

这里比较难理解的是第5步，此时的单形体存在两种情况:

1. 首次进入循环，单形体中只有一个初始的support点。如果投影小于0，说明沿着背离初始点的方向，无法找到一个能够跨越原点的support点了。也就是说，该点和初始点都在原点的同一侧；
2. 非首次进入循环，单形体中只有两个support点了，迭代方向是由步骤8生成的，该方向是垂直于单形体中剩余两个support点构成的直线。如果投影小于0，则说明单形体中仅剩的两点，已经是最接近原点两个support点了。同时这两点构成的线段，就是闵可夫斯基差集中最接近原点的边，该边是计算两个多边形最近距离的关键，下一章中会用到这条边。

需要注意一个特殊情况，步骤8中，如果原点恰好就在两个support点构成的直线上，说明原点就在闵可夫斯基差集的边界上。也就是说，两个多边形刚开始发生碰撞。

计算步骤的分解动图:

2.2 Support函数解析

这里用的support函数，与维基百科上给的support函数有差别，注意不要混淆了。

1. 沿着给定方向，在多边行中找到一个最远的点，也就是在该方向上投影最大的点；
2. 然后，沿着反方向，在另一多边形上也找到一个投影最大的点；
3. 最后，计算两个点的差值，作为support点。

Vector2 support(Vector2 dir)
{
    Vector2 a = getFarthestPointInDirection(shapeA, dir);
    Vector2 b = getFarthestPointInDirection(shapeA, -dir);
    return a - b;
}

Vector2 getFarthestPointInDirection(Shape shape, Vector2 dir)
{
    var vertices = shape.vertices;
    float maxDistance = float.MinValue;
    int maxIndex = 0;
    for(int i = 0; i < vertices.Count; ++i)
    {
        float distance = Vector2.Dot(vertices[i], dir);
        if(distance > maxDistance)
        {
            maxDistance = distance;
            maxIndex = i;
        }
    }
    return vertices[maxIndex];
}

2.3 首次方向选择

首次方向选择比较简单，可以随机选择一个方向；也可以选择两个多边形的中心连接起来的构成的向量；也可以在两个多边形上各取一个点，构成一个向量。但需要注意的是，如果取了两个点求方向，不要重合了，否则方向向量会是0向量。

2.4 后续迭代方向

如果单形体中存在3个support点，则使用最后一个点与前面两个点构成两条边，保留离原点更近的边的两个点，移除另外一个点。
将剩下的两个support点构成一条直线，计算原点到直线的垂线，并取垂线朝向原点的方向，作为下一次的迭代方向。

计算垂线的时候，可以先计算出垂足，然后垂足到原点的向量，就是下一次的迭代方向。但需要注意的是，原点可能就在直线上，则原点和垂足重合，无法计算出方向向量，不过这种情况下，说明已经开始发生碰撞了。

public Vector2 findNextDirection()
{
    if (simplex.count() == 2)
    {
    	// 计算原点到直线01的垂足
        Vector2 crossPoint = getPerpendicularToOrigin(simplex.get(0), simplex.get(1));
        // 取靠近原点方向的向量
        return Vector2.zero - crossPoint;
    }
    else if (simplex.count() == 3)
    {
	    // 计算原点到直线20的垂足
        Vector2 crossOnCA = getPerpendicularToOrigin(simplex.get(2), simplex.get(0));
        // 计算原点到直线21的垂足
        Vector2 crossOnCB = getPerpendicularToOrigin(simplex.get(2), simplex.get(1));
        // 保留距离原点近的，移除较远的那个点
        if (crossOnCA.sqrMagnitude < crossOnCB.sqrMagnitude)
        {
            simplex.remove(1);
            return Vector2.zero - crossOnCA;
        }
        else
        {
            simplex.remove(0);
            return Vector2.zero - crossOnCB;
        }
    }
    else
    {
        // 不应该执行到这里
        return new Vector2(0, 0);
    }
}

2.5 单形体包含检测

如果单形体中有3个support点，则用这三个点构成三角形，然后计算三角形是否包含原点。可以使用分离轴算法进行判断，只要有一条边能够将原点分离在三角形外侧，则说明三角形不包含原点。

2.6 计算闵可夫斯基差集

虽然碰撞检测不需要计算差集，但是为了调试方便，查看差集的形状，可以用最粗暴的方法进行计算。用多边形A的所有点，减去多边形B中所有的点得到的一个点集合，然后使用凸包算法，从差集中得到一个凸多边形。

3. 效果展示

本文Demo使用Unity3D引擎开发，使用了Unity的协程来做分步骤展示。工程已上传github: https://github.com/youlanhai/learn-physics/tree/master/Assets/03-gjk

4. 参考资料

• 维基百科: https://en.wikipedia.org/wiki/Gilbert%E2%80%93Johnson%E2%80%93Keerthi_distance_algorithm
• GJK算法论文: https://ieeexplore.ieee.org/document/2083?arnumber=2083
• dyn4j很详细的GJK算法教程(英文): http://www.dyn4j.org/2010/04/gjk-gilbert-johnson-keerthi/#gjk-iteration
• AndrewFan，对dyn4j教程的中文翻译: https://blog.csdn.net/AndrewFan/article/details/101694644
• Wyman，原始GJK详解: https://www.qiujiawei.com/collision-detection-2

本系列文章会和我的个人公众号同步更新，感兴趣的朋友可以关注下我的公众号：游戏引擎学习。扫下面的二维码加关注：