diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index 5feb5777..43c559b8 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -68,3 +68,8 @@
 - Andreas Mueller
   Example data set loader.
 
+- 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/__init__.py b/scikits/image/feature/__init__.py
new file mode 100644
index 00000000..eb6faa62
--- /dev/null
+++ b/scikits/image/feature/__init__.py
@@ -0,0 +1 @@
+from hog import hog
\ No newline at end of file
diff --git a/scikits/image/feature/hog.py b/scikits/image/feature/hog.py
new file mode 100644
index 00000000..8b71f8a3
--- /dev/null
+++ b/scikits/image/feature/hog.py
@@ -0,0 +1,193 @@
+import numpy as np
+from scipy import sqrt, pi, arctan2, cos, sin
+
+# XXX Replace with integral after merge
+from ..transform import sat_sum
+
+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.
+
+    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) normalising across blocks
+        4) flattening into a feature vector
+
+    Parameters
+    ----------
+    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)
+        Number of cells in each block.
+    visualise : bool, optional
+        Also return an image of the HOG.
+    normalise : bool, optional
+        Apply power law compression to normalise the image before
+        processing.
+
+    Returns
+    -------
+    newarr : ndarray
+        HOG for the image as a 1D (flattened) array.
+    hog_image : ndarray (if visualise=True)
+        A visualisation of the HOG image.
+
+    References
+    ----------
+    * http://en.wikipedia.org/wiki/Histogram_of_oriented_gradients
+
+    * 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)
+
+    """
+    The first stage applies an optional global image normalisation
+    equalisation that is designed to reduce the influence of illumination
+    effects. In practice we use gamma (power law) compression, either
+    computing the square root or the log of each colour channel.
+    Image texture strength is typically proportional to the local surface
+    illumination so this compression helps to reduce the effects of local
+    shadowing and illumination variations.
+    """
+
+    if image.ndim > 3:
+        raise ValueError("Currently only supports grey-level images")
+
+    if normalise:
+        image = sqrt(image)
+
+    """
+    The second stage computes first order image gradients. These capture
+    contour, silhouette and some texture information, while providing
+    further resistance to illumination variations. The locally dominant
+    colour channel is used, which provides colour invariance to a large
+    extent. Variant methods may also include second order image derivatives,
+    which act as primitive bar detectors - a useful feature for capturing,
+    e.g. bar like structures in bicycles and limbs in humans.
+    """
+
+    gx = np.zeros(image.shape)
+    gy = np.zeros(image.shape)
+    gx[:, :-1] = np.diff(image, n=1, axis=1)
+    gy[:-1, :] = np.diff(image, n=1, axis=0)
+
+    """
+    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.
+    """
+
+    magnitude = sqrt(gx ** 2 + gy ** 2)
+    orientation = arctan2(gy, (gx + 1e-15)) * (180 / pi) + 90
+
+    # compute orientations integral images
+    integral_images = []
+    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)
+
+        #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
+
+    n_cellsx = int(np.floor(sx // cx))  # number of cells in x
+    n_cellsy = int(np.floor(sy // cy))  # number of cells in y
+
+    # now for each cell, compute the histogram
+    orientation_histogram = np.zeros((n_cellsx, n_cellsy, orientations))
+
+    radius = min(cx, cy) // 2 - 1
+    hog_image = None
+    if visualise:
+        hog_image = np.zeros((sy, sx), dtype=float)
+
+    for x in range(n_cellsx):
+        for y in range(n_cellsy):
+            for o in range(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:
+        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) / 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
+    cells and contrast normalises their overall responses before passing
+    to next stage. Normalisation introduces better invariance to illumination,
+    shadowing, and edge contrast. It is performed by accumulating a measure
+    of local histogram "energy" over local groups of cells that we call
+    "blocks". The result is used to normalise each cell in the block.
+    Typically each individual cell is shared between several blocks, but
+    its normalisations are block dependent and thus different. The cell
+    thus appears several times in the final output vector with different
+    normalisations. This may seem redundant but it improves the performance.
+    We refer to the normalised block descriptors as Histogram of Oriented
+    Gradient (HOG) descriptors.
+    """
+
+    n_blocksx = (n_cellsx - bx) + 1
+    n_blocksy = (n_cellsy - by) + 1
+    normalised_blocks = np.zeros((n_blocksx, n_blocksy,
+                                  bx, by, orientations))
+
+    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.
+    """
+
+    if visualise:
+        return normalised_blocks.ravel(), hog_image
+    else:
+        return normalised_blocks.ravel()
diff --git a/scikits/image/feature/tests/test_hog.py b/scikits/image/feature/tests/test_hog.py
new file mode 100644
index 00000000..6f375cea
--- /dev/null
+++ b/scikits/image/feature/tests/test_hog.py
@@ -0,0 +1,17 @@
+import numpy as np
+import scipy
+
+from scikits.image.feature import hog 
+
+def test_histogram_of_oriented_gradients():
+    # Replace with scikits.image.data.lena() after merge
+    img = scipy.misc.lena().astype(np.int8) 
+    
+    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__':
+    from numpy.testing import run_module_suite
+    run_module_suite()
diff --git a/scikits/image/setup.py b/scikits/image/setup.py
index 2a55b2b9..1e68ccec 100644
--- a/scikits/image/setup.py
+++ b/scikits/image/setup.py
@@ -14,6 +14,8 @@ def configuration(parent_package='', top_path=None):
     config.add_subpackage('data')
     config.add_subpackage('util')
     config.add_subpackage('color')
+    config.add_subpackage('draw')
+    config.add_subpackage('feature')
 
     def add_test_directories(arg, dirname, fnames):
         if dirname.split(os.path.sep)[-1] == 'tests':