GJK碰撞检测算法原理与MATLAB实现详解

简介：这是一份基于MATLAB的GJK碰撞检测算法实现包，面向计算机图形学、物理模拟及机器人路径规划等领域的开发者与学习者。GJK算法通过支撑向量与Minkowski差快速判断三维物体是否相交，项目完整实现了支撑向量计算、Minkowski差构造、迭代求解最近距离等核心步骤，并附有示例脚本与形状数据生成代码，可帮助读者快速理解算法原理并应用于实际碰撞检测场景。压缩包共6个文件，包含4个.m函数与演示文件、1个README说明文档和1个.gitignore配置文件，整体仅6KB，结构精简便于阅读和二次开发。已有226人学习使用，适合具备一定MATLAB基础、希望深入掌握GJK算法及碰撞检测实现细节的研究者和工程师。通过分析源码，可以学习到如何将几何算法转化为可运行代码，并了解后续结合EPA精确计算碰撞点的扩展思路。

1. 为什么用GJK做碰撞检测：从Minkowski差到支撑点

做物理模拟或机器人路径规划时，宽阶段用AABB或OBB快速排除远距离物体，盒与盒一旦重叠，就要进入窄阶段精确回答“两个凸体到底碰没碰”。GJK算法是这个阶段最常用的方案之一，它不构造完整的Minkowski差多边形，也不枚举所有边面相交，而是利用支撑函数反复采样，把相交判断收敛到“原点是否落在Minkowski差内”的几何问题。这套MATLAB实现将主循环、形状数据和示例集中在一起，GJK.m负责迭代，SampleShapeData.m提供测试凸体，MAIN_example与兼容版本两个入口让旧版MATLAB也能直接跑。适合想读懂算法骨架、做图形学课程设计，或者准备移植到C++/Python的开发者。下面从最基本的定义开始，逐层拆到可在MATLAB里复现的实现。

2. GJK 的几何基础：Minkowski差、支撑函数与MATLAB形状格式

2.1 Minkowski差与支撑函数：用一组点积代替布尔运算

两个凸体 A 与 B 的 Minkowski差定义为 A ⊖ B = { a - b | a∈A, b∈B }。这个集合有一个漂亮的性质：原点属于 A ⊖ B，当且仅当 A 与 B 相交。所以“检测碰撞”变成了“判断原点是否在某个凸集内部”。直接构造 A ⊖ B 需要枚举所有顶点对，代价是 O(nm)，在实时物理引擎里不可接受。

GJK 用支撑函数绕开显式构造。支撑函数 s_T(d) 返回凸体 T 在方向 d 上最远的顶点，也就是 argmax_{v∈T} d·v。对离散顶点集，这一步只需要一次矩阵乘法：把所有顶点与 d 做点积，取最大值。MATLAB 里写成下面的函数，无论二维还是三维都适用。

function [v, idx] = supportPoint(verts, d) % verts: N x 2 或 N x 3 的顶点矩阵 % d : 搜索方向，列向量或行向量均可 score = verts * d(:); % 所有顶点在 d 上的投影长度 [~, idx] = max(score); v = verts(idx, :); end

代码说明：d(:)强制把方向向量转成列向量，这样verts * d(:)得到 N×1 的点积向量；max返回最大的投影值和对应行号。这里的idx在调试时很有用，可以高亮显示到底是哪个顶点充当了支撑点。需要提醒的是，d不要求是单位向量，但统一做归一化能让后续容差设置更稳定，尤其当两个物体尺度相差很大时。

Minkowski差上的支撑点可以从两个原始凸体直接算出：s_{A⊖B}(d) = s_A(d) - s_B(-d)。原因很简单，沿 d 方向最远的差值点，必然来自 A 上最靠 d 的点，减去 B 上最靠 -d 的点。这就是 GJK 每次只对两个形状各做一次supportPoint调用的基础。

2.2 单纯形与迭代推进：如何从离散点判断原点是否被包围

有了支撑点，GJK 维护一个称为单纯形（simplex）的极小点集。二维时单纯形是线段或三角形，三维时是线段、三角形或四面体。算法每一轮在新的搜索方向 d 上获得一个支撑点 p，然后判断当前单纯形是否已经“包围”原点；如果还没有，就收缩单纯形，并计算出离原点最近的子区域，作为下一轮搜索方向。

这里有一个容易混淆的双重终止条件。第一种是距离式终止：计算原点到单纯形最近点的距离，当距离小于容差时判定相交；第二种是方向式终止：当新支撑点 p 满足 p·d ≤ 0，即 p 没有超过当前单纯形沿 d 方向的最远点，说明 Minkowski差在 d 方向上无法继续前进，此时原点不可能被包含。两种写法在项目代码里都有迹可循，差别如下：

终止方式	判断条件	返回值	适用场景
方向衰减	p·d < tol	布尔碰撞标记	快速碰撞判定，可与 EPA 配合
距离收敛	原点到单纯形距离 < tol	碰撞标记 + 近似穿透距离	需要距离信息的窄阶段

项目里的GJK.m主要走方向式终止，因为它面向“是否碰撞”的布尔输出。若你需要在碰撞后求穿透深度，可以保留最后那个包含原点的单纯形，交给 EPA 继续扩展。

2.3 SampleShapeData 里的凸体表示与凹体预处理

SampleShapeData.m负责返回测试形状，最常见的组织方式是结构体数组，每个元素包含vertices、faces和name三个字段。下面是我在类似项目里常用的字段说明，你可以对照自己的数据做替换：

字段	含义	典型尺寸
vertices	凸体的全部顶点坐标	N×3
faces	三角面或多边形面的顶点索引	M×3
name	形状名称，用于日志输出	char

需要特别注意，GJK 的正确性依赖输入为凸体。凹多面体虽然在支撑函数阶段不会报错，但单纯形迭代可能收敛到错误的区域，导致漏报或误报。项目描述里强调“多边形”和“三维物体”，如果用于凹物体，常见做法是先用vHACD或CoACD做凸分解，再将多个凸包分别做 GJK 判断。另一个细节是顶点顺序对支撑函数没有影响，因为max只关心投影长度；但若之后做 EPA，面的法线方向必须一致，否则穿透方向会反向。

3. GJK.m 主循环实现：支撑点迭代与单纯形更新

3.1 将两个凸体的支撑点合成Minkowski差候选点

在GJK.m里，第一步不是直接判断相交，而是实现 Minkowski差支撑点。通常代码会拆成一个小函数，便于在每次迭代中反复调用：

function p = minkowskiSupport(A, B, d) [sa, ~] = supportPoint(A.vertices, d); [sb, ~] = supportPoint(B.vertices, -d); p = sa - sb; end

sa是 A 沿 d 方向的最远点，sb是 B 沿 -d 方向的最远点，两者相减得到的就是 A ⊖ B 沿 d 方向的支撑点。这个操作的时间复杂度是 O(N + M)，其中 N、M 是两个凸体的顶点数，不随迭代次数增长，所以 GJK 非常适合中等顶点数的凸体。如果输入顶点已经做过缓存，还可以把d变换到局部坐标系，从而复用预计算的凸包数据。

3.2 updateSimplex：线段与三角形的最近点判断

updateSimplex是整个算法里最容易出错的模块。它的输入是当前单纯形点集，输出是收缩后的单纯形和下一轮搜索方向。二维时，单纯形最多三个点，只需要处理两种情况。

function [simplex, d] = updateSimplex(simplex) if size(simplex, 1) == 2 d = edgeDirection(simplex); else [inside, d] = triangleDirection(simplex); if inside d = [0, 0, 0]; % 原点被三角形包围，碰撞成立 end end end function d = edgeDirection(s) A = s(1,:); B = s(2,:); AB = B - A; AO = -A; % AO 指向原点 if dot(AB, AO) > 0 d = cross(cross(AB, AO), AB); % 双叉积得到指向原点的垂线方向 d = d / norm(d); else d = A / norm(A); % 原点在 A 侧，改用顶点方向 end end function [inside, d] = triangleDirection(s) A = s(1,:); B = s(2,:); C = s(3,:); AB = B - A; AC = C - A; AO = -A; ABC = cross(AB, AC); if dot(ABC, AO) < 0 ABC = -ABC; % 保证法线朝向原点 end % 判断原点是否落在三角形内部，分别用边法线做两次叉积 if dot(cross(ABC, AC), AO) < 0 || dot(cross(AB, ABC), AO) < 0 inside = false; d = edgeDirection(s([1 2], :)); else inside = true; d = [0, 0, 0]; end end

逻辑说明：线段方向计算用的是“双叉积”技巧，cross(AB, AO)得到垂直纸面的向量，再与 AB 叉积就得到位于三角形平面上且指向原点的方向。三角形方向判断依赖法线向量ABC，通过对比每个边法线与 AO 的点积，可以判断原点是否在三角形内部。这个版本为了可读性简化了三维 Voronoi region 的完整分支，实际项目里要扩展成四面体的四个面逐一判断；验证时可以先从二维多边形开始，跑通单测再升到三维。

3.3 主循环与两种终止条件

有了支撑点和单纯形更新，主循环可以压缩成十几行。MAIN_example.m里演示的调用方式通常就是把这个循环包装成函数，返回碰撞布尔值和可选距离估计。

function [collision, dist] = GJK(shapeA, shapeB, tol) collision = false; dist = inf; d0 = shapeA.vertices(1,:) - shapeB.vertices(1,:); if norm(d0) < tol collision = true; dist = 0; return; end d = d0 / norm(d0); simplex = minkowskiSupport(shapeA, shapeB, d); d = -simplex(1,:); for iter = 1:50 p = minkowskiSupport(shapeA, shapeB, d); if dot(p, d) < tol * norm(d) dist = dot(simplex(1,:), d); % 保守近似距离 return; end simplex = [simplex; p]; %#ok<AGROW> [simplex, d] = updateSimplex(simplex); if norm(d) < tol collision = true; dist = 0; return; end end end

参数说明：tol控制浮点容差，工程上建议使用相对容差，例如1e-7 * max(shape大小)；50次迭代上限对绝大多数凸体足够，因为每次迭代单纯形包含的维度会递增，三维最多四次就能包围原点，余下是在细化搜索方向。dist的返回值需要小心：当dot(p, d) < tol时，dot(simplex(1,:), d)只是原点到一个支撑点在 d 方向投影的距离，不是严格意义的最短距离。如果需要精确最小距离，要额外跑一次距离求解，或者保留单纯形进入 EPA。

inside-out写法也常见，它与这里的outside-in方向相反：从原点位于某个初始单纯形内部开始，不断扩张寻找支撑点。两种写法最终收敛到同一几何结论，但outside-in实现更直观，调试时更容易用 plot 可视化，项目里的GJK.m更接近前者。

4. 跑通示例：MAIN_example 与 MATLAB 兼容性/数值容差排错

4.1 MAIN_example 与 MAIN_example_compatible 的差异判断

项目根目录给了两个入口文件：MAIN_example.m和MAIN_example_compatible.m。我对比后猜测，前者可能用了较新版本的语法特性，比如arguments块、string数组或隐式扩展；后者则回退成varargin和普通矩阵拼接，以便在 MATLAB R2016b 甚至更老的版本上运行。你在旧版 MATLAB 里如果遇到Function definitions are not supported in this context或Unrecognized function or variable之类的报错，可以直接尝试 compatible 版本。

运行前建议先执行一次clear all; close all;，因为两个示例脚本可能都定义了同名绘图变量。我一般先单步运行MAIN_example_compatible，确认 GJK 返回结果与patch显示相交区域一致。如果脚本里包含disp输出，留意它是打印顶点数还是迭代次数，这能帮助你定位是哪一步抛异常。

4.2 输入数据格式、单位与浮点容差调整

SampleShapeData返回的顶点矩阵可能是 N×3，也可能是 N×2。代码中的supportPoint需要兼容两种维度，d(:)与verts相乘时会自动处理维度。但有一个坑：当顶点数量很少且两个形状靠得很近时，点积最大值可能由多个顶点并列获得，max默认取第一个，导致支撑点在几次迭代中跳变。解决办法是在支撑函数返回后检查一下索引，或者给max结果的顶点坐标添加极小的高斯扰动。

浮点容差不能写死。两个单位立方体的包围盒约 2，而一个小齿轮的包围盒可能只有 0.02，同一tol会造成完全不同的判定行为。可以写一个自适应容差函数，在调用 GJK 之前算好。

function tol = autoTolerance(shapeA, shapeB) s1 = max(max(shapeA.vertices) - min(shapeA.vertices)); s2 = max(max(shapeB.vertices) - min(shapeB.vertices)); tol = 1e-9 * max(s1, s2); end

这段代码用max - min近似每个形状的包围盒边长，再取两者较大值，乘以 1e-9 得到相对容差。如果你遇到两个形状明明相交却返回 false，优先把系数改成 1e-8；如果两个形状相距很近但误报碰撞，改回 1e-10。这个调整比反复改主循环逻辑有效得多。

4.3 迭代过程可视化

GJK 的调试难点在于单纯形在迭代中不断变化，单看打印结果很难判断是支撑函数方向错了，还是单纯形收缩逻辑错了。我通常在updateSimplex入口加一个绘图断点：

plot3(simplex(:,1), simplex(:,2), simplex(:,3), 'o-', 'LineWidth', 1.5); hold on; quiver3(0, 0, 0, d(1), d(2), d(3), 0.2, 'r');

这段代码用plot3画当前单纯形的边，用quiver3从原点画出搜索方向d。如果单纯形像预期一样逐渐向原点收缩，方向箭头会越来越短；如果方向箭头来回震荡，多半是edgeDirection里双叉积的顺序错了，或者顶点顺序导致法线朝外。每次迭代暂停drawnow;可以观察实时动画，方便定位是哪一步开始偏离。

5. 进阶：GJK+EPA 求穿透深度，以及移植到 Python/NumPy

5.1 用 EPA 从单纯形中提取最小穿透向量

GJK 判定碰撞后，后续物理响应需要知道“重叠了多深、从哪个方向推开”。EPA 在 GJK 遗留的单纯形上继续扩展：反复挑选离原点最近的边或面，沿其法线方向取支撑点，如果支撑点接近当前面，就找到了最小穿透向量。MATLAB 里的骨架可以这样写：

while true [faceDist, faceIdx] = min(facesDistance); % 找到最近面 p = minkowskiSupport(shapeA, shapeB, faceNormal(faceIdx,:)); if abs(dot(p, faceNormal(faceIdx,:)) - faceDist) < tol penetration = p; break; end % 将 p 加入单纯形并更新面集合 end

这个循环的终止条件很关键，不要用固定的迭代次数，而要用支撑点投影距离与当前面距离的差是否小于容差。差越小，说明扩展后的面已经贴近真实边界。EPA 对凹体同样不适用，所以它和 GJK 一样，只负责凸包之间的穿透计算。

5.2 移植到 Python/NumPy 的关键改动

把项目代码搬到 Python 时，最常见的问题不是算法本身，而是 MATLAB 的索引从 1 开始且允许用end。核心支撑函数可以直接对照修改：

import numpy as np def support_point(verts, d): d = np.asarray(d).reshape(-1, 1) scores = verts @ d # 等价于 vertices * d(:) idx = np.argmax(scores.ravel()) return verts[idx], idx def minkowski_support(shape_a, shape_b, d): sa, _ = support_point(shape_a.vertices, d) sb, _ = support_point(shape_b.vertices, -d) return sa - sb

这里verts @ d利用了numpy的矩阵乘法，效果等同于 MATLAB 中的verts * d(:)。注意np.argmax返回第一个最大值，与 MATLABmax一致。另一个容易踩的坑是方向向量的形状：d若是一维 array，必须显式 reshape 成列向量，否则verts @ d会得到标量。移植后建议用项目里的SampleShapeData构造相同测试用例，逐轮对比单纯形顶点坐标，而不是只看最终布尔值。

5.3 性能优化与跨语言使用的实用技巧

在实际引擎里，GJK 很少单点调用，而是被物理步进循环反复执行。第一个优化是预计算凸包的局部坐标支撑点：由于物体在局部坐标系下顶点不变，可以提前把顶点矩阵转成只读数组，每次只需用旋转矩阵把方向d变换到局部系，再调用支撑函数。第二个优化是把 GJK 放在宽阶段之后，宽阶段用AABB或Sphere过滤掉绝大多数物体对，这样窄阶段每帧只处理真正可能接触的几对。第三个优化是迭代上限不要设置成固定的 50，而是根据凸包顶点数动态调整，常见凸体迭代典型值在 4 到 12 之间，超过 20 次仍未收敛就要检查输入是否是凹体或顶点包含 NaN。

如果项目后续要接入 Simulink，可以把GJK.m包装成 MATLAB Function 块，输入两个结构体会被自动映射成总线信号，但要注意 Simulink 的代码生成不支持动态增长数组，这时需要把simplex预分配为固定大小的4 x 3矩阵，用额外变量记录当前有效点数。这样改过之后，生成 C 代码时不会产生运行时内存分配，实时性能会稳定不少。

本文还有配套的精品资源，点击获取