My implementation of Bresenham's algorithm fails for lines at certain angles

Question

I've written an implementation of Bresenham's algorithm in Python (following the Wikipedia article), and it works correctly except for lines at certain angles. All lines that should extend between 45 and 90 degrees, or between 135 and 270 degrees, will instead extend along the line y = x.

Here's my code:

def bresenham(origin, dest):
    # debug code
    print origin
    print dest
    # end debug code
    x0 = origin[0]; y0 = origin[1]
    x1 = dest[0]; y1 = dest[1]
    steep = abs(y1 - y0) > abs(x1 - x0)
    backward = x0 > x1

    if steep:
        x0, y0 = y0, x0
        x1, y1 = y1, x1
    if backward:
        x0, x1 = x1, x0
        y0, y1 = y1, y0

    dx = x1 - x0
    dy = abs(y1 - y0)
    error = dx / 2
    y = y0

    if y0 < y1: ystep = 1 
    else: ystep = -1

    result = []
    #if x0 > x1: xstep = -1
    #else: xstep = 1
    # debug code
    print "x0 = %d" % (x0)
    print "x1 = %d" % (x1)
    print "y0 = %d" % (y0)
    print "y1 = %d" % (y1)
    for x in range(x0, x1):
        if steep: result.append((y,x))
        else: result.append((x,y))
        error -= dy
        if error < 0:
            y += ystep
            error += dx 
    # ensure the line extends from the starting point to the destination
    # and not vice-versa
    if backward: result.reverse()
    print result
    return result

Anyone see what I'm screwing up?

---

EDIT:

I added some printing code to the function.

(0,0) is at the top left of the display.

My test framework is pretty simple. It's a standalone function, so I just pass two points to it:

origin = (416, 384)
dest = (440, 347)
bresenham(origin, dest)
(416, 384)
(440, 347)
x0 = 384
x1 = 347
y0 = 416
y1 = 440
[]

Answer 1

I don't know why you're using an xstep variable. You don't really need one with the algorithm you're using.

@Gabe: xstep is needed because without it, if x0 > x1, then the for loop will terminate immediately, as the default step for a Python for loop is 1.

The reason you don't need an xstep variable is because, if it's going backwards, the coordinates were already switched (in the if backward: conditional at the beginning) so that the end-point is now the start-point and vice-versa, such that we now are still going left-to-right.

You just need this:

result = []

for x in range(x0, x1):
    if steep: result.append((y, x))
    else: result.append((x, y))
    error -= dy
    if error < 0:
        y += ystep
        error += dx

return result

If you want the list of coordinates in order from start to end-point, then you can do the check at the end:

if backward: return result.reverse()
else: return result

EDIT: The problem is that the backward boolean is being evaluated before it needs to be. If the steep conditional executes, then the values change, but by then your backward conditional is different. To fix this, instead of using a backward boolean, make it an explicit expression:

if x0 > x1:
    # swapping here

Then again, since you end up using the boolean later on, you could just define it before the conditional:

backward = x0 > x1

if backward:

Answer 2

The problem is that you are computing x0 > x1 before you swap x and y.

Instead of:

backward = x0 > x1 

if steep: 
    x0, y0 = y0, x0 
    x1, y1 = y1, x1 
if backward: 
    x0, x1 = x1, x0 
    y0, y1 = y1, y0 

You should have:

if steep: 
    x0, y0 = y0, x0 
    x1, y1 = y1, x1 

backward = x0 > x1 
if backward: 
    x0, x1 = x1, x0 
    y0, y1 = y1, y0 

Answer 3

Here is an implementation I am using, which is much faster than pure Python thanks to numba

from numba import njit
from typing import Tuple
import numpy as np


@njit
def bresenham_numba(x0, y0, x1, y1, w, h):
    dx = abs(x1 - x0)
    dy = abs(y1 - y0)
    sx = -1 if x0 > x1 else 1
    sy = -1 if y0 > y1 else 1
    err = dx - dy
    points = [(x0, y0)]
    while x0 != x1 or y0 != y1:
        e2 = err * 2
        if e2 > -dy:
            err -= dy
            x0 += sx
        if e2 < dx:
            err += dx
            y0 += sy
        points.append((x0, y0))

    points = [p for p in points if 0 <= p[0] < h and 0 <= p[1] < w]
    return points


def calculate_bresenham(start_point: Tuple[int, int], end_point: Tuple[int, int], resolution: Tuple[int, int]):
    w, h = resolution
    x0, y0 = start_point
    x1, y1 = end_point
    bresenham_list = bresenham_numba(x0, y0, x1, y1, w, h)
    return np.array(bresenham_list)