Ray intersection with oriented bounding box (OBB)

OBB_RAY_INTERSECT.js

const vec_sub = (a, b) => ({ x: a.x - b.x, y: a.y - b.y, z: a.z - b.z });
const vec_dot = (a, b) => a.x * b.x + a.y * b.y + a.z * b.z;
const vec_len = a => Math.sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
const vec_norm = a => {
	const norm = vec_len(a);
	return { x: a.x / norm, y: a.y / norm, z: a.z / norm };
}

/*
 * JS version of
 * https://github.com/opengl-tutorials/ogl/blob/master/misc05_picking/misc05_picking_custom.cpp#L83
 */
function testIntersectionRayOBB(
	ray_origin,        // Ray origin, in world space
	ray_direction,     // Ray direction (NOT target position!), in world space. Must be normalize()'d.
	aabb_min,          // Minimum X,Y,Z coords of the mesh when not transformed at all.
	aabb_max,          // Maximum X,Y,Z coords. Often aabb_min*-1 if your mesh is centered, but it's not always the case.
	matrix             // Transformation applied to the mesh (which will thus be also applied to its bounding box)
) {

	let tMin = 0.0;
	let tMax = Infinity; //100000.0;
	//const threshold = 0.001;
	const threshold = 0.0000000001;

	const pos = { x: matrix[3].x, y: matrix[3].y, z: matrix[3].z };
	const delta = vec_sub(pos, ray_origin);

	// Test intersection with the 2 planes perpendicular to the OBB's X axis
	{
		const xaxis = { x: matrix[0].x, y: matrix[0].y, z: matrix[0].z };
		const e = vec_dot(xaxis, delta);
		const f = vec_dot(ray_direction, xaxis);

		if (Math.abs(f) > threshold) { // Standard case

			let t1 = (e + aabb_min.x) / f; // Intersection with the "left" plane
			let t2 = (e + aabb_max.x) / f; // Intersection with the "right" plane
			// t1 and t2 now contain distances betwen ray origin and ray-plane intersections

			// We want t1 to represent the nearest intersection, 
			// so if it's not the case, invert t1 and t2
			if (t1 > t2) {
				const w = t1; t1 = t2; t2 = w; // swap t1 and t2
			}

			// tMax is the nearest "far" intersection (amongst the X,Y and Z planes pairs)
			if (t2 < tMax)
				tMax = t2;
			// tMin is the farthest "near" intersection (amongst the X,Y and Z planes pairs)
			if (t1 > tMin)
				tMin = t1;

			// And here's the trick :
			// If "far" is closer than "near", then there is NO intersection.
			// See the images in the tutorials for the visual explanation.
			if (tMax < tMin)
				return null;

		} else { // Rare case : the ray is almost parallel to the planes, so they don't have any "intersection"
			if (-e + aabb_min.x > 0.0 || -e + aabb_max.x < 0.0)
				return null;
		}
	}

	// Test intersection with the 2 planes perpendicular to the OBB's Y axis
	// Exactly the same thing than above.
	{
		const yaxis = { x: matrix[1].x, y: matrix[1].y, z: matrix[1].z };
		const e = vec_dot(yaxis, delta);
		const f = vec_dot(ray_direction, yaxis);

		if (Math.abs(f) > threshold) {

			let t1 = (e + aabb_min.y) / f;
			let t2 = (e + aabb_max.y) / f;

			if (t1 > t2) {
				const w = t1; t1 = t2; t2 = w;
			}

			if (t2 < tMax)
				tMax = t2;
			if (t1 > tMin)
				tMin = t1;
			if (tMin > tMax)
				return null;

		} else {
			if (-e + aabb_min.y > 0.0 || -e + aabb_max.y < 0.0)
				return null;
		}
	}

	// Test intersection with the 2 planes perpendicular to the OBB's Z axis
	// Exactly the same thing than above.
	{
		const zaxis = { x: matrix[2].x, y: matrix[2].y, z: matrix[2].z };
		const e = vec_dot(zaxis, delta);
		const f = vec_dot(ray_direction, zaxis);

		if (Math.abs(f) > threshold) {

			let t1 = (e + aabb_min.z) / f;
			let t2 = (e + aabb_max.z) / f;

			if (t1 > t2) {
				const w = t1; t1 = t2; t2 = w;
			}

			if (t2 < tMax)
				tMax = t2;
			if (t1 > tMin)
				tMin = t1;
			if (tMin > tMax)
				return null;

		} else {
			if (-e + aabb_min.z > 0.0 || -e + aabb_max.z < 0.0)
				return null;
		}
	}

	// intersection_distance: distance between ray_origin and the intersection with the OBB
	return tMin;
}

I translated your code to a GLSL version nonetheless it doesn't seems to work for me. I wrote it as follows:

float test_intersection_ray_OBB(
vec3 ray_origin, // Ray origin, in world space
vec3 ray_direction, // Ray direction (NOT target position!), in world space. Must be normalize()'d.
vec3 aabb_min, // Minimum X,Y,Z coords of the mesh when not transformed at all.
vec3 aabb_max, // Maximum X,Y,Z coords. Often aabb_min*-1 if your mesh is centered, but it's not always the case.
mat4 matrix // Transformation applied to the mesh (which will thus be also applied to its bounding box)
) {

float tMin = 0.0;
float tMax = 100000.0; //100000.0;
//const threshold = 0.001;
const float threshold = 0.0000000001;

const vec3 pos = vec3(matrix[3].x, matrix[3].y, matrix[3].z);
const vec3 delta = pos - ray_origin;

// Test intersection with the 2 planes perpendicular to the OBB's X axis
{
	const vec3 xaxis = vec3(matrix[0].x, matrix[0].y, matrix[0].z);
	const float e = dot(xaxis, delta);
	const float f = dot(ray_direction, xaxis);

	if (abs(f) > threshold) { // Standard case

		float t1 = (e + aabb_min.x) / f; // Intersection with the "left" plane
		float t2 = (e + aabb_max.x) / f; // Intersection with the "right" plane
		// t1 and t2 now contain distances betwen ray origin and ray-plane intersections

		// We want t1 to represent the nearest intersection, 
		// so if it's not the case, invert t1 and t2
		if (t1 > t2) {
			const float w = t1; t1 = t2; t2 = w; // swap t1 and t2
		}

		// tMax is the nearest "far" intersection (amongst the X,Y and Z planes pairs)
		if (t2 < tMax)
			tMax = t2;
		// tMin is the farthest "near" intersection (amongst the X,Y and Z planes pairs)
		if (t1 > tMin)
			tMin = t1;

		// And here's the trick :
		// If "far" is closer than "near", then there is NO intersection.
		// See the images in the tutorials for the visual explanation.
		if (tMax < tMin)
			return -1;

	} else { // Rare case : the ray is almost parallel to the planes, so they don't have any "intersection"
		if (-e + aabb_min.x > 0.0 || -e + aabb_max.x < 0.0)
			return -1;
	}
}

// Test intersection with the 2 planes perpendicular to the OBB's Y axis
// Exactly the same thing than above.
{
	const vec3 yaxis = vec3(matrix[1].x, matrix[1].y, matrix[1].z);
	const float e = dot(yaxis, delta);
	const float f = dot(ray_direction, yaxis);

	if (abs(f) > threshold) {

		float t1 = (e + aabb_min.y) / f;
		float t2 = (e + aabb_max.y) / f;

		if (t1 > t2) {
			const float w = t1; t1 = t2; t2 = w;
		}

		if (t2 < tMax)
			tMax = t2;
		if (t1 > tMin)
			tMin = t1;
		if (tMin > tMax)
			return -1;

	} else {
		if (-e + aabb_min.y > 0.0 || -e + aabb_max.y < 0.0)
			return -1;
	}
}

// Test intersection with the 2 planes perpendicular to the OBB's Z axis
// Exactly the same thing than above.
{
	const vec3 zaxis = vec3(matrix[2].x, matrix[2].y, matrix[2].z);
	const float e = dot(zaxis, delta);
	const float f = dot(ray_direction, zaxis);

	if (abs(f) > threshold) {

		float t1 = (e + aabb_min.z) / f;
		float t2 = (e + aabb_max.z) / f;

		if (t1 > t2) {
			const float w = t1; t1 = t2; t2 = w;
		}

		if (t2 < tMax)
			tMax = t2;
		if (t1 > tMin)
			tMin = t1;
		if (tMin > tMax)
			return -1;

	} else {
		if (-e + aabb_min.z > 0.0 || -e + aabb_max.z < 0.0)
			return -1;
	}
}

// intersection_distance: distance between ray_origin and the intersection with the OBB
return tMin;

}

I call it as:
main {

....

Ray ray;
ray.origin = gl_WorldRayOriginEXT;
ray.direction = gl_WorldRayDirectionEXT;
ray.direction = normalize(ray.direction);
Aabb aabb;
aabb.minimum = voxel_center - vec3(voxel.side);
aabb.maximum = voxel_center + vec3(voxel.side);

float tHit = test_intersection_ray_OBB(ray.origin, normalize(ray.direction),
aabb.minimum, aabb.maximum,
transformation);
}

Where for instance mat4 transformation data is:
{_11: 1.00000000, _12: 0.918377697, _13: 0.00000000, _14: 10.7665958,
_21: -0.918377697, _22: 1.00000000, _23: 0.00000000, _24: 0.283537865,
_31: 0.00000000, _32: 0.00000000, _33: 1.00000000, _34: -9.78049088,
_41: 0.00000000, _42: 0.00000000, _43: 0.00000000 _44: 1.00000000}

Greetings,
Jaime.