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:
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
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:
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:
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...
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