# Voxel Ray Casting in Zig

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:

- Take
*t*as the smallest component from \(t_\mathrm{candidate}\) - Add the corresponding component of \(t_\Delta\) to the selected component of \(t_\mathrm{candidate}\)
- IF \(t < 0\) THEN return to step 1
- 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:

- Voxel space is from (0,0,0) -> (16, 16, 16)
- The start position is within the voxel space
- 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.
$ ␣
```