Ray casting refers to tracing rays from a point to find the closest object that the ray will intersect.

A simple approach to this is to enumerate over candidate objects to find their intersections and then the object with the closest intersection point to the camera/user is the winner.

There are well-known, and often fairly obvious, algorithms for many shapes. Two important ones are for Axis-aligned bounding boxes (AABBs) and triangles.

Voxel data can be ray cast with either of these. Each voxel is an AABB, and is typically rendered using two triangles for each face. Unfortunately, this will be inefficient, because we will need to consider every voxel for every ray.

It is more efficient to use an approach that takes advantage of the structure of the voxel grid. We can “walk” the ray through the grid, stepping only far enough to touch the next grid cell.

Inspiration

This code was indirectly inspired by this 1987 paper by Amanatides and Woo. I only discovered this paper while writing this article though!

I have seen many implementations of this floating around, but rather than blindly implement, I would rather understand how something works. This then represents a summary of my understanding.

How this works

The path of the ray from a position can be expressed as \( pos^t = pos^0 + t \cdot ray\). We want to know the value of t when any component of \( pos^t \) has an integer value.

Looking at the same equation broken out into components

$$ \begin{align*} pos_x^t = pos_x^0 + t \cdot ray_x \\ pos_y^t = pos_y^0 + t \cdot ray_y \\ pos_z^t = pos_z^0 + t \cdot ray_z \\ \end{align*} $$

For now, let's just consider the x-component. and rearrange for t

$$ t = \left(x - x_0\right) / r_x $$

If we plug in the values \(a\) and \( a + 1 \)

$$ \begin{align*} t_a &= \left(a - x_0\right)/r_x \\ t_{a+1} &= \left(a + 1 - x_0\right)/r_x \\ &= \left(a - x_0\right)/r_x + 1/r_x \\ t_{a+1} &= t_a + 1/r_x \end{align*} $$

From this we can see that the increase in t needed to make a given component c increase by 1 is \(1/r_c\).

This allows us to compute \( t_\Delta = 1/ray \) or, more correctly, because the equations above assumed that we want to increase the value of a, but because the ray vector may have negative components we need the absolute value:

$$ t_\Delta = \left|1/ray\right| $$

\(t_\Delta\) is represented in the code below as t_delta.

This is fine for moving between integer values of the components, but our starting point is unlikely to be integer-aligned. It turns out not to matter that much. If we choose integer coordinates that are less than 1 away from the starting point — Here I'm using the floor for this, but in the code below I select the floor and ceiling values for each component based on the sign of the corresponding ray component — and compute the t values needed to intersect each of the component planes, then those t values will be close enough.

$$ \begin{align*} cell &= \lfloor pos \rfloor \\ t_{\mathrm{candidate}} &= \left(cell - pos\right) / ray \end{align*} $$

The algorithm is then quite simple. For each iteration:

  1. Take t as the smallest component from \(t_\mathrm{candidate}\)
  2. Add the corresponding component of \(t_\Delta\) to the selected component of \(t_\mathrm{candidate}\)
  3. IF \(t < 0\) THEN return to step 1
  4. Return the position with the current value for t. \(\left(pos + t \cdot ray\right)\)

Code

For my purposes, I did this as a zig iterator, and for convenience, define a Vec3f type for a vector of 3 f32.

Zig is convenient for code like this because it has native support for vectors and vector operations. However, I'm specifically using zig because it is the language that I'm using for a larger project.

const Vec3f = @Vector(3, f32);
const Self = @This();

The iterator struct has the following members

t_delta_x: Vec3f,
t_delta_y: Vec3f,
t_delta_z: Vec3f,
t_candidate: Vec3f,
pos: Vec3f,
ray: Vec3f,
first: bool = true,
t_prev: f32 = 0,

const ZERO: Vec3f = @splat(0);
const ONE: Vec3f = @splat(1);

The iterator is initialized with the position from which the ray will be cast and the direction of the ray.

This code has the following assumptions:

  1. Voxel space is from (0,0,0) -> (16, 16, 16)
  2. The start position is within the voxel space
  3. The ray vector is normalized

None of these are strictly necessary, the ray vector could be normalized in the initializer, and the restrictions on the position being within the voxel space are arbitrary.

To use this in a more general purpose setting you would typically find the point at which the ray enters the AABB of the voxel space and use that — probably nudged into the voxel space to avoid fussy edge cases — as the start point.

The iterator also doesn't need to terminate. That is an implementation dependent decision, the calling code can just as easily make the decision to stop iterating.

pub fn init(pos: Vec3f, ray: Vec3f) Self {

As described above, t_delta is a vector of the number of times the ray needs to be added to increment (or decrement) each axis by one. I.e. the value is t_delta[0] will be the t required to increment the 0th component of the position by 1. If the corresponding component of the ray vector is negative, this will decrement the component by 1.

Because we only ever use t_delta in a component-wise fashion, we store it in separate components — this makes the iteration code a little cleaner.

  const t_delta = @abs(ONE / ray);
  return .{
    .t_delta_x = .{t_delta[0], 0, 0},
    .t_delta_y = .{0, t_delta[1], 0},
    .t_delta_z = .{0, 0, t_delta[2]},

t_candidate is a vector of candidate values of t, each component is the amount required to advance the current position to the next intersection for that component. This is initialized based on either the floor or the ceiling of the current position. I use @select to choose which value to use based on the sign of ray.

@select creates a vector by selecting component values from one of two input vectors (@floor(pos) and @ceil(pos)) based on the respective component value of a predicate vector (in this case ray <= ZERO)

    .t_candidate = (@select(f32, ray <= ZERO, @floor(pos), @ceil(pos)) - pos) / ray,

and we also store the start position and ray vector, which will be used to calculate the position for a given t.

    .pos = pos,
    .ray = ray,
  };
}

This is the structure used for the result returned by each call into the iterator. t is the distance traveled by the ray from the start position, cell is the integer voxel cell that is being entered by this point, and pos is the position of this result.

pub const Intersection = struct {
  t: f32,
  cell: Vec3f,
  pos: Vec3f,
};

The first call into the iterator always returns the starting point.

pub fn next(self: *Self) ?Intersection {
  if (self.first) {
    self.first = false;
    self.t_prev = 0;
    return .{
      .t = 0,
      .cell = @floor(self.pos),
      .pos = self.pos,
    };
  }

For each subsequent call we determine which of the t_candidate components is smallest — there are clever tricks to do this with vectors, but I don't think they really help here, I'll discuss why below.

We assign the chosen components candidate value of t to ct and advance that component by the corresponding component of t_delta.

  while (true) {
    var ct: f32 = undefined;

    if (isFiniteAndLte(self.t_candidate[0], self.t_candidate[1]) and
        isFiniteAndLte(self.t_candidate[0], self.t_candidate[2])) {
      ct = self.t_candidate[0];
      self.t_candidate += self.t_delta_x;
    } else if (isFiniteAndLte(self.t_candidate[1], self.t_candidate[2])) {
      ct = self.t_candidate[1];
      self.t_candidate += self.t_delta_y;
    } else {
      ct = self.t_candidate[2];
      self.t_candidate += self.t_delta_z;
    }

    if (!std.math.isFinite(ct)) {
      unreachable;
    }

If the current value of ct has almost exactly the same as the previous, then we go around the loop again, otherwise we may return the same result for two iterations.

There are some edge cases where this might cause the algorithm to skip over a corner or edge of a voxel. A different approach would be to occasionally allow duplicate values from the iteration. This is very much an implementation choice.

    if (std.math.approxEqAbs(f32, 0, ct - self.t_prev, 1e-5)) {
      continue;
    }
    self.t_prev = ct;

We now compute the voxel cell for this iteration by multiplying the ray by the candidate t value (plus a little nudge) and taking the floor of that value.

The nudge is used to ensure that the correct cell is selected when the ray vector has a negative component, otherwise the adjacent cell will be selected.

The value of cell is then checked to ensure that it is within the voxel space bounds, and, if so, the iteration is stopped by returning null.

@reduce transforms a vector into a scalar by performing a sequential reduction using the specified operator. Used here to find the min and max values of the cell vector.

@splat creates a vector from a scalar, needed here to multiply each component of the ray vector by ct + 1e-2

    const cell = @floor(self.pos + self.ray * @as(Vec3f, @splat(ct + 1e-2)));

    if (@reduce(.Min, cell) < 0 or @reduce(.Max, cell) >= 16) {
      return null;
    }

This iteration result is then returned.

    return .{
      .t = ct,
      .cell = cell,
      .pos = self.pos + self.ray * @as(Vec3f, @splat(ct)),
    };
  }
}

This helper function is used when finding the smallest component of the t_candidate vector. This returns true IF:

  • a is finite AND
    • EITHER \( a \leq b\) OR
    • b is NOT finite.

This is needed because the candidate components will be inf if the corresponding component of ray is zero.

This is also why some of the vector tricks that can be used for this are hard to use because they don't necessarily work correctly in these edge cases.

inline fn isFiniteAndLte(a: f32, b:f32) bool {
  return std.math.isFinite(a) and (a <= b or !std.math.isFinite(b));
}

An alternative solution is to create a safe ray vector by replacing any zero component in the ray with a very small value. e.g.:

  const EPSILON: Vec3f = @splat(1e-10);
  const safe_ray = @select(f32, ray == ZERO, EPSILON, ray);

This safe_ray is then used to calculate t_candidate and t_delta, but the regular value of ray is used elsewhere in the algorithm.

Tests

Testing this should probably use some test vectors from an alternate implementation, but to be satisfied that it works, I assert a few invariants for each iteration.

Test helper

Those assertions are built into a single test helper function that can be called with different positions and vectors.

In order to normalize the ray, and to keep this code self contained, there's a one-liner to normalize the ray at the top of this test function.

fn testIterator(pos: Vec3f, ray: Vec3f) !void {
  const norm_ray = ray / @as(Vec3f, @splat(@sqrt(@reduce(.Add, ray * ray))));

  var it = Self.init(pos, norm_ray);

The first part of the test is to assert that the first intersection returned by the iterator represents the start point passed to the init function.

  if (it.next()) |curr| {
    try std.testing.expectEqual(0, curr.t);
    try std.testing.expectEqual(@floor(pos), curr.cell);
    try std.testing.expectEqual(pos, curr.pos);
  }

For each of the following reported intersections we check three invariants.

Firstly we validate the cell:

  • We check that the combined absolute difference in the cell components (compared with the prior cell), is an integer.
  • That between 1 and 3 components have changed.
  • That the largest absolute difference is 1

These checks verify that the cell has progressed to an immediate neighbor, including diagonal neighbors.

  var last = @floor(pos);
  var t_last: f32 = 0;

  while (it.next()) |curr| {
    const cell_diff = @abs(last - curr.cell);
    const changes = @reduce(.Add, cell_diff);
    try std.testing.expectEqual(changes, @floor(changes));
    try std.testing.expect(changes > 0);
    try std.testing.expect(changes <= 3);
    try std.testing.expectEqual(1, @reduce(.Max, cell_diff));

Then we validate that the 't' value is within expected bounds:

  • That it has increased from the previous iteration.
  • That it has not increased by more than \( \sqrt{3} \) (the length of a cell diagonal)

    const t_diff = curr.t - t_last;
    try std.testing.expect(t_diff > 0);
    try std.testing.expect(t_diff < @sqrt(3.0) or
      std.math.approxEqRel(f32, t_diff, @sqrt(3.0), 1e-5));

Lastly, that the intersection position is on the face of a voxel.

This is verified by checking that at least one component is close to an integer, which can be verified by taking the fractional part of the components and checking that the smallest of these fractions is sufficiently close to zero, or, failing that, that the largest is sufficiently close to one — we do this by subtracting 0.5 from the fractional part and taking the absolute value. This way the if the largest component is sufficiently close to 0.5 it was either very close to 0.0 or to 1.0.


    const fract = curr.pos - @floor(curr.pos) - @as(Vec3f, @splat(0.5));
    try std.testing.expectApproxEqRel(0.5, @reduce(.Max, @abs(fract)), 1e-5);

    last = curr.cell;
    t_last = curr.t;
  }
}

Test cases

test "simple case" {
  try testIterator(.{10.3, 11.4, 12.5}, .{1, 2, 3});
}

test "diagonal case" {
  try testIterator(.{10, 11, 12}, .{1, 1, 1});
}

test "2d case" {
  try testIterator(.{10.3, 11.4, 12.5}, .{1, -2, 0});
}

test "1d case" {
  try testIterator(.{10.3, 11.4, 12.5}, .{0, -1, 0});
}

Code

All put together, this gives us voxel_iterator.zig

const std = @import("std");

<<iterator>>
<<tests>>

I can then take the tangled output from this article and run the tests (using zig 0.12.0)

$ zig test voxel_iterator.zig
All 4 tests passed.
$ ␣