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]))

@KyriaAnnwyn
The author of the original post explained that this line:

plane_distance = np.dot(plane_normal, v0)

Should be:

plane_distance = np.dot(plane_normal, triangle[0])

@KyriaAnnwyn The author of the original post explained that this line:

plane_distance = np.dot(plane_normal, v0)

Should be:

plane_distance = np.dot(plane_normal, triangle[0])

That's not correct, that would only work when box center is at position 0,0,0
That part calculates distance between triangle and the box center in the direction of triangle normal.
By supplying v0, you calculate distance between box center and the surface of the triangle (in the direction of triangle normal).
Which is then compared to box extents projected onto triangle normal.

v0 is position of the triangle point relative to the box center.
A line below that is calculation of box extents projected onto triangle normal.

If there's triangle[0], box extents couldn't be used, but instead box corners would have to be used.

There's different error, missing abs, the line should be:

plane_distance = abs(np.dot(plane_normal, v0))

Full C# version (Unity3d engine - Vector3, Mathf):

public static float Max(float a, float b, float c)
{
    return Mathf.Max(Mathf.Max(a, b), c);
}

public static float Min(float a, float b, float c)
{
    return Mathf.Min(Mathf.Min(a, b), c);
}

public static bool BoxIntersectsTriangle(Bounds bounds, Vector3 triangle0, Vector3 triangle1, Vector3 triangle2)
{
    var box_center = bounds.center;
    var box_extents = bounds.extents;

    var v0 = triangle0 - box_center;
    var v1 = triangle1 - box_center;
    var v2 = triangle2 - box_center;

    // Compute edge vectors for triangle
    var f0 = triangle1 - triangle0;
    var f1 = triangle2 - triangle1;
    var f2 = triangle0 - triangle2;

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

    // Test axis a00
    var a00 = new Vector3(0, -f0.z, f0.y);
    var p0 = Vector3.Dot(v0, a00);
    var p1 = Vector3.Dot(v1, a00);
    var p2 = Vector3.Dot(v2, a00);
    var r = box_extents.y * Mathf.Abs(f0.z) + box_extents.z * Mathf.Abs(f0.y);
    if (Mathf.Max(-Max(p0, p1, p2), Min(p0, p1, p2)) > r)
        return false;

    // Test axis a01
    var a01 = new Vector3(0, -f1.z, f1.y);
    p0 = Vector3.Dot(v0, a01);
    p1 = Vector3.Dot(v1, a01);
    p2 = Vector3.Dot(v2, a01);
    r = box_extents.y * Mathf.Abs(f1.z) + box_extents.z * Mathf.Abs(f1.y);
    if (Mathf.Max(-Max(p0, p1, p2), Min(p0, p1, p2)) > r)
        return false;

    // Test axis a02
    var a02 = new Vector3(0, -f2.z, f2.y);
    p0 = Vector3.Dot(v0, a02);
    p1 = Vector3.Dot(v1, a02);
    p2 = Vector3.Dot(v2, a02);
    r = box_extents.y * Mathf.Abs(f2.z) + box_extents.z * Mathf.Abs(f2.y);
    if (Mathf.Max(-Max(p0, p1, p2), Min(p0, p1, p2)) > r)
        return false;

    // Test axis a10
    var a10 = new Vector3(f0.z, 0, -f0.x);
    p0 = Vector3.Dot(v0, a10);
    p1 = Vector3.Dot(v1, a10);
    p2 = Vector3.Dot(v2, a10);
    r = box_extents.x * Mathf.Abs(f0.z) + box_extents.z * Mathf.Abs(f0.x);
    if (Mathf.Max(-Max(p0, p1, p2), Min(p0, p1, p2)) > r)
        return false;

    // Test axis a11
    var a11 = new Vector3(f1.z, 0, -f1.x);
    p0 = Vector3.Dot(v0, a11);
    p1 = Vector3.Dot(v1, a11);
    p2 = Vector3.Dot(v2, a11);
    r = box_extents.x * Mathf.Abs(f1.z) + box_extents.z * Mathf.Abs(f1.x);
    if (Mathf.Max(-Max(p0, p1, p2), Min(p0, p1, p2)) > r)
        return false;

    // Test axis a12
    var a12 = new Vector3(f2.z, 0, -f2.x);
    p0 = Vector3.Dot(v0, a12);
    p1 = Vector3.Dot(v1, a12);
    p2 = Vector3.Dot(v2, a12);
    r = box_extents.x * Mathf.Abs(f2.z) + box_extents.z * Mathf.Abs(f2.x);
    if (Mathf.Max(-Max(p0, p1, p2), Min(p0, p1, p2)) > r)
        return false;

    // Test axis a20
    var a20 = new Vector3(-f0.y, f0.x, 0);
    p0 = Vector3.Dot(v0, a20);
    p1 = Vector3.Dot(v1, a20);
    p2 = Vector3.Dot(v2, a20);
    r = box_extents.x * Mathf.Abs(f0.y) + box_extents.y * Mathf.Abs(f0.x);
    if (Mathf.Max(-Max(p0, p1, p2), Min(p0, p1, p2)) > r)
        return false;

    // Test axis a21
    var a21 = new Vector3(-f1.y, f1.x, 0);
    p0 = Vector3.Dot(v0, a21);
    p1 = Vector3.Dot(v1, a21);
    p2 = Vector3.Dot(v2, a21);
    r = box_extents.x * Mathf.Abs(f1.y) + box_extents.y * Mathf.Abs(f1.x);
    if (Mathf.Max(-Max(p0, p1, p2), Min(p0, p1, p2)) > r)
        return false;

    // Test axis a22
    var a22 = new Vector3(-f2.y, f2.x, 0);
    p0 = Vector3.Dot(v0, a22);
    p1 = Vector3.Dot(v1, a22);
    p2 = Vector3.Dot(v2, a22);
    r = box_extents.x * Mathf.Abs(f2.y) + box_extents.y * Mathf.Abs(f2.x);
    if (Mathf.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 || 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 || 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 || Min(v0.z, v1.z, v2.z) > box_extents.z)
        return false;

    //// endregion

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

    var plane_normal = Vector3.Cross(f0, f1);
    var plane_distance = Mathf.Abs(Vector3.Dot(plane_normal, v0));

    // Compute the projection interval radius of b onto L(t) = b.c + t * p.n
    r = box_extents.x * Mathf.Abs(plane_normal.x) + box_extents.y * Mathf.Abs(plane_normal.y) + box_extents.z * Mathf.Abs(plane_normal.z);

    // Intersection occurs when plane distance falls within [-r,+r] interval
    if (plane_distance > r)
        return false;

    //// endregion

    return true;
}

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: