diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index 5082c508..8258eada 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -98,3 +98,6 @@
 
 - Nicolas Poilvert
   Shape views: ``util.shape.view_as_windows`` and ``util.shape.view_as_blocks``
+
+- Johannes Schönberger
+  Polygon, circle and ellipse drawing functions
diff --git a/doc/examples/plot_shapes.py b/doc/examples/plot_shapes.py
new file mode 100644
index 00000000..8107fc80
--- /dev/null
+++ b/doc/examples/plot_shapes.py
@@ -0,0 +1,46 @@
+"""
+===========
+Fill shapes
+===========
+
+This example shows how to fill several different shapes:
+* line
+* polygon
+* circle
+* ellipse
+
+"""
+
+import matplotlib.pyplot as plt
+
+from skimage.draw import line, polygon, circle, ellipse
+import numpy as np
+
+
+img = np.zeros((500, 500, 3), 'uint8')
+
+#: draw line
+rr, cc = line(120, 123, 20, 400)
+img[rr,cc,0] = 255
+
+#: fill polygon
+poly = np.array((
+    (300, 300),
+    (480, 320),
+    (380, 430),
+    (220, 590),
+    (300, 300),
+))
+rr, cc = polygon(poly[:,0], poly[:,1], img.shape)
+img[rr,cc,1] = 255
+
+#: fill circle
+rr, cc = circle(200, 200, 100, img.shape)
+img[rr,cc,:] = (255, 255, 0)
+
+#: fill ellipse
+rr, cc = ellipse(300, 300, 100, 200, img.shape)
+img[rr,cc,2] = 255
+
+plt.imshow(img)
+plt.show()
\ No newline at end of file
diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx
index a586ae03..9dc1d307 100644
--- a/skimage/draw/_draw.pyx
+++ b/skimage/draw/_draw.pyx
@@ -1,16 +1,20 @@
 import numpy as np
+import math
+from libc.math cimport sqrt
 cimport numpy as np
 cimport cython
 
-cdef extern from "math.h":
-    int abs(int i)
+
+cdef extern from "../morphology/_pnpoly.h":
+     int pnpoly(int nr_verts, double *xp, double *yp,
+                double x, double y)
+
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def bresenham(int y, int x, int y2, int x2):
-    """
-    Generate line pixel coordinates.
-    
+def line(int y, int x, int y2, int x2):
+    """Generate line pixel coordinates.
+
     Parameters
     ----------
     y, x : int
@@ -46,7 +50,7 @@ def bresenham(int y, int x, int y2, int x2):
 
     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
@@ -64,3 +68,123 @@ def bresenham(int y, int x, int y2, int x2):
     cc[dx] = x2
 
     return rr, cc
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.nonecheck(False)
+def polygon(y, x, shape=None):
+    """Generate coordinates of pixels within polygon.
+
+    Parameters
+    ----------
+    y : (N,) ndarray
+        y coordinates of vertices of polygon
+    x : (N,) ndarray
+        x coordinates of vertices of polygon
+    shape : tuple, optional
+        image shape which is used to determine maximum extents of output pixel
+        coordinates. This is useful for polygons which exceed the image size.
+        By default the full extents of the polygon are used.
+
+    Returns
+    -------
+    rr, cc : ndarray of int
+        Pixel coordinates of polygon.
+        May be used to directly index into an array, e.g.
+        ``img[rr, cc] = 1``.
+    """
+    cdef int nr_verts = x.shape[0]
+    cdef int minr = <int>max(0, y.min())
+    cdef int maxr = <int>math.ceil(y.max())
+    cdef int minc = <int>max(0, x.min())
+    cdef int maxc = <int>math.ceil(x.max())
+
+    # make sure output coordinates do not exceed image size
+    if shape is not None:
+        maxr = min(shape[0]-1, maxr)
+        maxc = min(shape[1]-1, maxc)
+
+    cdef int r, c
+
+    #: make contigous arrays for r, c coordinates
+    cdef np.ndarray contiguous_rdata, contiguous_cdata
+    contiguous_rdata = np.ascontiguousarray(y, 'double')
+    contiguous_cdata = np.ascontiguousarray(x, 'double')
+    cdef np.double_t* rptr = <np.double_t*>contiguous_rdata.data
+    cdef np.double_t* cptr = <np.double_t*>contiguous_cdata.data
+
+    #: output coordinate arrays
+    cdef list rr = list()
+    cdef list cc = list()
+
+    for r in range(minr, maxr+1):
+        for c in range(minc, maxc+1):
+            if pnpoly(nr_verts, cptr, rptr, c, r):
+                rr.append(r)
+                cc.append(c)
+
+    return np.array(rr), np.array(cc)
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.nonecheck(False)
+@cython.cdivision(True)
+def ellipse(double cy, double cx, double b, double a, shape=None):
+    """Generate coordinates of pixels within ellipse.
+
+    Parameters
+    ----------
+    cy, cx : double
+        centre coordinate of ellipse
+    b, a: double
+        minor and major semi-axes. (x/a)**2 + (y/b)**2 = 1
+
+    Returns
+    -------
+    rr, cc : ndarray of int
+        Pixel coordinates of ellipse.
+        May be used to directly index into an array, e.g.
+        ``img[rr, cc] = 1``.
+    """
+    cdef int minr = <int>max(0, cy-b)
+    cdef int maxr = <int>math.ceil(cy+b)
+    cdef int minc = <int>max(0, cx-a)
+    cdef int maxc = <int>math.ceil(cx+a)
+
+    # make sure output coordinates do not exceed image size
+    if shape is not None:
+        maxr = min(shape[0]-1, maxr)
+        maxc = min(shape[1]-1, maxc)
+
+    cdef int r, c
+
+    #: output coordinate arrays
+    cdef list rr = list()
+    cdef list cc = list()
+
+    for r in range(minr, maxr+1):
+        for c in range(minc, maxc+1):
+            if sqrt(((r - cy)/b)**2 + ((c - cx)/a)**2) < 1:
+                rr.append(r)
+                cc.append(c)
+
+    return np.array(rr), np.array(cc)
+
+def circle(double cy, double cx, double radius, shape=None):
+    """Generate coordinates of pixels within circle.
+
+    Parameters
+    ----------
+    cy, cx : double
+        centre coordinate of circle
+    radius: double
+        radius of circle
+
+    Returns
+    -------
+    rr, cc : ndarray of int
+        Pixel coordinates of ellipse.
+        May be used to directly index into an array, e.g.
+        ``img[rr, cc] = 1``.
+    """
+    return ellipse(cy, cx, radius, radius, shape)
diff --git a/skimage/draw/draw.py b/skimage/draw/draw.py
index e66faf1c..b7bb5b88 100644
--- a/skimage/draw/draw.py
+++ b/skimage/draw/draw.py
@@ -3,4 +3,4 @@ Methods to draw on arrays.
 
 """
 
-from ._draw import bresenham
+from ._draw import line, polygon, ellipse, circle
diff --git a/skimage/draw/tests/test_draw.py b/skimage/draw/tests/test_draw.py
index ee323d83..dd1c81e2 100644
--- a/skimage/draw/tests/test_draw.py
+++ b/skimage/draw/tests/test_draw.py
@@ -1,12 +1,13 @@
 from numpy.testing import assert_array_equal
 import numpy as np
 
-from skimage.draw import bresenham
+from skimage.draw import line, polygon, circle, ellipse
 
-def test_bresenham_horizontal():
+
+def test_line_horizontal():
     img = np.zeros((10, 10))
 
-    rr, cc = bresenham(0, 0, 0, 9)
+    rr, cc = line(0, 0, 0, 9)
     img[rr, cc] = 1
 
     img_ = np.zeros((10, 10))
@@ -14,10 +15,10 @@ def test_bresenham_horizontal():
 
     assert_array_equal(img, img_)
 
-def test_bresenham_vertical():
+def test_line_vertical():
     img = np.zeros((10, 10))
 
-    rr, cc = bresenham(0, 0, 9, 0)
+    rr, cc = line(0, 0, 9, 0)
     img[rr, cc] = 1
 
     img_ = np.zeros((10, 10))
@@ -25,10 +26,10 @@ def test_bresenham_vertical():
 
     assert_array_equal(img, img_)
 
-def test_reverse():
+def test_line_reverse():
     img = np.zeros((10, 10))
 
-    rr, cc = bresenham(0, 9, 0, 0)
+    rr, cc = line(0, 9, 0, 0)
     img[rr, cc] = 1
 
     img_ = np.zeros((10, 10))
@@ -36,10 +37,10 @@ def test_reverse():
 
     assert_array_equal(img, img_)
 
-def test_diag():
+def test_line_diag():
     img = np.zeros((5, 5))
 
-    rr, cc = bresenham(0, 0, 4, 4)
+    rr, cc = line(0, 0, 4, 4)
     img[rr, cc] = 1
 
     img_ = np.eye(5)
@@ -47,6 +48,127 @@ def test_diag():
     assert_array_equal(img, img_)
 
 
+def test_polygon_rectangle():
+    img = np.zeros((10, 10), 'uint8')
+    poly = np.array(((1, 1), (4, 1), (4, 4), (1, 4), (1, 1)))
+
+    rr, cc = polygon(poly[:,0], poly[:,1])
+    img[rr,cc] = 1
+
+    img_ = np.zeros((10, 10))
+    img_[1:4,1:4] = 1
+
+    assert_array_equal(img, img_)
+
+def test_polygon_rectangle_angular():
+    img = np.zeros((10, 10), 'uint8')
+    poly = np.array(((0, 3), (4, 7), (7, 4), (3, 0), (0, 3)))
+
+    rr, cc = polygon(poly[:,0], poly[:,1])
+    img[rr,cc] = 1
+
+    img_ = np.array(
+        [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+         [0, 0, 1, 1, 0, 0, 0, 0, 0, 0],
+         [0, 1, 1, 1, 1, 0, 0, 0, 0, 0],
+         [1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
+         [0, 1, 1, 1, 1, 1, 1, 0, 0, 0],
+         [0, 0, 1, 1, 1, 1, 0, 0, 0, 0],
+         [0, 0, 0, 1, 1, 0, 0, 0, 0, 0],
+         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
+    )
+
+    assert_array_equal(img, img_)
+
+def test_polygon_parallelogram():
+    img = np.zeros((10, 10), 'uint8')
+    poly = np.array(((1, 1), (5, 1), (7, 6), (3, 6), (1, 1)))
+
+    rr, cc = polygon(poly[:,0], poly[:,1])
+    img[rr,cc] = 1
+
+    img_ = np.array(
+        [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+         [0, 1, 1, 1, 0, 0, 0, 0, 0, 0],
+         [0, 1, 1, 1, 1, 1, 0, 0, 0, 0],
+         [0, 1, 1, 1, 1, 1, 0, 0, 0, 0],
+         [0, 1, 1, 1, 1, 1, 0, 0, 0, 0],
+         [0, 0, 0, 0, 1, 1, 0, 0, 0, 0],
+         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
+    )
+
+    assert_array_equal(img, img_)
+
+def test_polygon_exceed():
+    img = np.zeros((10, 10), 'uint8')
+    poly = np.array(((1, -1), (100, -1), (100, 100), (1, 100), (1, 1)))
+
+    rr, cc = polygon(poly[:,0], poly[:,1], img.shape)
+    img[rr,cc] = 1
+
+    img_ = np.zeros((10, 10))
+    img_[1:,:] = 1
+
+    assert_array_equal(img, img_)
+
+def test_circle():
+    img = np.zeros((15, 15), 'uint8')
+
+    rr, cc = circle(7, 7, 6)
+    img[rr,cc] = 1
+
+    img_ = np.array(
+      [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
+       [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
+       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+       [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
+       [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
+       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
+    )
+
+    assert_array_equal(img, img_)
+
+def test_ellipse():
+    img = np.zeros((15, 15), 'uint8')
+
+    rr, cc = ellipse(7, 7, 3, 7)
+    img[rr,cc] = 1
+
+    img_ = np.array(
+      [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+       [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
+       [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
+       [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
+       [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
+    )
+
+    assert_array_equal(img, img_)
+
+
 if __name__ == "__main__":
     from numpy.testing import run_module_suite
-    
+    run_module_suite()