1. Convex Polytopes and Support Functions

A convex polytope is the intersection of a finite set of half-spaces — equivalently, the convex hull of a finite point set. Every convex-vs-convex collision algorithm in this tutorial (SAT, GJK, EPA) is built on one primitive: the support function. Given a direction d, support(d) returns the vertex of the shape that is farthest along d:

// Support function for a convex hull stored as a vertex array
function support(vertices, d) {
  let best = vertices[0];
  let bestDot = best.x * d.x + best.y * d.y + best.z * d.z;
  for (let i = 1; i < vertices.length; i++) {
    const v = vertices[i];
    const dot = v.x * d.x + v.y * d.y + v.z * d.z;
    if (dot > bestDot) { bestDot = dot; best = v; }
  }
  return best;
}

// Support of a Minkowski difference A ⊖ B: farthest point of A
// along d, minus farthest point of B along -d
function supportMinkowski(hullA, hullB, d) {
  const pA = support(hullA, d);
  const negD = { x: -d.x, y: -d.y, z: -d.z };
  const pB = support(hullB, negD);
  return {
    x: pA.x - pB.x,
    y: pA.y - pB.y,
    z: pA.z - pB.z,
    a: pA, b: pB // keep origin points for contact reconstruction
  };
}

The key insight both SAT and GJK exploit: two convex shapes A and B overlap if and only if the origin lies inside the Minkowski difference A ⊖ B = { a − b : a ∈ A, b ∈ B }. SAT tests this indirectly via axis projections; GJK tests it directly by building a simplex inside the difference set.

2. SAT: Theory and Axis Selection

The Separating Axis Theorem states that two convex shapes do not intersect if and only if there exists an axis onto which their projections do not overlap. For polytopes, it is enough to test a finite set of candidate axes:

Face normals of A — catches face-vertex and face-edge contact where A's face is the separator.
Face normals of B — same, for B's faces.
Cross products of edge pairs (A × B) — in 3D, edge-edge contact (e.g. two boxes touching along their edges) is not caught by face normals alone and requires testing normalize(edgeA × edgeB) for every pair of edges.

Projection interval of a point set S onto axis n̂: [min, max] = [ min_{p∈S}(p·n̂), max_{p∈S}(p·n̂) ] Overlap test for intervals [minA,maxA], [minB,maxB]: overlap = min(maxA, maxB) − max(minA, minB) overlap ≤ 0 ⟹ n̂ is a separating axis ⟹ no collision overlap > 0 for ALL candidate axes ⟹ shapes intersect Minimum Translation Vector (MTV): the axis with the SMALLEST positive overlap gives both the contact normal and the penetration depth.

A box has 3 unique face normals (3 pairs of parallel faces), so two OBBs need 3 + 3 = 6 face-normal axes plus 3×3 = 9 edge-cross axes: 15 axes total. This is the classic result behind OBB-OBB SAT used in countless physics engines and the original Nintendo64-era collision code.

3. SAT for Oriented Boxes (OBB vs OBB)

Here is a complete 3D OBB-vs-OBB SAT test, including the 9 edge-cross axes required for edge-edge contacts:

// box = { center, halfExtents: {x,y,z}, axes: [ex, ey, ez] } (unit vectors)
function obbSAT(A, B) {
  const axes = [];
  for (const a of A.axes) axes.push(a);
  for (const b of B.axes) axes.push(b);
  for (const a of A.axes)
    for (const b of B.axes) {
      const c = cross(a, b);
      if (length(c) > 1e-6) axes.push(normalize(c)); // skip parallel edges
    }

  let minOverlap = Infinity, mtvAxis = null;
  const d = sub(B.center, A.center);

  for (const axis of axes) {
    const rA = projectedRadius(A, axis);
    const rB = projectedRadius(B, axis);
    const dist = Math.abs(dot(d, axis));
    const overlap = rA + rB - dist;
    if (overlap <= 0) return null; // separating axis → no collision
    if (overlap < minOverlap) { minOverlap = overlap; mtvAxis = axis; }
  }
  // normal must point from A to B
  if (dot(d, mtvAxis) < 0) mtvAxis = negate(mtvAxis);
  return { normal: mtvAxis, depth: minOverlap };
}

// Projected "radius" of an OBB onto an arbitrary axis
function projectedRadius(box, axis) {
  return box.halfExtents.x * Math.abs(dot(box.axes[0], axis))
       + box.halfExtents.y * Math.abs(dot(box.axes[1], axis))
       + box.halfExtents.z * Math.abs(dot(box.axes[2], axis));
}

Why overlap-of-radii, not projected vertices: for a box, the projection of all 8 vertices onto an axis is always centered at the projected box center with radius Σ hᵢ·|axisᵢ·n̂|. Computing it this way avoids projecting all 8 vertices explicitly — an O(1) formula instead of O(vertices).

4. SAT for Arbitrary Convex Hulls

For general convex polyhedra (not just boxes), face normals come from each face and edge-cross axes come from every pair of non-parallel edges — this is O(F₁ + F₂ + E₁·E₂), which is fine for typical collision meshes (tens of faces) but becomes expensive for high-poly hulls. Most engines therefore reduce collision geometry to simplified convex hulls (≤ 32 vertices) generated offline with quickhull, and reserve GJK/EPA — which don't enumerate axes at all — for anything with many vertices.

function hullSAT(hullA, hullB) {
  const axes = [
    ...hullA.faceNormals,
    ...hullB.faceNormals,
    ...edgeCrossAxes(hullA.edges, hullB.edges)
  ];
  let best = { depth: Infinity, normal: null };
  for (const n of axes) {
    const [minA, maxA] = projectHull(hullA.vertices, n);
    const [minB, maxB] = projectHull(hullB.vertices, n);
    const overlap = Math.min(maxA, maxB) - Math.max(minA, minB);
    if (overlap <= 0) return null;
    if (overlap < best.depth) best = { depth: overlap, normal: n };
  }
  return best;
}

5. GJK: Distance Between Convex Shapes

The Gilbert-Johnson-Keerthi (GJK) algorithm answers "are these two convex shapes intersecting?" without enumerating any axes — it works for spheres, capsules, and arbitrary hulls alike, which is why it replaced SAT as the general-purpose narrow phase in Bullet, PhysX, and Box2D. GJK iteratively builds a simplex (point, line, triangle, tetrahedron) out of Minkowski-difference support points that gets closer to the origin each iteration:

function gjkIntersect(hullA, hullB) {
  let d = { x: 1, y: 0, z: 0 }; // arbitrary initial direction
  let simplex = [supportMinkowski(hullA, hullB, d)];
  d = negate(simplex[0]);

  for (let iter = 0; iter < 64; iter++) {
    const p = supportMinkowski(hullA, hullB, d);
    if (dot(p, d) < 0) return { hit: false }; // can't reach origin → no overlap
    simplex.push(p);

    const result = doSimplex(simplex, d); // mutates simplex + direction
    if (result.containsOrigin) return { hit: true, simplex };
    d = result.newDirection;
  }
  return { hit: false }; // exceeded iteration budget (near-degenerate case)
}

// doSimplex reduces the simplex to the smallest sub-simplex still
// closest to the origin, and returns a new search direction
// perpendicular to that sub-simplex, pointing toward the origin.
// Case "line": direction = triple-product toward origin
// Case "triangle": direction = face normal toward origin, or an edge
// Case "tetrahedron": if origin is inside all 4 faces → containsOrigin

Why this terminates: each new support point in direction d is guaranteed to make progress toward the origin along d (that's what "support" maximizes). If it doesn't cross the origin's half-space, no point of the Minkowski difference can, so the shapes cannot overlap — this is the early-out check dot(p, d) < 0.

6. EPA: Penetration Depth From a GJK Simplex

GJK alone only returns a boolean. When shapes overlap, the final tetrahedron GJK found encloses the origin but tells us nothing about how deep the penetration is — for that we need the Expanding Polytope Algorithm (EPA). EPA starts from GJK's terminating simplex and repeatedly expands it toward the Minkowski surface until it converges on the closest face to the origin, whose normal and distance give the contact normal and penetration depth:

function epa(hullA, hullB, simplex) {
  let polytope = tetrahedronToFaces(simplex); // 4 triangular faces, CCW winding

  for (let iter = 0; iter < 32; iter++) {
    // 1. Find the face closest to the origin
    const closest = findClosestFace(polytope); // { normal, distance, verts }

    // 2. Get a new support point in that face's normal direction
    const p = supportMinkowski(hullA, hullB, closest.normal);
    const d = dot(p, closest.normal);

    // 3. If the new point barely extends the polytope, we've converged
    if (d - closest.distance < 1e-5) {
      return { normal: closest.normal, depth: closest.distance,
               points: closest.verts };
    }

    // 4. Otherwise expand: remove all faces the new point can "see",
    //    collect the resulting hole's boundary edges, re-triangulate
    const { remainingFaces, boundaryEdges } = removeVisibleFaces(polytope, p);
    const newFaces = boundaryEdges.map(edge => makeFace(edge[0], edge[1], p));
    polytope = [...remainingFaces, ...newFaces];
  }
  return null; // exceeded iteration budget
}

Distance of a triangular face {v0,v1,v2} from the origin: n̂ = normalize((v1 - v0) \times (v2 - v0)) [outward-facing] dist = n̂ \cdot v0 [signed distance] Convergence criterion: support(closest.normal) \cdot n̂ - dist < ε \to the polytope surface is within ε of the true Minkowski boundary along this normal, so distance \approx penetration depth

EPA's cost is dominated by findClosestFace (linear scan over polytope faces, or a heap for large polytopes) and the face-count growth each iteration. In practice 5–15 iterations converge to sub-millimeter accuracy for typical game-scale meshes.

7. Contact Manifold Generation

A single normal + depth from SAT or EPA is enough to separate two shapes, but a stable stack of boxes needs multiple contact points per pair — otherwise a box resting on a face will slowly rotate off its supporting normal each frame. The standard technique is clipping: identify the reference face (on the body with the flattest face aligned to the normal) and the incident face (on the other body), then clip the incident face's edges against the reference face's side planes (Sutherland-Hodgman polygon clipping):

function clipManifold(referenceFace, incidentFace, normal) {
  let points = incidentFace.vertices;
  // Clip incident polygon against each side plane of the reference face
  for (const sidePlane of sidePlanesOf(referenceFace)) {
    points = clipPolygonAgainstPlane(points, sidePlane);
    if (points.length === 0) return []; // fully clipped away
  }
  // Keep only points still below the reference face (penetrating)
  return points
    .filter(p => dot(sub(p, referenceFace.point), normal) <= 0)
    .map(p => ({
      point: p,
      depth: -dot(sub(p, referenceFace.point), normal)
    }));
}

Manifold caching: production engines (Bullet's btPersistentManifold) keep contact points alive across frames while the pair stays close, matching new clip results to old points by proximity so warm-started impulses (see the physics engine tutorial) remain valid and stacks don't "pop" every frame.

8. Numerical Pitfalls and Robustness

SAT and EPA are geometrically simple but numerically fragile at the edges. The failure modes below account for most of the "jittery box" and "objects sink through the floor" bug reports in physics forums:

Near-parallel edge-cross axes: when two edges are nearly parallel, their cross product approaches zero length and normalizing it amplifies floating-point noise. Skip axes whose cross-product length is below a small epsilon (as in the OBB code above) rather than normalizing near-zero vectors.
Degenerate GJK simplices: if four support points end up nearly coplanar, the initial tetrahedron has ~zero volume and EPA's face-normal computation becomes unstable. Perturb the initial search directions slightly, or fall back to SAT for axis-aligned primitives where it's cheap and exact.
EPA non-termination: always cap iterations (e.g. 32) and treat exceeding the cap as "return the current best face" rather than looping — floating-point round-off can otherwise prevent the convergence epsilon from ever being satisfied.
Speculative contacts: instead of waiting for deep overlap (which EPA must then resolve), many modern engines (Rapier, Bullet's contact margin) run narrow phase against slightly inflated hulls and generate contacts before actual penetration, keeping the solver in the cheap, numerically stable "shallow SAT" regime and avoiding EPA altogether in the common case.
Persistent manifold IDs: re-run full SAT/GJK/EPA only when the broad-phase pair is new or has moved significantly; otherwise incrementally refresh depths using the cached normal — this is what makes 60 FPS stacks of hundreds of boxes tractable.

Rule of thumb: use SAT directly whenever both shapes are boxes or convex hulls with few faces — it is exact, branch-predictable, and gives the MTV in one pass. Reach for GJK+EPA when you need spheres, capsules, or generic convex hulls in the same code path, since GJK's support-function abstraction hides shape-specific geometry behind one interface.