diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index 689f552b..f8f78c52 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -132,7 +132,8 @@
   Dense DAISY feature description, circle perimeter drawing.
 
 - François Boulogne
-  Drawing: Andres Method for circle perimeter, ellipse perimeter drawing, Bezier curve.
+  Drawing: Andres Method for circle perimeter, ellipse perimeter drawing,
+  Bezier curve, anti-aliasing.
   Circular and elliptical Hough Transforms
   Various fixes
 
@@ -154,3 +155,6 @@
 
 - Riaan van den Dool
   skimage.io plugin: GDAL
+
+- Fedor Morozov
+  Drawing: Wu's anti-aliased circle
diff --git a/doc/examples/plot_shapes.py b/doc/examples/plot_shapes.py
index 507452ef..8e082579 100644
--- a/doc/examples/plot_shapes.py
+++ b/doc/examples/plot_shapes.py
@@ -1,11 +1,12 @@
 """
-===========
-Fill shapes
-===========
+======
+Shapes
+======
 
-This example shows how to fill several different shapes:
+This example shows how to draw several different shapes:
 
 * line
+* Bezier curve
 * polygon
 * circle
 * ellipse
@@ -14,16 +15,17 @@ This example shows how to fill several different shapes:
 import numpy as np
 import matplotlib.pyplot as plt
 
-from skimage.draw import line, polygon, circle, circle_perimeter, \
-                         ellipse, ellipse_perimeter
-import numpy as np
+from skimage.draw import (line, polygon, circle,
+                          circle_perimeter,
+                          ellipse, ellipse_perimeter,
+                          bezier_curve)
 import math
 
 img = np.zeros((500, 500, 3), dtype=np.uint8)
 
 # draw line
 rr, cc = line(120, 123, 20, 400)
-img[rr,cc,0] = 255
+img[rr, cc, 0] = 255
 
 # fill polygon
 poly = np.array((
@@ -33,21 +35,25 @@ poly = np.array((
     (220, 590),
     (300, 300),
 ))
-rr, cc = polygon(poly[:,0], poly[:,1], img.shape)
-img[rr,cc,1] = 255
+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)
+img[rr, cc, :] = (255, 255, 0)
 
 # fill ellipse
 rr, cc = ellipse(300, 300, 100, 200, img.shape)
-img[rr,cc,2] = 255
+img[rr, cc, 2] = 255
 
 # circle
 rr, cc = circle_perimeter(120, 400, 15)
 img[rr, cc, :] = (255, 0, 0)
 
+# Bezier curve
+rr, cc = bezier_curve(70, 100, 10, 10, 150, 100, 1)
+img[rr, cc, :] = (255, 0, 0)
+
 # ellipses
 rr, cc = ellipse_perimeter(120, 400, 60, 20, orientation=math.pi / 4.)
 img[rr, cc, :] = (255, 0, 255)
@@ -58,3 +64,32 @@ img[rr, cc, :] = (255, 255, 255)
 
 plt.imshow(img)
 plt.show()
+
+"""
+
+Anti-aliased drawing for:
+* line
+* circle
+
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage.draw import (line_aa,
+                          circle_perimeter_aa)
+
+img = np.zeros((100, 100), dtype=np.uint8)
+
+# anti-aliased line
+rr, cc, val = line_aa(12, 12, 20, 50)
+img[rr, cc] = val * 255
+
+# anti-aliased circle
+rr, cc, val = circle_perimeter_aa(60, 40, 30)
+img[rr, cc] = val * 255
+
+
+plt.imshow(img, cmap=plt.cm.gray, interpolation='nearest')
+plt.title('Anti-aliasing')
+plt.show()
diff --git a/skimage/draw/__init__.py b/skimage/draw/__init__.py
index 4fb222d1..5f788ea7 100644
--- a/skimage/draw/__init__.py
+++ b/skimage/draw/__init__.py
@@ -1,9 +1,12 @@
 from .draw import circle, ellipse, set_color
-from ._draw import line, polygon, ellipse_perimeter, circle_perimeter, \
-                   bezier_segment
 from .draw3d import ellipsoid, ellipsoid_stats
+from ._draw import (line, line_aa, polygon, ellipse_perimeter,
+                    circle_perimeter, circle_perimeter_aa,
+                    _bezier_segment, bezier_curve)
 
 __all__ = ['line',
+           'line_aa',
+           'bezier_curve',
            'polygon',
            'ellipse',
            'ellipse_perimeter',
@@ -11,4 +14,5 @@ __all__ = ['line',
            'ellipsoid_stats',
            'circle',
            'circle_perimeter',
+           'circle_perimeter_aa',
            'set_color']
diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx
index 132ba1d4..3d48ce86 100644
--- a/skimage/draw/_draw.pyx
+++ b/skimage/draw/_draw.pyx
@@ -6,7 +6,7 @@ import math
 import numpy as np
 
 cimport numpy as cnp
-from libc.math cimport sqrt, sin, cos, floor
+from libc.math cimport sqrt, sin, cos, floor, ceil
 from skimage._shared.geometry cimport point_in_polygon
 
 
@@ -27,6 +27,10 @@ def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2):
         May be used to directly index into an array, e.g.
         ``img[rr, cc] = 1``.
 
+    Notes
+    -----
+    Anti-aliased line generator is available with `line_aa`.
+
     Examples
     --------
     >>> from skimage.draw import line
@@ -89,6 +93,99 @@ def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2):
     return np.asarray(rr), np.asarray(cc)
 
 
+def line_aa(Py_ssize_t y1, Py_ssize_t x1, Py_ssize_t y2, Py_ssize_t x2):
+    """Generate anti-aliased line pixel coordinates.
+
+    Parameters
+    ----------
+    y1, x1 : int
+        Starting position (row, column).
+    y2, x2 : int
+        End position (row, column).
+
+    Returns
+    -------
+    rr, cc, val : (N,) ndarray (int, int, float)
+        Indices of pixels (`rr`, `cc`) and intensity values (`val`).
+        ``img[rr, cc] = val``.
+
+    References
+    ----------
+    .. [1] A Rasterizing Algorithm for Drawing Curves, A. Zingl, 2012
+           http://members.chello.at/easyfilter/Bresenham.pdf
+
+    Examples
+    --------
+    >>> from skimage.draw import line_aa
+    >>> img = np.zeros((10, 10), dtype=np.uint8)
+    >>> rr, cc, val = line_aa(1, 1, 8, 8)
+    >>> img[rr, cc] = val * 255
+    >>> img
+    array([[  0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
+           [  0, 255,  56,   0,   0,   0,   0,   0,   0,   0],
+           [  0,  56, 255,  56,   0,   0,   0,   0,   0,   0],
+           [  0,   0,  56, 255,  56,   0,   0,   0,   0,   0],
+           [  0,   0,   0,  56, 255,  56,   0,   0,   0,   0],
+           [  0,   0,   0,   0,  56, 255,  56,   0,   0,   0],
+           [  0,   0,   0,   0,   0,  56, 255,  56,   0,   0],
+           [  0,   0,   0,   0,   0,   0,  56, 255,  56,   0],
+           [  0,   0,   0,   0,   0,   0,   0,  56, 255,   0],
+           [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0]], dtype=uint8)
+    """
+    cdef list rr = list()
+    cdef list cc = list()
+    cdef list val = list()
+
+    cdef int dx = abs(x1 - x2)
+    cdef int dy = abs(y1 - y2)
+    cdef int err = dx - dy
+    cdef int x, y, e, ed, sign_x, sign_y
+
+    if x1 < x2:
+        sign_x = 1
+    else:
+        sign_x = -1
+
+    if y1 < y2:
+        sign_y = 1
+    else:
+        sign_y = -1
+
+    if dx + dy == 0:
+        ed = 1
+    else:
+        ed = <int>(sqrt(dx*dx + dy*dy))
+
+    x, y = x1, y1
+    while True:
+        cc.append(x)
+        rr.append(y)
+        val.append(1. * abs(err - dx + dy) / <float>(ed))
+        e = err
+        if 2 * e >= -dx:
+            if x == x2:
+                break
+            if e + dy < ed:
+                cc.append(x)
+                rr.append(y + sign_y)
+                val.append(1. * abs(e + dy) / <float>(ed))
+            err -= dy
+            x += sign_x
+        if 2 * e <= dy:
+            if y == y2:
+                break
+            if dx - e < ed:
+                cc.append(x)
+                rr.append(y)
+                val.append(abs(dx - e) / <float>(ed))
+            err += dx
+            y += sign_y
+
+    return (np.array(rr, dtype=np.intp),
+            np.array(cc, dtype=np.intp),
+            1. - np.array(val, dtype=np.float))
+
+
 def polygon(y, x, shape=None):
     """Generate coordinates of pixels within polygon.
 
@@ -134,9 +231,9 @@ def polygon(y, x, shape=None):
 
     cdef Py_ssize_t nr_verts = x.shape[0]
     cdef Py_ssize_t minr = int(max(0, y.min()))
-    cdef Py_ssize_t maxr = int(math.ceil(y.max()))
+    cdef Py_ssize_t maxr = int(ceil(y.max()))
     cdef Py_ssize_t minc = int(max(0, x.min()))
-    cdef Py_ssize_t maxc = int(math.ceil(x.max()))
+    cdef Py_ssize_t maxc = int(ceil(x.max()))
 
     # make sure output coordinates do not exceed image size
     if shape is not None:
@@ -182,6 +279,7 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius,
     Returns
     -------
     rr, cc : (N,) ndarray of int
+        Bresenham and Andres' method:
         Indices of pixels that belong to the circle perimeter.
         May be used to directly index into an array, e.g.
         ``img[rr, cc] = 1``.
@@ -192,13 +290,14 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius,
     circles create a disc whereas Bresenham can make holes. There
     is also less distortions when Andres circles are rotated.
     Bresenham method is also known as midpoint circle algorithm.
+    Anti-aliased circle generator is available with `circle_perimeter_aa`.
 
     References
     ----------
     .. [1] J.E. Bresenham, "Algorithm for computer control of a digital
-           plotter", 4 (1965) 25-30.
-    .. [2] E. Andres, "Discrete circles, rings and spheres",
-           18 (1994) 695-706.
+           plotter", IBM Systems journal, 4 (1965) 25-30.
+    .. [2] E. Andres, "Discrete circles, rings and spheres", Computers &
+           Graphics, 18 (1994) 695-706.
 
     Examples
     --------
@@ -226,6 +325,10 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius,
     cdef Py_ssize_t x = 0
     cdef Py_ssize_t y = radius
     cdef Py_ssize_t d = 0
+
+    cdef double dceil = 0
+    cdef double dceil_prev = 0
+
     cdef char cmethod
     if method == 'bresenham':
         d = 3 - 2 * radius
@@ -258,8 +361,84 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius,
                 d = d + 2 * (y - x - 1)
                 y = y - 1
                 x = x + 1
+    return (np.array(rr, dtype=np.intp) + cy,
+            np.array(cc, dtype=np.intp) + cx)
 
-    return np.array(rr, dtype=np.intp) + cy, np.array(cc, dtype=np.intp) + cx
+
+def circle_perimeter_aa(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius):
+    """Generate anti-aliased circle perimeter coordinates.
+
+    Parameters
+    ----------
+    cy, cx : int
+        Centre coordinate of circle.
+    radius: int
+        Radius of circle.
+
+    Returns
+    -------
+    rr, cc, val : (N,) ndarray (int, int, float)
+        Indices of pixels (`rr`, `cc`) and intensity values (`val`).
+        ``img[rr, cc] = val``.
+
+    Notes
+    -----
+    Wu's method draws anti-aliased circle. This implementation doesn't use
+    lookup table optimization.
+
+    References
+    ----------
+    .. [1] X. Wu, "An efficient antialiasing technique", In ACM SIGGRAPH
+           Computer Graphics, 25 (1991) 143-152.
+
+    Examples
+    --------
+    >>> from skimage.draw import circle_perimeter_aa
+    >>> img = np.zeros((10, 10), dtype=np.uint8)
+    >>> rr, cc, val = circle_perimeter_aa(4, 4, 3)
+    >>> img[rr, cc] = val * 255
+    >>> img
+    array([[  0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
+           [  0,   0,  60, 211, 255, 211,  60,   0,   0,   0],
+           [  0,  60, 194,  43,   0,  43, 194,  60,   0,   0],
+           [  0, 211,  43,   0,   0,   0,  43, 211,   0,   0],
+           [  0, 255,   0,   0,   0,   0,   0, 255,   0,   0],
+           [  0, 211,  43,   0,   0,   0,  43, 211,   0,   0],
+           [  0,  60, 194,  43,   0,  43, 194,  60,   0,   0],
+           [  0,   0,  60, 211, 255, 211,  60,   0,   0,   0],
+           [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
+           [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0]], dtype=uint8)
+    """
+
+    cdef Py_ssize_t x = 0
+    cdef Py_ssize_t y = radius
+    cdef Py_ssize_t d = 0
+
+    cdef double dceil = 0
+    cdef double dceil_prev = 0
+
+    cdef list rr = [y, x,  y,  x, -y, -x, -y, -x]
+    cdef list cc = [x, y, -x, -y,  x,  y, -x, -y]
+    cdef list val = [1] * 8
+
+    while y > x + 1:
+        x += 1
+        dceil = sqrt(radius**2 - x**2)
+        dceil = ceil(dceil) - dceil
+        if dceil < dceil_prev:
+            y -= 1
+        rr.extend([y, y - 1, x, x, y, y - 1, x, x])
+        cc.extend([x, x, y, y - 1, -x, -x, -y, 1 - y])
+
+        rr.extend([-y, 1 - y, -x, -x, -y, 1 - y, -x, -x])
+        cc.extend([x, x, y, y - 1, -x, -x, -y, 1 - y])
+
+        val.extend([1 - dceil, dceil] * 8)
+        dceil_prev = dceil
+
+    return (np.array(rr, dtype=np.intp) + cy,
+            np.array(cc, dtype=np.intp) + cx,
+            np.array(val, dtype=np.float))
 
 
 def ellipse_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t yradius,
@@ -382,38 +561,38 @@ def ellipse_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t yradius,
         iyd = int(floor(ya * w + 0.5))
 
         # Draw the 4 quadrants
-        rr, cc = bezier_segment(iy0 + iyd, ix0, iy0, ix0, iy0, ix0 + ixd, 1-w)
+        rr, cc = _bezier_segment(iy0 + iyd, ix0, iy0, ix0, iy0, ix0 + ixd, 1-w)
         py.extend(rr)
         px.extend(cc)
-        rr, cc = bezier_segment(iy0 + iyd, ix0, iy1, ix0, iy1, ix1 - ixd, w)
+        rr, cc = _bezier_segment(iy0 + iyd, ix0, iy1, ix0, iy1, ix1 - ixd, w)
         py.extend(rr)
         px.extend(cc)
-        rr, cc = bezier_segment(iy1 - iyd, ix1, iy1, ix1, iy1, ix1 - ixd, 1-w)
+        rr, cc = _bezier_segment(iy1 - iyd, ix1, iy1, ix1, iy1, ix1 - ixd, 1-w)
         py.extend(rr)
         px.extend(cc)
-        rr, cc = bezier_segment(iy1 - iyd, ix1, iy0, ix1, iy0, ix0 + ixd,  w)
+        rr, cc = _bezier_segment(iy1 - iyd, ix1, iy0, ix1, iy0, ix0 + ixd,  w)
         py.extend(rr)
         px.extend(cc)
 
     return np.array(py, dtype=np.intp), np.array(px, dtype=np.intp)
 
 
-def bezier_segment(Py_ssize_t y0, Py_ssize_t x0,
-                   Py_ssize_t y1, Py_ssize_t x1,
-                   Py_ssize_t y2, Py_ssize_t x2,
-                   double weight):
+def _bezier_segment(Py_ssize_t y0, Py_ssize_t x0,
+                    Py_ssize_t y1, Py_ssize_t x1,
+                    Py_ssize_t y2, Py_ssize_t x2,
+                    double weight):
     """Generate Bezier segment coordinates.
 
     Parameters
     ----------
     y0, x0 : int
-        Coordinates of the first point
+        Coordinates of the first control point.
     y1, x1 : int
-        Coordinates of the middle point
+        Coordinates of the middle control point.
     y2, x2 : int
-        Coordinates of the last point
+        Coordinates of the last control point.
     weight : double
-        Middle point weight, it describes the line tension.
+        Middle control point weight, it describes the line tension.
 
     Returns
     -------
@@ -425,7 +604,7 @@ def bezier_segment(Py_ssize_t y0, Py_ssize_t x0,
     Notes
     -----
     The algorithm is the rational quadratic algorithm presented in
-    reference [1].
+    reference [1]_.
 
     References
     ----------
@@ -492,8 +671,8 @@ def bezier_segment(Py_ssize_t y0, Py_ssize_t x0,
             sy = floor((y0 + 2 * weight * y1 + y2) * xy * 0.5 + 0.5)
             dx = floor((weight * x1 + x0) * xy + 0.5)
             dy = floor((y1 * weight + y0) * xy + 0.5)
-            return bezier_segment(y0, x0, <Py_ssize_t>(dy), <Py_ssize_t>(dx),
-                                  <Py_ssize_t>(sy), <Py_ssize_t>(sx), cur)
+            return _bezier_segment(y0, x0, <Py_ssize_t>(dy), <Py_ssize_t>(dx),
+                                   <Py_ssize_t>(sy), <Py_ssize_t>(sx), cur)
 
         err = dx + dy - xy
         while dy <= xy and dx >= xy:
@@ -526,6 +705,137 @@ def bezier_segment(Py_ssize_t y0, Py_ssize_t x0,
     return np.array(py, dtype=np.intp), np.array(px, dtype=np.intp)
 
 
+def bezier_curve(Py_ssize_t y0, Py_ssize_t x0,
+                 Py_ssize_t y1, Py_ssize_t x1,
+                 Py_ssize_t y2, Py_ssize_t x2,
+                 double weight):
+    """Generate Bezier curve coordinates.
+
+    Parameters
+    ----------
+    y0, x0 : int
+        Coordinates of the first control point.
+    y1, x1 : int
+        Coordinates of the middle control point.
+    y2, x2 : int
+        Coordinates of the last control point.
+    weight : double
+        Middle control point weight, it describes the line tension.
+
+    Returns
+    -------
+    rr, cc : (N,) ndarray of int
+        Indices of pixels that belong to the Bezier curve.
+        May be used to directly index into an array, e.g.
+        ``img[rr, cc] = 1``.
+
+    Notes
+    -----
+    The algorithm is the rational quadratic algorithm presented in
+    reference [1]_.
+
+    References
+    ----------
+    .. [1] A Rasterizing Algorithm for Drawing Curves, A. Zingl, 2012
+           http://members.chello.at/easyfilter/Bresenham.pdf
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from skimage.draw import bezier_curve
+    >>> img = np.zeros((10, 10), dtype=np.uint8)
+    >>> rr, cc = bezier_curve(1, 5, 5, -2, 8, 8, 2)
+    >>> img[rr, cc] = 1
+    >>> img
+    array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+           [0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
+           [0, 0, 0, 1, 1, 0, 0, 0, 0, 0],
+           [0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
+           [0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
+           [0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
+           [0, 0, 1, 1, 0, 0, 0, 0, 0, 0],
+           [0, 0, 0, 0, 1, 1, 1, 0, 0, 0],
+           [0, 0, 0, 0, 0, 0, 0, 1, 1, 0],
+           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
+    """
+    # Pixels
+    cdef list px = list()
+    cdef list py = list()
+
+    cdef int x, y
+    cdef double xx, yy, ww, t, q
+    x = x0 - 2 * x1 + x2
+    y = y0 - 2 * y1 + y2
+
+    xx = x0 - x1
+    yy = y0 - y1
+
+    if xx * (x2 - x1) > 0:
+        if yy * (y2 - y1):
+            if abs(xx * y) > abs(yy * x):
+                x0 = x2
+                x2 = <Py_ssize_t>(xx + x1)
+                y0 = y2
+                y2 = <Py_ssize_t>(yy + y1)
+        if (x0 == x2) or (weight == 1.):
+            t = <double>(x0 - x1) / x
+        else:
+            q = sqrt(4. * weight * weight * (x0 - x1) * (x2 - x1) + (x2 - x0) * floor(x2 - x0))
+            if (x1 < x0):
+                q = -q
+            t = (2. * weight * (x0 - x1) - x0 + x2 + q) / (2. * (1. - weight) * (x2 - x0))
+
+        q = 1. / (2. * t * (1. - t) * (weight - 1.) + 1.0)
+        xx = (t * t * (x0 - 2. * weight * x1 + x2) + 2. * t * (weight * x1 - x0) + x0) * q
+        yy = (t * t * (y0 - 2. * weight * y1 + y2) + 2. * t * (weight * y1 - y0) + y0) * q
+        ww = t * (weight - 1.) + 1.
+        ww *= ww * q
+        weight = ((1. - t) * (weight - 1.) + 1.) * sqrt(q)
+        x = <int>(xx + 0.5)
+        y = <int>(yy + 0.5)
+        yy = (xx - x0) * (y1 - y0) / (x1 - x0) + y0
+
+        rr, cc = _bezier_segment(y0, x0, <int>(yy + 0.5), x, y, x, ww)
+        px.extend(rr)
+        py.extend(cc)
+
+        yy = (xx - x2) * (y1 - y2) / (x1 - x2) + y2
+        y1 = <int>(yy + 0.5)
+        x0 = x1 = x
+        y0 = y
+    if (y0 - y1) * floor(y2 - y1) > 0:
+        if (y0 == y2) or (weight == 1):
+            t = (y0 - y1) / (y0 - 2. * y1 + y2)
+        else:
+            q = sqrt(4. * weight * weight * (y0 - y1) * (y2 - y1) + (y2 - y0) * floor(y2 - y0))
+            if y1 < y0:
+                q = -q
+            t = (2. * weight * (y0 - y1) - y0 + y2 + q) / (2. * (1. - weight) * (y2 - y0))
+        q = 1. / (2. * t * (1. - t) * (weight - 1.) + 1.)
+        xx = (t * t * (x0 - 2. * weight * x1 + x2) + 2. * t * (weight * x1 - x0) + x0) * q
+        yy = (t * t * (y0 - 2. * weight * y1 + y2) + 2. * t * (weight * y1 - y0) + y0) * q
+        ww = t * (weight - 1.) + 1.
+        ww *= ww * q
+        weight = ((1. - t) * (weight - 1.) + 1.) * sqrt(q)
+        x = <int>(xx + 0.5)
+        y = <int>(yy + 0.5)
+        xx = (x1 - x0) * (yy - y0) / (y1 - y0) + x0
+
+        rr, cc = _bezier_segment(y0, x0, y, <int>(xx + 0.5), y, x, ww)
+        px.extend(rr)
+        py.extend(cc)
+
+        xx = (x1 - x2) * (yy - y2) / (y1 - y2) + x2
+        x1 = <int>(xx + 0.5)
+        x0 = x
+        y0 = y1 = y
+
+    rr, cc = _bezier_segment(y0, x0, y1, x1, y2, x2, weight * weight)
+    px.extend(rr)
+    py.extend(cc)
+    return np.array(px, dtype=np.intp), np.array(py, dtype=np.intp)
+
+
 def set_color(img, coords, color):
     """Set pixel color in the image at the given coordinates.
 
diff --git a/skimage/draw/draw.py b/skimage/draw/draw.py
index d01bc2b0..cbf3ced2 100644
--- a/skimage/draw/draw.py
+++ b/skimage/draw/draw.py
@@ -52,7 +52,7 @@ def ellipse(cy, cx, yradius, xradius, shape=None):
     dc = 1 / float(xradius)
 
     r, c = np.ogrid[-1:1:dr, -1:1:dc]
-    rr, cc = np.nonzero(r ** 2  + c ** 2 < 1)
+    rr, cc = np.nonzero(r ** 2 + c ** 2 < 1)
 
     rr.flags.writeable = True
     cc.flags.writeable = True
diff --git a/skimage/draw/tests/test_draw.py b/skimage/draw/tests/test_draw.py
index f15fba01..d3cab811 100644
--- a/skimage/draw/tests/test_draw.py
+++ b/skimage/draw/tests/test_draw.py
@@ -1,8 +1,11 @@
-from numpy.testing import assert_array_equal
+from numpy.testing import assert_array_equal, assert_equal
 import numpy as np
 
-from skimage.draw import line, polygon, circle, circle_perimeter, \
-                         ellipse, ellipse_perimeter, bezier_segment
+from skimage.draw import (line, line_aa, polygon,
+                          circle, circle_perimeter, circle_perimeter_aa,
+                          ellipse, ellipse_perimeter,
+                          _bezier_segment, bezier_curve,
+                          )
 
 
 def test_line_horizontal():
@@ -52,6 +55,43 @@ def test_line_diag():
     assert_array_equal(img, img_)
 
 
+def test_line_aa_horizontal():
+    img = np.zeros((10, 10))
+
+    rr, cc, val = line_aa(0, 0, 0, 9)
+    img[rr, cc] = val
+
+    img_ = np.zeros((10, 10))
+    img_[0, :] = 1
+
+    assert_array_equal(img, img_)
+
+
+def test_line_aa_vertical():
+    img = np.zeros((10, 10))
+
+    rr, cc, val = line_aa(0, 0, 9, 0)
+    img[rr, cc] = val
+
+    img_ = np.zeros((10, 10))
+    img_[:, 0] = 1
+
+    assert_array_equal(img, img_)
+
+
+def test_line_aa_diagonal():
+    img = np.zeros((10, 10))
+
+    rr, cc, val = line_aa(0, 0, 9, 6)
+    img[rr, cc] = 1
+
+    # Check that each pixel belonging to line,
+    # also belongs to line_aa
+    r, c = line(0, 0, 9, 6)
+    for x, y in zip(r, c):
+        assert_equal(img[r, c], 1)
+
+
 def test_polygon_rectangle():
     img = np.zeros((10, 10), 'uint8')
     poly = np.array(((1, 1), (4, 1), (4, 4), (1, 4), (1, 1)))
@@ -215,6 +255,38 @@ def test_circle_perimeter_andres():
     assert_array_equal(img, img_)
 
 
+def test_circle_perimeter_aa():
+    img = np.zeros((15, 15), 'uint8')
+    rr, cc, val = circle_perimeter_aa(7, 7, 0)
+    img[rr, cc] = 1
+    assert(np.sum(img) == 1)
+    assert(img[7][7] == 1)
+
+    img = np.zeros((17, 17), 'uint8')
+    rr, cc, val = circle_perimeter_aa(8, 8, 7)
+    img[rr, cc] = val * 255
+    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,  82, 180, 236, 255, 236, 180,  82,   0,   0,   0,   0,   0],
+         [  0,   0,   0,   0, 189, 172,  74,  18,   0,  18,  74, 172, 189,   0,   0,   0,   0],
+         [  0,   0,   0, 229,  25,   0,   0,   0,   0,   0,   0,   0,  25, 229,   0,   0,   0],
+         [  0,   0, 189,  25,   0,   0,   0,   0,   0,   0,   0,   0,   0,  25, 189,   0,   0],
+         [  0,  82, 172,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0, 172,  82,   0],
+         [  0, 180,  74,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,  74, 180,   0],
+         [  0, 236,  18,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,  18, 236,   0],
+         [  0, 255,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0, 255,   0],
+         [  0, 236,  18,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,  18, 236,   0],
+         [  0, 180,  74,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,  74, 180,   0],
+         [  0,  82, 172,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0, 172,  82,   0],
+         [  0,   0, 189,  25,   0,   0,   0,   0,   0,   0,   0,   0,   0,  25, 189,   0,   0],
+         [  0,   0,   0, 229,  25,   0,   0,   0,   0,   0,   0,   0,  25, 229,   0,   0,   0],
+         [  0,   0,   0,   0, 189, 172,  74,  18,   0,  18,  74, 172, 189,   0,   0,   0,   0],
+         [  0,   0,   0,   0,   0,  82, 180, 236, 255, 236, 180,  82,   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')
 
@@ -360,18 +432,21 @@ def test_bezier_segment_straight():
     y1 = 50
     x2 = 150
     y2 = 150
-    rr, cc = bezier_segment(x0, y0, x1, y1, x2, y2, 0)
-    image [rr, cc] = 1
+    rr, cc = _bezier_segment(x0, y0, x1, y1, x2, y2, 0)
+    image[rr, cc] = 1
 
     image2 = np.zeros((200, 200), dtype=int)
     rr, cc = line(x0, y0, x2, y2)
-    image2 [rr, cc] = 1
+    image2[rr, cc] = 1
     assert_array_equal(image, image2)
 
 
 def test_bezier_segment_curved():
     img = np.zeros((25, 25), 'uint8')
-    rr, cc = bezier_segment(20, 20, 20, 2, 2, 2, 1)
+    x1, y1 = 20, 20
+    x2, y2 = 20, 2
+    x3, y3 = 2, 2
+    rr, cc = _bezier_segment(x1, y1, x2, y2, x3, y3, 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, 0, 0, 0, 0],
@@ -400,9 +475,101 @@ def test_bezier_segment_curved():
             [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_equal(img[x1, y1], 1)
+    assert_equal(img[x3, y3], 1)
     assert_array_equal(img, img_)
 
 
+def test_bezier_curve_straight():
+    image = np.zeros((200, 200), dtype=int)
+    x0 = 50
+    y0 = 50
+    x1 = 150
+    y1 = 50
+    x2 = 150
+    y2 = 150
+    rr, cc = bezier_curve(x0, y0, x1, y1, x2, y2, 0)
+    image [rr, cc] = 1
+
+    image2 = np.zeros((200, 200), dtype=int)
+    rr, cc = line(x0, y0, x2, y2)
+    image2 [rr, cc] = 1
+    assert_array_equal(image, image2)
+
+
+def test_bezier_curved_weight_eq_1():
+    img = np.zeros((23, 8), 'uint8')
+    x1, y1 = (1, 1)
+    x2, y2 = (11, 11)
+    x3, y3 = (21, 1)
+    rr, cc = bezier_curve(x1, y1, x2, y2, x3, y3, 1)
+    img[rr, cc] = 1
+    assert_equal(img[x1, y1], 1)
+    assert_equal(img[x3, y3], 1)
+    img_ = np.array(
+           [[0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 1, 0, 0, 0, 0, 0, 0],
+            [0, 0, 1, 0, 0, 0, 0, 0],
+            [0, 0, 0, 1, 0, 0, 0, 0],
+            [0, 0, 0, 0, 1, 0, 0, 0],
+            [0, 0, 0, 0, 1, 0, 0, 0],
+            [0, 0, 0, 0, 0, 1, 0, 0],
+            [0, 0, 0, 0, 0, 1, 0, 0],
+            [0, 0, 0, 0, 0, 0, 1, 0],
+            [0, 0, 0, 0, 0, 0, 1, 0],
+            [0, 0, 0, 0, 0, 0, 1, 0],
+            [0, 0, 0, 0, 0, 0, 1, 0],
+            [0, 0, 0, 0, 0, 0, 1, 0],
+            [0, 0, 0, 0, 0, 0, 1, 0],
+            [0, 0, 0, 0, 0, 0, 1, 0],
+            [0, 0, 0, 0, 0, 1, 0, 0],
+            [0, 0, 0, 0, 0, 1, 0, 0],
+            [0, 0, 0, 0, 1, 0, 0, 0],
+            [0, 0, 0, 0, 1, 0, 0, 0],
+            [0, 0, 0, 1, 0, 0, 0, 0],
+            [0, 0, 1, 0, 0, 0, 0, 0],
+            [0, 1, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0]]
+            )
+    assert_equal(img, img_)
+
+
+def test_bezier_curved_weight_neq_1():
+    img = np.zeros((23, 10), 'uint8')
+    x1, y1 = (1, 1)
+    x2, y2 = (11, 11)
+    x3, y3 = (21, 1)
+    rr, cc = bezier_curve(x1, y1, x2, y2, x3, y3, 2)
+    img[rr, cc] = 1
+    assert_equal(img[x1, y1], 1)
+    assert_equal(img[x3, y3], 1)
+    img_ = np.array(
+            [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+             [0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
+             [0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
+             [0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
+             [0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
+             [0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
+             [0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
+             [0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
+             [0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
+             [0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+             [0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+             [0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+             [0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+             [0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+             [0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
+             [0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
+             [0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
+             [0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
+             [0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
+             [0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
+             [0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
+             [0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
+             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
+            )
+    assert_equal(img, img_)
+
 if __name__ == "__main__":
     from numpy.testing import run_module_suite
     run_module_suite()