GJK and EPA

Deep dive

How it works

In: two convex shapes, each given only by its support function · Out: overlap flag or closest points and distance; with EPA, penetration depth and normal · State: none (cached simplex optional) · Cost: a handful of iterations, two support calls each; EPA O(iterations × faces)

GJK (boolean)

1. s = [support(d0)], d = −s0.
2. Loop: p = support(d). If p·d < 0, the farthest point in direction d does not reach the origin: separated.
3. Append p and reduce the simplex to the feature (vertex, edge, face) whose Voronoi region contains the origin; set d perpendicular to that feature toward the origin. Because p is newest and passed the origin, only regions touching p need testing.
4. If the simplex is a tetrahedron enclosing the origin: overlap.

GJK (distance): same loop, but compute the point of the simplex closest to the origin (barycentric coordinates); stop when |v|² − v·p ≤ ε|v|². That point maps back to closest points on A and B through the barycentric weights.

EPA

1. Start from GJK's tetrahedron, faces oriented outward.
2. Take the face closest to the origin, normal n, distance δ. p = support(n).
3. If p·n − δ < ε, that face lies on the boundary of A − B: depth δ, normal n.
4. Otherwise delete every face that p can see; the edges used by only one deleted face form the horizon loop; connect each horizon edge to p. Repeat.

 A − B (Minkowski difference)        EPA: nearest face → support along its normal
     ___________                       → delete faces visible from p
    /           \                      → stitch horizon edges to p
   /   ·origin   \   overlap ⇔         stop when p lies on that face
   \_____________/   origin inside

Code

use std::ops::{Add, Mul, Neg, Sub};
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct V3(pub f64, pub f64, pub f64);
impl Add for V3 { type Output = V3; fn add(self, o: V3) -> V3 { V3(self.0 + o.0, self.1 + o.1, self.2 + o.2) } }
impl Sub for V3 { type Output = V3; fn sub(self, o: V3) -> V3 { V3(self.0 - o.0, self.1 - o.1, self.2 - o.2) } }
impl Neg for V3 { type Output = V3; fn neg(self) -> V3 { V3(-self.0, -self.1, -self.2) } }
impl Mul<f64> for V3 { type Output = V3; fn mul(self, s: f64) -> V3 { V3(self.0 * s, self.1 * s, self.2 * s) } }
fn dot(a: V3, b: V3) -> f64 { a.0 * b.0 + a.1 * b.1 + a.2 * b.2 }
fn cross(a: V3, b: V3) -> V3 { V3(a.1 * b.2 - a.2 * b.1, a.2 * b.0 - a.0 * b.2, a.0 * b.1 - a.1 * b.0) }

/// A convex shape is fully described by its support function: farthest point along d.
pub trait Convex { fn support(&self, d: V3) -> V3; }

/// Support point of the Minkowski difference A − B.
fn support(a: &dyn Convex, b: &dyn Convex, d: V3) -> V3 { a.support(d) - b.support(-d) }

/// Direction from segment (a → b) toward the origin; any perpendicular if the origin is on it.
fn perp_toward_origin(ab: V3, ao: V3) -> V3 {
    let d = cross(cross(ab, ao), ab);
    if dot(d, d) > 1e-24 { return d; }
    let axis = if ab.0.abs() < 0.57 { V3(1.0, 0.0, 0.0) } else { V3(0.0, 1.0, 0.0) };
    cross(ab, axis)
}

/// Boolean GJK. On overlap returns a tetrahedron of A − B that contains the origin.
pub fn gjk(a: &dyn Convex, b: &dyn Convex) -> Option<Vec<V3>> {
    let mut s = vec![support(a, b, V3(1.0, 0.0, 0.0))];
    let mut d = -s[0];
    for _ in 0..64 {
        if dot(d, d) < 1e-24 { return None; } // origin exactly on a vertex: touching
        let p = support(a, b, d);
        if dot(p, d) < 0.0 { return None; }   // farthest point doesn't pass the origin
        s.push(p);
        if next_simplex(&mut s, &mut d) { return Some(s); }
    }
    None
}

/// Newest point is last. Keep the sub-simplex nearest the origin and aim d at the origin.
/// Triangle invariant: for [s0, s1, s2], (s1 − s2) × (s0 − s2) points toward d.
fn next_simplex(s: &mut Vec<V3>, d: &mut V3) -> bool {
    let a = s[s.len() - 1];
    let ao = -a;
    match s.len() {
        2 => {
            let ab = s[0] - a;
            if dot(ab, ao) > 0.0 { *d = perp_toward_origin(ab, ao); } else { *s = vec![a]; *d = ao; }
            false
        }
        3 => {
            let (b, c) = (s[1], s[0]);
            let (ab, ac) = (b - a, c - a);
            let abc = cross(ab, ac);
            if dot(cross(abc, ac), ao) > 0.0 {
                if dot(ac, ao) > 0.0 { *s = vec![c, a]; *d = perp_toward_origin(ac, ao); false }
                else { *s = vec![b, a]; next_simplex(s, d) }
            } else if dot(cross(ab, abc), ao) > 0.0 {
                *s = vec![b, a]; next_simplex(s, d)
            } else if dot(abc, ao) > 0.0 { *d = abc; false }
            else { *s = vec![b, c, a]; *d = -abc; false } // origin below: flip winding
        }
        _ => {
            let (b, c, e) = (s[2], s[1], s[0]);
            let (ab, ac, ae) = (b - a, c - a, e - a);
            // the three faces through the newest point, with outward normals
            for (tri, n) in [([c, b, a], cross(ab, ac)), ([e, c, a], cross(ac, ae)),
                             ([b, e, a], cross(ae, ab))] {
                if dot(n, ao) > 0.0 { *s = tri.to_vec(); return next_simplex(s, d); }
            }
            true // origin is behind all faces: inside the tetrahedron
        }
    }
}

/// EPA: grow the GJK tetrahedron toward the boundary of A − B.
/// Returns (unit contact normal pointing from A to B, penetration depth).
pub fn epa(a: &dyn Convex, b: &dyn Convex, tet: &[V3]) -> (V3, f64) {
    let mut v: Vec<V3> = tet.to_vec();
    let mut faces: Vec<[usize; 3]> = vec![[0, 1, 2], [0, 3, 1], [0, 2, 3], [1, 3, 2]];
    let normal = |v: &[V3], f: [usize; 3]| {
        let n = cross(v[f[1]] - v[f[0]], v[f[2]] - v[f[0]]);
        n * (1.0 / dot(n, n).sqrt())
    };
    let centroid = (v[0] + v[1] + v[2] + v[3]) * 0.25;
    for f in faces.iter_mut() { // orient outward (not via the origin: it may lie on a face)
        if dot(normal(&v, *f), v[f[0]] - centroid) < 0.0 { f.swap(1, 2); }
    }
    let (mut best_n, mut best_d) = (V3(0.0, 0.0, 0.0), f64::INFINITY);
    for _ in 0..64 { // curved shapes converge slowly: see Pitfalls
        best_d = f64::INFINITY;
        for &f in &faces {
            let n = normal(&v, f);
            let dist = dot(n, v[f[0]]); // distance of the face plane from the origin
            if dist < best_d { best_d = dist; best_n = n; }
        }
        let p = support(a, b, best_n);
        if dot(p, best_n) - best_d < 1e-6 { break; } // face lies on the boundary: done
        // Remove faces p can see; edges used by exactly one removed face form the horizon.
        let mut horizon: Vec<(usize, usize)> = Vec::new();
        faces.retain(|&f| {
            if dot(normal(&v, f), p - v[f[0]]) <= 1e-10 { return true; } // eps: no coplanar flips
            for (i, j) in [(f[0], f[1]), (f[1], f[2]), (f[2], f[0])] {
                match horizon.iter().position(|&e| e == (j, i)) {
                    Some(k) => { horizon.swap_remove(k); } // shared with another removed face
                    None => horizon.push((i, j)),
                }
            }
            false
        });
        v.push(p);
        let pi = v.len() - 1;
        faces.extend(horizon.into_iter().map(|(i, j)| [i, j, pi])); // keeps outward winding
    }
    (best_n, best_d)
}

A sphere of radius 1 at the origin and one of radius 1 at (1.5, 0, 0) yield normal ≈ (1, 0, 0) and depth ≈ 0.5. Boxes are exact: a unit-half-size box and a 0.5-half-size box at (0.3, 1.3, 0.2) give normal (0, 1, 0), depth 0.2.

Complexity and constants

• Support function cost: O(1) for spheres, boxes, capsules; O(n) brute force for an n-vertex hull, much less with hill climbing over vertex adjacency. Cache the last support vertex per shape for warm-started hill climbing.
• GJK: iterations are bounded for polytopes (no simplex repeats); for typical game shapes a handful.
• EPA: each iteration adds one vertex; the face scan here is O(F) per iteration, F = 2V − 4 for a closed triangle mesh. A heap of faces keyed by distance makes it O(log F).
• Curved shapes: in tests of the code above on random sphere pairs, 64 iterations left depth errors up to order of 0.07 units for radii ~1, while 400 iterations brought them down to order of 10⁻⁵.

Efficiency

• Cost: a handful of iterations, two support calls each; bound: the support function, order of 100 ns per pair on boxes (not measured here).
• Levers:
  • warm start from last step's simplex — 1-2 iterations instead of a handful — coherent motion
  • accelerated GJK — 5x to 15x faster than Bullet, FCL and libccd in the Coal authors' benchmarks (Coal README; Montaut 2022) — robotics queries
• Naive → standard → state of the art:
  • explicit Minkowski difference · n·m points, then a hull
  • GJK with warm start and margins, EPA when deep · a handful of support calls
  • accelerated GJK, signed-volume sub-algorithm (Montanari 2017) · 5x to 15x vs older libraries

Implementation notes

• Margins: shrink each shape by a radius (Bullet's collision margin, Box2D's polygon radius); GJK on the cores stays separated for shallow contacts, depth = radii − core distance.
• Termination guards: cap iterations (64 is common), stop when the new point duplicates a simplex vertex, and use a relative tolerance for large-scale worlds.
• Signed-volume subalgorithm (Montanari et al.) replaces Johnson's distance subalgorithm for the distance variant; more robust for degenerate simplices.
• Work in A's local frame and transform B's support direction by the relative rotation; avoids transforming every vertex.
• Return the witness points: store support_A and support_B with each simplex/polytope vertex; the barycentric weights of the closest point on the final face give contact points on each body.
• f32 vs f64: f32 works with margins and tolerances near 1e-4 relative; the EPA horizon test needs an epsilon either way.

Pitfalls

• Origin exactly on a segment or triangle: the triple product (ab × ao) × ab is zero and naive GJK reports "no overlap" for a deep hit.
• Orienting EPA's initial faces by the origin: if the origin lies on a face (touching or symmetric configurations) the orientation is a coin flip. Orient by the tetrahedron centroid.
• No epsilon in EPA's visibility test: a nearly coplanar support point breaks the horizon, the polytope turns non-convex and the depth is wrong.
• Flat GJK tetrahedron handed to EPA (shapes merely touch): check its volume, report depth 0.
• Support functions that return interior points (e.g. a box support with a sign bug): GJK silently returns wrong answers. Unit-test supports alone.

Variants and alternatives

In production

Games, VR and effects: Bullet (btGjkPairDetector, btGjkEpa2), Jolt Physics (GJKClosestPoint, EPAPenetrationDepth), Box2D (b2Distance), ODE/libccd. Robotics: MuJoCo's default convex collision pipeline is GJK/EPA (docs); Coal has its own GJK/EPA and serves Pinocchio and the Humanoid Path Planner (README).

Implementations

Sources

• Elmer G. Gilbert, Daniel W. Johnson, S. Sathiya Keerthi, A Fast Procedure for Computing the Distance Between Complex Objects in Three-Dimensional Space, IEEE Journal on Robotics and Automation, 1988. paper
• Gino van den Bergen, Proximity Queries and Penetration Depth Computation on 3D Game Objects, GDC 2001.
• Gino van den Bergen, Collision Detection in Interactive 3D Environments, Morgan Kaufmann, 2003.
• Christer Ericson, Real-Time Collision Detection, Morgan Kaufmann, 2005 (ch. 9).
• Mattia Montanari, Nik Petrinic, Ettore Barbieri, Improving the GJK Algorithm for Faster and More Reliable Distance Queries Between Convex Objects, ACM Transactions on Graphics, 2017. paper
• Casey Muratori, Implementing GJK, video, 2006.
• Louis Montaut, Quentin Le Lidec, Vladimir Petrik, Josef Sivic, Justin Carpentier, Collision Detection Accelerated: An Optimization Perspective, RSS 2022. paper · code