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

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)

Example:

+------+------+------+------+
|      |      |      |      |
|      |      | 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! :)

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.

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