Raytracing on a grid

Gnast’e Villon fires his death ray wildly as Hirsutia, intrepid heroine, dives for cover. Will she escape his palace or leave a nasty stain on his tapestry?

The last two articles have been about computing everything that’s visible from a single viewpoint. Sometimes, though, all we need is a line-of-sight check between two specific points (i.e. a gun and a target). We want to identify all grid squares that intersect the line segment. If any square is solid, the line of sight is blocked.

The Bresenham line-drawing algorithm is fairly well known and might seem to be a reasonable candidate for this task. However, in its usual form it does not generate all squares that intersect the line. Fortunately there is a simple algorithm that will.

Here is a sample problem:



The numbers at top and right are the grid coordinates; the line segment goes from (0.5, 0.5) to (2.5, 3.5). We want to identify all the squares marked in pink.

In pseudocode the algorithm is as follows:

Start at the square containing the line segment’s starting endpoint.
For the number of intersected squares:
    If the next intersection is with a vertical grid line:
        Move one square horizontally toward the other endpoint
    Else:
        Move one square vertically toward the other endpoint

The number of squares intersected by the line is fairly easy to compute; it’s one (for the starting square) plus the number of horizontal and vertical grid line crossings.

As we move along the line, we can track where we are as a fraction of the line’s total length. Call this variable t. In the example, we’d start with t = 0 in the lower left corner. The first horizontal grid intersection is at 1/4, while the first vertical intersection is at 1/6. Since the vertical one is closer, we step up one square and advance t to 1/6. Now the next vertical intersection is at 3/6, while the next horizontal intersection is still at 1/4, so we move horizontally, advancing t to 1/4, and so forth.

The fraction of the line between grid intersections in the horizontal direction is simply one over the line’s horizontal extent. Likewise, the line fraction between vertical grid intersections is one over the line’s vertical extent. These values go to infinity when the line is exactly horizontal or exactly vertical. Fortunately, IEEE 754 floating-point (used by all modern computers) can represent positive or negative infinity exactly and compute correctly with them, so the divisions by zero work fine.

The only gotcha (which I missed until it was pointed out in the comments) is that zero times infinity is undefined; the result is the dreaded NaN (Not a Number). So we still have to check for zero horizontal or vertical movement to avoid trying to multiply by dt_dx or dt_dy.

Here’s code to do this:

void raytrace(double x0, double y0, double x1, double y1)
{
    double dx = fabs(x1 - x0);
    double dy = fabs(y1 - y0);

    int x = int(floor(x0));
    int y = int(floor(y0));

    double dt_dx = 1.0 / dx;
    double dt_dy = 1.0 / dy;

    double t = 0;

    int n = 1;
    int x_inc, y_inc;
    double t_next_vertical, t_next_horizontal;

    if (dx == 0)
    {
        x_inc = 0;
        t_next_horizontal = dt_dx; // infinity
    }
    else if (x1 > x0)
    {
        x_inc = 1;
        n += int(floor(x1)) - x;
        t_next_horizontal = (floor(x0) + 1 - x0) * dt_dx;
    }
    else
    {
        x_inc = -1;
        n += x - int(floor(x1));
        t_next_horizontal = (x0 - floor(x0)) * dt_dx;
    }

    if (dy == 0)
    {
        y_inc = 0;
        t_next_vertical = dt_dy; // infinity
    }
    else if (y1 > y0)
    {
        y_inc = 1;
        n += int(floor(y1)) - y;
        t_next_vertical = (floor(y0) + 1 - y0) * dt_dy;
    }
    else
    {
        y_inc = -1;
        n += y - int(floor(y1));
        t_next_vertical = (y0 - floor(y0)) * dt_dy;
    }

    for (; n > 0; --n)
    {
        visit(x, y);

        if (t_next_vertical < t_next_horizontal)
        {
            y += y_inc;
            t = t_next_vertical;
            t_next_vertical += dt_dy;
        }
        else
        {
            x += x_inc;
            t = t_next_horizontal;
            t_next_horizontal += dt_dx;
        }
    }
}

Simplifying, Round One

The above algorithm generalizes nicely to three dimensions, if you ever find yourself needing to traverse a voxel grid. However, the 2D case can be simplified a bit.

First of all, we can avoid the two divisions by scaling our t value by the product of the two denominators. Instead of going from 0 to 1, it goes from 0 to dx*dy. This makes dt_dx = dy, and dt_dy = dx. Note that now, when a line is horizontal, dt_dx is zero, instead of dt_dy being infinity. (Conversely for vertical lines.) This means that the next horizontal crossing is always at 0, and thus t never advances. Since the next vertical crossing is greater than zero, the algorithm ends up always moving horizontally. This works fine as long as we don’t rely on t reaching the end of the line to indicate completion.

Secondly, we can combine the t, t_next_horizontal, and t_next_vertical variables into one “error” term. The “error” is the difference between t_next_horizontal and t_next_vertical; if it’s positive, we know one is closer; if it’s negative, the other is closer. We add or subtract dt_dx and dt_dy as appropriate when we move.

Here’s the simplified code:

#include <limits> // for infinity
void raytrace(double x0, double y0, double x1, double y1)
{
    double dx = fabs(x1 - x0);
    double dy = fabs(y1 - y0);

    int x = int(floor(x0));
    int y = int(floor(y0));

    int n = 1;
    int x_inc, y_inc;
    double error;

    if (dx == 0)
    {
        x_inc = 0;
        error = std::numeric_limits<double>::infinity();
    }
    else if (x1 > x0)
    {
        x_inc = 1;
        n += int(floor(x1)) - x;
        error = (floor(x0) + 1 - x0) * dy;
    }
    else
    {
        x_inc = -1;
        n += x - int(floor(x1));
        error = (x0 - floor(x0)) * dy;
    }

    if (dy == 0)
    {
        y_inc = 0;
        error -= std::numeric_limits<double>::infinity();
    }
    else if (y1 > y0)
    {
        y_inc = 1;
        n += int(floor(y1)) - y;
        error -= (floor(y0) + 1 - y0) * dx;
    }
    else
    {
        y_inc = -1;
        n += y - int(floor(y1));
        error -= (y0 - floor(y0)) * dx;
    }

    for (; n > 0; --n)
    {
        visit(x, y);

        if (error > 0)
        {
            y += y_inc;
            error -= dx;
        }
        else
        {
            x += x_inc;
            error += dy;
        }
    }
}

All-Integer Math

The above algorithm works with any input coordinates. If we restrict ourselves to integral input coordinates, or known fractions thereof, the algorithm can be converted to operate entirely with integer math. As an example let us assume that, if the (integral) input coordinates are (x0, y0) and (x1, y1), the line should be traced from (x0 + 0.5, y0 + 0.5) to (x1 + 0.5, y1 + 0.5). This is what’s going on in the example diagram, and is a reasonable choice for line-of-sight of actors who are standing in the centers of the grid squares.

As can be seen in the previous code, if the input coordinates are an integral distance apart, the only fractional number is the error term. We can once again scale everything up to make this number an integer as well. With endpoints offset from the grid by 0.5, it suffices to multiply by 2. The numerous floor() and ceil() operations are no longer needed, and the resulting code is quite simple:

void raytrace(int x0, int y0, int x1, int y1)
{
    int dx = abs(x1 - x0);
    int dy = abs(y1 - y0);
    int x = x0;
    int y = y0;
    int n = 1 + dx + dy;
    int x_inc = (x1 > x0) ? 1 : -1;
    int y_inc = (y1 > y0) ? 1 : -1;
    int error = dx - dy;
    dx *= 2;
    dy *= 2;

    for (; n > 0; --n)
    {
        visit(x, y);

        if (error > 0)
        {
            x += x_inc;
            error -= dy;
        }
        else
        {
            y += y_inc;
            error += dx;
        }
    }
}

Finally: The code above does not always return the same set of squares if you swap the endpoints. When error is zero, the line is passing through a vertical grid line and a horizontal grid line simultaneously. In this case, the code currently will always move vertically (the else clause), then horizontally. If this is undesirable, you could make the if statement break ties differently when moving up vs. down; or you could have a third clause for error == 0 which considers both moves (horizontal-then-vertical and vertical-then-horizontal).