Download as Notebook

Bresenham Algorithm

This notebook presents the Bresenham algorithm for lines. This algorithm is used to rasterize a line into a raster and is largely based on the observation that there is always a fast direction. That is, while rasterizing a line, it is sufficient to consider one step into the fast direction and from time to time a step into the slow direction. With simple symmetries, it suffices to implement this algorithm in a single octant and we do this for the case that increasing X is the fast direction and increasing Y is the slow direction.

All other cases are simple reflections of the input and output coordinates.

Setup

In the following snippet, some libraries are loaded, a matrix of size 70x200 pixel is generated and two coordinates p and q are stored between which we want to draw a line segment.

## Setup
import numpy as np;
from  matplotlib import pyplot as plt;
%matplotlib inline
raster = np.zeros(200*70).reshape(70,200)

p = [10,10]
q = [150.0,40.0]

Simple Line Rasterization

The algorithm now proceeds by starting with the Y coordinate of P, looping over all X and tracking the error (using and error of zero at the beginning and adding deltaerr, which is the Y-error per step in X)

If the error grows over 0.5, we do a step and reduce the error by one. Given that we model the pixel center, an error between -0.5 and +0.5 lies in the current pixel, if we exceed it in one direction, we need to step in the slow direction and correct the error.



def line_simple(raster,p,q):
    error=0
    y = p[1]
    deltaerr = float(q[1]-p[1]) / (q[0]-p[0])
    for x in range(int(p[0]),int(q[0])):
        raster[y][x] = 1
        error += deltaerr
        if error >= 0.5:
            y = y + 1
            error = error -1
        #print(x,y)
    return raster

Anti-Aliasing

Now, these lines can have artifacts and look a bit ugly, because the observed thickness of the line for a human eye depends a lot on the angle. A horizontal line is exactly one pixel wide, a diagonal line looks very ragged, because it is sometimes $\sqrt 2$ wide followed by a zero width.

Therefore, we need to apply a technique called anti-aliasing in which all pixels are set with a varying intensity depending on how much of the ideal line of given thickness overlaps the pixel at hand. A simple observation for lines of width one is that only the pixel above and below the naive Bresenham are relevant.

A simplistic (non-aesthetic) option to perform anti-aliasing is given in the following snippet. Find better versions by calculating different values for fragmente!

def line_simple_antialias(raster,p,q):
    error=0
    y = p[1]
    deltaerr = float(q[1]-p[1]) / (q[0]-p[0])
    for x in range(int(p[0]),int(q[0])):
        fragmente = np.array([max(0,-error),abs(1-error),max(0,error)])
        fragmente = fragmente / np.sum(fragmente)
        raster[(y-1):(y+2),x] = fragmente
        
        error += deltaerr
        if error >= 0.5:
            y = y + 1
            error = error -1
    return raster

Circles

For the sake of completeness, we also give a circle rasterization method, which is based on the observation that circles have a valid distinction of slow and fast direction in each octant. Therefore, we actually render a single octant and put pixels by reflections.

Note that this is good for matrix displays, for a vector display where you want to walk around the circle, one would have to loop through all quadrants.

def circle(raster, p, r):
    x = r
    y = 0
    err = r
    raster[p[1]+y,p[0]+x]=2
    while y < x:
        dy = 2*y +1
        y = y+1
        err  = err -dy
        if err< 0:
            dx = 1-2*x
            x = x-1
            err = err -dx
        raster[p[1]+y,p[0]+x]=2
        raster[p[1]-y,p[0]+x]=2
        raster[p[1]+y,p[0]-x]=2
        raster[p[1]-y,p[0]-x]=2
        raster[p[1]+x,p[0]+y]=2
        raster[p[1]-x,p[0]+y]=2
        raster[p[1]+x,p[0]-y]=2
        raster[p[1]-x,p[0]-y]=2
        
    return raster

The Main Program

The following snippet draws a line from p to q (change to enable anti-aliasing if you whish) and a circle and presents the result.


raster = line_simple(raster,p,q)
raster = circle(raster, [30,30],15)

%matplotlib inline
plt.imshow(raster[10:30,10:30], cmap="gray")
plt.show()
plt.imshow(raster, cmap="gray")
plt.show()

png

png