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

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.

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:

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

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:

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:

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

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.

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!

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

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

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

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?

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.

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.

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.

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:

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.

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 expectationsbool 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;

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.

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!

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.

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?

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.

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.

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:

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

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 !

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.

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.