diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index fd3ff235..9f3f2aa7 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -126,3 +126,6 @@
 
 - Steven Silvester, Karel Zuiderveld
   Adaptive Histogram Equalization
+
+- Anders Boesen Lindbo Larsen
+  Dense DAISY feature description, circle perimeter drawing.
diff --git a/doc/examples/plot_daisy.py b/doc/examples/plot_daisy.py
new file mode 100644
index 00000000..af78103d
--- /dev/null
+++ b/doc/examples/plot_daisy.py
@@ -0,0 +1,28 @@
+"""
+===============================
+Dense DAISY feature description
+===============================
+
+The DAISY local image descriptor is based on gradient orientation histograms
+similar to the SIFT descriptor. It is formulated in a way that allows for fast
+dense extraction which is useful for e.g. bag-of-features image
+representations.
+
+In this example a limited number of DAISY descriptors are extracted at a large
+scale for illustrative purposes.
+"""
+
+from skimage.feature import daisy
+from skimage import data
+import matplotlib.pyplot as plt
+
+
+img = data.camera()
+descs, descs_img = daisy(img, step=180, radius=58, rings=2, histograms=6,
+                         orientations=8, visualize=True)
+
+plt.axis('off')
+plt.imshow(descs_img)
+descs_num = descs.shape[0] * descs.shape[1]
+plt.title('%i DAISY descriptors extracted:' % descs_num)
+plt.show()
diff --git a/doc/examples/plot_shapes.py b/doc/examples/plot_shapes.py
index 9e586209..25a48ebe 100644
--- a/doc/examples/plot_shapes.py
+++ b/doc/examples/plot_shapes.py
@@ -13,7 +13,7 @@ This example shows how to fill several different shapes:
 
 import matplotlib.pyplot as plt
 
-from skimage.draw import line, polygon, circle, ellipse
+from skimage.draw import line, polygon, circle, circle_perimeter, ellipse
 import numpy as np
 
 
@@ -42,5 +42,9 @@ img[rr,cc,:] = (255, 255, 0)
 rr, cc = ellipse(300, 300, 100, 200, img.shape)
 img[rr,cc,2] = 255
 
+# circle
+rr, cc = circle_perimeter(120, 400, 50)
+img[rr, cc, :] = (255, 0, 255)
+
 plt.imshow(img)
-plt.show()
\ No newline at end of file
+plt.show()
diff --git a/skimage/draw/__init__.py b/skimage/draw/__init__.py
index 4eac9b63..4b30ea18 100644
--- a/skimage/draw/__init__.py
+++ b/skimage/draw/__init__.py
@@ -1,2 +1,2 @@
-from ._draw import line, polygon, ellipse, circle
+from ._draw import line, polygon, ellipse, circle, circle_perimeter, set_color
 bresenham = line
diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx
index ca0c502a..509796bb 100644
--- a/skimage/draw/_draw.pyx
+++ b/skimage/draw/_draw.pyx
@@ -187,3 +187,69 @@ def circle(double cy, double cx, double radius, shape=None):
         ``img[rr, cc] = 1``.
     """
     return ellipse(cy, cx, radius, radius, shape)
+
+
+def circle_perimeter(int cy, int cx, int radius):
+    """Generate circle perimeter coordinates.
+
+    Parameters
+    ----------
+    cy, cx : int
+        Centre coordinate of circle.
+    radius: int
+        Radius of circle.
+
+    Returns
+    -------
+    rr, cc : (N,) ndarray of int
+        Indices of pixels that belong to the circle perimeter.
+        May be used to directly index into an array, e.g.
+        ``img[rr, cc] = 1``.
+
+    """
+
+    cdef list rr = list()
+    cdef list cc = list()
+
+    cdef int x = 0
+    cdef int y = radius
+    cdef int d = 3 - 2 * radius
+
+    while y >= x:
+        rr.extend([y, -y, y, -y, x, -x, x, -x])
+        cc.extend([x, x, -x, -x, y, y, -y, -y])
+        if d < 0:
+            d += 4 * x + 6
+        else:
+            d += 4 * (x - y) + 10
+            y -= 1
+        x += 1
+
+    return np.array(rr) + cy, np.array(cc) + cx
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def set_color(img, coords, color):
+    """Set pixel color in the image at the given coordiantes. Coordinates that
+    exceeed the shape of the image will be ignored.
+
+    Parameters
+    ----------
+    img : (M, N, D) ndarray
+        Image
+    coords : ((P,) ndarray, (P,) ndarray)
+        Coordinates of pixels to be colored.
+    color : (D,) ndarray
+        Color to be assigned to coordinates in the image.
+
+    Returns
+    -------
+    img : (M, N, D) ndarray
+        The updated image.
+    """
+    rr, cc = coords
+    rr_inside = np.logical_and(rr >= 0, rr < img.shape[0])
+    cc_inside = np.logical_and(cc >= 0, cc < img.shape[1])
+    inside = np.logical_and(rr_inside, cc_inside)
+    img[rr[inside], cc[inside]] = color
diff --git a/skimage/draw/tests/test_draw.py b/skimage/draw/tests/test_draw.py
index 9e6ca8e8..04cbae92 100644
--- a/skimage/draw/tests/test_draw.py
+++ b/skimage/draw/tests/test_draw.py
@@ -1,7 +1,7 @@
 from numpy.testing import assert_array_equal
 import numpy as np
 
-from skimage.draw import line, polygon, circle, ellipse
+from skimage.draw import line, polygon, circle, circle_perimeter, ellipse
 
 
 def test_line_horizontal():
@@ -150,6 +150,37 @@ def test_circle():
     assert_array_equal(img, img_)
 
 
+def test_circle_perimeter():
+    img = np.zeros((15, 15), 'uint8')
+    rr, cc = circle_perimeter(7, 7, 0)
+    img[rr, cc] = 1
+    assert(np.sum(img) == 1)
+
+    img = np.zeros((17, 15), 'uint8')
+    rr, cc = circle_perimeter(7, 7, 7)
+    img[rr, cc] = 1
+    img_ = np.array(
+        [[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, 1, 1, 0, 0, 0],
+         [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
+         [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+         [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+         [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+         [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+         [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+         [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+         [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+         [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+         [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+         [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
+         [0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0],
+         [0, 0, 0, 0, 0, 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]]
+    )
+    assert_array_equal(img, img_)
+
+
 def test_ellipse():
     img = np.zeros((15, 15), 'uint8')
 
diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py
index e0756bfe..9f1e5f92 100644
--- a/skimage/feature/__init__.py
+++ b/skimage/feature/__init__.py
@@ -1,3 +1,4 @@
+from ._daisy import daisy
 from ._hog import hog
 from .texture import greycomatrix, greycoprops, local_binary_pattern
 from .peak import peak_local_max
diff --git a/skimage/feature/_daisy.py b/skimage/feature/_daisy.py
new file mode 100644
index 00000000..f83605b5
--- /dev/null
+++ b/skimage/feature/_daisy.py
@@ -0,0 +1,217 @@
+import numpy as np
+from scipy import sqrt, pi, arctan2, cos, sin, exp
+from scipy.ndimage import gaussian_filter
+import skimage.color
+from skimage import img_as_float, draw
+
+
+def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
+          normalization='l1', sigmas=None, ring_radii=None, visualize=False):
+    '''Extract DAISY feature descriptors densely for the given image.
+
+    DAISY is a feature descriptor similar to SIFT formulated in a way that
+    allows for fast dense extraction. Typically, this is practical for
+    bag-of-features image representations.
+
+    The implementation follows Tola et al. [1]_ but deviate on the following
+    points:
+
+      * Histogram bin contribution are smoothed with a circular Gaussian
+        window over the tonal range (the angular range).
+      * The sigma values of the spatial Gaussian smoothing in this code do not
+        match the sigma values in the original code by Tola et al. [2]_. In
+        their code, spatial smoothing is applied to both the input image and
+        the center histogram. However, this smoothing is not documented in [1]_
+        and, therefore, it is omitted.
+
+    Parameters
+    ----------
+    img : (M, N) array
+        Input image (greyscale).
+    step : int, optional
+        Distance between descriptor sampling points.
+    radius : int, optional
+        Radius (in pixels) of the outermost ring.
+    rings : int, optional
+        Number of rings.
+    histograms  : int, optional
+        Number of histograms sampled per ring.
+    orientations : int, optional
+        Number of orientations (bins) per histogram.
+    normalization : [ 'l1' | 'l2' | 'daisy' | 'off' ], optional
+        How to normalize the descriptors
+
+          * 'l1': L1-normalization of each descriptor.
+          * 'l2': L2-normalization of each descriptor.
+          * 'daisy': L2-normalization of individual histograms.
+          * 'off': Disable normalization.
+
+    sigmas : 1D array of float, optional
+        Standard deviation of spatial Gaussian smoothing for the center
+        histogram and for each ring of histograms. The array of sigmas should
+        be sorted from the center and out. I.e. the first sigma value defines
+        the spatial smoothing of the center histogram and the last sigma value
+        defines the spatial smoothing of the outermost ring. Specifying sigmas
+        overrides the following parameter.
+           ``rings = len(sigmas)-1``
+
+    ring_radii : 1D array of int, optional
+        Radius (in pixels) for each ring. Specifying ring_radii overrides the
+        following two parameters.
+          | ``rings = len(ring_radii)``
+          | ``radius = ring_radii[-1]``
+
+        If both sigmas and ring_radii are given, they must satisfy the
+        following predicate since no radius is needed for the center
+        histogram.
+            ``len(ring_radii) == len(sigmas)+1``
+    visualize : bool, optional
+        Generate a visualization of the DAISY descriptors
+
+
+    Returns
+    -------
+    descs : array
+        Grid of DAISY descriptors for the given image as an array
+        dimensionality  (P, Q, R) where
+          | ``P = ceil((M-radius*2)/step)``
+          | ``Q = ceil((N-radius*2)/step)``
+          | ``R = (rings*histograms + 1)*orientations``
+    descs_img : (M, N, 3) array (only if visualize==True)
+        Visualization of the DAISY descriptors.
+
+    References
+    ----------
+    .. [1] Tola et al. "Daisy: An efficient dense descriptor applied to wide-
+           baseline stereo." Pattern Analysis and Machine Intelligence, IEEE
+           Transactions on 32.5 (2010): 815-830.
+    .. [2] http://cvlab.epfl.ch/alumni/tola/daisy.html
+    '''
+
+    # Validate image format.
+    if img.ndim > 2:
+        raise ValueError('Only grey-level images are supported.')
+    if img.dtype.kind != 'f':
+        img = img_as_float(img)
+
+    # Validate parameters.
+    if sigmas is not None and ring_radii is not None \
+            and len(sigmas) - 1 != len(ring_radii):
+        raise ValueError('len(sigmas)-1 != len(ring_radii)')
+    if ring_radii is not None:
+        rings = len(ring_radii)
+        radius = ring_radii[-1]
+    if sigmas is not None:
+        rings = len(sigmas) - 1
+    if sigmas is None:
+        sigmas = [radius * (i + 1) / float(2 * rings) for i in range(rings)]
+    if ring_radii is None:
+        ring_radii = [radius * (i + 1) / float(rings) for i in range(rings)]
+    if normalization not in ['l1', 'l2', 'daisy', 'off']:
+        raise ValueError('Invalid normalization method.')
+
+    # Compute image derivatives.
+    dx = np.zeros(img.shape)
+    dy = np.zeros(img.shape)
+    dx[:, :-1] = np.diff(img, n=1, axis=1)
+    dy[:-1, :] = np.diff(img, n=1, axis=0)
+
+    # Compute gradient orientation and magnitude and their contribution
+    # to the histograms.
+    grad_mag = sqrt(dx ** 2 + dy ** 2)
+    grad_ori = arctan2(dy, dx)
+    orientation_kappa = orientations / pi
+    orientation_angles = [2 * o * pi / orientations - pi
+                          for o in range(orientations)]
+    hist = np.empty((orientations,) + img.shape, dtype=float)
+    for i, o in enumerate(orientation_angles):
+        # Weigh bin contribution by the circular normal distribution
+        hist[i, :, :] = exp(orientation_kappa * cos(grad_ori - o))
+        # Weigh bin contribution by the gradient magnitude
+        hist[i, :, :] = np.multiply(hist[i, :, :], grad_mag)
+
+    # Smooth orientation histograms for the center and all rings.
+    sigmas = [sigmas[0]] + sigmas
+    hist_smooth = np.empty((rings + 1,) + hist.shape, dtype=float)
+    for i in range(rings + 1):
+        for j in range(orientations):
+            hist_smooth[i, j, :, :] = gaussian_filter(hist[j, :, :],
+                                                      sigma=sigmas[i])
+
+    # Assemble descriptor grid.
+    theta = [2 * pi * j / histograms for j in range(histograms)]
+    desc_dims = (rings * histograms + 1) * orientations
+    descs = np.empty((desc_dims, img.shape[0] - 2 * radius,
+                      img.shape[1] - 2 * radius))
+    descs[:orientations, :, :] = hist_smooth[0, :, radius:-radius,
+                                             radius:-radius]
+    idx = orientations
+    for i in range(rings):
+        for j in range(histograms):
+            y_min = radius + int(round(ring_radii[i] * sin(theta[j])))
+            y_max = descs.shape[1] + y_min
+            x_min = radius + int(round(ring_radii[i] * cos(theta[j])))
+            x_max = descs.shape[2] + x_min
+            descs[idx:idx + orientations, :, :] = hist_smooth[i + 1, :,
+                                                              y_min:y_max,
+                                                              x_min:x_max]
+            idx += orientations
+    descs = descs[:, ::step, ::step]
+    descs = descs.swapaxes(0, 1).swapaxes(1, 2)
+
+    # Normalize descriptors.
+    if normalization != 'off':
+        descs += 1e-10
+        if normalization == 'l1':
+            descs /= np.sum(descs, axis=2)[:, :, np.newaxis]
+        elif normalization == 'l2':
+            descs /= sqrt(np.sum(descs ** 2, axis=2))[:, :, np.newaxis]
+        elif normalization == 'daisy':
+            for i in range(0, desc_dims, orientations):
+                norms = sqrt(np.sum(descs[:, :, i:i + orientations] ** 2,
+                                    axis=2))
+                descs[:, :, i:i + orientations] /= norms[:, :, np.newaxis]
+
+    if visualize:
+        descs_img = skimage.color.gray2rgb(img)
+        for i in range(descs.shape[0]):
+            for j in range(descs.shape[1]):
+                # Draw center histogram sigma
+                color = (1, 0, 0)
+                desc_y = i * step + radius
+                desc_x = j * step + radius
+                coords = draw.circle_perimeter(desc_y, desc_x, sigmas[0])
+                draw.set_color(descs_img, coords, color)
+                max_bin = np.max(descs[i, j, :])
+                for o_num, o in enumerate(orientation_angles):
+                    # Draw center histogram bins
+                    bin_size = descs[i, j, o_num] / max_bin
+                    dy = sigmas[0] * bin_size * sin(o)
+                    dx = sigmas[0] * bin_size * cos(o)
+                    coords = draw.line(desc_y, desc_x, desc_y + dy,
+                                       desc_x + dx)
+                    draw.set_color(descs_img, coords, color)
+                for r_num, r in enumerate(ring_radii):
+                    color_offset = float(1 + r_num) / rings
+                    color = (1 - color_offset, 1, color_offset)
+                    for t_num, t in enumerate(theta):
+                        # Draw ring histogram sigmas
+                        hist_y = desc_y + int(round(r * sin(t)))
+                        hist_x = desc_x + int(round(r * cos(t)))
+                        coords = draw.circle_perimeter(hist_y, hist_x,
+                                                       sigmas[r_num + 1])
+                        draw.set_color(descs_img, coords, color)
+                        for o_num, o in enumerate(orientation_angles):
+                            # Draw histogram bins
+                            bin_size = descs[i, j, orientations + r_num *
+                                             histograms * orientations +
+                                             t_num * orientations + o_num]
+                            bin_size /= max_bin
+                            dy = sigmas[r_num + 1] * bin_size * sin(o)
+                            dx = sigmas[r_num + 1] * bin_size * cos(o)
+                            coords = draw.line(hist_y, hist_x, hist_y + dy,
+                                               hist_x + dx)
+                            draw.set_color(descs_img, coords, color)
+        return descs, descs_img
+    else:
+        return descs
diff --git a/skimage/feature/tests/test_daisy.py b/skimage/feature/tests/test_daisy.py
new file mode 100644
index 00000000..40781a64
--- /dev/null
+++ b/skimage/feature/tests/test_daisy.py
@@ -0,0 +1,95 @@
+import numpy as np
+from numpy.testing import assert_raises, assert_almost_equal
+from numpy import sqrt, ceil
+
+from skimage import data
+from skimage import img_as_float
+from skimage.feature import daisy
+
+
+def test_daisy_color_image_unsupported_error():
+    img = np.zeros((20, 20, 3))
+    assert_raises(ValueError, daisy, img)
+
+
+def test_daisy_desc_dims():
+    img = img_as_float(data.lena()[:128, :128].mean(axis=2))
+    rings = 2
+    histograms = 4
+    orientations = 3
+    descs = daisy(img, rings=rings, histograms=histograms,
+                  orientations=orientations)
+    assert(descs.shape[2] == (rings * histograms + 1) * orientations)
+
+    rings = 4
+    histograms = 5
+    orientations = 13
+    descs = daisy(img, rings=rings, histograms=histograms,
+                  orientations=orientations)
+    assert(descs.shape[2] == (rings * histograms + 1) * orientations)
+
+
+def test_descs_shape():
+    img = img_as_float(data.lena()[:256, :256].mean(axis=2))
+    radius = 20
+    step = 8
+    descs = daisy(img, radius=radius, step=step)
+    assert(descs.shape[0] == ceil((img.shape[0] - radius * 2) / float(step)))
+    assert(descs.shape[1] == ceil((img.shape[1] - radius * 2) / float(step)))
+
+    img = img[:-1, :-2]
+    radius = 5
+    step = 3
+    descs = daisy(img, radius=radius, step=step)
+    assert(descs.shape[0] == ceil((img.shape[0] - radius * 2) / float(step)))
+    assert(descs.shape[1] == ceil((img.shape[1] - radius * 2) / float(step)))
+
+
+def test_daisy_incompatible_sigmas_and_radii():
+    img = img_as_float(data.lena()[:128, :128].mean(axis=2))
+    sigmas = [1, 2]
+    radii = [1, 2]
+    assert_raises(ValueError, daisy, img, sigmas=sigmas, ring_radii=radii)
+
+
+def test_daisy_normalization():
+    img = img_as_float(data.lena()[:64, :64].mean(axis=2))
+
+    descs = daisy(img, normalization='l1')
+    for i in range(descs.shape[0]):
+        for j in range(descs.shape[1]):
+            assert_almost_equal(np.sum(descs[i, j, :]), 1)
+    descs_ = daisy(img)
+    assert_almost_equal(descs, descs_)
+
+    descs = daisy(img, normalization='l2')
+    for i in range(descs.shape[0]):
+        for j in range(descs.shape[1]):
+            assert_almost_equal(sqrt(np.sum(descs[i, j, :] ** 2)), 1)
+
+    orientations = 8
+    descs = daisy(img, orientations=orientations, normalization='daisy')
+    desc_dims = descs.shape[2]
+    for i in range(descs.shape[0]):
+        for j in range(descs.shape[1]):
+            for k in range(0, desc_dims, orientations):
+                assert_almost_equal(sqrt(np.sum(
+                    descs[i, j, k:k + orientations] ** 2)), 1)
+
+    img = np.zeros((50, 50))
+    descs = daisy(img, normalization='off')
+    for i in range(descs.shape[0]):
+        for j in range(descs.shape[1]):
+            assert_almost_equal(np.sum(descs[i, j, :]), 0)
+
+    assert_raises(ValueError, daisy, img, normalization='does_not_exist')
+
+
+def test_daisy_visualization():
+    img = img_as_float(data.lena()[:128, :128].mean(axis=2))
+    descs, descs_img = daisy(img, visualize=True)
+    assert(descs_img.shape == (128, 128, 3))
+
+if __name__ == '__main__':
+    from numpy import testing
+    testing.run_module_suite()