# How do I find voxels along a ray?

I am making a voxel engine, and I want to make it possible to create and destroy voxels with the mouse. I use C++ and OpenGL, so C++ examples would be best.

So far, I have the camera position and appropriate angle as Vector3 values. I need to find voxels along a ray extending from the camera position along the angle. I tried to solve it myself, and essentially came up with a slow, 3D version of Bresenham's line. This would be okay, but it skips diagonals, meaning if I point toward the edge of a voxel, it usually misses.

The "bad" line skips over where the line crosses diagonals, while the good line hits them. I currently have the bad, but I want the good:

I have tried reading around this, but most of what I find relates to sparse octree rendering.

• You could divide the lines in your image into "spans", in this case horizontal spans. To find where a span starts, you solve for the intersection of the ray with an edge/line (horizontal lines in your case). The span ends at the next intersection, and the subsequent span starts at that point. In 3D it would be intersections with planes. This is just me postulating, though, I don't know how well it would work in practice. Mar 18, 2014 at 22:32
• The difference between your Bad and Good example is that when two hit voxels aren't adjacent (by your adjacency criteria), the Good fills the space in. Have you tried extending Bresenham thus?
– Anko
Mar 18, 2014 at 22:36
• You may be interested in A Fast Voxel Traversal Algorithm for Ray Tracing. Jun 24, 2021 at 14:08

I actually found anothe stackexchange answer that answers this. Cast ray to select block in voxel game

It turns out I was searching for the wrong stuff, and I was looking for C++ specific code. The link is in JS, but I ported it to C++. Here is my whole voxel function. In the answer, he also has some extra functions (mod, signum, and intbound), but those are easy enough to implement yourself.

void World::mouseVoxel(GLFWwindow* window, int button, int action, int mods) {

if (action == GLFW_PRESS && (button == GLFW_MOUSE_BUTTON_1 || button == GLFW_MOUSE_BUTTON_2)) { //make voxel
//calculate angle vector of the mouse
double cursorX, cursorY;
glfwGetCursorPos(window, &cursorX, &cursorY);
//std::cout << cursorX << ", " << cursorY << std::endl;
float x = (2.0f * cursorX) / 1280.0f - 1.0f;
float y = 1.0f - (2.0f * cursorY) / 720.0f;
vec4 rayClip = vec4(x, y, -1.0, 1.0);
vec4 rayEye = inverse(projectionMatrix) * rayClip;
rayEye = vec4 (rayEye.v[0], rayEye.v[1], -1.0f, 0.0f);
vec3 rayWOR = vec4(inverse(camera->getMatrix()) * rayEye);
rayWOR = normalise(rayWOR);

float range = 64.0f; //max range to check (in voxels)
vec3 camPos = camera->getPosition() * (1 / voxelSize);
float xPos = floor(camPos.v[0]);
float yPos = floor(camPos.v[1]);
float zPos = floor(camPos.v[2]);
int stepX = signum(rayWOR.v[0]);
int stepY = signum(rayWOR.v[1]);
int stepZ = signum(rayWOR.v[2]);
vec3 tMax(intbound(camPos.v[0], rayWOR.v[0]), intbound(camPos.v[1], rayWOR.v[1]), intbound(camPos.v[2], rayWOR.v[2]));
vec3 tDelta((float)stepX / rayWOR.v[0], (float)stepY / rayWOR.v[1], (float)stepZ / rayWOR.v[2]);
float faceX;
float faceY;
float faceZ;

print(camPos);

do {
if (isVoxelSolid(xPos,yPos,zPos)) {
std::cout << "boom";
if (button == GLFW_MOUSE_BUTTON_2) setVoxelType(xPos,yPos,zPos, 0, true);
else setVoxelType(xPos + faceX, yPos + faceY, zPos + faceZ, 1, true);
break;
}
if (tMax.v[0] < tMax.v[1]) {
if (tMax.v[0] < tMax.v[2]) {
if (tMax.v[0] > range) break;

xPos += stepX;
tMax.v[0] += tDelta.v[0];

faceX = -stepX;
faceY = 0;
faceZ = 0;
} else {
if (tMax.v[2] > range) break;
zPos += stepZ;
tMax.v[2] += tDelta.v[2];
faceX = 0;
faceY = 0;
faceZ = -stepZ;
}
} else {
if (tMax.v[1] < tMax.v[2]) {
if (tMax.v[1] > range) break;
yPos += stepY;
tMax.v[1] += tDelta.v[1];
faceX = 0;
faceY = -stepY;
faceZ = 0;
} else {
if (tMax.v[2] > range) break;
zPos += stepZ;
tMax.v[2] += tDelta.v[2];
faceX = 0;
faceY = 0;
faceZ = -stepZ;
}
}
} while (true);
}
}


I'm not sure how efficient this is. In the future when I run it every frame to show an indicator, I'm not sure how it will hold up.

• How's it holding up if I may ask? I'm curious about voxel rendering and came upon a similar prob with Bresenham.
– user77245
Jan 9, 2018 at 8:28

A good source for voxel based questions is PolyVox. For this particular case, just check out this: link. There is also a lot of documentation included there.

For the sake of simplicity lets stick with two dimensions: (x,y)

The problem:

We have some start position xs :: vec2(float).

We have a (normalized) direction vector dv :: vec2(float).

Given the above, we want to march in the direction of dv from xs, iterating through every voxel cell we cross.

The naïve solution

Advance the position along by some small epsilon:

xs += dv * 0.1f


At each step, I would floor my xs vector to get a vec2(int) - which represents the voxel cell I was currently in.

This solution works well but fails when casting near the corners of voxels. If your epsilon is too large, you could end up entirely skipping a voxel that your line barely intersects. This can be remedied by decreasing the size of epsilon, but this comes with the cost of requiring us to perform more iterations to march the same distance.

The efficient solution

At a high level, this solution is similar to above, but instead of marching a fixed epsilon - its dynamic: we work out the maximum distance we can move before we hit the next voxel cell and move that amount.

Define a voxel cursor cur :: vec2(int) - which represents the voxel cell we are currently living in.

Initialize the cursor to the cell we begin in:

cur = floor(xs)


Now lets enumerate through all faces of a square:

TOP,
LEFT,
RIGHT,
BOTTOM


For each face, face, let us grab the normal:

normal = get_normal(face)


Now lets grab the component of dv that travels according to normal:

proj_dv = dot(dv, normal)


If proj_dv <=0 then there is no motion in the direction of the normal and we can skip.

Next step is to work out how far we would have to travel in the normal's direction to hit a voxel boundary. This can be accomplished as follows:

// Will be 1 if the normal has a positive component and -1 if its negative
normal_sign = get_normal_sign(normal)

// Extract the component of the position that is parallel to the normal.
curr_pos = get_component(xs)

// Travelling in the direction of the normal, this is the next voxel
// cell we would hit (for this component)
next_pos = (normal_sign > 0 ? floor : ceil)(curr_pos + normal_sign)

delta_pos = abs(next_pos - curr_pos)


Finally, we must calculate how far along the dv vector we would need to travel to affect this change:

m = delta_pos / proj_dv


We can now loop through each face, returning the minimum m (and the face that corresponds to it).

We now have an m (min_m) such that xs + dv * min_m will take us to next nearest voxel boundary. Thus we update:

cursor += normal(min_face)
xs += dv * min_m


We can emit the cursor and continue...