Monday, March 26, 2007

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).

56 comments:

Anonymous said...

This article was an absolute lifesafer. I'm developing a player-centric overhead tile RPG and I have been really struggling trying to figure out how to implement LOS. thank you so much!

Anonymous said...

Thankyou sooooooooo much :-) This problem has been doing my head in for days. After so many failed attempts your solution worked perfectly.

Unknown said...

Thanks, was a great find, I really needed this algorithm for LOS on a grid that did not 'peek' through corners. :)

Anonymous said...

This was instrumental in putting me on the right track for my game - thanks so much!

http://vitaro.wordpress.com/heroquest/

Anonymous said...

I'm not certain, but i believe this line is bugged:

t_next_horizontal = (ceil(x0) - x0) * dt_dx;

since (ceil(x0) - x0) will return 0 when x0 is an exact integer. Instead, you can change that to:

t_next_horizontal = (floor(x0) + 1 - x0) * dt_dx;

Same for ceil(y0)

James McNeill said...

Anonymous, you are absolutely correct. Thanks for catching that! I'm sorry that was incorrect for so long. I've changed the code as you indicated.

SyntaxError said...

I think there is still a bug in case of a straight line when using an integral starting point.

In case of a vertical line dx will be 0 and dt_dx will be positive infinity. If x0 is a fractional number, t_next_horizontal will be positive infinity too, because multiplying a positive number with positive infinity yields positive infinity. If x0 is a integral number, however, t_next_horizontal will be NaN, because multiplying zero with infinity is not defined.

The if comparison inside the loop will be always return false, because any comparison with NaN returns false. Instead of the vertical line a horizontal line will be drawn.

Since it happens to be correct for horizontal lines, you could simply change the if-statement in the loop to

if (t_next_vertical < t_next_horizontal || dx == 0)

to fix this. However, I would suggest to handle this case in a clean way:

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 t = 0, dt_dx, dt_dy;
..int n = 1;
..int x_inc, y_inc;
..double t_next_vertical, t_next_horizontal;
..
..if (dx == 0)
..{
....x_inc = 0;
....dt_dx = 0;
....t_next_horizontal = POSITIVE_INFINITY;
..}
..else if (x1 > x0)
..{
....x_inc = 1;
....dt_dx = 1.0 / dx;
....n += int(floor(x1)) - x;
....t_next_horizontal = (floor(x0) + 1 - x0) * dt_dx;
..}
..else
..{
....x_inc = -1;
....dt_dx = 1.0 / dx;
....n += x - int(floor(x1));
....t_next_horizontal = (x0 - floor(x0)) * dt_dx;
..}
..
..if (dy == 0)
..{
....y_inc = 0;
....dt_dy = 0;
....t_next_vertical = POSITIVE_INFINITY;
..}
..else if (y1 > y0)
..{
....y_inc = 1;
....dt_dy = 1.0 / dy;
....n += int(floor(y1)) - y;
....t_next_vertical = (floor(y0) + 1 - y0) * dt_dy;
..}
..else
..{
....y_inc = -1;
....dt_dy = 1.0 / dy;
....n += y - int(floor(y1));
....t_next_vertical = (y0 - floor(y0)) * dt_dy;
..}
..
..while (true)
..{
....visit(x, y);
....
....if (--n == 0) // saves us an unnecessary update step
......break;
....
....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;
....}
..}
}

By the way, how is this algorithm called?

SyntaxError said...

Actually, if you use the code I've posted, you don't have to rely on IEEE infinity arithmetic rules at all. Since t will always be <= 1, you can use any number > 1 instead of POSITIVE_INFINITY.

SyntaxError said...

I just took a look at your simplified version which does not use infinity arithmetic. Unfortunately, it suffers from a similar problem in the same case. If you draw a vertical line starting at an integral position, dx and error will be zero. Therefore the loop will make a horizontal step first and you end up with a L-shape with ends one pixel too early.

Now you can't simply replace > with a >= inside the loop. Then you would have the same problem with horizontal lines. But you can use the same fix I proposed for the elaborate version by replacing the if-statement in the loop with:

if (error > 0 || v.x == 0)

SyntaxError said...

And analog to above, here the clean version of the optimized algorithm. It saves you the additional comparison inside of the loop.

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 = [some large value, e.g. DOUBLE_MAXVAL or POSITIVE_INFINTIY]
..}
..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 -= [some large value, e.g. DOUBLE_MAXVAL or POSITIVE_INFINTIY]
..}
..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;
..}
..
..while (true)
..{
....visit(x, y);
....
....if (--n == 0)
......break;
....
....if (error > 0)
....{
......y += y_inc;
......error -= dx;
....}
....else
....{
......x += x_inc;
......error += dy;
....}
..}
}

James McNeill said...

Thank you for the fixes! I've tested them and incorporated them into the post above.

Unknown said...

Thanks a ton for the article, but I think I must be missing something obvious. I'm not a C++ programmer, but I know C# so I feel like I have at least some idea of what your code is doing, but I can't figure out for the life of me, how your simplified raytrace function does anything except traverse the grid.

What does the "visit" method do (I don't see it defined)? Also, your method returns void, when I would expect it to return at least a boolean (or perhaps the cell where the ray impacted).

Any light you can shed is much appreciated.

James McNeill said...

Kai,

You're right, visit() isn't defined here. It's just a stub for whatever you want to do. It basically says "the ray goes through this square". So if you wanted to know whether the ray hits a solid square, say (it sounds like that's what you're interested in), you'd replace visit(x, y) with something like:

if (solid(x, y)) return true;

...and then return false at the end of the function. Obviously solid() is a function you'd need to define.

Duke said...

Hey, I know this is old, but I found it and am using some of the code to program a line-of-sight component. I was getting innacurate results with the INT version until I made a slight modification. When the "line" crossed directly over the vertex, it would visit an adjacent tile perpendicular to the "line." Instead, I wanted it to pass through a vertex with no visit. All I did was add a 3rd case if (error==0) then I advance both x and y and then n---. Works great!

for (; n > 0; --n)
{
//trace(n)
//visit(xTmp, yTmp); tile2d_array[xTmp][yTmp].turnOn(val);

if (error > 0)
{
xTmp += x_inc;
error -= dy;
}
else if (error < 0)
{
yTmp += y_inc;
error += dx;
}else if(error == 0){
xTmp += x_inc;
error -= dy;
yTmp += y_inc;
error += dx;
n--
}
}

James McNeill said...

Duke: Yeah, that's a valid way to handle the case where the line goes over a vertex exactly. It all depends on what you need from the line tracer. Thanks for the code!

Unknown said...

2Duke: In your solution is little mistake. Each time you go straight through a vertex, you increase both x and y and thus skip one step of the FOR cycle. Then the algorithm goes as far behind the end point as many vertices are on the way to the end point. Therefore use WHILE cycle in this case ;)

James McNeill said...

I think the code is correct (although I haven't run it to verify). The (error == 0) condition decrements n; this combines with the for loop's decrement to ensure the correct number of steps are taken.

Unknown said...

Oh, I am sorry, I've missed that line.
It's just not standard approach to change loop control variable in FOR cycle. It's dangerous and people should use WHILE in that case.

Unknown said...

> The above algorithm generalizes nicely
> to three dimensions, if you ever find
> yourself needing to traverse a voxel
> grid.

Hi, I tried to generalize your code to three dimensions, but I'm unsure whether I did it correctly. Also, I'm wondering about how to do simplifications for integer values only, as you have done for 2 dimensions.

Here is my code (in Java):

public void raytrace(double x0, double y0, double z0, double x1, double y1, double z1)
{
double dx = Math.abs(x1 - x0);
double dy = Math.abs(y1 - y0);
double dz = Math.abs(z1 - z0);

int x = (int) (Math.floor(x0));
int y = (int) (Math.floor(y0));
int z = (int) (Math.floor(z0));

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

double t = 0;

int n = 1;
int x_inc, y_inc, z_inc;
double t_next_y, t_next_x, t_next_z;

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

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

if (dz == 0)
{
z_inc = 0;
t_next_z = dt_dz; // infinity
}
else if (z1 > z0)
{
z_inc = 1;
n += (int) (Math.floor(z1)) - z;
t_next_z = (Math.floor(z0) + 1 - z0) * dt_dz;
}
else
{
z_inc = -1;
n += z - (int) (Math.floor(z1));
t_next_z = (z0 - Math.floor(z0)) * dt_dz;
}

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

if (t_next_x <= t_next_y && t_next_x <= t_next_z) // t_next_x is smallest
{
x += x_inc;
t = t_next_x;
t_next_x += dt_dx;
}
else if (t_next_y <= t_next_x && t_next_y <= t_next_z) // t_next_y is smallest
{
y += y_inc;
t = t_next_y;
t_next_y += dt_dy;
}
else // t_next_z is smallest
{
z += z_inc;
t = t_next_z;
t_next_z += dt_dz;
}
}
}

James McNeill said...

I'm sorry, skrjablin! Blogger flagged your comment as spam for some reason, and it took me a few days to figure out what had happened and un-flag it.

Your code looks correct to me (I haven't tried running it, though).

Simplifying to integer math is not going to change the structure of the code, I think; it's probably simplest to do it the way you have there. You'd figure out a common denominator such that you can represent the size of a step in each of the X, Y, and Z dimensions as integral amounts. For instance, if you're moving 2 units in X, 3 units in Y, and 4 units in Z you could multiply all those together to get a denominator of 24. A step in X would move you 12/24 of the way along the line, a step in Y would move you 8/24, and a step in Z would move you 6/24. (In this case there's a smaller common denominator due to the common factors of 2 and 4, but I wouldn't worry about finding the smallest common denominator. Any common denominator will work.)

Unknown said...

Thank you for the response! It wasn't that long a wait for it after all. If the code is correct, I think I will keep it like that for now and only optimize for integers if I need to speed up things.

emily said...

Thanks so much for your codes, it is very clearly to me. I'm trying to apply your ideas to my project.

But I have a few questions about the accuracy issue of your codes.

As you see, there is another simple way to get the intersections of the example using line equation y=ax+b. By iterating the vertical grid lines and horizontal grid lines, the intersections coordinates (x,y) can be easily obtained and then the intersected grid cells are got.

But,for float endpoint cases, when the line segment passes the vertex of a grid cell, your codes and the y=ax+b method will produce different results because of different floating point operations, right?

Another thing, when the grid is really large, your codes will produce propagation errors because of additive floating operations, right?

I am wondering if there is literature related to your codes, could I call it as extended Bresenham's algorithm?

Thanks a lot!

James McNeill said...

I'd have to think about the floating-point precision issues a bit to be able to give a good answer.

I didn't come up with these algorithms. For instance, Daniel Cohen has an article in Graphics Gems 4: "Voxel Traversal along a 3D Line" which covers this, including code for the integer version of 3D traversal. I suppose it is fairly related to the Bresenham line-drawing algorithm.

Q said...

Awesome, thanks man!
Works perfectly.
I wasted hours on a much more complicated algorithm. arghhhhh

Anonymous said...
This comment has been removed by the author.
Anonymous said...
This comment has been removed by the author.
Anonymous said...

EDIT: well no html tags in Blogger comments? Ah well my graphs are squished but I'll post anyway:



I see this post is from 2007 but I just stumbled across it. Saved me a ton of work and pushed my knowledge along a bit as well...

Algorithm works great, but there is a further issue with the particular game I am modding. It is square (tile) based game and the obstructions to LOS use a graphic that does not fill the entire tile, it is a round object that graphically does not fill the entire tile like a square building the same size of the tile would. The code is correctly blocking LOS but intuitively (from the player's perspective) the LOS should be clear. In layman's terms the player's unit should be able to see (shoot) across the corners of a blocked tile but not across the center.

I need to mod this code to break up each game tile into 9 virtual tiles, only the center one containing either the firing unit, obstruction or target.

In this example the entire (virtual) 3X3 grid represents ONE game tile where XX is either firing unit, obstruction or target.

!----!----!----!
! ! ! !
!----!----!----!
! ! XX ! !
!----!----!----!
! ! ! !
!----!----!----!


Now in this second example EACH SQUARE represents a game tile. In each case, the firing unit needs to be able to see (shoot) at it's target. using the current algrithm all of these LOS would be blocked. I need to be able to calculate if a unit can see (fire) through the *corners* of a game tile.

UU = firing unit, XX = obstruction, TT = target tile

!----!----!----!----!----!----!
! XX ! TT ! ! ! XX ! TT !
!----!----!----!----!----!----!
! UU ! ! ! UU ! ! !
!----!----!----!----!----!----!
! ! ! ! ! ! !
!----!----!----!----!----!----!
! ! ! ! ! XX ! TT !
!----!----!----!----!----!----!
! ! ! UU ! ! ! !
!----!----!----!----!----!----!
! ! ! ! ! ! !
!----!----!----!----!----!----!


Hope this all makes sense, I am currently working on a solution. Will post here when complete (barring brain explosion). Any input appreciated.

-c

James McNeill said...

Having some sort of higher-resolution interior for squares in terms of how they block visibility seems totally reasonable. I did something like that for my visibility computation in the previous two posts:

http://playtechs.blogspot.com/2007/03/2d-portal-visibility-part-1.html

That's more of an area approach, though. If you just want to cast a single ray from one spot to another on the map and see if/where it is obstructed, another approach might be to build a small 2D BSP for each tile type that represents its solid and open sub-regions. Then, as you run the raytrace across the grid (via the method in this post), look up each interesected tile's BSP and clip the line against that to see if it's obstructed by something inside the tile.

Anonymous said...
This comment has been removed by the author.
Anonymous said...

This is my solution. I'm converting tiles to 3X3 and mapping start and target tiles on this grid, checking for tiles intersected by the ray using your modified Bresenham algorithm, and checking only the center square of each 3X3 grid that is intersected for an obstruction to the ray. This works nicely in this particular game so that the result meets the expectation from the player's perspective, as the game graphic for a Peak (mountain) is visually a cone-shaped object in the middle of the tile. Thanks a million for your blog posts, all of your blogs are EXTREMELY helpful and useful, especially for a rookie.

// this is a 2D line of sight check based on Bresenham algorithm interpretation by James McNeill
// with a twist because the Peak graphic in Civ4 does not fill the entire tile and
// it looks counterintuitive to the player to not be able to shoot through the corner of a Peak tile
// we are going to expand the tiles x3 onto a virtual grid so that each tile is now 3X3 or 9 tiles
// the firing unit peak and target always occupy the center virtual tile in the 9 virtual tiles
// we will then check only the center of these 9 virtual tiles for a Peak obstacle
// so that we can shoot past peaks but not over or through them and we meet player expectations
bool CvPlot::checkForPeak(int startX, int startY, int targetX, int targetY) const // ChazMod Ranged Fire
{
// init init
CvPlot* pCheckPlot;
int vstartX;
int vstartY;
int vtargetX;
int vtargetY;
int testX = 0;
int testY = 0;
int vobsX;
int vobsY;
int obsX;
int obsY;

// convert startXY and targetXY to place them on a virtual grid
if (startX <= targetX)
{
vstartX = 1;
vtargetX = (((targetX - startX)*3)+vstartX);
}
else
{
vtargetX = 1;
vstartX = (((startX - targetX)*3)+vtargetX);
}


if (startY <= targetY)
{
vstartY = 1;
vtargetY = (((targetY - startY)*3)+vstartY);
}
else
{
vtargetY = 1;
vstartY = (((startY - targetY)*3)+vtargetY);
}

//init based on virtual grid
int iDDX = abs(vtargetX - vstartX);
int iDDY = abs(vtargetY - vstartY);
int iX = vstartX;
int iY = vstartY;
int iN = 1 + iDDX + iDDY;
int iX_inc = (vtargetX > vstartX) ? 1 : -1;
int iY_inc = (vtargetY > vstartY) ? 1 : -1;
int error = iDDX - iDDY;
iDDX *= 2;
iDDY *= 2;

// now we have vstartXY and vtargetXY on an expanded virtual grid and can check for obstacles
for (; iN > 0; --iN)
{
// note that we are only checking the center virtual tile in each 3X3 grid for an obstacle
if (iX > vstartX)
{
if (((iX - vstartX) % 3) == 0)
{
testX = 1;
}
}
else if (iX < vstartX)
{
if (((vstartX - iX) % 3) == 0)
{
testX = 1;
}
}
else
{
testX = 1;
}

if (iY > vstartY)
{
if (((iY - vstartY) % 3) == 0)
{
testY = 1;
}
}
else if (iY < vstartY)
{
if (((vstartY - iY) % 3)== 0)
{
testY = 1;
}
}
else
{
testY = 1;
}

if (testX == 1 && testY == 1)
{
// convert the check tile from virtualXY to gameXY so that we can check isPeak
vobsX = iX;
vobsY = iY;

if (vobsX > vstartX)
{
obsX = (startX + ((vobsX - vstartX)/3));
}
else if (vstartX > vobsX)
{
obsX = (startX - ((vstartX - vobsX)/3));
}
else
{
obsX = startX;
}

if (vobsY > vstartY)
{
obsY = (startY + ((vobsY - vstartY)/3));
}
else if (vstartY > vobsY)
{
obsY = (startY - ((vstartY - vobsY)/3));
}
else
{
obsY = startY;
}

pCheckPlot = GC.getMapINLINE().plotINLINE(obsX, obsY);

if (pCheckPlot->isPeak())
{
return true;
}
}

testX = 0;
testY = 0;

if (error > 0)
{
iX += iX_inc;
error -= iDDY;
}
else
{
iY += iY_inc;
error += iDDX;
}
}

return false;
}

Sirius the Great Dane said...

This article really helped me! I have been searching for days for a similar article. I am using this in a 3-dimensional voxel space and it works great.

Thanks again!

Conor said...

Hi James

This is a great piece of code - just what i was looking for.

I'm trying to alter it slightly, so it will operate with rectangular grids rather than 1x1 squares but I am not having any luck.

I know the widths and lengths of the grid rectangles, so i can edit the x_inc and y_inc accordingly - but im not sure how to change it so the algorithm knows how to calculate the error value correctly.

This piece;

error = (floor(x0) + 1 - x0) * dy;

I'm using the first simplified one for 2D space.

Any suggestions?

Conor said...

Think I got it cracked - I just need to re-introduce the t variable into the code.

Unknown said...

I have implemented a similar algorithm in the past, but your approach is much nicer - handles the corner cases very cleanly. Thank you for sharing.

Are your code snippets on this blog under any particular license or condition of use?

James McNeill said...

The code snippets are free to use as you see fit and no attribution is necessary; I hadn't thought as far as a license. I developed and wrote them myself (the ones in the article, at least) so they shouldn't be encumbered by anyone else's IP. I'm glad it's useful!

Unknown said...

Cheers!

I've used a 3D variant of the first technique in a modified version of the JBullet physics engine to support voxel worlds. The voxels are handled as a special shape that takes advantage of the grid nature of the world.

dmplusge said...

In the first code example, if both dx and dy are <1.0, the final t value is outside of the range [0,1] and not suitable for calculating the length of the line travelled. dt_dx and dt_dy are very large in this case due to the 1.0/dx calculation, and the final t value will be either 0.0 (if visit() breaks on the first iteration) or at least at large as min(dt_dx, dt_dy).

Do you have any suggestions for computing the fraction of the line travelled when dx and dy are <1.0?

James McNeill said...

Looking at that first example I'm not sure why I even have the variable t; it isn't used. Nevertheless, the value of t when visit() is called represents the fraction of the line segment traversed prior to entry into the square being visited.

After all the squares have been visited you'd expect all of the line segment to have been traversed, right? Which would correspond to t=1. What I'm saying is I don't think it's meaningful to look at t after the loop has finished. In this case, as you point out, it represents the value of t when the line (not the line segment) exits the last visited square, which will generally be greater than one, not just when dx and dy are less than one.

dmplusge said...

I was referring to the t value in the case where visit() caused a break and the line was only partially travelled.

Either way, there was a mistake in my t_next comparisons when extended to 3 dimensions that caused the problem to appear (comparing infinities) when the small line segment crossed at least 2 cube boundaries.

Thanks for the code.

Anonymous said...

Please give me a little push...

I have implemented the algorithm, but I need it to determine where my 'for of war' needs to be.

https://gist.github.com/4192017

FYI, wallsGFX[][] contains the wall data. If wallsGFX■[y] == 0, you can see through it.

rays[][] contains the current ray, and it is re-initialised to 0 at the start of each iteration of the loop.

lineOfSight[][] is where I want to store a '1' for tiles that can be seen, and a '0' for tiles that cannot.

How would you go about transferring this information?

Have I been clear?

I will try to post a little example

wallsGFX[][]

0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 1 1 1 1 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0

rays[][]

0 0 0 1 0 0 0 0
0 0 0 1 0 0 0 0
0 0 0 1 0 0 0 0
0 0 0 0 1 0 0 0
0 0 0 0 1 0 0 0
0 0 0 0 1 0 0 0
0 0 0 0 0 1 0 0

what lineOfSight[][] should show after this one iteration

0 0 0 1 0 0 0 0
0 0 0 1 0 0 0 0
0 0 0 1 0 0 0 0
0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0

James McNeill said...

honeyspider,

I took a look at your code on github. One way to do it might be to get rid of your rays[] array and write directly to lineOfSight[] instead, stopping the ray-tracing loop when you hit a 1 in wallsGFX[].

If you're implementing line of sight, though, you might take a look at my article on that; it has a fairly simple recursive solution that avoids the aliasing problems you'll get with the ray-tracing approach:

http://playtechs.blogspot.com/2007/03/2d-portal-visibility-part-1.html
http://playtechs.blogspot.com/2007/03/2d-portal-visibility-part-2.html

Good luck!

Unknown said...

Hmm, using the second of the three sets of code, I am finding it to be, sometimes off by an entire row of tiles that the line is actually over.

James McNeill said...

Can you give an example set of inputs that would illustrate the problem?

Bhanu Chander V said...

wow, great and simple logic, thank u so much!

Buttonius said...

Very useful example code. Thank your for that. Important question near the end.

I've adapted the code for a project where I need to accumulate values that correspond to the cells in a grid, accompanied by the path length in each of those cells.

As a usage example for such code consider a weather map with a grid that has one value indicating the thickness of the snow cover for each cell. How much snow must be displaced to clear a (one meter wide) path from A to B?

This took a little extra code to correctly compute the path length to the first grid boundary and then to each subsequent grid line boundary and - at the end - to the end point.

Question. Are there any limitations on the use of your example code?

If you like, I can make my Java code available. It is rather large to post here (319 lines, including demo and test code).

James McNeill said...

No, no limitations on the use of the code. Glad it was useful!

Unknown said...
This comment has been removed by the author.
Unknown said...

Hey,

Thanks for this tutorial that was really helpfull. I adapted your algorithm in JavaScript to create a 2D grid games line of sight simulation and it works like a charm. Thanks also to Duke for your contribution, I've tried to do the same but I forgot the n-- statement !

Unknown said...

Nit: the first version has a variable "float t" which is assigned to but never used.

James McNeill said...

Good point! That could be removed from the code.

goiken said...

Wondering: How would I keep track of t while traversing. Instead of `visit(x,y)`, I’d like to have `visit(x,y,dt)` where dt is the length of the line in the pixel…

Buttonius said...

Commenting on Golken's comment...

I would implement that as a two-phase operation.
Phase 1 would build a list of visited grid values (or cell coordinates) with the length traveled in each visited cell.
Phase 2 would determine travel time in each element in the list built in phase 1 and accumulate whatever property you want to accumulate over the traveled path.

The logic of the two phases is simpler when you do one thing at a time. If the travel speed in a grid cell depends on the values of that cell in a way that differs for different travelers, you still have to perform phase 1 once.

James McNeill said...

Golken:

If the dt you want to pass to visit() is the line fraction between entry and exit for the square, I think the first version of the code ought to work. It's computing the t value for the next edge the line will cross; subtracting that from the previous t value will give the fraction of the line inside the square.

goiken said...

So I’d declare `double t_old = 0;` before the for-loop, I `visit(x,y,t - t_old);` and I set `t_old=t;` at the end of the for-loop?

goiken said...

Nah… That doesn’t work. What does seem to work is to move the visit function into the beginning of each case in the for-loop and to `visit(x,y,t_next_vertical -t)` or …_horizontal respectively.

Unknown said...

If you're looking for a lower-rent line drawing algorithm, a sliding point like this might be a good alternative.

void raytrace(double sx, double sy, double dx, double dy) {
const int steps = std::max(abs(sx-dx), abs(sy-dy));

for (int step = 0; step <= steps; step++) {
double x{ sx + (dx-sx) * step/steps };
double y{ sy + (dy-sy) * step/steps };
visit(x, y);
}
}

If you don't want to cover the starting point, switch int step = 1, or if you don't want to cover the last point, switch the condition to step < steps. ☺