Octree raycasting/raytracing - best ray/leaf intersection without recursion

Question

Could anyone provide a short & sweet explanation (or suggest a good tutorial) on how to cast a ray against a voxel octree without recursion?

I have a complex model baked into an octree, and I need to find the best/closest leaf that intersects a ray. A standard drill-down iterative tree walk:

1. Grab the root node
2. Check for intersection
3. No? Exit
4. Yes? Find child that intersects the ray that is closest to the ray's origin
5. Loop until I reach a leaf or exit the tree

Always returns a leaf, but in instances where the tree stores, say, terrain, the closest node to the ray's origin doesn't necessarily contain the leaf that's the best match. This isn't suprising - taller objects in farther nodes won't get tested using this approach.

I can do this recursively by finding all of the intersecting leaves in the tree, sorting by distance and picking the closest one to the ray's position. However, this is slow and requires recursion.

I've read a little about using the Bresenham line algorithm to walk the tree, which seems to require that each node contain pointers to adjacent neighbors, but I'm unclear on how to implement this in a useful way.

Any suggestions? I can fake a stack in HLSL using a fixed-length array or a struct with an element for each potential stack entry, but the memory requirements for that can become crippling with a sufficiently large tree.

Help.

Answer 1

I've managed to get this mostly working using a 3D DDA algorithm and neighbor node pointers.

I'm still working out a few bugs, but here's a C# version that appears to work. This one stops when it reaches the first leaf, but that's not entirely necessary.

/// <param name="ray"></param>
public OctreeNode DDATraverse(Ray ray)
{
    float tmin;
    float tmax;


    /// make sure the ray hits the bounding box of the root octree node
    if (!RayCasting.HitsBox(ray, root.BoundingBox.Min, root.BoundingBox.Max, out tmin, out tmax))
        return null;


    /// move the ray position to the point of intersection with the bounding volume.
    ray.Position += ray.Direction * MathHelper.Min(tmin, tmax);// intersectionDistance.Value;

    /// get integer cell coordinates for the given position
    ///     leafSize is a Vector3 containing the dimensions of a leaf node in world-space coordinates
    ///     cellCount is a Vector3 containng the number of cells in each direction, or the size of the tree root divided by leafSize.

    var cell = Vector3.Min(((ray.Position - boundingBox.Min) / leafSize).Truncate(), cellCount - Vector3.One);

    /// get the Vector3 where of the intersection point relative to the tree root.
    var pos = ray.Position - boundingBox.Min;

    /// get the bounds of the starting cell - leaf size offset by "pos"
    var cellBounds = GetCellBounds(cell);

    /// calculate initial t values for each axis based on the sign of the ray.
    /// See any good 3D DDA tutorial for an explanation of t, but it basically tells us the 
    /// distance we have to move from ray.Position along ray.Direction to reach the next cell boundary
    /// This calculates t values for both positive and negative ray directions.
    var tMaxNeg = (cellBounds.Min - ray.Position) / ray.Direction;
    var tMaxPos = (cellBounds.Max - ray.Position) / ray.Direction;

    /// calculate t values within the cell along the ray direction.
    /// This may be buggy as it seems odd to mix and match ray directions
    var tMax = new Vector3(
        ray.Direction.X < 0 ? tMaxNeg.X : tMaxPos.X
        ,
        ray.Direction.Y < 0 ? tMaxNeg.Y : tMaxPos.Y
        ,
        ray.Direction.Z < 0 ? tMaxNeg.Z : tMaxPos.Z
        );

    /// get cell coordinate step directions
    /// .Sign() is an extension method that returns a Vector3 with each component set to +/- 1
    var step = ray.Direction.Sign();

    /// calculate distance along the ray direction to move to advance from one cell boundary 
    /// to the next on each axis. Assumes ray.Direction is normalized.
    /// Takes the absolute value of each ray component since this value is in units along the
    /// ray direction, which makes sure the sign is correct.
    var tDelta = (leafSize / ray.Direction).Abs();

    /// neighbor node indices to use when exiting cells
    /// GridDirection.East = Vector3.Right
    /// GridDirection.West = Vector3.Left
    /// GridDirection.North = Vector3.Forward
    /// GridDirection.South = Vector4.Back
    /// GridDirection.Up = Vector3.Up
    /// GridDirection.Down = Vector3.Down
    var neighborDirections = new[] { 
        (step.X < 0) ? GridDirection.West : GridDirection.East
        ,
        (step.Y < 0) ? GridDirection.Down : GridDirection.Up
        ,
        (step.Z < 0) ? GridDirection.North : GridDirection.South
    };



    OctreeNode node=root;


    /// step across the bounding volume, generating a marker entity at each
    /// cell that we touch. Extension methods GreaterThanOrEEqual and LessThan
    /// ensure that we stay within the bounding volume.
    while (node!=null)
    {
        /// if the current node isn't a leaf, find one.
        /// this version exits when it encounters the first leaf.
        if (!node.Leaf)
            for (var i = 0; i < OctreeNode.ChildCount; i++)
            {
                var child = node.Children[i];
                if (child != null && child.Contains(cell))
                {
                    //SetNode(ref node, child, visitedNodes);
                    node = child;
                    i = -1;

                    if (node.Leaf)
                        return node;
                }
            }

        /// index into the node's Neighbor[] array to move
        int dir = 0;

        /// This is off-the-shelf DDA.
        if (tMax.X < tMax.Y)
        {
            if (tMax.X < tMax.Z)
            {
                tMax.X += tDelta.X;
                cell.X += step.X;
                dir = 0;

            }
            else
            {
                tMax.Z += tDelta.Z;
                cell.Z += step.Z;
                dir = 2;

            }
        }
        else
        {
            if (tMax.Y < tMax.Z)
            {
                tMax.Y += tDelta.Y;
                cell.Y += step.Y;
                dir = 1;
            }
            else
            {
                tMax.Z += tDelta.Z;
                cell.Z += step.Z;
                dir = 2;
            }
        }

        /// see if the new cell coordinates fall within the current node.
        /// this is important when moving from a leaf into empty space within 
        /// the tree.
        if (!node.Contains(cell))
        {
            /// if we stepped out of this node, grab the appropriate neighbor. 
            var neighborDir = neighborDirections[dir];
            node = node.GetNeighbor(neighborDir);
        }
        else if (node.Leaf && stopAtFirstLeaf)
            return node;
    }

    return null;

}

Feel free to point out any bugs. I'll post the HLSL version if there's any demand.

Here's another version that just steps through the tree in leaf-sized steps without intersection checking. This is useful as a 3D DDA demonstration:

/// <summary>
/// draw a 3D DDA "line" in units of leaf size where the ray intersects the
/// tree's bounding volume/
/// </summary>
/// <param name="ray"></param>
public IEnumerable<Vector3> DDA(Ray ray)
{

    float tmin;
    float tmax;


    if (!RayCasting.HitsBox(ray, root.BoundingBox.Min, root.BoundingBox.Max, out tmin, out tmax))
        yield break;

    /// move the ray position to the point of intersection with the bounding volume.
    ray.Position += ray.Direction * tmin;

    /// get integer cell coordinates for the given position
    var cell = Vector3.Min(((ray.Position - boundingBox.Min) / leafSize).Truncate(), cellCount - Vector3.One);

    /// get the bounds of the starting cell.
    var cellBounds = GetCellBounds(cell);

    /// calculate initial t values for each axis based on the sign of the ray.
    var tMaxNeg = (cellBounds.Min - ray.Position) / ray.Direction;
    var tMaxPos = (cellBounds.Max - ray.Position) / ray.Direction;

    /// calculate t values within the cell along the ray direction.
    var tMax = new Vector3(
        ray.Direction.X < 0 ? tMaxNeg.X : tMaxPos.X
        ,
        ray.Direction.Y < 0 ? tMaxNeg.Y : tMaxPos.Y
        ,
        ray.Direction.Z < 0 ? tMaxNeg.Z : tMaxPos.Z
        );

    /// get cell coordinate step directions
    var step = ray.Direction.Sign();

    /// calculate distance along the ray direction to move to advance from one cell boundary 
    /// to the next on each axis. Assumes ray.Direction is normalized.
    var tDelta = (leafSize / ray.Direction).Abs();

    /// step across the bounding volume, generating a marker entity at each
    /// cell that we touch. Extension methods GreaterThanOrEEqual and LessThan
    /// ensure that we stay within the bounding volume.
    while (cell.GreaterThanOrEqual(Vector3.Zero) && cell.LessThan(cellCount))
    {
        yield return boundingBox.Min + cell * leafSize;
        ///create a cube at the given cell coordinates, and add it to the draw list.
        if (tMax.X < tMax.Y)
        {
            if (tMax.X < tMax.Z)
            {
                tMax.X += tDelta.X;
                cell.X += step.X;
            }
            else
            {
                tMax.Z += tDelta.Z;
                cell.Z += step.Z;
            }
        }
        else
        {
            if (tMax.Y < tMax.Z)
            {
                tMax.Y += tDelta.Y;
                cell.Y += step.Y;

            }
            else
            {
                tMax.Z += tDelta.Z;
                cell.Z += step.Z;
            }
        }
    }

}

And an HLSL version that just stores the tree in a Texture3D, without neighbors or any "sparseness" to the tree.

This is still buggy. The first test with hitbox() works correctly, but the ray winds up getting refracted within the tree. This looks very cool, but isn't correct.

Here's what it looks like when I stop at the root bounds, without using the DDA to traverse the tree:

/*
find which leaf, if any, the ray intersects.
Returns transparency (Color(0,0,0,0)) if no intersection was found.

TestValue is a shader constant parameter passed from the caller which is used to dynamically adjust the number of loops the shader code will execute. Useful for debugging.

intrinsics:
step(y,x) : (x >= y) ? 1 : 0


*/
float4 DDATraverse(Ray ray)
{
    float3 bounds_min = OctreeCenter-OctreeObjectSize/2;
    float3 bounds_max = OctreeCenter+OctreeObjectSize/2;
    float4 cellsPerSide = float4(trunc((bounds_max-bounds_min)/CellSize),1);
    float3 vector3_one = float3(1,1,1);

    float tmin;
    float tmax;

    if(hitbox(ray,bounds_min,bounds_max,tmin,tmax))
    {
        ray.Position+=ray.Direction*tmin;

        float4 cell = float4((ray.Position-bounds_min)/CellSize,1); 


        float3 tMaxNeg = (bounds_min-ray.Position)/ray.Direction;
        float3 tMaxPos = (bounds_max-ray.Position)/ray.Direction;

        float3 tmax = float3(
            ray.Direction.x < 0 ? tMaxNeg.x : tMaxPos.x
            ,
            ray.Direction.y < 0 ? tMaxNeg.y : tMaxPos.y
            ,
            ray.Direction.z < 0 ? tMaxNeg.z : tMaxPos.z
            );


        float3 tstep = sign(ray.Direction);
        float3 dt = abs(CellSize/ray.Direction);
        float4 texel;

        float4 color;

        for(int i=0;i<TestValue;i++)
        {
            texel=smoothstep(float4(0,0,0,0),cellsPerSide,cell);
            if (color.a < 0.9)
                color = tex3Dlod(octreeSampler,texel);

            if (tmax.x < tmax.y)
            {
                if (tmax.x < tmax.z)
                {
                    tmax.x+=dt.x;
                    cell.x+=tstep.x;
                }
                else
                {
                    tmax.z+=dt.z;
                    cell.z+=tstep.z;
                }
            }
            else
            {
                if (tmax.y < tmax.z)
                {
                    tmax.y+=dt.y;
                    cell.y+=tstep.y;
                }
                else
                {
                    tmax.z+=dt.z;
                    cell.z+=tstep.z;
                }

            }
        }

        return color;


    }
    else
        return float4(1,0,0,1);

}

update

Found a very good volume rendering tutorial!

http://graphicsrunner.blogspot.com/search?updated-max=2009-08-27T02%3A45%3A00-04%3A00&max-results=10