Rendering

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 \in \R^n$ as $R(t) = \overrightarrow{u} + t\overrightarrow{v}$ parameterized on $t \in [0, inf)$ with $\overrightarrow{u}$ being the origin position and $\overrightarrow{v}$ being a vector representing direction of movement in $\mathbb{R}^n$ . Looking first at the x axis, we can easily find the parameterizations $t_{x_i}$ of $R(t)$ where $R(t_{x_i})_x \equiv 0 \pmod 1$ - that is, where $R$ crosses from one voxel to another along the x-axis. This is given by:

t_{x_i}= t_{x_{i-1}} + \Delta_x = t_{x_{i-1}} + \left\lvert \frac{1}{\overrightarrow{v}_x} \right\rvert

t_{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 $n$ . We can then find the list of parameterizations $t_i$ of $R(t)$ where $R(t_i)_n \equiv 0 \pmod 1$ , that is, where $R$ crosses from one voxel to another along any axis by merging all of $t_{x_i}, t_{y_i}, ... t_{n_i}$ into a single list. This list can be iterated through by finding the next smallest value $t_{n_i}$ in that set:

t_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 $t_i$ until we either find a voxel at $R(t_i)$ or escape the bounds of the world.

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 $\overrightarrow{tmax} = \langle t_{x_0}, t_{y_0}, t_{z_0} \rangle$

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 $t_{n_i}= t_{n_{i-1}} + \Delta_n$ , we can define $\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(!isVoxelInABlock && i < limitVoxels) {
    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) {
      voxel.y += step.y;
      tmax.y += delta.y;
    } else {
      voxel.z += step.z;
      tmax.z += delta.z;
    }
    isVoxelInABlock = checkVoxelInBlock(voxel);
    if(voxel.z > 20.0f) {
      return vec4(1);
    }
    i++;
  }