A Stupidly Simple, Fast Octree Traversal Algorithm for Ray Intersection

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;
    }
}