[SciPy-User] Bresenham algorithm?

[Sebastian Haase]

On Sun, Sep 20, 2009 at 5:30 PM, David Warde-Farley <dwf at cs.toronto.edu> wrote:
> On 20-Sep-09, at 12:27 PM, Xavier Gnata wrote:
>
>> Does scipy provide the Bresenham algorithm ?
>> For sure it is simple to implement this algo in python but it is
>> slow :(
>
> It doesn't, that I know of, but it would be easy enough to use Cython
> to speed it up. This is precisely the sort of thing Cython is good
> for: iterative algorithms that cannot easily be vectorized. It'll be
> particularly fast since Bresenham's algorithm is all in integer
> arithmetic (note that you're still dealing with ints in Python).
>
> Here's an implementation I found here: http://mail.python.org/pipermail/python-list/1999-July/007163.html
>
> def bresenham(x,y,x2,y2):
>     """Brensenham line algorithm"""
>     steep = 0
>     coords = []
>     dx = abs(x2 - x)
>     if (x2 - x) > 0: sx = 1
>     else: sx = -1
>     dy = abs(y2 - y)
>     if (y2 - y) > 0: sy = 1
>     else: sy = -1
>     if dy > dx:
>         steep = 1
>         x,y = y,x
>         dx,dy = dy,dx
>         sx,sy = sy,sx
>     d = (2 * dy) - dx
>     for i in range(0,dx):
>         if steep: coords.append((y,x))
>         else: coords.append((x,y))
>         while d >= 0:
>             y = y + sy
>             d = d - (2 * dx)
>         x = x + sx
>         d = d + (2 * dy)
>     return coords #added by me
>
> Now here's the version I've marked up with Cython:
>
> cdef extern from "math.h":
>     int abs(int i)
>
> def bresenham(int x, int y, int x2, int y2):
>     cdef int steep = 0
>     cdef int dx = abs(x2 - x)
>     cdef int dy = abs(y2 - y)
>     cdef int sx, sy, d, i
>     coords = []
>     if (x2 - x) > 0: sx = 1
>     else: sx = -1
>     if (y2 - y) > 0: sy = 1
>     else: sy = -1
>     if dy > dx:
>         steep = 1
>         x,y = y,x
>         dx,dy = dy,dx
>         sx,sy = sy,sx
>     d = (2 * dy) - dx
>     for i in range(0,dx):
>         if steep:
>             coords.append((y,x))
>         else:
>             coords.append((x,y))
>         while d >= 0:
>             y = y + sy
>             d = d - (2 * dx)
>         x = x + sx
>         d = d + (2 * dy)
>     return coords
>
>
> And here is the speed comparison:
>
> In [1]: from bresenham import bresenham as bresenham_py
>
> In [2]: from bresenham_cython import bresenham as bresenham_cy
>
> In [3]: a = bresenham_py(0, 0, 12900, 10500)
>
> In [4]: b = bresenham_cy(0, 0, 12900, 10500)
>
> In [5]: a == b # Check that they produce the same results
> Out[5]: True
>
> In [6]: timeit bresenham_py(0, 0, 12900, 10500)
> 100 loops, best of 3: 12.6 ms per loop
>
> In [7]: timeit bresenham_cy(0, 0, 12900, 10500) # python was already
> pretty fast
> 100 loops, best of 3: 2.27 ms per loop
>
> So, with that minimal effort and the compilation step we've already
> got a 5x speedup. Note that my Cython'ed version still uses a Python
> list. It could be easily modified so that you could pass in your array
> and do the summing inside this function, bypassing the need for any
> Python API calls at all, and gaining further speed ups.
>
> If I comment out the lines involving 'coords' (the list/tuples
> manipulation)
>
> In [7]: timeit bresenham_cy(0, 0, 12900, 10050)
> 10000 loops, best of 3: 28.7 us per loop
>
> which is almost an additional hundred-fold speedup.
>
> For information on how to pass in ndarray arguments into a Cython
> function, see:
>
>        http://wiki.cython.org/tutorials/numpy
>
> David

If you comment out the list/tuples handling - assuming you want to
plug in numpy arrays here, you have to know the number of points up
front.
(To be able to pre-allocate the array correctly)
Is this easy ?  (rounding errors ?)

Cheers,

- Sebastian Haase