Triangle AABB Intersection Test

intersection.hpp

// reference: https://gdbooks.gitbooks.io/3dcollisions/content/Chapter4/aabb-triangle.html
//            https://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/code/tribox2.txt
#pragma once

#include <glm/glm.hpp>

namespace glm {
template<int I, typename T, qualifier Q>
inline __device__ __host__ glm::vec<I, T, Q> min(const glm::vec<I, T, Q>& v1, const glm::vec<I, T, Q>& v2, const glm::vec<I, T, Q>& v3) {
	return glm::min(v1, glm::min(v2, v3));
}
template<int I, typename T, qualifier Q>
inline __device__ __host__ glm::vec<I, T, Q> max(const glm::vec<I, T, Q>& v1, const glm::vec<I, T, Q>& v2, const glm::vec<I, T, Q>& v3) {
	return glm::max(v1, glm::max(v2, v3));
}
}

struct BoundingBox { // AABB
	glm::vec3 min{ FLT_MAX, FLT_MAX, FLT_MAX };
	glm::vec3 max{ -FLT_MAX, -FLT_MAX, -FLT_MAX };
	glm::vec3 center;
	glm::vec3 halfExtent;
};

struct Triangle {
	glm::vec3 vert0;
	glm::vec3 vert1;
	glm::vec3 vert2;
	BoundingBox bbox;
};

inline bool boxBoxIntersect(const BoundingBox& a, const BoundingBox& b) {
	const glm::vec3 amin{ a.min - b.center };
	const glm::vec3 amax{ a.max - b.center };
	if (amin.x > b.halfExtent.x || amax.x < -b.halfExtent.x ||
		amin.y > b.halfExtent.y || amax.y < -b.halfExtent.y ||
		amin.z > b.halfExtent.z || amax.z < -b.halfExtent.z)
		return false;

	return true;
}

inline bool axisProjectTest(
	const glm::vec3& v0, const glm::vec3& v1, const glm::vec3& v2, 
	const glm::vec3& u0, const glm::vec3& u1, const glm::vec3& u2, 
	const glm::vec3& e, const glm::vec3& axis) {
	float p0 = glm::dot(v0, axis);
	float p1 = glm::dot(v1, axis);
	float p2 = glm::dot(v2, axis);
	float r = e.x * fabsf(glm::dot(u0, axis)) +
		e.y * fabsf(glm::dot(u1, axis)) +
		e.z * fabsf(glm::dot(u2, axis));

	// This means BOTH of the points of the projected triangle
	// are outside the projected half-length of the AABB
	// Therefore the axis is seperating and we can exit
	return fmaxf(-fmaxf(p0, p1, p2), fminf(p0, p1, p2)) > r;
}

inline bool tirgBoxIntersect(const Triangle& t, const BoundingBox& b) {
	glm::vec3 v0 = t.vert0 - b.center;
	glm::vec3 v1 = t.vert1 - b.center;
	glm::vec3 v2 = t.vert2 - b.center;
	glm::vec3 extent = b.max - b.min;

	const glm::vec3 e0{ v1 - v0 };
	const glm::vec3 e1{ v2 - v1 };
	const glm::vec3 e2{ v0 - v2 };
	const glm::vec3 n{ glm::cross(e0, e1) };

	const glm::vec3 u0{ 1.0f, 0.0f, 0.0f };
	const glm::vec3 u1{ 0.0f, 1.0f, 0.0f };
	const glm::vec3 u2{ 0.0f, 0.0f, 1.0f };

	// We first test against 9 axis, these axis are given by
	// cross product combinations of the edges of the triangle
	// and the edges of the AABB.

	const glm::vec3 axis_u0_e0 = glm::cross(u0, e0);
	const glm::vec3 axis_u0_e1 = glm::cross(u0, e1);
	const glm::vec3 axis_u0_e2 = glm::cross(u0, e2);

	const glm::vec3 axis_u1_e0 = glm::cross(u1, e0);
	const glm::vec3 axis_u1_e1 = glm::cross(u1, e1);
	const glm::vec3 axis_u1_e2 = glm::cross(u2, e2);

	const glm::vec3 axis_u2_e0 = glm::cross(u2, e0);
	const glm::vec3 axis_u2_e1 = glm::cross(u2, e1);
	const glm::vec3 axis_u2_e2 = glm::cross(u2, e2);

	// test 9 axis
	if (axisProjectTest(v0, v1, v2, u0, u1, u2, extent, axis_u0_e0)) return false;
	if (axisProjectTest(v0, v1, v2, u0, u1, u2, extent, axis_u0_e1)) return false;
	if (axisProjectTest(v0, v1, v2, u0, u1, u2, extent, axis_u0_e2)) return false;
	if (axisProjectTest(v0, v1, v2, u0, u1, u2, extent, axis_u1_e0)) return false;
	if (axisProjectTest(v0, v1, v2, u0, u1, u2, extent, axis_u1_e1)) return false;
	if (axisProjectTest(v0, v1, v2, u0, u1, u2, extent, axis_u1_e2)) return false;
	if (axisProjectTest(v0, v1, v2, u0, u1, u2, extent, axis_u2_e0)) return false;
	if (axisProjectTest(v0, v1, v2, u0, u1, u2, extent, axis_u2_e1)) return false;
	if (axisProjectTest(v0, v1, v2, u0, u1, u2, extent, axis_u2_e2)) return false;

	if (boxBoxIntersect(t.bbox, b)) { // test 3 box normals. equivalent to test two aabb
		// test if the box intersects the plane of the triangle
		// compute plane equation of triangle: normal*x+d=0
		float d = -glm::dot(n, v0);
		glm::vec3 vmin, vmax;

		for (int q = 0; q <= 2; q++) {
			if (n[q] > 0.f) {
				vmin[q] = -b.halfExtent[q];
				vmax[q] = b.halfExtent[q];
			}
			else {
				vmin[q] = b.halfExtent[q];
				vmax[q] = -b.halfExtent[q];
			}
		}
		if (glm::dot(n, vmax) + d >= 0.f) return true;
	}

	return false;
}