3 replies to this topic

[taby]

I've got some ray-AABB intersection code that I was hoping someone could take a look at to see if they see any major errors. As far as I know, it works, but I thought I would run it past the community. The code is fair game for anyone who would like to use it as well.

void ray_aabb_in(const vector_3 &start, const vector_3 &normalized_dir, const vector_3 &aabb_min, const vector_3 &aabb_max, vector_3 &entry_point)
{
	// Axis aligned bounding box face normals.
	static const vector_3 x_plus_normal(1, 0, 0);
	static const vector_3 x_minus_normal(-1, 0, 0);
	static const vector_3 y_plus_normal(0, 1, 0);
	static const vector_3 y_minus_normal(0, -1, 0);
	static const vector_3 z_plus_normal(0, 0, 1);
	static const vector_3 z_minus_normal(0, 0, -1);

	// If within the box, the entry point is the ray origin.
	if(within(start.x, aabb_min.x, aabb_max.x) && within(start.y, aabb_min.y, aabb_max.y) && within(start.z, aabb_min.z, aabb_max.z))
	{
		entry_point = start;
		return;
	}

	// Find which plane normal points opposite to normalized_dir.

	// Test against +x plane.
	if(start.x > aabb_max.x && 0 > normalized_dir.dot(x_plus_normal))
	{
		// Extend the unit line segment to whatever length that meets the +x plane.
		vector_3 final_dir = normalized_dir;
		final_dir *= d(start.x, aabb_max.x) / d(start.x, start.x + normalized_dir.x);

		vector_3 potential_point = start;
		potential_point += final_dir;

		// Is the end of the extended line segment within the volume defined by the infinitely extended yz plane?
		if(within(potential_point.y, aabb_min.y, aabb_max.y) && within(potential_point.z, aabb_min.z, aabb_max.z))
		{
			entry_point = potential_point;
			return;
		}
	}

	// Test again -x plane.
	else if(start.x < aabb_min.x && 0 > normalized_dir.dot(x_minus_normal))
	{
		vector_3 final_dir = normalized_dir;
		final_dir *= d(start.x, aabb_min.x) / d(start.x, start.x + normalized_dir.x);

		vector_3 potential_point = start;
		potential_point += final_dir;

		if(within(potential_point.y, aabb_min.y, aabb_max.y) && within(potential_point.z, aabb_min.z, aabb_max.z))
		{
			entry_point = potential_point;
			return;
		}
	}

	// Test against +y plane.
	if(start.y > aabb_max.y && 0 > normalized_dir.dot(y_plus_normal))
	{
		vector_3 final_dir = normalized_dir;
		final_dir *= d(start.y, aabb_max.y) / d(start.y, start.y + normalized_dir.y);

		vector_3 potential_point = start;
		potential_point += final_dir;

		if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.z, aabb_min.z, aabb_max.z))
		{
			entry_point = potential_point;
			return;
		}
	}

	// Test against -y plane.
	else if(start.y < aabb_min.y && 0 > normalized_dir.dot(y_minus_normal))
	{
		vector_3 final_dir = normalized_dir;
		final_dir *= d(start.y, aabb_min.y) / d(start.y, start.y + normalized_dir.y);

		vector_3 potential_point = start;
		potential_point += final_dir;

		if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.z, aabb_min.z, aabb_max.z))
		{
			entry_point = potential_point;
			return;
		}
	}

	// Test against +z plane.
	if(start.z > aabb_max.z && 0 > normalized_dir.dot(z_plus_normal))
	{
		vector_3 final_dir = normalized_dir;
		final_dir *= d(start.z, aabb_max.z) / d(start.z, start.z + normalized_dir.z);

		vector_3 potential_point = start;
		potential_point += final_dir;

		if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.y, aabb_min.y, aabb_max.y))
		{
			entry_point = potential_point;
			return;
		}
	}

	// Test against -z plane.
	else if(start.z < aabb_min.z && 0 > normalized_dir.dot(z_minus_normal))
	{
		vector_3 final_dir = normalized_dir;
		final_dir *= d(start.z, aabb_min.z) / d(start.z, start.z + normalized_dir.z);

		vector_3 potential_point = start;
		potential_point += final_dir;

		if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.y, aabb_min.y, aabb_max.y))
		{
			entry_point = potential_point;
			return;
		}
	}

	// If here, no match found, and entry_point is invalidated.
	entry_point.x = numeric_limits<double>::infinity();
}

void ray_aabb_out(const vector_3 &start, const vector_3 &normalized_dir, const vector_3 &aabb_min, const vector_3 &aabb_max, vector_3 &exit_point)
{
	// Axis aligned bounding box face normals.
	static const vector_3 x_plus_normal(1, 0, 0);
	static const vector_3 x_minus_normal(-1, 0, 0);
	static const vector_3 y_plus_normal(0, 1, 0);
	static const vector_3 y_minus_normal(0, -1, 0);
	static const vector_3 z_plus_normal(0, 0, 1);
	static const vector_3 z_minus_normal(0, 0, -1);

	// Find which plane normal points in the same direction as normalized_dir.

	// Test against +x plane.
	if(start.x < aabb_max.x && 0 < normalized_dir.dot(x_plus_normal))
	{
		// Extend the unit line segment to whatever length that meets the +x plane.
		vector_3 final_dir = normalized_dir;
		final_dir *= d(start.x, aabb_max.x) / d(start.x, start.x + normalized_dir.x);

		vector_3 potential_point = start;
		potential_point += final_dir;

		// Is the end of the extended line segment within the volume defined by the infinitely extended yz plane?
		if(within(potential_point.y, aabb_min.y, aabb_max.y) && within(potential_point.z, aabb_min.z, aabb_max.z))
		{
			exit_point = potential_point;
			return;
		}
	}

	// Test against -x plane.
	else if(start.x > aabb_min.x && 0 < normalized_dir.dot(x_minus_normal))
	{
		vector_3 final_dir = normalized_dir;
		final_dir *= d(start.x, aabb_min.x) / d(start.x, start.x + normalized_dir.x);

		vector_3 potential_point = start;
		potential_point += final_dir;

		if(within(potential_point.y, aabb_min.y, aabb_max.y) && within(potential_point.z, aabb_min.z, aabb_max.z))
		{
			exit_point = potential_point;
			return;
		}
	}

	// Test against +y plane.
	if(start.y < aabb_max.y && 0 < normalized_dir.dot(y_plus_normal))
	{
		vector_3 final_dir = normalized_dir;
		final_dir *= d(start.y, aabb_max.y) / d(start.y, start.y + normalized_dir.y);

		vector_3 potential_point = start;
		potential_point += final_dir;

		if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.z, aabb_min.z, aabb_max.z))
		{
			exit_point = potential_point;
			return;
		}
	}

	// Test against -y plane.
	else if(start.y > aabb_min.y && 0 < normalized_dir.dot(y_minus_normal))
	{
		vector_3 final_dir = normalized_dir;
		final_dir *= d(start.y, aabb_min.y) / d(start.y, start.y + normalized_dir.y);

		vector_3 potential_point = start;
		potential_point += final_dir;

		if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.z, aabb_min.z, aabb_max.z))
		{
			exit_point = potential_point;
			return;
		}
	}

	// Test against +z plane.
	if(start.z < aabb_max.z && 0 < normalized_dir.dot(z_plus_normal))
	{
		vector_3 final_dir = normalized_dir;
		final_dir *= d(start.z, aabb_max.z) / d(start.z, start.z + normalized_dir.z);

		vector_3 potential_point = start;
		potential_point += final_dir;

		if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.y, aabb_min.y, aabb_max.y))
		{
			exit_point = potential_point;
			return;
		}
	}

	// Test against -z plane.
	else if(start.z > aabb_min.z && 0 < normalized_dir.dot(z_minus_normal))
	{
		vector_3 final_dir = normalized_dir;
		final_dir *= d(start.z, aabb_min.z) / d(start.z, start.z + normalized_dir.z);

		vector_3 potential_point = start;
		potential_point += final_dir;

		if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.y, aabb_min.y, aabb_max.y))
		{
			exit_point = potential_point;
			return;
		}
	}

	// If here, no match found, and exit_point is invalidated.
	exit_point.x = numeric_limits<double>::infinity();
}

bool ray_aabb_intersection(const vector_3 &start, const vector_3 &normalized_dir, const vector_3 &aabb_min, const vector_3 &aabb_max, vector_3 &entry_point, vector_3 &exit_point)
{
	// Find potentially valid entry point.
	ray_aabb_in(start, normalized_dir, aabb_min, aabb_max, entry_point);

	// Is entry point valid?
	if(numeric_limits<double>::infinity() == entry_point.x)
		return false;

	// If so, find exit point.
	ray_aabb_out(start, normalized_dir, aabb_min, aabb_max, exit_point);

	return true;
}

[taby]

Here's a visual of many rays (top) intersecting with a series of axis-aligned bounding boxes (bottom).

[JuNC]

I haven't waded through the entire code so maybe you need something else, but the following is about as simple as you can get:

/*

 * Ray-box intersection using IEEE numerical properties to ensure that the

 * test is both robust and efficient, as described in:

 *

 *      Amy Williams, Steve Barrus, R. Keith Morley, and Peter Shirley

 *      "An Efficient and Robust Ray-Box Intersection Algorithm"

 *      Journal of graphics tools, 10(1):49-54, 2005

 *

 */

static bool intersect_bbox(const Vec3f bbox[2],const Ray& ray, float t0, float t1)

{

  float tmin = (bbox[ray.sign[0]][0] - ray.org[0]) * ray.inv_dir[0];

  float tmax = (bbox[1-ray.sign[0]][0] - ray.org[0]) * ray.inv_dir[0];



  float tymin = (bbox[ray.sign[1]][1] - ray.org[1]) * ray.inv_dir[1];

  float tymax = (bbox[1-ray.sign[1]][1] - ray.org[1]) * ray.inv_dir[1];



  if ( (tmin > tymax) || (tymin > tmax) )

    return false;

  if (tymin > tmin)

    tmin = tymin;

  if (tymax < tmax)

    tmax = tymax;



  float tzmin = (bbox[ray.sign[2]][2] - ray.org[2]) * ray.inv_dir[2];

  float tzmax = (bbox[1-ray.sign[2]][2] - ray.org[2]) * ray.inv_dir[2];



  if ( (tmin > tzmax) || (tzmin > tmax) )

    return false;

  if (tzmin > tmin)

    tmin = tzmin;

  if (tzmax < tmax)

    tmax = tzmax;

  return ( (tmin < t1) && (tmax > t0) );



}

You can get an intersection point by grabbing tmin (if return is true).

ray.inv_dir[i] is 1.0 / ray.dir[i] and ray.sign[i] = (ray.inv_dir[i] < 0);

[taby]

Hey thank you very much. I will try this code and see how it compares (hopefully it's identical!).