Linear-Time Voxel Walking for Octrees

Jim Arvo - 1988

Here is a new way to attack the problem of "voxel walking" in octrees (at least I think it's new). By voxel walking I mean identifying the successive voxels along the path of a ray. This is more for theoretical interest than anything else, though the algorithm described below may actually be practical in some situations. I make no claims about the practicality, however, and stick to theoretical time complexity for the most part.

For this discussion assume that we have recursively subdivided a cubical volume of space into a collection of equal-sized voxels using a BSP tree -- i.e. each level imposes a single axis-orthogonal partitioning plane. The algorithm is much easier to describe using BSP trees, and from the point of view of computational complexity, there is basically no difference between BSP trees and octrees. Also, assuming that the subdivision has been carried out to uniform depth throughout simplifies the discussion, but is by no means a prerequisite. This would defeat the whole purpose because we all know how to efficiently walk the voxels along a ray in the context of uniform subdivision -- i.e. use a 3DDDA.

Assuming that the leaf nodes form an NxNxN array of voxels, any given ray will pierce at most O(N) voxels. The actual bound is something like 3N, but the point is that it's linear in N. Now, suppose that we use a "re-traversal" technique to move from one voxel to the next along the ray. That is, we create a point that is guaranteed to lie within the next voxel and then traverse the hierarchy from the root node until we find the leaf node, or voxel, containing this point. This requires O( log N ) operations. In real life this is quite insignificant, but here we are talking about the actual time complexity. Therefore, in the worst case situation of following a ray through O( N ) voxels, the "re-traversal" scheme requires O( N log N ) operations just to do the "voxel walking." Assuming that there is an upper bound on the number of objects in any voxel (i.e. independent of N), this is also the worst case time complexity for intersecting a ray with the environment.

In this note I propose a new "voxel walking" algorithm for octrees (call it the "partitioning" algorithm for now) which has a worst case time complexity of O( N ) under the conditions outlined above. In the best case of finding a hit "right away" (after O(1) voxels), both "re-traversal" and "partitioning" have a time complexity of O( log N ). That is:

                    BEST CASE: O(1) voxels    WORST CASE: O(N) voxels
                    searched before a hit.    searched before a hit.
                  +---------------------------------------------------+
                  |                                                   |
  Re-traveral     |     O( log N )               O( N Log N )         |
                  |                                                   |
  Partitioning    |     O( log N )               O( N )               |
                  |                                                   |
                  +---------------------------------------------------+

The new algorithm proceeds by recursively partitioning the ray into little line segments which intersect the leaf voxels. The top-down nature of the recursive search ensures that partition nodes are only considered ONCE PER RAY. In addition, the voxels will be visited in the correct order, thereby retaining the O( log N ) best case behavior.

Below is a pseudo code description of the "partitioning" algorithm. It is the routine for intersecting a ray with an environment which has been subdivided using a BSP tree. Little things like checking to make sure the intersection is within the appropriate interval have been omitted. The input arguments to this routine are:

Node - A BSP tree node which comes in two flavors -- a partition node or a leaf node. A partition node defines a plane and points to two child nodes which further partition the "positive" and "negative" half-spaces. A leaf node points to a list of candidate objects.
P - The ray origin. Actually, think of this as an endpoint of a 3D line segment, since we are constraining the "ray" to be of finite length.
D - A unit vector indicating the ray direction.
len - The length of the "ray" -- or, more appropriately, the line segment. This is measured from the origin, P, along the direction vector, D.

The function "Intersect" is initially passed the root node of the BSP tree, the origin and direction of the ray, and a length, "len", indicating the maximum distance to intersections which are to be considered. This starts out being the distance to the far end of the original bounding cube.

FUNCTION Intersect( Node, P, D, len ) RETURNING "results of intersection"

    IF Node is NIL THEN RETURN( "no intersection" )

    IF Node is a leaf THEN BEGIN

        intersect ray (P,D) with objects in the candidate list

        RETURN( "the closest resulting intersection" )

        END IF

    dist := signed distance along ray (P,D) to plane defined by Node

    near := child of Node in half-space which contains P

    IF 0 < dist < len THEN BEGIN  /* the interval intersects the plane */

        hit_data := Intersect( near, P, D, dist )

        IF hit_data <> "no intersection" THEN RETURN( hit_data )

        Q := P + dist * D   /* 3D coords of point of intersection */

        far := child of Node in half-space which does NOT contain P

        RETURN( Intersect( far, Q, D, len - dist ) )

        END IF

    ELSE RETURN( Intersect( near, P, D, len ) ) 

    END

As the BSP tree is traversed, the line segments are chopped up by the partitioning nodes. The "shrinking" of the line segments is critical to ensure that only relevent branches of the tree will be traversed.

The actual encodings of the intersection data, the partitioning planes, and the nodes of the tree are all irrelevant to this discussion. These are "constant time" details. Granted, they become exceedingly important when considering whether the algorithm is really practial. Let's save this for later.

A naive (and incorrect) proof of the claim that the time complexity of this algorithm is O(N) would go something like this:

The voxel walking that we perform on behalf of a single ray is really just a search of a binary tree with voxels at the leaves. Since each node is only processed once, and since a binary tree with k leaves has k - 1 internal nodes, the total number of nodes which are processed in the entire operation must be of the same order as the number of leaves. We know that there are O( N ) leaves. Therefore, the time complexity is O( N ). But wait! The tree that we search is not truly binary since many of the internal nodes have one NIL branch. This happens when we discover that the entire current line segment is on one side of a partitioning plane and we prune off the branch on the other side. This is essential because there are really N**3 leaves and we need to discard branches leading to all but O( N ) of them. Thus, k leaves does not imply that there are only k - 1 internal nodes. The quention is, "Can there be more than O( k ) internal nodes?".

Suppose we were to pick N random voxels from the N**3 possible choices, then walk up the BSP tree marking all the nodes in the tree which eventually lead to these N leaves. Let's call this the subtree "generated" by the original N voxels. Clearly this is a tree and it's uniquely determined by the leaves. A very simple argument shows that the generated subtree can have as many as 2 * ( N - 1 ) * log N nodes. This puts us right back where we started from, with a time complexity of O( N log N ), even if we visit these nodes only once. This makes sense, because the "re-traversal" method, which is also O( N log N ), treats the nodes as though they were unrelated. That is, it does not take advantage of the fact that paths leading to neighboring voxels are likely to be almost identical, diverging only very near the leaves. Therefore, if the "partitioning" scheme really does visit only O( N ) nodes, it does so because the voxels along a ray are far from random. It must implicitly take advantage of the fact that the voxels are much more likely to be brothers than distant cousins.

This is in fact the case. To prove it I found that all I needed to assume about the voxels was connectedness -- provided I made some assumptions about the "niceness" of the BSP tree. To give a careful proof of this is very tedious, so I'll just outline the strategy (which I *think* is correct). But first let's define a couple of convenient terms:

Two voxels are "connected" (actually "26-connected") if they meet at a face, an edge, or a corner. We will say that a collection of voxels is connected if there is a path of connected voxels between any two of them.
A "regular" BSP tree is one in which each axis-orthogonal partition divides the parent volume in half, and the partitions cycle: X, Y, Z, X, Y, Z, etc. (Actually, we can weaken both of these requirements considerably and still make the proof work. If we're dealing with "standard" octrees, the regularity is automatic.)

Here is a sequence of little theorems which leads to the main result:

A ray pierces O(N) voxels.
The voxels pierced by a ray form a connected set.
Given a collection of voxels defined by a "regular" BSP tree, any connected subset of K voxels generates a unique subtree with O( K ) nodes.
The "partitioning" algorithm visits exactly the nodes of the subtree generated by the voxels pierced by a ray. Furthermore, each of these nodes is visited exaclty once per ray.
The "partitioning" algorithm has a worst case complexity of O( N ) for walking the voxels pierced by a ray.

Theorems 1 and 2 are trivial. With the exception of the "uniqueness" part, theorem 3 is a little tricky to prove. I found that if I completely removed either of the "regularity" properties of the BSP tree (as opposed to just weakening them), I could construct a counterexample. I think that theorem 3 is true as stated, but I don't like my "proof" yet. I'm looking for an easy and intuitive proof. Theorem 4 is not hard to prove at all. All the facts become fairly clear if you see what the algorithm is doing. Finally, theorem 5, the main result, follows immediately from theorems 1 through 4.

Some Practical Matters

Since log N is typically going to be very small -- bounded by 10, say -- this whole discussion may be purely academic. However, just for the heck of it, I'll mention some things which could make this a (maybe) competative algorithm for real-life situations (in as much as ray tracing can ever be considered to be "real life").

First of all, it would probably be advisable to avoid recursive procedure calls in the "inner loop" of a voxel walker. This means maintaining an explicit stack. At the very least one should "longjump" out of the recursion once an intersection is found.

The calculation of "dist" is very simple for axis-orthogonal planes, consisting of a subtract and a multiply (assuming that the reciprocals of the direction components are computed once up front, before the recursion begins).

A nice thing which falls out for free is that arbitrary partitioning planes can be used if desired. The only penalty is a more costly distance calculation. The rest of the algorithm works without modification. There may be some situations in which this extra cost is justified.

Sigh. This turned out to be much longer than I had planned...

Here is a slightly improved version of the algorithm in my previous mail. It turns out that you never need to explicitly compute the points of intersection with the partitioning planes. This makes it a little more attractive.

-- Jim

FUNCTION BSP_Intersect( Ray, Node, min, max ) RETURNING "intersection results"

BEGIN

    IF Node is NIL THEN RETURN( "no intersection" )

    IF Node is a leaf THEN BEGIN  /* Do the real intersection checking */

        intersect Ray with each object in the candidate 
        list discarding those farther away than "max." 

        RETURN( "the closest resulting intersection" )

        END IF

    dist := signed distance along Ray to plane defined by Node

    near := child of Node for half-space containing the origin of Ray

    far  := the "other" child of Node -- i.e. not equal to near.

    IF dist > max OR dist < 0 THEN  /* Whole interval is on near side. */

        RETURN( BSP_Intersect( Ray, near, min, max ) )

    ELSE IF dist < min THEN  /* Whole interval is on far side. */

        RETURN( BSP_Intersect( Ray, far , min, max ) )

    ELSE BEGIN  /* the interval intersects the plane */

        hit_data := BSP_Intersect( Ray, near, min, dist ) /* Test near side */

        IF hit_data indicates that there was a hit THEN RETURN( hit_data )

        RETURN( BSP_Intersect( Ray, far, dist, max ) )  /* Test far side. */

        END IF

    END

It should be noted that the recursive traversal algorithm introduced independently by Frederik Jansen in 1986 is a different implementation of the same idea.