Cris’ Image Analysis Blog

J.E. Bresenham published an algorithm in 1965, which is used to draw 2D digital lines. Many people have extended the algorithm to work in 3D, and I’ve even found a website with code for 4D, 5D and 6D versions of the algorithm. However, all the code I come across is way more complex that it needs to be, giving the impression that this is a complicated algorithm. Nothing is further from the truth! Drawing digital lines is the most simple and straight-forward task. I’m going to deviate from the classical thinking that an efficient algorithm should use only integer values. This will make our task trivial.

The first thing we need to realize is that, when walking from the start pixel to the end pixel, there is always one dimension along which we take integer steps. Along the other dimensions sometimes we’ll have to step, ans sometimes we won’t. This is, to me, the key finding in the Bresenham algorithm. This tells us how many pixels will compose the line. If we use more pixels, the line will be thickened at some points, or we’ll use some pixels more than once. If we use fewer pixels, the line will not appear continuous. Based on this 1-pixel step size in one dimension, we can easily calculate the sub-pixel step sizes we have to take in the other dimensions. We will then take these steps, computing floating-point coordinates, until we reach the end point. Each of the points we compute we’ll round to integer pixel coordinates for output. In the original algorithm, instead of computing a floating-point coordinate, an integer coordinate is computed, and an additional integer value accumulates the fractional error. These two approaches are equivalent, except that, when using floating-point values, everything becomes easier to understand.

Let’s write some MATLAB code. We’ll start with two random points to draw a line between (I’m using 3D points here, but you can use any number of dimensions):

p1 = round(rand(1,3)*20);
p2 = round(rand(1,3)*20);

Next, we take a floating-point vector variable, p, which will contain the sub-pixel coordinates of the line. We also compute the step size, s, which is the total distance travelled in each dimension divided by the longest distance. Thus, the vector s contains one element with value 1 (or -1), and other values that are smaller.

p = p1;
d = p2-p1;
N = max(abs(d));
s = d/N;

We need to take N steps of size s, then round each of the floating-point positions found:

disp(p1)
for ii=1:N
   p = p+s;
   disp(round(p))
end

    14     6    19
    13     6    18
    12     6    17
    11     7    16
    10     7    16
     9     7    15
     8     7    14
     7     8    13
     6     8    12
     5     8    11
     4     8    11
     3     9    10
     2     9     9
     1     9     8

There it is! Bresenham’s algorithm, for any number of dimensions, in 9 lines of code.