碰撞检测

发表于 更新于 本文字数： 4.6k 阅读时长 ≈ 17 分钟

粗筛阶段（Broad Phase）

SAP（扫描线算法）

不论碰撞体本身是圆形、OBB、多边形，都可以获得顶点集合的minPoint: (minX, minZ)和maxPoint: (maxX, maxZ)，得到最大外接矩形也就是AABB，将所有待检测对象的AABB投影到对应的轴上。每个对象在每个轴上都会得到一个线段，两个不同对象的线段如果在某一个轴上没有交集，就证明一定不存在碰撞，反之，在所有轴上都存在交集，则这两个对象的AABB一定发生了碰撞，如果对象本身就是AABB那就检测完毕，否则需要进一步判断两个对象是否实际碰撞。

算法

先从X轴开始，遍历每个对象，投影可以产生两个端点minValue, maxValue，将所有对象的投影端点都放入一个数组内，并按坐标从小到大进行排序。维护两个列表：活跃列表[]和重叠对列表[]。扫描线从左到右扫描这些端点。

  • 当遇到对象A的minValue时，将其加入活跃列表[A]，表示它的区间开始了。
  • 当遇到对象B的minValue时，将其加入活跃列表[A, B]，在当前轴，B与活跃列表的所有对象发生了重叠，因此将所有重叠对加入重叠对列表[(A, B)]。
  • 当遇到对象C的minValue时，将其加入活跃列表[A, B, C]，在当前轴，C与活跃列表的所有对象发生了重叠，因此将所有重叠对加入重叠对列表[(A, B), (A, C), (B, C)]。
  • 当遇到A的maxValue时，将A移除活跃列表[B, C]
  • 当遇到B的maxValue时，将B移除活跃列表[C]
  • 当遇到C的maxValue时，将C移除活跃列表[]

遍历上一步得到的X轴重叠对列表，按照先前相同的思路进行Z轴的重叠判定，可以得到XZ轴重叠列表，如果是二维碰撞，此时已经完成了SAP判断；否则遍历XZ轴重叠列表，进行Y轴的重叠判定，最终得到XYZ轴重叠列表。

如果物体本就是AABB，那最终的XYZ轴重叠列表就是碰撞对，否则需要对这些碰撞对进行更精确的碰撞检测。

时空相关性

游戏里，物体每一帧之间的位置变化通常非常小，意味着上一帧排好的顺序，在这一帧几乎还是对的，只需要用插入排序即可以接近O(N)的复杂度重新恢复顺序，只有物体“互相超车”端点顺序发生交换、新物体加入、旧物体删除时才会进行调整

补充：对于逆序对较少的列表，插入排序的速度极快，这也是希尔排序的中心思想

策略

之前的例子是从X轴开始，这肯定不是必须的，最好的方法是从物体分布最稀疏的轴开始

精细阶段（Narrow Phase）

下面介绍的两个算法都是基于凸多边形的，如果是凹多边形，需要将其划分成多个凸多边形来进行处理（所有凹面多边形都可以拆成凸面多边形）。

  • 凸多边形：多边形任意两点的连线都在多边形内部，凹多边形总能找到反例

1. SAT

从分离的角度判断物体间的碰撞

SAT非常简单，核心思路就几点：

  • 取两个多边形所有的边
  • 对每条边求一个垂直轴（法线、分离轴）
    • 把两个多边形的每个边依次投影到该轴上，求得投影的$(min, max)$区间
  • 如果存在某个轴，求得的两个区间不重叠，代表不碰撞；反之，不存在这样的轴，则表明碰撞
  • 且可以通过重叠区域在对应轴的投影长度，求得最小平移向量（MTV, Minimum Translation Vector）
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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
public static class Sat
{
    public static bool Intersects(ConvexPolygon polyA, ConvexPolygon polyB, out Vector2 mtv)
    {
        mtv = default;
        var minSeparateDist = float.MaxValue;
        Vector2 smallestAxis = default;

        if (!TestPolygonAxes(polyA, polyB, ref minSeparateDist, ref smallestAxis))
            return false;

        if (!TestPolygonAxes(polyB, polyA, ref minSeparateDist, ref smallestAxis))
            return false;
        
        var centerA = ComputeCenter(polyA);
        var centerB = ComputeCenter(polyB);
        // 方向修正：保证从A指向B
        if (Vector2.Dot(centerB-centerA, smallestAxis) < 0) smallestAxis *= -1f;
        mtv = smallestAxis * minSeparateDist;
        return true;
    }

    private static bool TestPolygonAxes(ConvexPolygon polyA, ConvexPolygon polyB, ref float minSeparateDist, ref Vector2 smallestAxis)
    {
        var count = polyA.Vertices.Length;
        
        
        for (int i = 0; i < count; i++)
        {
            var p1 = polyA.Vertices[i];
            var p2 = polyA.Vertices[(i + 1) % count];

            var edge = p2 - p1;
            var axis = new Vector2(-edge.y, edge.x).normalized; // 2D垂直向量

            Project(polyA, axis, out var minA, out var maxA);
            Project(polyB, axis, out var minB, out var maxB);

            if (maxA < minB || maxB < minA) return false; // 找到分离轴

            var separateDist = Mathf.Min(maxA - minB, maxB - minA); // 脱离距离：在当前轴完成脱离的最小距离
            // var overlap = Mathf.Min(maxA, maxB) - Mathf.Max(minA, minB); // overlap在包含的情况下失去意义
            if (separateDist >= minSeparateDist) continue;
            
            minSeparateDist = separateDist;
            smallestAxis = axis;
        }
        return true;
    }

    private static void Project(ConvexPolygon poly, Vector2 axis, out float min, out float max)
    {
        var projection = Vector2.Dot(poly.Vertices[0], axis);
        min = projection;
        max = projection;

        for (var i = 1; i < poly.Vertices.Length; i++)
        {
            projection = Vector2.Dot(poly.Vertices[i], axis);

            if (projection < min) min = projection;
            if (projection > max) max = projection;
        }
    }

    private static Vector2 ComputeCenter(ConvexPolygon poly)
    {
        var x = 0f;
        var y = 0f;
        var length = poly.Vertices.Length;
        for (var i = 0; i < poly.Vertices.Length; i++)
        {
            x += poly.Vertices[i].x;
            y += poly.Vertices[i].y;
        }
        return new Vector2(x / length, y / length);
    }
}

代码里不需要MTV可以去掉，同时也可以去掉向量的标准化，优化效率。

1.1 OBB优化

OBB是平行四边形，每个图形只需要取两个边的轴进行判断，节省一半运算。此外根据OBB的实现还可以进一步优化：

  • OBB包含顺时针或逆时针顺序的顶点相对坐标：可以取连续两个顶点，进行投影取其中最大值，这就是该OBB的“半径”。然后计算两个OBB的中心距离，两个半径减去中心距离即可得出沿当前轴的分离距离，如果不相交，则分离距离小于0
  • OBB包含分离的单位向量axisX、axisY和两个轴的伸展extend：我们可以分别计算OBB的两个朝向的投影长度，乘以对应轴的extend，然后将结果相加，此时得到的值也是“半径”，后续操作和上述一致
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
46
47
48
49
50
51
52
53
public class SatObb
{
    // mtv: A推开B（B远离A）
    public static bool Intersects(Obb obbA, Obb obbB, out Vector2 mtv)
    {
        mtv = default;
        var minSeparateDist = float.MaxValue;
        Vector2 smallestAxis = default;

        if (!TestObbAxes(obbA, obbB, ref minSeparateDist, ref smallestAxis))
            return false;

        if (!TestObbAxes(obbB, obbA, ref minSeparateDist, ref smallestAxis))
            return false;
    
        // 方向修正：保证从A指向B
        if (Vector3.Dot(obbB.Center - obbA.Center, smallestAxis) < 0) smallestAxis *= -1f;
        mtv = smallestAxis * minSeparateDist;
        return true;
    }

    private static bool TestObbAxes(Obb obbA, Obb obbB, ref float minSeparateDist, ref Vector2 smallestAxis)
    {
        var separateDist = CalcSeparateDistInAxis(obbA, obbB, obbA.AxisX);
        if (separateDist <= 0) return false;
        if (separateDist < minSeparateDist)
        {
            minSeparateDist = separateDist;
            smallestAxis = obbA.AxisX;
        }
        separateDist = CalcSeparateDistInAxis(obbA, obbB, obbA.AxisY);
        if (separateDist <= 0) return false;
        if (separateDist < minSeparateDist)
        {
            minSeparateDist = separateDist;
            smallestAxis = obbA.AxisY;
        }
        return true;
    }

    private static float CalcSeparateDistInAxis(Obb obbA, Obb obbB, Vector2 axis)
    {
        var dist = MathF.Abs(Vector2.Dot(obbA.Center - obbB.Center, axis));
        var aXProj = MathF.Abs(Vector2.Dot(obbA.AxisX , axis)) * obbA.Extend.x;
        var aYProj = MathF.Abs(Vector2.Dot(obbA.AxisY , axis)) * obbA.Extend.y;
        var aRadius = aXProj + aYProj;
    
        var bXProj = MathF.Abs(Vector2.Dot(obbB.AxisX, axis)) * obbB.Extend.x;
        var bYProj = MathF.Abs(Vector2.Dot(obbB.AxisY, axis)) * obbB.Extend.y;
        var bRadius = bXProj + bYProj;
        return aRadius + bRadius - dist;
    }
}

2. GJK

从重叠的角度判断物体间的碰撞
现在我们要判断两个图形相交，我们不再考虑图形的边界线，而是将图形视为一组无限点集。这种情况下，如何判断两个图形相交呢？本能的，只要确定了两个点集有一个共同点就可以判定他们相交：$(A_{x}, A_{y})=(B_{x}, B_{y})$，我们将他们挪到等号一侧：$(A_{x}, A_{y})-(B_{x}, B_{y})=0$，而在数学上，存在两个点集的运算：闵可夫斯基和/差。

2.1 闵可夫斯基和/差

两个图形的每个点都当成从原点出发的向量，然后两两相加，得到新的一堆坐标点。换句话说，就像把一个图形在另一个图形边上扫描一圈，形成的一个膨胀区域。下面就是闵可夫斯基和的数学公式：
$$
A⊕B={{a+b|a∈A,b∈B}}
$$
而在刚才推导的等式形式上则是闵可夫斯基差，说白了，其实就是其中一个图形与另一个图形取反（原点中心对称）后的闵可夫斯基和：
$$
A⊖B={{a+(-b)|a∈A,b∈B}}
$$
GJK中两条关键特性：

  • 两个凸多边形的闵可夫斯基和/差，结果依旧是凸多边形
    • A和B是凸集 -> $A⊖B convex$
  • 如果两个形状相交，那么他们的闵可夫斯基差集合一定包含原点（意味着他们至少有一个共同点）
    • A和B相交 -> $(0,0)∈A⊖B$

至此，判断两个图形相交被转化为判断两个图形的闵可夫斯基差是否包含原点。但问题来了，这样无限点集的定义，如何用计算机来计算呢？实际上，我们压根不需要去求得完整的闵可夫斯基差集，我们只需要判断所得出的新图形包含原点就行了。进一步，由于图形是凸的，他们的闵可夫斯基差也是凸的，内部任意三个点形成的三角形区域必然在差集内，只要存在某个构建的三角形包含了原点，那就证明了两个图形相交。

在数学上，这个“三角形”被称为“单纯形”，指代能推广到任意维度的广义三角形概念，简单来说，单纯形就是能包住某个维度点的最基础形状，所以二维空间的单纯形就是三角形，三维空间的单纯形就是四面体，这块不是重点，在本文，理解为三角形就行了。

现在问题就是如何从无限点集内，合理快速的找到这个包围原点的三角形呢？在考虑该问题前，需要介绍凸集的另一个特性：形状上的每个顶点，都存在一个方向，使它成为凸多边形最远的点。换句话说，只要我挨个检查所有方向，找出凸多边形在该方向最远的点，就能得到图形上所有的顶点。

现在，咱们可以把某个方向，对应到凸多边形上的顶点了：$S_A(\vec{d})$->凸多边形某个顶点。我们将“输入某个方向向量，返回最远顶点”的函数称为支撑函数，对应的那个顶点就叫做支撑点。用原始形状的支撑点，就能算出闵可夫斯基和与差的支撑点：
$$
S_A(\vec{d}) + S_B(\vec{d}) = S_C(\vec{d})
$$
计算闵可夫斯基差很简单：随便选一个方向$\vec{d}$，找到第一个图形在$\vec{d}$的支撑点，然后找到第二个图形在反方向$-\vec{d}$的支撑点，然后将这两个点相减，即可得到最终闵可夫斯基差集合中$\vec{d}$方向最远的顶点。简单来说：我们其实是在两个图形上分别用正反方向的最远点相减，来找出闵可夫斯基差集的边界点。
$$
S_A(\vec{d}) - S_B(-\vec{d}) = S_C(\vec{d})
$$
现在我们不用完整计算闵可夫斯基差集就能直接找到差集图形边界上的顶点了！这样一来，我们就可以按方向来选点，等会儿要在闵可夫斯基差集边界上找包围原点的三角形时，这个技巧就可以派上大用场。

2.2 获取支撑点

支撑函数输入某个方向向量$\vec{d}$，返回凸集$B$中最远的顶点$v$，数学上就是遍历顶点，找到某个顶点，它与$\vec{d}$的点积最大。$\vec{d}$是单位向量，点积的结果就是投影长度，也就是该凸集在该方向的“最远点”。
$$
S_B(\vec{d}) = v = \underset{v \in B}{\argmax} \ v^T\vec{d}
$$

  • 凸多边形：只需要计算每个顶点与$\vec{d}$的点积，找到最大值的那一个顶点
  • 圆形：直接用方向向量乘以半径加上圆心偏移就能得到支撑点$S_B(\vec{d}) = C + r\vec{d}$
  • 其他几何图形都可以自定义计算支撑点的方法，因此GJK压根不关心具体的形状（要凸）

2.3 具体步骤

具体步骤参考文献内Youtuber Reducible的视频讲解的非常好。思想就是朝向原点方向迭代构建三角形尝试去包围原点，涉及到操作有：点积判断是否越过原点、向量三重积求取线段的有向垂线、根据构建得出的三角形沃罗诺伊区域（Voronoi）判断点在三角形内以及后续的操作。视频结合动画实在效果一流，远非博客文字和图片所能代替，因此这里不再赘述。

2.4 代码实现

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
public static class Gjk
{
    private static readonly List<Vector2> Simplex = new();
    
    public static bool Intersects(ConvexPolygon a, ConvexPolygon b)
    {
        Simplex.Clear();
        var d = Vector2.right;

        Simplex.Add(Support(a, b, d));
        d = -Simplex[0];

        while (true)
        {
            var p = Support(a, b, d);

            if (Vector2.Dot(p, d) <= 0)
                return false; // 不相交

            Simplex.Add(p);

            if (HandleSimplex(Simplex, ref d))
                return true; // 相交
        }
    }

    public static Vector2 Support(ConvexPolygon a, ConvexPolygon b, Vector2 d)
    {
        return a.Support(d) - b.Support(-d);
    }

    // 2D 只有两种情况：
    // 当前 simplex 是线段
    // 当前 simplex 是三角形
    public static bool HandleSimplex(List<Vector2> simplex, ref Vector2 d)
    {
        if (simplex.Count == 2)
        {
            // simplex: [B, A]（注意：最后加入的是 A）
            var a = simplex[1];
            var b = simplex[0];

            var ab = b - a;
            var ao = -a;

            // AB 垂直方向，指向原点那边
            d = TripleProduct(ab, ao, ab);
            
            // 防止出现共线情况（两个正方形紧邻，交集是一条线）
            if (d.sqrMagnitude < 1e-6f)
            {
                d = new Vector2(-ab.y, ab.x);
                // 这里主要是防止浮点数误差，对于数学上严格共线的情况下，点积必然为0，法向量方向不用在意
                if (Vector2.Dot(ao, d) < 0) d = -d;
            }

            return false;
        }
        else
        {
            // simplex: [C, B, A]（注意：最后加入的是 A）
            var a = simplex[2];
            var b = simplex[1];
            var c = simplex[0];

            var ab = b - a;
            var ac = c - a;
            var ao = -a;

            var abPerp = TripleProduct(ac, ab, ab);
            var acPerp = TripleProduct(ab, ac, ac);

            // 原点在 AB 外侧
            if (Vector2.Dot(abPerp, ao) > 0)
            {
                simplex.RemoveAt(0); // 移除 C
                d = abPerp;
                return false;
            }

            // 原点在 AC 外侧
            if (Vector2.Dot(acPerp, ao) > 0)
            {
                simplex.RemoveAt(1); // 移除 B
                d = acPerp;
                return false;
            }

            // 原点在三角形内部
            return true;
        }
    }

    // 向量三重积公式来求取向量A向原点的法向量
    // aXbXc结果是c的垂线（与(b-a)点积大于0的方向）
    public static Vector2 TripleProduct(Vector2 a, Vector2 b, Vector2 c)
    {
        var ac = Vector2.Dot(a, c);
        var bc = Vector2.Dot(b, c);
        return b * ac - a * bc;
    }
}

在求解点是否在三角形内时没有使用传统的循环变叉积（inside test）同向方法。因为我们没有“完整三角形”，这些点也不是固定顺序，不知道逆时针/顺时针，simplex也是动态变化的。此外不仅仅只判断点在三角形内，还需要知道在哪个边外，从而丢掉离远点最远的点并且更新方向。实际上由一个三角形，可以划分七个区域：沃罗诺伊区域（Voronoi），而原点只会存在于3个区域。

此外，求某个线的垂线为什么用三重积，因为可以推广到三维，此外三重积自动确认方向，如果写$(-y,x)$只能确保垂直，无法确认方向，还需要额外做点积。

2.5 EPA

全称Expanding Polytope Algorithm，就是用来在GJK碰撞确认后，求穿透深度和穿透向量的算法。理解完GJK后，该算法非常简单。思路就是：

  • 在GJK求取的包含原点的单纯形内，继续扩展。
  • 不断地去找到离原点最近的边然后向原点的反方向作垂向量，然后得到新的闵可夫斯基差顶点，将它加入定点列表
  • 如果下一个多边形顶点加入后，其正好是该方向的最远点，就代表找到了
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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
public static class Epa
{
    private const float Epsilon = 1e-3f; // 太小，float会误判
    private const int MaxIterations = 32;
    private static readonly List<Vector2> Simplex = new();

    public static bool Intersects(ConvexPolygon a, ConvexPolygon b, out Vector2 mtv)
    {
        mtv = default;
        Simplex.Clear();

        if (!GjkInternal(a, b, Simplex)) return false;

        // 进入 EPA
        mtv = EPA(a, b, Simplex);
        return true;
    }

    #region GJK
    private static bool GjkInternal(ConvexPolygon a, ConvexPolygon b, List<Vector2> simplex)
    {
        var d = Vector2.right;

        simplex.Add(Gjk.Support(a, b, d));
        d = -simplex[0];

        while (true)
        {
            var p = Gjk.Support(a, b, d);

            if (Vector2.Dot(p, d) <= 0)
                return false;

            simplex.Add(p);

            if (Gjk.HandleSimplex(simplex, ref d))
                return true;
        }
    }
    #endregion

    #region EPA
    private static Vector2 EPA(ConvexPolygon A, ConvexPolygon B, List<Vector2> simplex)
    {
        // 包含三个顶点
        var polytope = new List<Vector2>(simplex);
        for (var iter = 0; iter < MaxIterations; iter++)
        {
            var closestEdgeIndex = -1;
            var minDist = float.MaxValue;
            var closestNormal = Vector2.zero;

            // 找离原点最近的边
            for (var i = 0; i < polytope.Count; i++)
            {
                var j = (i + 1) % polytope.Count;
                var a = polytope[i];
                var b = polytope[j];
                var edge = b - a;
                var toOrigin = -a;

                var normal = Gjk.TripleProduct(edge, toOrigin, edge).normalized; // 指向原点
                if (normal.sqrMagnitude < Epsilon)
                {
                    normal = new Vector2(a.y, -a.x);
                    if (Vector2.Dot(normal, toOrigin) < 0) normal = -normal;
                    normal = normal.normalized;
                }
                
                var dist = Mathf.Abs(Vector2.Dot(normal, a));

                if (dist >= minDist) continue;
                minDist = dist;
                closestEdgeIndex = j;
                closestNormal = -normal; // 背向原点
            }

            // 沿背向原点的法线做 support
            var p = Gjk.Support(A, B, closestNormal);
            var d = Vector2.Dot(p, closestNormal);

            if (d - minDist < Epsilon) return closestNormal * d;

            // 插入新点
            polytope.Insert(closestEdgeIndex, p);
        }
        return Vector2.zero;
    }
    #endregion
}

EPA实现起来坑很多：

  • 浮点数误差，Epsilon不要搞太小
  • 初始的单纯形的边可能会非常接近甚至经过原点，此时求出的法线是接近零向量，用$(-y,x)$坐标直接获取法向量并做点积，尝试拿到原点反向的法向量，方向不对就换成$(y,-x)$。但是对于经过原点的边，就只能默认$(-y,x)$，后面收集的顶点列表里可能会重复，甚至就不是一个纯粹的图形了，但是最终结果确实对的。（比如两个2x2正方形，都在原点，此时就会遇到这个情况）
  • 里面技巧的一点是polytope.Insert(closestEdgeIndex, p)，正常情况插入的新顶点是在本轮循环找到的最近边的两个顶点之间的，这样依旧能保证顶点是顺序的（虽然上一条发现，不是顺序的也能求）

profiler

在进行OBB检测，每帧1k次

算法 时间
SAT-OBB 10ms
SAT 20ms
GJK 15ms
EPA 50ms

参考文献