Elegant/Clean (special case) Straight-line Grid Traversal Algorithm?

Question

I'm dusting off an old project of mine. One of the things it had to do was -- given a Cartesian grid system, and two squares on the grid, find a list of all squares that a line joining the center of those two squares would pass through.

The special case here is that all start and end points are confined to the exact center of squares/cells.

Here are some examples -- with pairs of sample starting and ending points. The shaded squares are the ones that should be returned by the respective function call

removed dead ImageShack link - example

The starting and end points are referred to by the squares they are in. In the above picture, assuming that the bottom left is [1,1], the line on the bottom right would be identified as [6,2] to [9,5].

That is, from the (center of the) square on the sixth column from the left, on the second row from the bottom to the (center of the) square on the ninth column from the left, on the fifth row from the bottom

Which really doesn't seem that complicated. However, I somehow seemed to have found some complex algorithm online and implemented it.

I do recall that it was very, very fast. Like, optimized-for-a-hundreds-or-thousands-of-times-per-frames fast.

Basically, it jumped from border to border of the squares, along the line (the points where the line crosses the grid lines). It knew where the next crossing point was by seeing which crossing point was closer -- a horizontal one or a vertical one -- and moved to that next one.

Which is sort of okay in concept, but the actual implementation turned out to be pretty not-so-pretty, and I'm afraid that the level of optimization might be way too high for what I practically need (I'm calling this traversal algorithm maybe five or six times a minute).

Is there a simple, easy-to-understand, transparent straight-line grid traversal algorithm?

In programmatic terms:

def traverse(start_point,end_point)
  # returns a list of all squares that this line would pass through
end

where the given coordinates identify the squares themselves.

Some examples:

traverse([0,0],[0,4])
# => [0,0], [0,1], [0,2], [0,3], [0,4]
traverse([0,0],[3,2])
# => [0,0], [0,1], [1,1], [2,1], [2,2], [3,2]
traverse([0,0],[3,3])
# => [0,0], [1,1], [2,2], [3,3]

Note that lines that move directly through corners should not include squares on the "wing" of the line.

(Good ol' Bresenham's might work here, but it's a bit backwards from what I want. As far as I know, in order to use it, I'd basically have to apply it to the line and then scan every single square on the grid for true or false. Infeasible -- or at least inelegant -- for large grids)

(I am re-looking into Bresenham, and Bresenham-based algorithms, due to a misunderstanding of mine)

---

For clarification, one possible application of this would be, if I was storing all of my objects in a game inside zones (a grid), and I have a ray, and want to see which objects the ray touches. Using this algorithm, I could test the ray to only the objects that are within the given zones, instead of every object on the map.

The actual use of this in my application is that every tile has an effect associated with it, and an object moves through a given line segment every turn. At every turn, it is necessary to check to see which squares the object has traversed through, and therefore, which effects to apply to the object.

Note that, at this point, the current implementation I have does work. This question is mostly for curiosity's purpose. There has to be a simpler way...somehow...for such a simple problem.

---

What am I looking for exactly? Something conceptually/neat and clean. Also, I've realized that due to what I am exactly specifying, all start and end points are always going to be in the center of squares/cells; so perhaps something that takes advantage of that would be neat as well.

Solution 1

What you want is a particular case of a supercover, which is all of the pixels intersected by a geometric object. The object may be a line or a triangle, and there are generalizations to higher dimensions.

Anyways, here is one implementation for line segments. That page also compares the supercover with the result of Bresenham's algorithm - they are different.
_{(source: free.fr)}

I don't know whether you consider the algorithm there as elegant/clean, but it certainly seems simple enough to adapt the code and move on to other parts of your project.

By the way, your question implies that Bresenham's algorithm is not efficient for large grids. That's not true - it generates only the pixels on the line. You don't have to do a true/false test for every pixel on the grid.

Update 1: I noticed that in the picture there are two 'extra' blue squares that I believe the line does not pass through. One of them is adjacent to the 'h' in 'This algorithm'. I don't know if that reflects an error in the algorithm or the diagram (but see @kikito's comment below).

In general, the 'hard' cases are probably when the line passes exactly through a grid point. I speculate that if you are using floating point that float point error may mess you up in these cases. In other words, algorithms should probably stick to integer arithmetic.

Update 2: Another implementation.

Solution 2

A paper on this topic can be found here. This is in regards to raytracing, but that seems quite relevant to what you are after anyway, and I assume you will be able to work with it a bit.

There is also another paper here, dealing with something similar.

Both of these papers are linked in part 4 of Jakko Bikker's excellent tutorials on raytracing (which also include his source code, so you can browse / examine his implementations).

Solution 3

There's a very easy algorithm for Your problem that runs in linear time:

1. Given two points A and B, determine the intersection points of the line (A, B) with every vertical line of Your grid, that lies within this interval.
2. Insert two special intersection points inside the cells containing A and B at the start/end of the list from point 1.
3. Interpret every two sequent intersection points as the min and max vectors of a axis aligned rectangle and mark all grid cells that lie inside this rectangle (this is very easy (intersection of two axis aligned rects), especially considering, that the rectangle has a width of 1 and therefore occupies only 1 column of Your grid)

+------+------+------+------+
|      |      |      |      |
|      |      | B    *      |
|      |      |/     |      |
+------+------*------+------+
|      |     /|      |      |
|      |    / |      |      |
|      |   /  |      |      |
+------+--/---+------+------+
|      | /    |      |      |
|      |/     |      |      |
|      *      |      |      |
+-----/+------+------+------+
|    / |      |      |      |
*   A  |      |      |      |
|      |      |      |      |
+------+------+------+------+ 

"A" and "B" are the points terminating the line represented by "/". "*" marks intersection points of the line with the grid. The two special intersection points are needed in order to mark the cells that contain A & B and to handle special cases like A.x == B.x

An optimized implementation needs Θ(|B.x - A.x| + |B.y - A.y|) time for a line (A, B). Further, one can write this algorithm to determine intersection points with horizontal grid lines, if that is easier for the implementer.

Update: Border cases

As brainjam correctly points out in his answer, the hard cases are the ones, when a line goes exactly through a grid point. Lets assume such a case occurs and the floating point arithmetic operations correctly return a intersection point with integral coordinates. In this case, the proposed algorithm marks only the correct cells (as specified by the image provided by the OP).

However, floating point errors will occur sooner or later and yield incorrect results. From my opinion even using double won't suffice and one should switch to a Decimal number type. An optimized implementation will perform Θ(|max.x - min.x|) additions on that data type each one taking Θ(log max.y) time. That means in the worst case (the line ((0, 0), (N, N)) with huge N (> 10⁶) the algorithm degrades to a O(N log N) worst case runtime. Even switching between vertical/horizontal grid line intersection detection depending on the slope of the line (A, B) doesn't help in this worst case, but it sure does in the average case - I would only consider to implement such a switch if a profiler yields the Decimal operations to be the bottle neck.

Lastly, I can imagine that some clever people might be able to come up with a O(N) solution that correctly deals with this border cases. All Your suggestions are welcome!

Correction

brainjam pointed out, that a decimal data type is not satisfying even if it can represent arbitrary precision floating point numbers, since, for example, ¹/₃ can't be represented correctly. Therefore one should use a fraction data type, which should be able to handle border cases correctly. Thx brainjam! :)

Comments

[a_m0d]

Where do you want the ends of the line to be - corners of the squares (if so, which ones) or in the centers?

[Justin L.]

@a_m0d thank you for pointing out the ambiguity; I have added a picture for clarification. I mean the centers.

[a_m0d]

thanks for the clarification and the image - makes it a lot clearer

[Drew Hall]

+1 particularly for the observation about Bresenham. Bresenham is the right tool for the job (possibly including the antialiased version).

[Justin L.]

Unfortunately, it appears that the first paper is exactly the algorithm I am already using. It does everything fine; I'd just like to find maybe one that was perhaps a bit more elegant. I couldn't find anything of relevance on the second paper =/ but thanks for the suggestions

[Justin L.]

Thank you for pointing that about about Bresenham. My understanding of it was flawed (I thought it was a way to test if a pixel was a part of a given line). But I'm looking into the alternative algorithm now :)

[Justin L.]

I looked through it; I honestly don't believe that the algorithm in the paper has much elegance compared to the one I already had (see a_m0d's answer), but I'll look into supercover algorithms.

[a_m0d]

@Justin - sorry bout that :( - what was the chances of coming up with exactly the same paper?

[Justin L.]

haha probably small. but it only means that you and i thought along the same lines.

[brainjam]

+1. Nice! I think all of these algorithms may have to be a little careful with floating point errors to avoid trouble when an intersection is exactly on a grid point (but the errors make it slightly off the grid point).

[brainjam]

@Justin L: I agree that the algorithm looks a little messy and lacks a simple explanation. But it avoids floating point (which may mess up cases where the line goes exactly through a grid point).

[Dave O.]

@Justin L. LOL, I just noticed, that the algorithm I proposed is the same, that a_m0d links to in his answer (I skipped the first "here") and that it is the one that You already use. What do You dislike about Your algorithm? I can imagine a beautiful implementation of it - unfortunately not everything can be a one-liner. Maybe some abstraction could help^^ or starting a code golf xD.

[Justin L.]

This is actually conceptually different from the paper in a_m0d's answer; it steps through the squares taking note of both horizontal and vertical axes intersections. But I like the generalization that it must intersect all squares in between the y intersections.

[Kirk Broadhurst]

I would expect that the cases where a line goes exactly through a grid point would have to form a Pythagorean triple. It should then be easy to determine if a line has a chance of crossing a grid point by determining whether it matches a known Pythagorean triple, where the hypotenuse value is no greater than the length of the line. I think this problem is made more obscure by the fact that the lines begin between grid points rather than one grid points.

[brainjam]

@Dave, if you're referring to the Python Decimal type, I don't think it buys you much. 1/3 can't be represented exactly. However, the Python Fraction type should do the trick.

[Dave O.]

@brainjam I'm referring to a pseudo decimal number type which can represent arbitrary precision floating point numbers like java's BigDecimal for example. But yeah, You're right, a fraction type is much better suited for that task as it should be able to handle the border cases.

[kikito]

"I noticed that in the picture there are two 'extra' blue squares that I believe the line does not pass through" I think this is addressed by the author at the end of the page, where he says "Here we have supposed that if the line passes through a corner, both squares are drawn. If you want to remove this, you can simply remove the else part dealing with the corner".

[Karussell]

As @DrewHall already pointed out: a modified version of the antialiased version could be used too: en.wikipedia.org/wiki/Xiaolin_Wu%27s_line_algorithm