diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index d5a71f60..fda2d646 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -63,3 +63,6 @@
 
 - Brian Holt
   Histograms of Oriented Gradients
+
+- David-Warde Farley, Sturla Molden
+  Bresenheim line drawing, from snippets on numpy-discussion.
diff --git a/TASKS.txt b/TASKS.txt
index 8c99a40a..1d32b373 100644
--- a/TASKS.txt
+++ b/TASKS.txt
@@ -122,3 +122,7 @@ Implement Algorithms
 - Graph cut segmentation
 - Probabilistic Hough transform
 
+Drawing
+```````
+- Wu's algorithm for lines and circles
+
diff --git a/scikits/image/draw/__init__.py b/scikits/image/draw/__init__.py
new file mode 100644
index 00000000..907a1ba7
--- /dev/null
+++ b/scikits/image/draw/__init__.py
@@ -0,0 +1 @@
+from draw import *
diff --git a/scikits/image/draw/_draw.pyx b/scikits/image/draw/_draw.pyx
new file mode 100644
index 00000000..a586ae03
--- /dev/null
+++ b/scikits/image/draw/_draw.pyx
@@ -0,0 +1,66 @@
+import numpy as np
+cimport numpy as np
+cimport cython
+
+cdef extern from "math.h":
+    int abs(int i)
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def bresenham(int y, int x, int y2, int x2):
+    """
+    Generate line pixel coordinates.
+    
+    Parameters
+    ----------
+    y, x : int
+        Starting position (row, column).
+    y2, x2 : int
+        End position (row, column).
+
+    Returns
+    -------
+    rr, cc : (N,) ndarray of int
+        Indices of pixels that belong to the line.
+        May be used to directly index into an array, e.g.
+        ``img[rr, cc] = 1``.
+
+    """
+    cdef np.ndarray[np.int32_t, ndim=1, mode="c"] rr, cc
+
+    cdef int steep = 0
+    cdef int dx = abs(x2 - x)
+    cdef int dy = abs(y2 - y)
+    cdef int sx, sy, d, i
+
+    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
+
+    rr = np.zeros(int(dx) + 1, dtype=np.int32)
+    cc = np.zeros(int(dx) + 1, dtype=np.int32)
+    
+    for i in range(dx):
+        if steep:
+            rr[i] = x
+            cc[i] = y
+        else:
+            rr[i] = y
+            cc[i] = x
+        while d >= 0:
+            y = y + sy
+            d = d - (2 * dx)
+        x = x + sx
+        d = d + (2 * dy)
+
+    rr[dx] = y2
+    cc[dx] = x2
+
+    return rr, cc
diff --git a/scikits/image/draw/draw.py b/scikits/image/draw/draw.py
new file mode 100644
index 00000000..9360623a
--- /dev/null
+++ b/scikits/image/draw/draw.py
@@ -0,0 +1,6 @@
+"""
+Methods to draw on arrays.
+
+"""
+
+from _draw import bresenham
diff --git a/scikits/image/draw/setup.py b/scikits/image/draw/setup.py
new file mode 100644
index 00000000..3e0cd1a4
--- /dev/null
+++ b/scikits/image/draw/setup.py
@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+
+import os
+from scikits.image._build import cython
+
+base_path = os.path.abspath(os.path.dirname(__file__))
+
+def configuration(parent_package='', top_path=None):
+    from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
+
+    config = Configuration('draw', parent_package, top_path)
+    config.add_data_dir('tests')
+
+    cython(['_draw.pyx'], working_path=base_path)
+
+    config.add_extension('_draw', sources=['_draw.c'],
+                         include_dirs=[get_numpy_include_dirs()])
+
+    return config
+
+if __name__ == '__main__':
+    from numpy.distutils.core import setup
+    setup(maintainer = 'Scikits-image developers',
+          author = 'Scikits-image developers',
+          maintainer_email = 'scikits-image@googlegroups.com',
+          description = 'Drawing',
+          url = 'https://github.com/scikits-image/scikits.image',
+          license = 'SciPy License (BSD Style)',
+          **(configuration(top_path='').todict())
+          )
diff --git a/scikits/image/draw/tests/test_draw.py b/scikits/image/draw/tests/test_draw.py
new file mode 100644
index 00000000..06ebaac0
--- /dev/null
+++ b/scikits/image/draw/tests/test_draw.py
@@ -0,0 +1,52 @@
+from numpy.testing import assert_array_equal
+import numpy as np
+
+from scikits.image.draw import bresenham
+
+def test_bresenham_horizontal():
+    img = np.zeros((10, 10))
+
+    rr, cc = bresenham(0, 0, 0, 9)
+    img[rr, cc] = 1
+
+    img_ = np.zeros((10, 10))
+    img_[0, :] = 1
+
+    assert_array_equal(img, img_)
+
+def test_bresenham_vertical():
+    img = np.zeros((10, 10))
+
+    rr, cc = bresenham(0, 0, 9, 0)
+    img[rr, cc] = 1
+
+    img_ = np.zeros((10, 10))
+    img_[:, 0] = 1
+
+    assert_array_equal(img, img_)
+
+def test_reverse():
+    img = np.zeros((10, 10))
+
+    rr, cc = bresenham(0, 9, 0, 0)
+    img[rr, cc] = 1
+
+    img_ = np.zeros((10, 10))
+    img_[0, :] = 1
+
+    assert_array_equal(img, img_)
+
+def test_diag():
+    img = np.zeros((5, 5))
+
+    rr, cc = bresenham(0, 0, 4, 4)
+    img[rr, cc] = 1
+
+    img_ = np.eye(5)
+
+    assert_array_equal(img, img_)
+
+
+if __name__ == "__main__":
+    from numpy.testing import run_module_suite
+    
diff --git a/scikits/image/feature/hog.py b/scikits/image/feature/hog.py
index 8826ddba..53fafa40 100644
--- a/scikits/image/feature/hog.py
+++ b/scikits/image/feature/hog.py
@@ -1,51 +1,43 @@
-"""
-:author: Brian Holt, 2011
-:license: modified BSD
-"""
-
 import numpy as np
 from scipy import sqrt, pi, arctan2, cos, sin
+from scipy.ndimage import uniform_filter
+
+# XXX Replace with integral after merge
 from ..transform import sat_sum
 
-def hog(image, n_orientations=9, pixels_per_cell=(8, 8),
+def hog(image, orientations=9, pixels_per_cell=(8, 8),
         cells_per_block=(3, 3), visualise=False, normalise=False):
-    """  Extract Histogram of Oriented Gradients (HOG) for a given image.
+    """Extract Histogram of Oriented Gradients (HOG) for a given image.
 
     Compute a Histogram of Oriented Gradients (HOG) by
         1) (optional) global image normalisation
         2) computing the gradient image in x and y
         3) computing gradient histograms
-        3) normalise across blocks
-        4) flatten into a feature vector
+        3) normalising across blocks
+        4) flattening into a feature vector
 
     Parameters
     ----------
-    image: ndarray, 2D
-        2D image (greyscale)
-
-    n_orientations  : int
-        number of orientation bins
-
-    pixels_per_cell  : 2 tuple (int,int)
-        pixels per cell, size in pixels of a cell
-
+    image : (M, N) ndarray
+        Input image (greyscale).
+    orientations : int
+        Number of orientation bins.
+    pixels_per_cell : 2 tuple (int, int)
+        Size (in pixels) of a cell.
     cells_per_block  : 2 tuple (int,int)
-        cells per block, number of cells in each block
-
+        Number of cells in each block.
     visualise : bool, optional
-        return an image of the HOG
-
+        Also return an image of the HOG.
     normalise : bool, optional
-        apply power law compression to normalise the image before
-        processing
+        Apply power law compression to normalise the image before
+        processing.
 
     Returns
     -------
     newarr : ndarray
         HOG for the image as a 1D (flattened) array.
-
-    hog_image : PIL Image (if visualise=True)
-        A visualisation of the HOG image
+    hog_image : ndarray (if visualise=True)
+        A visualisation of the HOG image.
 
     References
     ----------
@@ -54,8 +46,8 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8),
     * Dalal, N and Triggs, B, Histograms of Oriented Gradients for
       Human Detection, IEEE Computer Society Conference on Computer
       Vision and Pattern Recognition 2005 San Diego, CA, USA
-    """
 
+    """
     image = np.atleast_2d(image)
 
     """
@@ -68,9 +60,9 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8),
     shadowing and illumination variations.
     """
 
-    if image.ndim == 3:
-        # replace RGB with locally dominant colour channel
-        pass  # TODO
+    if image.ndim > 3:
+        raise ValueError("Currently only supports grey-level images")
+
     if normalise:
         image = sqrt(image)
 
@@ -91,40 +83,22 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8),
 
     """
     The third stage aims to produce an encoding that is sensitive to
-    local image content while remaining resistant to small changes in pose
-    or appearance. The adopted method pools gradient orientation information
-    locally in the same way as the SIFT [Lowe 2004] feature. The image window
-    is divided into small spatial regions, called "cells". For each cell we
-    accumulate a local 1-D histogram of gradient or edge orientations over
-    all the pixels in the cell. This combined cell-level 1-D histogram
-    forms the basic "orientation histogram" representation. Each orientation
-    histogram divides the gradient angle range into a fixed number of
-    predetermined bins. The gradient magnitudes of the pixels in the cell
-    are used to vote into the orientation histogram.
+    local image content while remaining resistant to small changes in
+    pose or appearance. The adopted method pools gradient orientation
+    information locally in the same way as the SIFT [Lowe 2004]
+    feature. The image window is divided into small spatial regions,
+    called "cells". For each cell we accumulate a local 1-D histogram
+    of gradient or edge orientations over all the pixels in the
+    cell. This combined cell-level 1-D histogram forms the basic
+    "orientation histogram" representation. Each orientation histogram
+    divides the gradient angle range into a fixed number of
+    predetermined bins. The gradient magnitudes of the pixels in the
+    cell are used to vote into the orientation histogram.
     """
 
     magnitude = sqrt(gx ** 2 + gy ** 2)
     orientation = arctan2(gy, (gx + 1e-15)) * (180 / pi) + 90
 
-    # compute n_orientations integral images
-    integral_images = []
-    for i in range(0, n_orientations):
-        #create new integral image for this orientation
-        # isolate orientations in this range
-
-        temp_ori = np.where(orientation < 180 / n_orientations * (i + 1),
-                            orientation, 0)
-        temp_ori = np.where(orientation >= 180 / n_orientations * i,
-                            temp_ori, 0)
-        # select magnitudes for those orientations
-        cond2 = temp_ori > 0
-        temp_mag = np.where(cond2, magnitude, 0)
-
-        #compute integral image
-        integral = np.cumsum(np.cumsum(temp_mag, axis=0, dtype=float),
-                axis=1, dtype=float)
-        integral_images.append(integral)
-
     sx, sy = image.shape
     cx, cy = pixels_per_cell
     bx, by = cells_per_block
@@ -132,34 +106,43 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8),
     n_cellsx = int(np.floor(sx // cx))  # number of cells in x
     n_cellsy = int(np.floor(sy // cy))  # number of cells in y
 
+    # compute orientations integral images
+    orientation_histogram = np.zeros((n_cellsx, n_cellsy, orientations))
+    for i in range(orientations):
+        #create new integral image for this orientation
+        # isolate orientations in this range
+
+        temp_ori = np.where(orientation < 180 / orientations * (i + 1),
+                            orientation, 0)
+        temp_ori = np.where(orientation >= 180 / orientations * i,
+                            temp_ori, 0)
+        # select magnitudes for those orientations
+        cond2 = temp_ori > 0
+        temp_mag = np.where(cond2, magnitude, 0)
+
+        orientation_histogram[:,:,i] = uniform_filter(temp_mag, size=(cx, cy))[cx/2::cx, cy/2::cy].T
+
+
     # now for each cell, compute the histogram
-    orientation_histogram = np.zeros((n_cellsx, n_cellsy, n_orientations))
+    #orientation_histogram = np.zeros((n_cellsx, n_cellsy, orientations))
 
     radius = min(cx, cy) // 2 - 1
     hog_image = None
     if visualise:
-        import Image
-        import ImageDraw
-        hog_image = Image.new("F", (sy, sx))
-        draw = ImageDraw.Draw(hog_image)
+        hog_image = np.zeros((sy, sx), dtype=float)
 
-    for x in range(0, n_cellsx):
-        for y in range(0, n_cellsy):
-            for o in range(0, n_orientations):
-                # compute the histogram from integral image
-                orientation_histogram[x, y, o] = sat_sum(integral_images[o],
-                                                         y * cy,
-                                                         x * cx,
-                                                         (y + 1) * cy - 1,
-                                                         (x + 1) * cx - 1)
-
-                if visualise:
+    if visualise:
+        from scikits.image import draw
+        
+        for x in range(n_cellsx):
+            for y in range(n_cellsy):
+                for o in range(orientations):
                     centre = tuple([y * cy + cy // 2, x * cx + cx // 2])
-                    dx = radius * cos(float(o) / n_orientations * np.pi)
-                    dy = radius * sin(float(o) / n_orientations * np.pi)
-                    draw.line([(centre[0] - dx, centre[1] - dy),
-                               (centre[0] + dx, centre[1] + dy)],
-                              fill=orientation_histogram[x, y, o])
+                    dx = radius * cos(float(o) / orientations * np.pi)
+                    dy = radius * sin(float(o) / orientations * np.pi)
+                    rr, cc = draw.bresenham(centre[0] - dx, centre[1] - dy,
+                                            centre[0] + dx, centre[1] + dy)
+                    hog_image[rr, cc] += orientation_histogram[x, y, o]
 
     """
     The fourth stage computes normalisation, which takes local groups of
@@ -179,18 +162,18 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8),
     n_blocksx = (n_cellsx - bx) + 1
     n_blocksy = (n_cellsy - by) + 1
     normalised_blocks = np.zeros((n_blocksx, n_blocksy,
-                                  bx, by, n_orientations))
+                                  bx, by, orientations))
 
-    for x in range(0, n_blocksx):
-        for y in range(0, n_blocksy):
+    for x in range(n_blocksx):
+        for y in range(n_blocksy):
             block = orientation_histogram[x:x + bx, y:y + by, :]
             eps = 1e-5
             normalised_blocks[x, y, :] = block / sqrt(block.sum() ** 2 + eps)
 
     """
-     The final step collects the HOG descriptors from all blocks of a dense
-     overlapping grid of blocks covering the detection window into a combined
-     feature vector for use in the window classifier
+    The final step collects the HOG descriptors from all blocks of a dense
+    overlapping grid of blocks covering the detection window into a combined
+    feature vector for use in the window classifier.
     """
 
     if visualise:
diff --git a/scikits/image/feature/tests/test_hog.py b/scikits/image/feature/tests/test_hog.py
index b37e0e98..6f375cea 100644
--- a/scikits/image/feature/tests/test_hog.py
+++ b/scikits/image/feature/tests/test_hog.py
@@ -4,10 +4,12 @@ import scipy
 from scikits.image.feature import hog 
 
 def test_histogram_of_oriented_gradients():
-    img = scipy.lena().astype(np.int8) 
+    # Replace with scikits.image.data.lena() after merge
+    img = scipy.misc.lena().astype(np.int8) 
     
-    fd = hog(img, n_orientations=9, pixels_per_cell=(8, 8), 
+    fd = hog(img, orientations=9, pixels_per_cell=(8, 8), 
              cells_per_block=(1, 1))
+
     assert len(fd) == 9 * (512//8) ** 2
     
 if __name__ == '__main__':