Genius Bresenham Line Algorithm Impl. - How Does It Work?

[rwosdev]

I found this line drawing algorithm (among others) on GitHub: https://gist.github.com/bert/1085538#file-plot_line-c

I've modified it slightly for my own reading preferences and understanding and with some questions:

Code: Select all

		int sx; int sy; // Is S sign or slope? Either way it's the increment or decrement.
		if (x0 < x1) sx = 1; else sx  = -1;
		if (y0 < y1) sy = 1; else sy  = -1;
		
		int dx = iabs(x1 - x0);
		int dy = -iabs(y1 - y0); // Negative absolute distance. It's negative for the algorithm.
		
		int err = dx + dy;
		int e2; // Error * 2
		
		while (true) {
			GFXPlot(gfx, x0, y0, colour);
			if (x0 == x1 && y0 == y1) break;
			
			e2 = (err << 1);
			// This is where the selection happens. If e2 is larger than the negative value in Y delta, increment X.
			// If e2 is smaller than the positive value in X delta, increment Y.
			if (e2 >= dy) { err += dy; x0 += sx; } // err += dy actually subtracts EACH TIME due to Y delta being negative.
			if (e2 <= dx) { err += dx; y0 += sy; }
		}

So I've been trying to wrap my brain around the basics Bresenham's line algorithm, I can understand the floating point variations I've seen (which involve the use of a slope variable), but not this. This implementation is so succinct and simple and works so perfectly, compared to the other ones I've seen, I think it's genius, but I can't fully understand how it determines moving next X, and/or next Y. I kind of get it but not fully.

Could someone with more of a mathematical brain explain it in plain English without formulas/mathematical notation, and possibly the general thinking process behind something like this?

[alexfru]

So, you want a line whose X-axis projection is DX (=x1-x0) pixels, and Y-axis projection is DY (=y1-y0) pixels.
(Cases with DX=0 and DY=0 are trivial. And so are DX=+/-DY.)

Your line's formula is (excluding the trivial cases with division by 0):
y = y0 + (x - x0) * DY/DX.
DY/DX (or DX/DY) may not be an integer, obviously, so it may seem that this kind of formula can't be computed using only integers.
But think about it this way...

Suppose DY<DX (and both are positive).
This means y changes slower than x as you move from (x0,y0) to (x1,y1). IOW, advancing x by 1 results in y advancing by less than 1.
Yet IOW, it looks like there are a bunch of horizontal segments that make up the line, in each of which y is constant in two or more adjacent pixels (constant because your screen is discrete and assumes no fractional coordinates). This is nothing new to you, right?

So, as you render each segment advancing x by 1 at every step, you also advance y by DY/DX. But you don't see a change in y on the screen until you accumulate enough DY/DX's in y (you round/truncate y to an integer when rendering a pixel). After some N pixels in a segment you accumulate N*DY/DX in y and if that value becomes 1 or greater, the change in y becomes visible on the screen. It doesn't matter what N is, you just keep adding to a running sum: DY/DX + DY/DX + ... And when that sum becomes 1 or more, that's when you start a new line segment at y+1 (and at that moment you subtract 1 from the running sum, to restart the process and maintain the line's slope).

So, you kinda need to do DY/DX + DY/DX + ... and test if this sum is >= 1. But you can't do it with integers.
But you can do DY + DY + ... with integers and test if this sum is >= DX. See the trick?

If the line is more vertical than horizontal, you swap the roles of x/DX and y/DY. And you trivially take care of the signs of DX and DY.
There's nothing more to it.

[rwosdev]

Thank you Alex for your reply.

I absolutely see what you're talking about mostly, and I see the accumulation aspect more clearly.

I wrote a floating point version of this that was as basic as it gets (I have little experience of graphing algorithms) but it was always off by a pixel or 2 depending on the Y endpoint.

There's no slope value here though, so what I really don't understand is how the exact sloping is implicit in the use of the delta values and ends up in the correct place every time? Like how would you KNOW to calculate it like that and what fundamentals are there for reference?

[alexfru]

There's no slope value here though, so what I really don't understand is how the exact sloping is implicit in the use of the delta values and ends up in the correct place every time? Like how would you KNOW to calculate it like that and what fundamentals are there for reference?

In the end you have DX steps on the x axis, each adding DY/DX to y. DX*(DY/DX) = DY. So, as your x advances by DX, your y advances by DY. Of course, if you have infinite precision. You don't lose anything here because you aren't doing any imprecise division (DY/DX) or addition (DY/DX + DY/DX + ...). You're doing an equivalent computation with integers and it is exact. (In principle, you can have your ints overflow given appropriate end points, but that's a separate problem)

[simeonz]

Apparently I have forgotten some of Bresenham's algorithm. I cannot remember it being that smart, but I may have forgotten. Alexfru's explanation about the use of exact integer arithmetic for fractional slopes is spot on, but I want to add a little explanation for the error part. Unfortunately, it does involve some math.

First, the algorithm as you can see uses sx and sy to transform the result to the necessary "quadrant". I am going to assume that dx and dy are positive as well. I am going to treat only the case for the line starting at the origin, because the starting point has no significant impact on the computations.

You may want to skip parts of the following if you feel that it is unhelpful. Anyway. We are at a pixel with centerpoint (x0, y0). If this is not the final pixel, adjacent pixel must be selected in the direction of the line slope. There are three choices - with centerpoints (x0+1,y0), (x0,y0+1), (x0+1,y0+1). It is senseless to explain without illustration the argument, but let's assume for the moment that if (x0+0.5,y0+1) (midpoint between (x0,y0+1) and (x0+1,y0+1)) lies between the line and (x0+1,y0+1), then we must choose (x0,y0+1). This geometrically means that the line falls away from (x0+1,y0+1) and closer to (x0,y0+1) and is expressed in the slope (x0+0.5)/(y0+1) being greater then dx/dy. Similarly, when (x0+1,y0+0.5) (midpoint between (x0+1,y0) and (x0+1,y0+1)) lies between the line and the (x0+1,y0+1), then we must choose (x0+1,y0). Again, this means that the line falls away from (x0+1,y0+1) and closer to (x0+1,y0), and is the same as the slope (y0+0.5)/(x0+1) being greater then dy/dx. If neither of these conditions are satisfied, we must choose (x0+1,y0+1). The verbiage can be summarized as follows:

Code: Select all

(x0+0.5)/(y0+1) > dx/dy choose (x0,y0+1)
(y0+0.5)/(x0+1) > dy/dx choose (x0+1,y0)
otherwise choose (x0+1,y0+1)

We can decouple the choice for x0 and y0 from above:

Code: Select all

(x0+0.5)/(y0+1) <= dx/dy choose x0+1
(y0+0.5)/(x0+1) <= dy/dx choose y0+1

We can multiply each side by the denominator of the other side in both inequalities:

Code: Select all

(x0+0.5) * dy <= (y0+1) * dx
(y0+0.5) * dx <= (x0+1) * dy

We can add and subtract 0.5 * dy and 0.5 * dx respectively:

Code: Select all

(x0+1) * dy - 0.5 * dy <= (y0+1) * dx
(y0+1) * dx - 0.5 * dx <= (x0+1) * dy

We can move terms:

Code: Select all

(y0+1) * dx - (x0+1) * dy >= - 0.5 * dy
(y0+1) * dx - (x0+1) * dy <= 0.5 * dx

Multiply by two both sides:

Code: Select all

2 * ((y0+1) * dx - (x0+1) * dy) >= -dy (choose x0+1)
2 * ((y0+1) * dx - (x0+1) * dy) <= dx (choose y0+1)

We need efficient way to compute (y0+1) * dx - (x0+1) * dy. We start with dx - dy for (0, 0), and add dx when incrementing y0, subtract dy when incrementing x0.

[Sik]

The main problem with Bresenham is that often explanations are pretty awful.

Let's assume the line goes from top left to bottom right, being mostly horizontal, the gist of the algorithm would be as follows: first figure out the distance between the two endpoints in each axis

Code: Select all

xdist = x2 - x1
ydist = y2 - y1

then in the loop do this

Code: Select all

x += 1
error += ydist
if error >= xdist
    y += 1
    error -= xdist
end

Yeah, it's just error accumulation, no need for division or anything. Learn the above and you can already figure out the rest yourself for anything that needs interpolation.

--

The issue with Bresenham lines is that there are eight possible combinations (excluding perfectly horizontal or vertical lines which can be handled which a much simpler loop):

• Top left to bottom right, mostly horizontal
• Top left to bottom right, mostly vertical
• Top right to bottom left, mostly horizontal
• Top right to bottom left, mostly vertical
• Bottom left to top right, mostly horizontal
• Bottom left to top right, mostly vertical
• Bottom right to top left, mostly horizontal
• Bottom right to top left, mostly vertical

If you swap the endpoints to ensure it always goes in one direction, you halve the amount of combinations. Getting clever enough with variable use you can reduce even further to just two combinations (mostly horizontal vs mostly vertical) or even just one. The problem is that everybody prefers to teach showing loops for all eight combinations and that makes it look a lot more complex than needed. That may have made sense a couple of decades ago (with x86 only having six usable registers from a compiler's viewpoint), but nowadays the cache is your biggest performance problem and x86 has more registers so throwing in a few more variables is not gonna hurt.

There's also the issue of line clipping at the bitmap boundaries. That's a different topic though.

[simeonz]

In restrospect, what I wrote above doesn't explain how one can come with the idea of using a common term "(y0+1) * dx - (x0+1) * dy" for the decisions to move horizontally and vertically. Honestly, I don't know. You need to sit under the right apple tree I would say. However, if you already know that the choice for any of the three possible pixels can be broken into two choices - for horizontal and vertical movement (which is not trivial), and you know that the latter is based on the horizontal and vertical separation between the diagonal choice (x0+1,y0+1) and the line (in relation to half a pixel width/height), then you can actually reason that transforming the space can help you. Namely, by stretching the plane anisotropically - by factor -dy in the x direction and factor dx in the y direction, the drawn line actually will become after the transformation a bisector of the right angle between the -x and y axes. That is, in the second quadrant. For such a line, the horizontal and vertical separation of any arbitrary point to it are equal to each other. Hence, the two choices - for vertical and horizontal motion of the algorithm - can reuse most of the computation in the new space, as long as you don't compare to half a pixel width and height anymore, but to half a -dy and half a dx, which is the same, but in the transformed plane.

P.S.: I suspect that the error term can be used for computing an alpha value to perform approximate anti-aliasing. Note that Bresenham's algorithm does not produce lines of correct thickness, which can be determined by noticing that it draws max(dx, dy) pixels for any line, which is less than the length of the line. And the area covered by a line of thickness 1 pixel is proportional to its length, not its horizontal or vertical projection. Thus, it is intended for rough and fast display.

[BenLunt]

Code: Select all

		while (true) {
			GFXPlot(gfx, x0, y0, colour);
			if (x0 == x1 && y0 == y1) break;

(Sorry for the late reply, I have been busy)
I don't have a comment about the technique of the line drawing, but actually, the ending of the technique. Checking if current pixel location == ending is not going to work for certain lengths and angles. For example, imagine that the following is a nearly horz line, S being the starting pixel, E being the ending position, 'x' being the pixels plotted between them.

Code: Select all

  Sxxxxx
        xxxxxx
              xxxxxxE
                    xxxxxx

Notice that the 'Y' coord is incremented just at the ending point. You might want to check the surrounding pixel(s) as the END pixel.

- Ben

[Sik]

Or just check the axis you're iterating through (which is a single comparison), since you'll always go through every value in it exactly once.

[simeonz]

Notice that the 'Y' coord is incremented just at the ending point. You might want to check the surrounding pixel(s) as the END pixel.

That's a good catch. I'm thinking however that it is a bug essentially. Drawing a line from A to B is expected to have both A and B included conceptually (even if B is not drawn for usage in polylines). I believe that making the error inequalities strict (that is "e2 > dy" and "e2 < dx") should fix that issue.

P.S. I am starting to think that even with the current inequalities this shouldn't happen at all. For one, it makes no sense geometrically, as the line we are drawing is still way above the midline separating the two raster rows that we are considering (the one of the ending point and the one after it). Also, note that at the starting and ending point "err" should hold "dx+dy" exactly, so at the point horizontally to the left of the ending point, err should hold "dy" exactly. Meaning that (e2 >= dy), which is the same as (dy * 2 >= dy) should not pass as dy is negative.

As another footnote, I believe that x << y is undefined behavior for negative x, so err << 1 should be replaced with err * 2 if portability is imperative, and the compiler will most likely optimize it to shift anyway (; unless there is a compiler optimization bug that is).

[simeonz]

Note - this is the actual output from the original code without changes to the inequalities (generated from python admittedly). It works as expected, considering that the algorithm aligns the starting and ending coordinates on pixel centerpoints and not verteces. Hence the rather atypical pattern.

The code does not really use a fixed stride. It evaluates where pixel edge midpoints reside relative to the drawn line. But the fixed stride picture is a good intuition for how integer arithmetic can be used, despite the involvement of fractional slope.

[Schol-R-LEA]

I seem to recall that the now-ancient first edition of Foley and van Dam I studied ages ago had a really great explanation, but since I haven't read it since, oh, 1993 or so, I can't recall it very well. There is a new 3rd edition, however, from just three years ago, and I am guessing that part won't really have changed - though the prices for it are rather steep.

I should also mention that there are other line drawing algorithms which actually perform better on average than the original version of Bresenham (though most of these are based on it), and others (such as Wu's) which do better in terms of appearance.

While algorithms such as Wu's algorithm are also frequently used in modern computer graphics because they can support antialiasing, the speed and simplicity of Bresenham's line algorithm means that it is still important. The algorithm is used in hardware such as plotters and in the graphics chips of modern graphics cards. It can also be found in many software graphics libraries. Because the algorithm is very simple, it is often implemented in either the firmware or the graphics hardware of modern graphics cards.

The Bresenham algorithm has held on because, a) it was the first really efficient line algorithm, so everyone ended up hearing about it, b) it is indeed very, very fast, and had minimal memory costs, so it worked well on early home computers in the 1970s and 1980s, c) it is a simple enough algorithm that even someone who doesn't completely understand it can muddle through implementing it, d) it is well-suited to hardware implementations, as well as ones in ROM or some other type of firmware, and e) it can be generalized to work with circles, ellipses, and several other types of curves.

The Wu algorithm is apparently a bit more complex, though this explanation might help if you want to look at it, too.

The Wicked-Pedo article on line-drawing algorithms in general also mentions the Gupta-Sproull algorithm, but apparently there isn't a separate article on that. It is easy enough to find a few pages on it via a web search, however, such as this one and this one, and this paper discusses it's efficiency, but my impression is that for line drawing, it isn't any better than Wu's, and possibly worse. However, from what I've read, the line algorithm is actually just a specific application of a general anti-aliasing algorithm, so if you are doing more than just drawing lines, it might make more sense to use that algorithm (or some similar one) regardless of performance.

Dunno if any of this helps you understand the Bresenham algo, though.

I will point out that, if your goal is to implement for your OS, you might be better off focusing on your GPU driver instead - as the Wicked-Pedo article quoted earlier indicates, almost all of the major GPUs in widespread use now (which AFAICT means both discrete GPUs such as Nvidia and AMD, and integrated ones such as Intel's or the various SoC ones) implement a version of the algorithm natively, and are likely to be a lot faster with it than would be possible at all using CPU-bound code.