August 27, 2018

I’ve been doing some game dev stuff lately and I needed to intersect a ray with an octree of triangles, for collision detection. I first implemented a naive algorithm that simply checked if the AABB of each octant intersected the ray, then found the closest point. This was devastatingly slow, as you might expect. I then implemented the algorithm described by Revelles et al which is a nice algorithm, but limited (all octants must be half the size of their parents, for instance; this means it can work only on true octrees and not “loose octrees” or k-d trees) and fairly complicated.

Today I had a random thought while doing day-job work: what if I treat the octree divisions as splitting planes and essentially do a binary search? By knowing which plane my ray is closest to at a given step, I know which nodes I need to search. To my surprise – and slight horror, because it’s never a good sign – this worked perfectly on the first try and is crazy fast. I can’t see any performance drop whatsoever when I enable physics in my game anymore, which is always a good sign.

Given how simple this is, I find it exceptionally hard to believe that this is actually a novel algorithm, but I haven’t found this documented anywhere. Anyone have an idea?

The algorithm #

Here’s a simple description of the algorithm as a recursive process. It’s possible to make this iterative, of course.

If ray origin is outside the octree, either find its intersection point and begin there or (easier but slower) fall back to a naive traversal.

For each step (maximum of three per node): Determine which side of the X, Y, and Z planes the ray origin lies on Generate the child octant number by ORing together the corresponding bits for each plane; +X is 4, +Y is 2, +Z is 1 Test the corresponding child octant (run these steps anew) If hit, return – this is always the nearest intersection Find the nearest plane intersection If no intersection No hit in this octant Set the ray origin to the previously discovered intersection point Flip the side value corresponding to the plane that was intersected If it’s the X plane, for instance, you would XOR the child octant number by 4 Repeat



The code #

Here’s my implementation of this in C# along with the interfaces it expects.

// Planes are in normal+distance form class Plane { Vector3 Normal; float Distance; float RayDistance(Vector3 origin, Vector3 direction); // Gives the distance along the ray (+direction) where intersection with the plane occurs } interface AABB { Vector3 Center; // Center point of the bounding box Plane[] MidPlanes = new[] { // These are the planes which split a bounding box in half in each direction new Plane(Vector3.UnitZ, Center.Z), new Plane(Vector3.UnitY, Center.Y), new Plane(Vector3.UnitX, Center.X) }; bool Contains(Vector3 point); // True if the point lies within the bounding box } class Octree { bool Empty; // Whether this octant contains nothing Octree[] Nodes; // This is an array of 8 child octants, laid out where 0 bits correspond to the negative side of the plane, 1 bits are positive // X == 4, Y == 2, Z == 1 AABB BoundingBox; // Axis-aligned bounding box for this octant Mesh Leaf; // Triangles for a given leaf octant (Triangle, Vector3)? RayIntersection(Vector3 _origin, Vector3 direction) { if(Empty) return null; var origin = _origin; if(Leaf != null) return Leaf.RayIntersection(origin, direction); var planes = BoundingBox.MidPlanes; var side = ( X: Vector3.Dot(origin, planes[2].Normal) - planes[2].Distance >= 0, Y: Vector3.Dot(origin, planes[1].Normal) - planes[1].Distance >= 0, Z: Vector3.Dot(origin, planes[0].Normal) - planes[0].Distance >= 0 ); var xDist = side.X == direction.X < 0 ? planes[2].RayDistance(origin, direction) : float.PositiveInfinity; var yDist = side.Y == direction.Y < 0 ? planes[1].RayDistance(origin, direction) : float.PositiveInfinity; var zDist = side.Z == direction.Z < 0 ? planes[0].RayDistance(origin, direction) : float.PositiveInfinity; for(var i = 0; i < 3; ++i) { var idx = (side.Z ? 1 : 0) | (side.Y ? 2 : 0) | (side.X ? 4 : 0); var ret = Nodes[idx].RayIntersection(origin, direction); if(ret != null) return ret; var minDist = MathF.Min(MathF.Min(xDist, yDist), zDist); if(float.IsInfinity(minDist)) return null; origin = _origin + direction * minDist; if(!BoundingBox.Contains(origin)) return null; if(minDist == xDist) { side.X = !side.X; xDist = float.PositiveInfinity; } else if(minDist == yDist) { side.Y = !side.Y; yDist = float.PositiveInfinity; } else if(minDist == zDist) { side.Z = !side.Z; zDist = float.PositiveInfinity; } } return null; } }

169 Kudos