-
-
Save zvonicek/fe73ba9903f49d57314cf7e8e0f05dcf to your computer and use it in GitHub Desktop.
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);
}
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:
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]))