# Cast ray to select block in voxel game

I am developing a game with a Minecraft-like terrain made out of blocks. Since basic rendering and chunk loading is done now, I want to implement block selecting.

Therefore I need to find out what block the first person camera is facing. I already heard of unprojecting the whole scene but I decided against that because it sounds hacky and isn't accurate. Maybe I could somehow cast a ray in view direction but I do not know how to check the collision with a block in my voxel data. Of course this calculations must be done on the CPU since I need the results to perform game logic operations.

So how could I find out which block is in front of the camera? If it is preferable, how could I cast a ray and check collisions?

• I've never done it myself. But couldn't you just have a "ray" ( linesegment in this case ) from the camera plane, a normal vector, with a certain length( you only want it to be within a radius) and see if it intersects with one of the blocks. I assume partial spacing and clipping is implemented as well. So knowing which blocks to test with shouldn't be that much of an issue...i think? – Sidar Jan 14 '13 at 12:17

When I had this problem while working on my Cubes, I found the paper "A Fast Voxel Traversal Algorithm for Ray Tracing" by John Amanatides and Andrew Woo, 1987 which describes an algorithm which can be applied to this task; it is accurate and needs only one loop iteration per voxel intersected.

I have written an implementation of the relevant parts of the paper's algorithm in JavaScript. My implementation adds two features: it allows specifying a limit on the distance of the raycast (useful for avoiding performance issues as well as defining a limited 'reach'), and also computes which face of each voxel the ray entered.

The input origin vector must be scaled such that the side length of a voxel is 1. The length of the direction vector is not significant but may affect the numerical accuracy of the algorithm.

The algorithm operates by using a parameterized representation of the ray, origin + t * direction. For each coordinate axis, we keep track of the t value which we would have if we took a step sufficient to cross a voxel boundary along that axis (i.e. change the integer part of the coordinate) in the variables tMaxX, tMaxY, and tMaxZ. Then, we take a step (using the step and tDelta variables) along whichever axis has the least tMax — i.e. whichever voxel-boundary is closest.

/**
* Call the callback with (x,y,z,value,face) of all blocks along the line
* segment from point 'origin' in vector direction 'direction' of length
*
* 'face' is the normal vector of the face of that block that was entered.
* It should not be used after the callback returns.
*
* If the callback returns a true value, the traversal will be stopped.
*/
function raycast(origin, direction, radius, callback) {
// From "A Fast Voxel Traversal Algorithm for Ray Tracing"
// by John Amanatides and Andrew Woo, 1987
// <http://www.cse.yorku.ca/~amana/research/grid.pdf>
// <http://citeseer.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.3443>
// Extensions to the described algorithm:
//   • Imposed a distance limit.
//   • The face passed through to reach the current cube is provided to
//     the callback.

// The foundation of this algorithm is a parameterized representation of
// the provided ray,
//                    origin + t * direction,
// except that t is not actually stored; rather, at any given point in the
// traversal, we keep track of the *greater* t values which we would have
// if we took a step sufficient to cross a cube boundary along that axis
// (i.e. change the integer part of the coordinate) in the variables
// tMaxX, tMaxY, and tMaxZ.

// Cube containing origin point.
var x = Math.floor(origin);
var y = Math.floor(origin);
var z = Math.floor(origin);
// Break out direction vector.
var dx = direction;
var dy = direction;
var dz = direction;
// Direction to increment x,y,z when stepping.
var stepX = signum(dx);
var stepY = signum(dy);
var stepZ = signum(dz);
// See description above. The initial values depend on the fractional
// part of the origin.
var tMaxX = intbound(origin, dx);
var tMaxY = intbound(origin, dy);
var tMaxZ = intbound(origin, dz);
// The change in t when taking a step (always positive).
var tDeltaX = stepX/dx;
var tDeltaY = stepY/dy;
var tDeltaZ = stepZ/dz;
// Buffer for reporting faces to the callback.
var face = vec3.create();

// Avoids an infinite loop.
if (dx === 0 && dy === 0 && dz === 0)
throw new RangeError("Raycast in zero direction!");

// Rescale from units of 1 cube-edge to units of 'direction' so we can
// compare with 't'.

while (/* ray has not gone past bounds of world */
(stepX > 0 ? x < wx : x >= 0) &&
(stepY > 0 ? y < wy : y >= 0) &&
(stepZ > 0 ? z < wz : z >= 0)) {

// Invoke the callback, unless we are not *yet* within the bounds of the
// world.
if (!(x < 0 || y < 0 || z < 0 || x >= wx || y >= wy || z >= wz))
if (callback(x, y, z, blocks[x*wy*wz + y*wz + z], face))
break;

// tMaxX stores the t-value at which we cross a cube boundary along the
// X axis, and similarly for Y and Z. Therefore, choosing the least tMax
// chooses the closest cube boundary. Only the first case of the four
// has been commented in detail.
if (tMaxX < tMaxY) {
if (tMaxX < tMaxZ) {
// Update which cube we are now in.
x += stepX;
// Adjust tMaxX to the next X-oriented boundary crossing.
tMaxX += tDeltaX;
// Record the normal vector of the cube face we entered.
face = -stepX;
face = 0;
face = 0;
} else {
z += stepZ;
tMaxZ += tDeltaZ;
face = 0;
face = 0;
face = -stepZ;
}
} else {
if (tMaxY < tMaxZ) {
y += stepY;
tMaxY += tDeltaY;
face = 0;
face = -stepY;
face = 0;
} else {
// Identical to the second case, repeated for simplicity in
// the conditionals.
z += stepZ;
tMaxZ += tDeltaZ;
face = 0;
face = 0;
face = -stepZ;
}
}
}
}

function intbound(s, ds) {
// Find the smallest positive t such that s+t*ds is an integer.
if (ds < 0) {
return intbound(-s, -ds);
} else {
s = mod(s, 1);
// problem is now s+t*ds = 1
return (1-s)/ds;
}
}

function signum(x) {
return x > 0 ? 1 : x < 0 ? -1 : 0;
}

function mod(value, modulus) {
return (value % modulus + modulus) % modulus;
}

• Does this algorithm also work for negative number space? I implemented the algorithm only just and generally I am impressed. It works great for positive coordinates. But for some reason I get strange results if negative coordinates are involved sometimes. – danijar Mar 11 '13 at 23:15
• @danijar I couldn't get the intbounds/mod stuff to work with negative space, so I use this: function intbounds(s,ds) { return (ds > 0? Math.ceil(s)-s: s-Math.floor(s)) / Math.abs(ds); }. As Infinity is greater than all numbers, I don't think you need to guard against ds being 0 there either. – Will Jul 5 '13 at 21:00
• @BotskoNet That sounds like you have a problem with unprojecting to find your ray. I had problems like that early on. Suggestion: Draw a line from origin to origin+direction, in world space. If that line is not under the cursor, or if it does not appear as a point (since projected X and Y should be equal) then you have a problem in the unprojection (not part of this answer's code). If it's reliably a point under the cursor then the problem is in the raycast. If you still have a problem, please ask a separate question instead of extending this thread. – Kevin Reid Dec 11 '13 at 23:18
• The edge case is where a coordinate of the ray origin is an integer value, and the corresponding part of the ray direction is negative. The initial tMax value for that axis should be zero, since the origin is already at the bottom edge of its cell, but it is instead 1/ds causing one of the other axes to be incremented instead. The fix is to write intfloor to check if both ds is negative and s is an integer value (mod returns 0), and return 0.0 in that case. – codewarrior Dec 24 '14 at 12:00
• Here is my port to Unity: gist.github.com/dogfuntom/cc881c8fc86ad43d55d8 . Though, with some additional changes: integrated Will's and codewarrior's contributions and made possible to cast in an unlimited world. – Maxim Kamalov Sep 12 '15 at 1:28

Perhaps look into Bresenham's line algorithm, particularly if you're working with unit-blocks (as most minecraftish games tend to).

Basically this takes any two points, and traces an unbroken line between them. If you cast a vector from the player to their maximum picking distance, you can use this, and the players positions as points.

I have a 3D implementation in python here: bresenham3d.py.

• A Bresenham-type algorithm will miss some blocks, though. It doesn't consider every block the ray passes through; it'll skip some in which the ray doesn't get close enough to the block center. You can see this clearly from the diagram on Wikipedia. The block 3rd down and 3rd right from the top-left corner is an example: the line passes through it (barely) but Bresenham's algorithm doesn't hit it. – Nathan Reed Jan 16 '13 at 7:02

To find the first block in front of the camera, create a for loop that loops from 0 to some maximum distance. Then, multiply the camera's forward vector by the counter and check if the block at that position is solid. If it is, then store the position of the block for later use and stop looping.

If you also want to be able to place blocks, face-picking is no harder. Simply loop back from the block and find the first empty block.

• Wouldn't work, with an angled forward vector it would be very possible to have a point before one part of a block, and the subsequent point after, missing the block. The only solution with this would be to reduce the size of the increment, but you'd have to get it so small as to make other algorithms far more effective. – Phil Jan 16 '13 at 3:07
• This works pretty well with my engine; I use an interval of 0.1. – untitled Jan 17 '13 at 23:27
• Like @Phil pointed out, the algorithm would miss blocks where only a small edge is seen. Furthermore looping backwards for placing blocks wouldn't work. We would have to loop forward as well and decrement the result by one. – danijar Feb 17 '13 at 14:53

I made a post on Reddit with my implementation, which uses Bresenham's Line Algorithm. Here's an example of how you would use it:

// A plotter with 0, 0, 0 as the origin and blocks that are 1x1x1.
PlotCell3f plotter = new PlotCell3f(0, 0, 0, 1, 1, 1);
// From the center of the camera and its direction...
plotter.plot( camera.position, camera.direction, 100);
// Find the first non-air block
while ( plotter.next() ) {
Vec3i v = plotter.get();
Block b = map.getBlock(v);
if (b != null && !b.isAir()) {
plotter.end();
// set selected block to v
}
}


Here is the implementation itself:

public interface Plot<T>
{
public boolean next();
public void reset();
public void end();
public T get();
}

public class PlotCell3f implements Plot<Vec3i>
{

private final Vec3f size = new Vec3f();
private final Vec3f off = new Vec3f();
private final Vec3f pos = new Vec3f();
private final Vec3f dir = new Vec3f();

private final Vec3i index = new Vec3i();

private final Vec3f delta = new Vec3f();
private final Vec3i sign = new Vec3i();
private final Vec3f max = new Vec3f();

private int limit;
private int plotted;

public PlotCell3f(float offx, float offy, float offz, float width, float height, float depth)
{
off.set( offx, offy, offz );
size.set( width, height, depth );
}

public void plot(Vec3f position, Vec3f direction, int cells)
{
limit = cells;

pos.set( position );
dir.norm( direction );

delta.set( size );
delta.div( dir );

sign.x = (dir.x > 0) ? 1 : (dir.x < 0 ? -1 : 0);
sign.y = (dir.y > 0) ? 1 : (dir.y < 0 ? -1 : 0);
sign.z = (dir.z > 0) ? 1 : (dir.z < 0 ? -1 : 0);

reset();
}

@Override
public boolean next()
{
if (plotted++ > 0)
{
float mx = sign.x * max.x;
float my = sign.y * max.y;
float mz = sign.z * max.z;

if (mx < my && mx < mz)
{
max.x += delta.x;
index.x += sign.x;
}
else if (mz < my && mz < mx)
{
max.z += delta.z;
index.z += sign.z;
}
else
{
max.y += delta.y;
index.y += sign.y;
}
}
return (plotted <= limit);
}

@Override
public void reset()
{
plotted = 0;

index.x = (int)Math.floor((pos.x - off.x) / size.x);
index.y = (int)Math.floor((pos.y - off.y) / size.y);
index.z = (int)Math.floor((pos.z - off.z) / size.z);

float ax = index.x * size.x + off.x;
float ay = index.y * size.y + off.y;
float az = index.z * size.z + off.z;

max.x = (sign.x > 0) ? ax + size.x - pos.x : pos.x - ax;
max.y = (sign.y > 0) ? ay + size.y - pos.y : pos.y - ay;
max.z = (sign.z > 0) ? az + size.z - pos.z : pos.z - az;
max.div( dir );
}

@Override
public void end()
{
plotted = limit + 1;
}

@Override
public Vec3i get()
{
return index;
}

public Vec3f actual() {
return new Vec3f(index.x * size.x + off.x,
index.y * size.y + off.y,
index.z * size.z + off.z);
}

public Vec3f size() {
return size;
}

public void size(float w, float h, float d) {
size.set(w, h, d);
}

public Vec3f offset() {
return off;
}

public void offset(float x, float y, float z) {
off.set(x, y, z);
}

public Vec3f position() {
return pos;
}

public Vec3f direction() {
return dir;
}

public Vec3i sign() {
return sign;
}

public Vec3f delta() {
return delta;
}

public Vec3f max() {
return max;
}

public int limit() {
return limit;
}

public int plotted() {
return plotted;
}

}

• As someone in the comments noticed, your code is undocumented. While the code may be helpful, it doesn't quite answer the question. – Anko Jan 16 '13 at 7:03

My approach is somewhat different, yet heavily relies on answers in this thread (especially @Maxim Kamalov's comment with his gist - was great starting point, basically main logic is taken from there and improved).

The problem with existing code for me was different edge-cases (coordinate delta equal to zero ("straight lines"), hitting exactly the edge (where tMax values are equal), negative traversal, not accounting the fact that blocks in my data structure are .5f-based, not 0-based (center is in the very center of the cube), floating point precision, and a few others).

Information upon this is not so spread, so I decided to share my traversal approach. You can think of it as variation of Amanatides and Woo algorithm. Main difference is that I'm not using ray, but rather 2 points - but it shouldn't be too hard to adjust.

Note that It might (and probably will) be slower, but it is more precise than face-to-face traversal.

Also, I'm not quite fond of using local functions, but that made developing the logic for me easier.

Hope it helps to whoever it may concern :)

using System.Collections.Generic;

public static class Foo
{
public static IEnumerable<Vector3> GetIntersectedWorldBlocksD(Vector3 pointA, Vector3 pointB)
{
var line = new HashSet<Vector3>
{
new Vector3(Mathf.Round(pointA.X), Mathf.Round(pointA.Y), Mathf.Round(pointA.Z))
};

var x = (float) Mathf.Round(pointA.X);
var y = (float) Mathf.Round(pointA.Y);
var z = (float) Mathf.Round(pointA.Z);

var delta = pointB - pointA;
var deltaNormal = delta.Normal;
var dx = deltaNormal.X;
var dy = deltaNormal.Y;
var dz = deltaNormal.Z;

var stepX = Signum(dx);
var stepY = Signum(dy);
var stepZ = Signum(dz);

var tMaxX = Intbound(pointA.X - 0.5f, deltaNormal.x);
var tMaxY = Intbound(pointA.Y - 0.5f, deltaNormal.y);
var tMaxZ = Intbound(pointA.Z - 0.5f, deltaNormal.z);

var tDeltaX = stepX / dx;
var tDeltaY = stepY / dy;
var tDeltaZ = stepZ / dz;

var lineLength = delta.Length;

bool StepX()
{
if (tMaxX >= lineLength) return true;
x += stepX;
tMaxX += tDeltaX;
return false;
}

bool StepY()
{
if (tMaxY >= lineLength) return true;
y += stepY;
tMaxY += tDeltaY;
return false;
}

bool StepZ()
{
if (tMaxZ >= lineLength) return true;
z += stepZ;
tMaxZ += tDeltaZ;
return false;
}

var res = true;

if (stepX != 0 && stepY != 0 && stepZ != 0)
while (true)
{
var tX = tMaxX.Precision5();
var tY = tMaxY.Precision5();
var tZ = tMaxZ.Precision5();
if (tX < tY)
{
if (tX < tZ)
res = StepX();
else if (tZ < tX)
res = StepZ();
else
res = StepZ() || StepX();
}
else if (tY < tX)
{
if (tY < tZ)
res = StepY();
else if (tZ < tY)
res = StepZ();
else
res = StepY() || StepZ();
}
else
{
if (tY < tZ)
res = StepX() || StepY();
else if (tZ < tY)
res = StepZ();
else
res = StepX() || StepY() || StepZ();
}

if (res) break;
}
else if (stepX == 0 && stepY != 0 && stepZ != 0)
while (true)
{
var tY = tMaxY.Precision5();
var tZ = tMaxZ.Precision5();
if (tY < tZ)
res = StepY();
else if (tZ < tY)
res = StepZ();
else
res = StepY() || StepZ();
if (res) break;
}
else if (stepX != 0 && stepY == 0 && stepZ != 0)
while (true)
{
var tX = tMaxX.Precision5();
var tZ = tMaxZ.Precision5();
if (tX < tZ)
res = StepX();
else if (tZ < tX)
res = StepZ();
else
res = StepZ() || StepX();

if (res)
{
break;
}
}
else if (stepX != 0 && stepY != 0 && stepZ == 0)
while (true)
{
var tY = tMaxY.Precision5();
var tX = tMaxZ.Precision5();
if (tY < tX)
res = StepY();
else if (tX < tY)
res = StepX();
else
res = StepY() || StepX();

if (res) break;
}
else if (stepX != 0 && stepY == 0 && stepZ == 0)
while (true)
{
res = StepX();
if (res) break;
}
else if (stepX == 0 && stepY != 0 && stepZ == 0)
while (true)
{
res = StepY();
if (res) break;
}
else if (stepX == 0 && stepY == 0 && stepZ != 0)
while (true)
{
res = StepZ();
if (res) break;
}

return line;
}

private static int Signum(float x)
{
return x > 0 ? 1 : x < 0 ? -1 : 0;
}

private static float Ceil(float s)
{
return s == 0f ? 1f : Mathf.Ceil(s);
}

private static float Intbound(float s, float ds)
{
if (ds < 0 && Mathf.Round(s) == s) return 0;
s = Mod(s, 1);
return (ds > 0 ? Ceil(s) - s : s - Mathf.Floor(s)) / Mathf.Abs(ds);
}

private static float Mod(float value, float modulus)
{
return (value % modulus + modulus) % modulus;
}

public static float Precision5(this float val)
{
return  ((int)(100000 * val)) / 100000.0f;
}
}
$$$$
`