Triangle/AABB (Triangle-Box) Intersection Test using the "Separating Axis Theorem" in Python

TriangleBoxIntersect.py

import numpy as np

# based on http://www.gamedev.net/topic/534655-aabb-triangleplane-intersection--distance-to-plane-is-incorrect-i-have-solved-it/

def intersects_box(triangle, box_center, box_extents):
    X, Y, Z = 0, 1, 2

    # Translate triangle as conceptually moving AABB to origin
    v0 = triangle[0] - box_center
    v1 = triangle[1] - box_center
    v2 = triangle[2] - box_center

    # Compute edge vectors for triangle
    f0 = triangle[1] - triangle[0]
    f1 = triangle[2] - triangle[1]
    f2 = triangle[0] - triangle[2]

    ## region Test axes a00..a22 (category 3)

    # Test axis a00
    a00 = np.array([0, -f0[Z], f0[Y]])
    p0 = np.dot(v0, a00)
    p1 = np.dot(v1, a00)
    p2 = np.dot(v2, a00)
    r = box_extents[Y] * abs(f0[Z]) + box_extents[Z] * abs(f0[Y])
    if (max(-max(p0, p1, p2), min(p0, p1, p2))) > r:
        return False

    # Test axis a01
    a01 = np.array([0, -f1[Z], f1[Y]])
    p0 = np.dot(v0, a01)
    p1 = np.dot(v1, a01)
    p2 = np.dot(v2, a01)
    r = box_extents[Y] * abs(f1[Z]) + box_extents[Z] * abs(f1[Y])
    if (max(-max(p0, p1, p2), min(p0, p1, p2))) > r:
        return False

    # Test axis a02
    a02 = np.array([0, -f2[Z], f2[Y]])
    p0 = np.dot(v0, a02)
    p1 = np.dot(v1, a02)
    p2 = np.dot(v2, a02)
    r = box_extents[Y] * abs(f2[Z]) + box_extents[Z] * abs(f2[Y])
    if (max(-max(p0, p1, p2), min(p0, p1, p2))) > r:
        return False

    # Test axis a10
    a10 = np.array([f0[Z], 0, -f0[X]])
    p0 = np.dot(v0, a10)
    p1 = np.dot(v1, a10)
    p2 = np.dot(v2, a10)
    r = box_extents[X] * abs(f0[Z]) + box_extents[Z] * abs(f0[X])
    if (max(-max(p0, p1, p2), min(p0, p1, p2))) > r:
        return False

    # Test axis a11
    a11 = np.array([f1[Z], 0, -f1[X]])
    p0 = np.dot(v0, a11)
    p1 = np.dot(v1, a11)
    p2 = np.dot(v2, a11)
    r = box_extents[X] * abs(f1[Z]) + box_extents[Z] * abs(f1[X])
    if (max(-max(p0, p1, p2), min(p0, p1, p2))) > r:
        return False

    # Test axis a12
    a11 = np.array([f2[Z], 0, -f2[X]])
    p0 = np.dot(v0, a11)
    p1 = np.dot(v1, a11)
    p2 = np.dot(v2, a11)
    r = box_extents[X] * abs(f2[Z]) + box_extents[Z] * abs(f2[X])
    if (max(-max(p0, p1, p2), min(p0, p1, p2))) > r:
        return False

    # Test axis a20
    a20 = np.array([-f0[Y], f0[X], 0])
    p0 = np.dot(v0, a20)
    p1 = np.dot(v1, a20)
    p2 = np.dot(v2, a20)
    r = box_extents[X] * abs(f0[Y]) + box_extents[Y] * abs(f0[X])
    if (max(-max(p0, p1, p2), min(p0, p1, p2))) > r:
        return False

    # Test axis a21
    a21 = np.array([-f1[Y], f1[X], 0])
    p0 = np.dot(v0, a21)
    p1 = np.dot(v1, a21)
    p2 = np.dot(v2, a21)
    r = box_extents[X] * abs(f1[Y]) + box_extents[Y] * abs(f1[X])
    if (max(-max(p0, p1, p2), min(p0, p1, p2))) > r:
        return False

    # Test axis a22
    a22 = np.array([-f2[Y], f2[X], 0])
    p0 = np.dot(v0, a22)
    p1 = np.dot(v1, a22)
    p2 = np.dot(v2, a22)
    r = box_extents[X] * abs(f2[Y]) + box_extents[Y] * abs(f2[X])
    if (max(-max(p0, p1, p2), min(p0, p1, p2))) > r:
        return False

    ## endregion

    ## region Test the three axes corresponding to the face normals of AABB b (category 1)

    # Exit if...
    # ... [-extents.X, extents.X] and [min(v0.X,v1.X,v2.X), max(v0.X,v1.X,v2.X)] do not overlap
    if max(v0[X], v1[X], v2[X]) < -box_extents[X] or min(v0[X], v1[X], v2[X]) > box_extents[X]:
        return False

    # ... [-extents.Y, extents.Y] and [min(v0.Y,v1.Y,v2.Y), max(v0.Y,v1.Y,v2.Y)] do not overlap
    if max(v0[Y], v1[Y], v2[Y]) < -box_extents[Y] or min(v0[Y], v1[Y], v2[Y]) > box_extents[Y]:
        return False

    # ... [-extents.Z, extents.Z] and [min(v0.Z,v1.Z,v2.Z), max(v0.Z,v1.Z,v2.Z)] do not overlap
    if max(v0[Z], v1[Z], v2[Z]) < -box_extents[Z] or min(v0[Z], v1[Z], v2[Z]) > box_extents[Z]:
        return False

    ## endregion

    ## region Test separating axis corresponding to triangle face normal (category 2)

    plane_normal = np.cross(f0, f1)
    plane_distance = np.dot(plane_normal, v0)

    # Compute the projection interval radius of b onto L(t) = b.c + t * p.n
    r = box_extents[X] * abs(plane_normal[X]) + box_extents[Y] * abs(plane_normal[Y]) + box_extents[Z] * abs(plane_normal[Z])

    # Intersection occurs when plane distance falls within [-r,+r] interval
    if plane_distance > r:
        return False

    ## endregion

    return True


# example:
# intersects_box(np.array([[20,9,44], [7,7,7], [8,8,8]]), np.array([0,1,3]), np.array([8,8,4]))

In a note called "Fast 3D Triangle-Box Overlap Testing" of Tomas Akenine-Moller it is written
r = hx|ax| + hy|ay| + hz|az| = hy|ay| + hz|az|
On line 25 it is written
r = box_extents[Y] * abs(f0[Z]) + box_extents[Z] * abs(f0[Y])
that does not match with the line above. Why?

The various solutions posted here appears to yield incorrect results as some triangles are missed on some of the geometry I have been testing e.g spheres and the monkey head from Blender (Suzanne). Bronson Zgeb has a solution that appears to work well: https://bronsonzgeb.com/index.php/2021/05/29/gpu-mesh-voxelizer-part-2/

I am using Godot C#, but converting it to Unity (or otherwise) should be non-trivial. Please note the way in which the AABB box extents and box centre are calculated:

// Drop into your gameplay loop or rename to whatever you wish
void CheckLoop(){
  Vector3 box_extents = (aabbMax - aabbMin)/2;
  Vector3 box_centre = aabbMax - box_extents;

  if(AABB_Tri_Intersect(triangle.v0, triangle.v1, triangle.v2, box_centre, box_extents)){
    // Do stuff if intersecting
  }
}

bool AABB_Tri_Intersect(Vector3 v0, Vector3 v1, Vector3 v2, Vector3 aabbCentre, Vector3 aabbExtents){
  v0 -= aabbCentre;
  v1 -= aabbCentre;
  v2 -= aabbCentre;
  
  Vector3 ab = (v1 - v0).Normalized();
  Vector3 bc = (v2 - v1).Normalized();
  Vector3 ca = (v0 - v2).Normalized();
  
  //Cross ab, bc, and ca with (1, 0, 0)
  Vector3 a00 = new (0.0f, -ab.Z, ab.Y);
  Vector3 a01 = new (0.0f, -bc.Z, bc.Y);
  Vector3 a02 = new (0.0f, -ca.Z, ca.Y);
  
  //Cross ab, bc, and ca with (0, 1, 0)
  Vector3 a10 = new (ab.Z, 0.0f, -ab.X);
  Vector3 a11 = new (bc.Z, 0.0f, -bc.X);
  Vector3 a12 = new (ca.Z, 0.0f, -ca.X);
  
  //Cross ab, bc, and ca with (0, 0, 1)
  Vector3 a20 = new (-ab.Y, ab.X, 0.0f);
  Vector3 a21 = new (-bc.Y, bc.X, 0.0f);
  Vector3 a22 = new (-ca.Y, ca.X, 0.0f);
  
  if (!AABB_Tri_SAT(v0, v1, v2, aabbExtents, a00) ||
    !AABB_Tri_SAT(v0, v1, v2, aabbExtents, a01) ||
    !AABB_Tri_SAT(v0, v1, v2, aabbExtents, a02) ||
    !AABB_Tri_SAT(v0, v1, v2, aabbExtents, a10) ||
    !AABB_Tri_SAT(v0, v1, v2, aabbExtents, a11) ||
    !AABB_Tri_SAT(v0, v1, v2, aabbExtents, a12) ||
    !AABB_Tri_SAT(v0, v1, v2, aabbExtents, a20) ||
    !AABB_Tri_SAT(v0, v1, v2, aabbExtents, a21) ||
    !AABB_Tri_SAT(v0, v1, v2, aabbExtents, a22) ||
    !AABB_Tri_SAT(v0, v1, v2, aabbExtents, new Vector3(1, 0, 0)) ||
    !AABB_Tri_SAT(v0, v1, v2, aabbExtents, new Vector3(0, 1, 0)) ||
    !AABB_Tri_SAT(v0, v1, v2, aabbExtents, new Vector3(0, 0, 1)) ||
    !AABB_Tri_SAT(v0, v1, v2, aabbExtents, ab.Cross(bc))) 
  {
    return false;
  }

  return true;
}

bool AABB_Tri_SAT(Vector3 v0, Vector3 v1, Vector3 v2, Vector3 aabbExtents, Vector3 axis){
  float p0 = v0.Dot(axis);
  float p1 = v1.Dot(axis);
  float p2 = v2.Dot(axis);
  
  float r = aabbExtents.X * Mathf.Abs(new Vector3(1, 0, 0).Dot(axis)) +
	  aabbExtents.Y * Mathf.Abs(new Vector3(0, 1, 0).Dot(axis)) +
	  aabbExtents.Z * Mathf.Abs(new Vector3(0, 0, 1).Dot(axis));
  
  float maxP = Mathf.Max(p0, Mathf.Max(p1, p2));
  float minP = Mathf.Min(p0, Mathf.Min(p1, p2));
  
  return !(Mathf.Max(-maxP, minP) > r);
}

@KyriaAnnwyn

In some cases this works incorrectly, try these triangle and cube: intersects_box(np.array([[-0.236396998,0.0570373014,-0.0280157998], [-0.232334003,0.0409222990,0.0128640002], [0.197703004,-0.190671995,-0.312638015]]), np.array([0.0159800090,-0.0884050056,-0.200349063]), np.array([0.01,0.01,0.01]))

You are wrong, as your blender screenshot (showing its not intersecting) does not match your written numbers.
If i enter your written numbers into Blender myself i get this, showing that your Triangle/Box combo indeed does intersect, just like the SAT-Function sees it: