Rendering with Voxel Raytracing

Conventional rendering has us send a lot of triangles to the gpu, and if we want to render a few kilometers of meter cubes, that turns out on one hand to be a lot of triangles, and on another hand it is not obvious how to update the triangles. This is presumably a solved problem in graphics programing, but I have barely any knowledge on graphics programming. Sending all the triangles is too many triangles, but just sending the visible ones means that an update might be an insert in the middle of a triangle strip, and that seems difficult to handle memory-localitywise.

On the other hand, a voxel landscape does seem like something that should be relatively straightforward to render. After poking around a bit I found this nice paper on "A Fast Voxel Traversal Algorithm for Ray Tracing", which outlines an algorithm by which we can perform ray tracing on a voxel environment using minimal computation. It should allow us to send the raw voxel data to the gpu, so updates are just changes to the existing matrix. It is also pretty straightforward to implement!

The Math

Define a ray R∈RnR \in \R^n as R(t)=u→+tv→R(t) = \overrightarrow{u} + t\overrightarrow{v} parameterized on t∈[0,inf)t \in [0, inf) with u→\overrightarrow{u} being the origin position and v→\overrightarrow{v} being a vector representing direction of movement in Rn\mathbb{R}^n. Looking first at the x axis, we can easily find the parameterizations txit_{x_i} of R(t)R(t) where R(txi)x≡0(mod1)R(t_{x_i})_x \equiv 0 \pmod 1 - that is, where RR crosses from one voxel to another along the x-axis. This is given by:

txi=txi−1+Δx=txi−1+∣1v→x∣t_{x_i}= t_{x_{i-1}} + \Delta_x = t_{x_{i-1}} + \left\lvert \frac{1}{\overrightarrow{v}_x} \right\rvert tx0={⌊u→x⌋−u→x+1v→xif v→x>0⌊u→x⌋−u→xv→xif v→x<0t_{x_0} = \begin{cases} \frac{\lfloor \overrightarrow{u}_x \rfloor - \overrightarrow{u}_x + 1}{\overrightarrow{v}_x} & \text{if } \overrightarrow{v}_x > 0 \\ \frac{\lfloor \overrightarrow{u}_x \rfloor - \overrightarrow{u}_x}{\overrightarrow{v}_x} & \text{if } \overrightarrow{v}_x < 0 \end{cases}

This math works for all dimensions nn. We can then find the list of parameterizations tit_i of R(t)R(t) where R(ti)n≡0(mod1)R(t_i)_n \equiv 0 \pmod 1, that is, where RR crosses from one voxel to another along any axis by merging all of txi,tyi,...tnit_{x_i}, t_{y_i}, ... t_{n_i} into a single list. This list can be iterated through by finding the next smallest value tnit_{n_i} in that set:

ti=min⁡n,jtnjsuch that tnj>ti−1t_i = \min\limits_{n,j} t_{n_j} \text{such that } t_{n_j} > t_{i - 1}

And with that, we have all the tools necessary to raytrace through a voxel world, by iterating over tit_i until we either find a voxel at R(ti)R(t_i) or escape the bounds of the world.

Illustration of meaning of variables on a graph in R^2

The Implementation

Lets start to implement the above, defining a function vec4 raytrace(vec3 u, vec3 v) that returns the color of a pixel based on the path a ray of light would travel from a given location in a given direction.

First, we define tmax→=⟨tx0,ty0,tz0⟩\overrightarrow{tmax} = \langle t_{x_0}, t_{y_0}, t_{z_0} \rangle, which will always hold the next-smallest tnjt_{n_j} value for each of the three dimensions

float tMaxOf(
  /** position coordinate */
  float u,
  /** direction vector coordinate */
  float v,
  /** floored coordinate */
  float c
) {
  return (c - u + (v > 0.0f ? 1.0f : 0.0f)) / v;
}

vec4 raytrace(vec3 u, vec3 v) {
  vec3 voxel = floor(u);
  vec3 tmax = vec3(
    tMaxOf(u.x, direction.x, voxel.x),
    tMaxOf(u.y, direction.y, voxel.y), 
    tMaxOf(u.z, direction.z, voxel.z));
  ...
}

Then, continuing to solve for tni=tni−1+Δnt_{n_i}= t_{n_{i-1}} + \Delta_n, we can define Δn\Delta_n:

  vec3 delta = sign(v) / v;

We also need to know what direction each integer change of voxel coordinates will go in, using vec3 step = sign(v), and then have all the tools necessary to infinitely iterate over the ray:

vec4 raytrace(vec3 origin, vec3 direction) {
  vec3 voxel = floor(origin);
  vec3 step = sign(direction);
  vec3 delta = sign(direction) / direction;
  vec3 tmax = vec3(
    tMaxOf(u.x, direction.x, voxel.x),
    tMaxOf(u.y, direction.y, voxel.y), 
    tMaxOf(u.z, direction.z, voxel.z));
  
  while(true) {
    if(tmax.x <= tmax.y && tmax.x <= tmax.z) {
      // if t_{x_i} is less than t_{y_j} and t_{z_k}, move in the x direction
      voxel.x += step.x;
      tmax.x += delta.x;
    } else if(tmax.y <= tmax.z) {
       // if t_{y_i} is less than t_{z_j} move in the y direction
      voxel.y += step.y;
      tmax.y += delta.y;
    } else {
      voxel.z += step.z;
      tmax.z += delta.z;
    }
  }
}

Then we just need to return early if we start in a voxel, cap the number of iterations to prevent infinite loops, keep track of which voxel face we are traversing through at each step, and return the background color (call it white for now, although we will get into the sky coloring later) if no voxels are intersected with or a distinct color based on the face that was interacted with. In this first example, we can just distribute voxels in a linear matrix every 8 blocks - and by capping the iterations to 1000, this should be rendering about 250 cubes in each direction, or a total of 15.625 million.

TODO

References