Bresenham's line algorithm

9 Aug 2023

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^[1] is the famous one. As Wikipedia notes:

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^[2] and also on roguebasin^[3] (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)?

I decided to look at Bresenham's original 1965 paper^[4] (mirror) 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.

1 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^[5], 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?

Blue vs Orange

Do lerp and dda match?

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:

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?

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

2 Symmetry#

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.

3 Numeric#

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.

I think DDA in particular will behave differently.

TODO: draw a line with 0,0 origin and also 1000000,1000000 origin and see if they take the same steps. Also try -1000000,-1000000

4 Performance#

In my brief tests, I found optimized linear interpolation (called DDA^[6]) 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^[7]; I can't get their benchmark code to run properly though.

1 Output#

2 Symmetry#

3 Numeric#

4 Performance#

Endnotes