A common approach to ray-tracing voxels is a method described in the 2015 paper by Thjis van Wingerden titled “Real-time Ray tracing and Editing of Large Voxel Scenes”.

Ignoring details like how the voxels are stored, how LoD is handled or how the data to be streamed is determined, the paper describes two simple data structres, namely “brickgrid"s and “brickmap"s and the traversal thereof.

This post is supposed to be a crash course on said structures and the traversal, and not on the afformentioned specifics which would have to be tailored on the type of program you, the reader, want to write.

The so-called “3D-DDA”#

The goal of this algorithm is to determine the next point on a grid corresponding to the closest “box” from a given point and on a given direction.

It doesn’t really matter to the explanation but I’d like to mention that I believe the name to be a misnomer in this context as the lines drawn by the titular DDA algorithm do not match the lines that need to be drawn to trace rays.

It can be seen that the box to be chosen has to be one of the four boxes directly adjacent to the one within which the starting point resides. This is because, to reach one of the corner boxes, we have to pass through the adjacent ones anyway:

And if we assume that the ray is north-eastwardly, the candidate boxes can be further whittled down:

At this point, we need to recall the implicit function for a ray: $\mathbf{x} = \mathbf{o} + t\hat{\mathbf{d}}$ where a point in space along the ray (denoted $\mathbf{x}$) is characterised by its distance $t$ from the starting point of the ray (denoted $\mathbf{o}$) in terms of the direction $\hat{\mathbf{d}}$.

And with that, it can be seen that the problem is finding the closer plane to the ray in terms of $t$:

Recall that for a plane described by an origin $\mathbf{p}_0$ and a normal at that point $\hat{\mathbf{n}}$ can be described by the implicit function $(\mathbf{x} - \mathbf{p}\_0) \cdot \hat{\mathbf{n}} = 0$. We can then use the implicit function for a ray $\mathbf{x} = \mathbf{o} + t\hat{\mathbf{d}}$ to solve for $\mathbf{x}$:

$$(\mathbf{o} + t\hat{\mathbf{d}} - \mathbf{p}_0) \cdot \hat{\mathbf{n}} = 0$$ $$\mathbf{o} \cdot \hat{\mathbf{n}} + t\hat{\mathbf{d}} \cdot \hat{\mathbf{n}} - \mathbf{p}_0 \cdot \hat{\mathbf{n}} = 0$$ $$t(\hat{\mathbf{d}} \cdot \hat{\mathbf{n}}) = \mathbf{p}_0 \cdot \hat{\mathbf{n}} - \mathbf{o} \cdot \hat{\mathbf{n}}$$ $$t = {(\mathbf{p}_0 - \mathbf{o}) \cdot \hat{\mathbf{n}} \over \hat{\mathbf{d}} \cdot \hat{\mathbf{n}}}$$

And when $\hat{\mathbf{n}}$ has a single component, computing a dot product between it and another vector amounts to ‘selecting’ an axis of said vector. So when we want to know the $t$ value for the y-axis-plane (where the normal is $\begin{Bmatrix} -1 \\ 0 \\ 0 \end{Bmatrix}$) in the image above, we would do:

$$t = {-(\mathbf{p}_0 - \mathbf{o})_x \over -\hat{\mathbf{d}}_x}$$

This result actually has a geometric interpretation as well: In determining the $t$ value of a ray to a plane describing an axis, we can project down the $\hat{\mathbf{d}}$ to a plane orthogonal to it, and it would make no difference to the $t$ despite the reduced magnitude:

And it can be seen that this computation can be done on all three axes in one fell swoop time by an elementwise division of the vector $\mathbf{p}_0 - \mathbf{o}$ (represented by $\Delta$) by the vector $\hat{\mathbf{d}}$:

And after said division, the smallest axis of the resultant vector determines the iteration direction and the value therein is our $t$.

The so-called 3D-DDA algorithm is a repeated application of these steps.

“Why not just DDA it up?”#

Consider the following:

Now imagine this scenario in a 4096x4096x4096 grid.

And by “scenario”, haha, well. Let’s just say. A worst case scenario of $n\sqrt{2}$ iterations for an nxn grid.

So what’s the solution? While there are many solutions (namely distance fields, SVOs, and so on and so forth), the one proposed in the paper is much simpler: Store voxels in groups of 8x8x8 called “brickmaps” and don’t bother with storing the empty ones at all. Have an additional and infinite layer on top called the “brickgrid” that references these groups. Traverse the brickgrid first and only if the position in the brickgrid contains a brickmap, traverse the brickmap.

As is suggested in the image, the brickgrid has no particular size. It is the brickmaps that must be statically sized.

My Gripes#

I feel that the paper stops too short. But oh, what’s that? They acknowledge it?

Oh yeah, they do.

I’ll be talking about my take on the paper and kinds of fat-trees you can add with relative ease to this scheme on the next post, whenever that might come.