TODO:
[X]
add tests[X]
bresenham from roguebasin[X]
bresenham from rosettacode[X]
bresenham from wikipedia[X]
lerp line from my site[X]
optimized lerp line from my site[X]
dda from my site[X]
dda from https://en.wikipedia.org/wiki/Digital_differential_analyzer_(graphics_algorithm)[1][ ]
dda from libtcod https://github.com/libtcod/libtcod-fov/blob/wip/src/libtcod-fov/dda.c[2][ ]
bresenham from libtcod[ ]
bresenham from michael abrash (where to find? here[3])
[X]
write a script to run all the tests[X]
write a program to parse the output files[X]
update the non-bresenham programs to handle translation[X]
update all the tests to test the reverse line (for symmetry test)[-]
one-algorithm tests[X]
test correct order of output, A→B should output A first and B last[ ]
test symmetry, A→B and B→A should be reversed output → need the data files to have more data[ ]
test rotation, A→B should be the same as when x,y are rotated 90°, 180°, 270°[ ]
test mirrored, A→B should be the same as when x,y are flipped[ ]
test translation (for non-bresenham), A→B should be the same as (A+P)→(B+P) - P, but this will mainly affect dda for long lines
[ ]
two-algorithm tests[ ]
outputs match, maybe count mismatch steps (or outputs)[ ]
summarize mismatches by making a tiny plot, one pixel per test point[ ]
arrange the tiny plots into a splom-style plot?
[ ]
generate variants of some algorithm with < vs <=, and put them into two-algorithm tests
On the main page I describe line-drawing using linear interpolation. It's simple. It's fast. It generalizes to 3d, hexagons, colors, and other types of spaces. But Bresenham's line algorithm[4] is the famous one. The original paper says:
The algorithm may be programmed without multiplication or division instructions and is efficient with respect to speed of execution and memory utilization.
and Wikipedia says "it uses only integer addition, subtraction, and bit shifting, all of which are very cheap operations in historically common computer architectures".
The label "Bresenham" is used today for a family of algorithms extending or modifying Bresenham's original algorithm.
This bothered me a little bit because there's the original algorithm with separate code for each of eight octants, and there are also simplified versions that collapse some of the octants together. There's a collection of these simplified versions on rosettacode[5] and also on roguebasin[6] (a site for roguelike-games). I have some questions:
- Does Bresenham's algorithm produce the same output as linear interpolation?
- Do the simplified versions produce the same output as original Bresenham's?
- Is Bresenham's line algorithm the fastest algorithm?
- Which of these algorithms are symmetric (line from P to Q the same as Q to P)?
- Which of these algorithms are rotationally consistent (line from P to Q the same as if P,Q are rotated 90°)?
- Which of these algorithms are translationally consistent (line from P to Q the same shape as P+D to Q+D)?
I decided to look at Bresenham's 1965 paper[7] (mirror). The description at the top is:
and … there isn't separate code for the eight octants. I've been wrong all this time. In fact there isn't any code at all. But there's a table of eight octants and how to use that to set up the equations, and then from the equations how to construct an algorithm. So that makes me think it was intended to be eight separate loops. Some of those loops might be combinable. Also, the original paper wasn't about graphics on a screen. It was about how to draw lines on a pen plotter. I think it would be a few more years before people would be drawing raster graphics on screens, although I'm not able to find a reference right now.
Also, Bresenham's 1965 paper says:
This paper is based on “An incremental algorithm for digital plotting,” presented by the author at the ACM National Conference at Denver, Colorado, on August 30, 1963.
However I can't find that 1963 paper.
1 Summary of findings#
- All the line drawing algorithms produce the same output except for rounding at positions exactly halfway between two pixels. One algorithm chooses one way and the other chooses the other way. (This is good)
- The different versions of Bresenham's Algorithm do not all produce the same output.
- Bresenham's probably is as fast as DDA (optimized lerp), but I am not sure. Try downloading line-drawing-benchmark.cpp and run it yourself with c++ -O3.
- [unconfirmed] The only symmetric algorithms are the ones that swap the inputs before running.
- [untested] Don't know the answer yet
- DDA (optimized lerp) has major problems with translational consistency :-( and lerp occasionally has issues.
2 Output#
Let's start with the linear interpolation algorithm from the main page.
function diagonal_distance(p0, p1) { let dx = p1.x - p0.x, dy = p1.y - p0.y; return Math.max(Math.abs(dx), Math.abs(dy)); } function lerp(start, end, t) { return start * (1.0 - t) + t * end; } function lerp_line(p0, p1) { let points = []; let N = diagonal_distance(p0, p1); for (let step = 0; step <= N; step++) { let t = N === 0? 0.0 : step / N; let x = Math.round(lerp(p0.x, p1.x, t)); let y = Math.round(lerp(p0.y, p1.y, t)); points.push({x, y}); } return points; }
And let's compare to DDA line drawing[8], which is a version of the above code transformed with inlining and unrolling optimizations:
function dda_line(p0, p1) { let points = []; let dx = p1.x - p0.x; let dy = p1.y - p0.y; let N = Math.max(Math.abs(dx), Math.abs(dy)); dx /= N; dy /= N; let {x, y} = p0; for (let step = 0; step <= N; step++) { points.push({x: Math.round(x), y: Math.round(y)}); x += dx; y += dy; } return points; }
Since DDA is an optimized version of linear interpolation lines, we should get the same results. Do we?
The answer is … no. That's interesting. They mostly match but DDA uses repeated addition, and that means the floating point error accumulates more than using the original multiplication.
Let's compare the two implementations of Bresenham's on Wikipedia[9]:
function bresenham_wikipedia_1_line(p0, p1) { let points = []; let dx = p1.x - p0.x; let dy = p1.y - p0.y; function plotLineLow(q0, q1) { let dx = q1.x - q0.x; let dy = q1.y - q0.y; let yi = dy < 0 ? -1 : +1; dy = Math.abs(dy); let D = (2 * dy) - dx; let y = q0.y; for (let x = q0.x; x <= q1.x; x++) { points.push({x, y}); if (D > 0) { y += yi; D += 2 * (dy - dx); } else { D += 2 * dy; } } } function plotLineHigh(q0, q1) { let dx = q1.x - q0.x; let dy = q1.y - q0.y; let xi = dx < 0 ? -1 : +1; dx = Math.abs(dx); let D = (2 * dx) - dy; let x = q0.x; for (let y = q0.y; y <= q1.y; y++) { points.push({x, y}); if (D > 0) { x += xi; D += 2 * (dx - dy); } else { D += 2 * dx; } } } if (Math.abs(dy) < Math.abs(dx)) { if (p0.x > p1.x) { plotLineLow(p1, p0); points.reverse(); } else { plotLineLow(p0, p1); } } else { if (p0.y > p1.y) { plotLineHigh(p1, p0); points.reverse(); } else { plotLineHigh(p0, p1); } } return points; } function bresenham_wikipedia_2_line(p0, p1) { let points = []; let dx = Math.abs(p1.x - p0.x); let sx = p0.x < p1.x ? +1 : -1; let dy = -Math.abs(p1.y - p0.y); let sy = p0.y < p1.y ? +1 : -1; let error = dx + dy; let {x, y} = p0; while (true) { points.push({x, y}); if (x === p1.x && y === p1.y) break; let e2 = 2 * error; if (e2 >= dy) { if (x === p1.x) break; error += dy; x += sx; } if (e2 <= dx) { if (y === p1.y) break; error += dx; y += sy; } } return points; }
They don't match either!
How about others? RosettaCode[10] and RogueBasin[11]
function bresenham_rosettacode({x: x0, y: y0}, {x: x1, y: y1}) { // https://rosettacode.org/wiki/Bitmap/Bresenham's_line_algorithm#JavaScript let points = []; var dx = Math.abs(x1 - x0), sx = x0 < x1 ? 1 : -1; var dy = Math.abs(y1 - y0), sy = y0 < y1 ? 1 : -1; var err = (dx>dy ? dx : -dy)/2; while (true) { points.push({x: x0, y: y0}); if (x0 === x1 && y0 === y1) break; var e2 = err; if (e2 > -dx) { err -= dy; x0 += sx; } if (e2 < dy) { err += dx; y0 += sy; } } return points; } function bresenham_roguebasin({x: x0, y: y0}, {x: x1, y: y1}) { // https://roguebasin.com/index.php/Bresenham%27s_Line_Algorithm#JavaScript let points = []; var tmp; var steep = Math.abs(y1-y0) > Math.abs(x1-x0); if(steep){ //swap x0,y0 tmp=x0; x0=y0; y0=tmp; //swap x1,y1 tmp=x1; x1=y1; y1=tmp; } var sign = 1; if(x0>x1){ sign = -1; x0 *= -1; x1 *= -1; } var dx = x1-x0; var dy = Math.abs(y1-y0); var err = ((dx/2)); var ystep = y0 < y1 ? 1:-1; var y = y0; for(var x=x0;x<=x1;x++){ points.push(steep ? {x: y, y: sign*x} : {x: sign*x, y: y}); err = (err - dy); if(err < 0){ y+=ystep; err+=dx; } } return points; } function bresenham_zingl({x: x0, y: y0}, {x: x1, y: y1}) { // https://zingl.github.io/bresenham.html let points = []; let dx = Math.abs(x1-x0), sx = x0<x1 ? 1 : -1; let dy = -Math.abs(y1-y0), sy = y0<y1 ? 1 : -1; let err = dx+dy, e2; /* error value e_xy */ for(;;){ /* loop */ points.push({x: x0, y: y0}); if (x0===x1 && y0===y1) break; e2 = 2*err; if (e2 >= dy) { err += dy; x0 += sx; } /* e_xy+e_x > 0 */ if (e2 <= dx) { err += dx; y0 += sy; } /* e_xy+e_y < 0 */ } return points; } /* Michael Abrash's code * Draws a line in octant 0 or 3 ( |DeltaX| >= DeltaY ). */ function abrash_Octant0(X0, Y0, DeltaX, DeltaY, XDirection, points) { let DeltaYx2; let DeltaYx2MinusDeltaXx2; let ErrorTerm; /* Set up initial error term and values used inside drawing loop */ DeltaYx2 = DeltaY * 2; DeltaYx2MinusDeltaXx2 = DeltaYx2 - ( DeltaX * 2 ); ErrorTerm = DeltaYx2 - DeltaX; /* Draw the line */ points.push({x: X0, y: Y0}); /* draw the first pixel */ while ( DeltaX-- ) { /* See if it’s time to advance the Y coordinate */ if ( ErrorTerm >= 0 ) { /* Advance the Y coordinate & adjust the error term back down */ Y0++; ErrorTerm += DeltaYx2MinusDeltaXx2; } else { /* Add to the error term */ ErrorTerm += DeltaYx2; } X0 += XDirection; /* advance the X coordinate */ points.push({x: X0, y: Y0}); /* draw a pixel */ } } /* Michael Abrash's code * Draws a line in octant 1 or 2 ( |DeltaX| < DeltaY ). */ function abrash_Octant1(X0, Y0, DeltaX, DeltaY, XDirection, points) { let DeltaXx2; let DeltaXx2MinusDeltaYx2; let ErrorTerm; /* Set up initial error term and values used inside drawing loop */ DeltaXx2 = DeltaX * 2; DeltaXx2MinusDeltaYx2 = DeltaXx2 - ( DeltaY * 2 ); ErrorTerm = DeltaXx2 - DeltaY; points.push({x: X0, y: Y0}); /* draw the first pixel */ while ( DeltaY-- ) { /* See if it’s time to advance the X coordinate */ if ( ErrorTerm >= 0 ) { /* Advance the X coordinate & adjust the error term back down */ X0 += XDirection; ErrorTerm += DeltaXx2MinusDeltaYx2; } else { /* Add to the error term */ ErrorTerm += DeltaXx2; } Y0++; /* advance the Y coordinate */ points.push({x: X0, y: Y0}); /* draw a pixel */ } } // Michael Abrash's code function bresenham_abrash({x: X0, y: Y0}, {x: X1, y: Y1}) { // http://www.phatcode.net/res/224/files/html/ch35/35-03.html let points = []; let DeltaX, DeltaY; let Temp; let points_reversed = false; // amitp modification /* Save half the line-drawing cases by swapping Y0 with Y1 and X0 with X1 if Y0 is greater than Y1. As a result, DeltaY is always > 0, and only the octant 0-3 cases need to be handled. */ if ( Y0 > Y1 ) { Temp = Y0; Y0 = Y1; Y1 = Temp; Temp = X0; X0 = X1; X1 = Temp; points_reversed = true; } /* Handle as four separate cases, for the four octants in which Y1 is greater than Y0 */ DeltaX = X1 - X0; /* calculate the length of the line in each coordinate */ DeltaY = Y1 - Y0; if ( DeltaX > 0 ) { if ( DeltaX > DeltaY ) { abrash_Octant0(X0, Y0, DeltaX, DeltaY, 1, points); } else { abrash_Octant1(X0, Y0, DeltaX, DeltaY, 1, points); } } else { DeltaX = -DeltaX; /* absolute value of DeltaX */ if ( DeltaX > DeltaY ) { abrash_Octant0(X0, Y0, DeltaX, DeltaY, -1, points); } else { abrash_Octant1(X0, Y0, DeltaX, DeltaY, -1, points); } } if (points_reversed) points.reverse(); return points; }
TODO: variants of these that use >= instead of > or vice versa
TODO: splom to show all comparisons at once, maybe just a summary of the count; probably should generate this offline
TODO: visualization that compares each one algorithm to the majority output from all the algorithms? not sure this is interesting
3 Symmetry & Rotation#
Take a line from P to Q and calculate the line from Q to P, and see if they hit the same points. Flag the ones that don't. Count how many mismatches there are for each algorithm.
The algorithms that flip P and Q naturally won't have any mismatches.
I'm curious whether adding an ε value like 1e-6 will make a difference here. I should test that.
An extension of left/right symmetry is that the line should be rotateable 90°, and then left/right symmetry is 180°
How about reflection?
4 Translation#
I've assumed that the line drawing algorithms behave identically if both points are translated by a fixed amount. This may not be the case. I can take some lines and translate them and see if they produce the same output.
Notes:
- Linear interpolation has minor changes at different offsets.
- DDA has a lot of changes at higher offsets. It's losing precision.
- Bresenham has no changes at any offsets. This is good.
5 Performance#
In my brief tests, I found optimized linear interpolation (called DDA[12]) to be as fast as Bresenham's line algorithm, but benchmarking is tricky and I don't have high confidence in this. I can believe it to be as fast though, since Bresenham's algorithm was designed for the computers of the 1960s. Floating point was slow and branches were fast. Today's CPUs are different. Floating point is fast and branches are slow. The linear interpolation algorithm has a much tighter inner loop, with no branches. There are other algorithms that claim to be faster than Bresenham's[13]; I can't get their benchmark code to run properly though.
Bresenham also made a second line drawing algorithm, described by demofox[14] and tomforsyth[15].