From 2eb03c0fc95c8edd61c4e14fea7aee26919223af Mon Sep 17 00:00:00 2001
From: Anders Boesen Lindbo Larsen <abæll@imm.dtu.dk>
Date: Sat, 8 Dec 2012 18:26:13 +0100
Subject: [PATCH 01/10] Add dense DAISY feature description.

---
 skimage/feature/__init__.py         |   1 +
 skimage/feature/_daisy.py           | 167 ++++++++++++++++++++++++++++
 skimage/feature/tests/test_daisy.py |  88 +++++++++++++++
 3 files changed, 256 insertions(+)
 create mode 100644 skimage/feature/_daisy.py
 create mode 100644 skimage/feature/tests/test_daisy.py

diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py
index 1597e9a8..b1dabad3 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..7fb0364c
--- /dev/null
+++ b/skimage/feature/_daisy.py
@@ -0,0 +1,167 @@
+import numpy as np
+import scipy as sp
+from numpy.linalg import norm
+from scipy import sqrt, pi, arctan2, cos, sin, exp
+from scipy.ndimage import gaussian_filter
+from scipy.special import iv
+
+def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
+        normalization='l1', sigmas=None, ring_radii=None):
+    '''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 specifies the spatial smoothing of the center histogram
+        and the last sigma value specifies 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
+          len(ring_radii) == len(sigmas)+1
+        since no radius is needed for the center histogram.
+
+    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
+
+    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 == 'u':
+        img = img.astype(float)
+        img = img/255.
+
+    # Validate parameters.
+    if sigmas != None and ring_radii != None \
+            and len(sigmas)-1 != len(ring_radii):
+        raise ValueError('len(sigmas)-1 != len(ring_radii)')
+    if ring_radii != None:
+        rings = len(ring_radii)
+        radius = ring_radii[-1]
+    if sigmas != None:
+        rings = len(sigmas)-1
+    if sigmas == None:
+        sigmas = [radius*(i+1)/float(2*rings) for i in range(rings)]
+    if ring_radii == 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)
+    hist_sigma = pi/orientations
+    kappa = 1./hist_sigma;
+    bessel = iv(0, kappa)
+    hist = np.empty((orientations,) + img.shape, dtype=float)
+    for i in range(orientations):
+        mu = 2*i*pi/orientations-pi
+        # Weigh bin contribution by the circular normal distribution
+        hist[i,:,:] = exp(kappa*cos(grad_ori-mu))/(2*pi*bessel)
+        # 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]
+
+    return descs
+
diff --git a/skimage/feature/tests/test_daisy.py b/skimage/feature/tests/test_daisy.py
new file mode 100644
index 00000000..9e4c663a
--- /dev/null
+++ b/skimage/feature/tests/test_daisy.py
@@ -0,0 +1,88 @@
+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')
+
+
+if __name__ == '__main__':
+    from numpy import testing
+    testing.run_module_suite()

From 979a15cd9a417b030ef42b386646b12b8db218f5 Mon Sep 17 00:00:00 2001
From: Anders Boesen Lindbo Larsen <abæll@imm.dtu.dk>
Date: Sat, 8 Dec 2012 22:37:36 +0100
Subject: [PATCH 02/10] PEP8 style corrections.

---
 skimage/feature/_daisy.py           | 150 ++++++++++++++--------------
 skimage/feature/tests/test_daisy.py |  28 +++---
 2 files changed, 89 insertions(+), 89 deletions(-)

diff --git a/skimage/feature/_daisy.py b/skimage/feature/_daisy.py
index 7fb0364c..afe059be 100644
--- a/skimage/feature/_daisy.py
+++ b/skimage/feature/_daisy.py
@@ -1,27 +1,26 @@
 import numpy as np
-import scipy as sp
-from numpy.linalg import norm
 from scipy import sqrt, pi, arctan2, cos, sin, exp
 from scipy.ndimage import gaussian_filter
 from scipy.special import iv
 
+
 def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
-        normalization='l1', sigmas=None, ring_radii=None):
+          normalization='l1', sigmas=None, ring_radii=None):
     '''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.
+    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.
+    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
     ----------
@@ -39,26 +38,25 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
         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.
+            * '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 specifies the spatial smoothing of the center histogram
-        and the last sigma value specifies the spatial smoothing of the
-        outermost ring.
-        Specifying sigmas overrides the following parameter:
-          rings = len(sigmas)-1
+        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]
+        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
-          len(ring_radii) == len(sigmas)+1
+            len(ring_radii) == len(sigmas)+1
         since no radius is needed for the center histogram.
 
     Returns
@@ -66,16 +64,16 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
     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
+            P = ceil((M-radius*2)/step)
+            Q = ceil((N-radius*2)/step)
+            R = (rings*histograms + 1)*orientations
 
     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
+    [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.
@@ -83,21 +81,21 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
         raise ValueError('Only grey-level images are supported.')
     if img.dtype.kind == 'u':
         img = img.astype(float)
-        img = img/255.
+        img = img / 255.
 
     # Validate parameters.
-    if sigmas != None and ring_radii != None \
-            and len(sigmas)-1 != len(ring_radii):
+    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 != None:
+    if ring_radii is not None:
         rings = len(ring_radii)
         radius = ring_radii[-1]
-    if sigmas != None:
-        rings = len(sigmas)-1
-    if sigmas == None:
-        sigmas = [radius*(i+1)/float(2*rings) for i in range(rings)]
-    if ring_radii == None:
-        ring_radii = [radius*(i+1)/float(rings) for i in range(rings)]
+    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.')
 
@@ -109,59 +107,59 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
 
     # Compute gradient orientation and magnitude and their contribution
     # to the histograms.
-    grad_mag = sqrt(dx**2 + dy**2)
+    grad_mag = sqrt(dx ** 2 + dy ** 2)
     grad_ori = arctan2(dy, dx)
-    hist_sigma = pi/orientations
-    kappa = 1./hist_sigma;
+    hist_sigma = pi / orientations
+    kappa = 1. / hist_sigma
     bessel = iv(0, kappa)
     hist = np.empty((orientations,) + img.shape, dtype=float)
     for i in range(orientations):
-        mu = 2*i*pi/orientations-pi
+        mu = 2 * i * pi / orientations - pi
         # Weigh bin contribution by the circular normal distribution
-        hist[i,:,:] = exp(kappa*cos(grad_ori-mu))/(2*pi*bessel)
+        hist[i, :, :] = exp(kappa * cos(grad_ori - mu)) / (2 * pi * bessel)
         # Weigh bin contribution by the gradient magnitude
-        hist[i,:,:] = np.multiply(hist[i,:,:], grad_mag)
+        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):
+    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])
+            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]
+    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_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_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]
+            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)
+    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]
+            descs /= np.sum(descs, axis=2)[:, :, np.newaxis]
         elif normalization == 'l2':
-            descs /=  sqrt(np.sum(descs**2, axis=2))[:,:,np.newaxis]
+            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]
+                norms = sqrt(np.sum(descs[:, :, i:i + orientations] ** 2,
+                                    axis=2))
+                descs[:, :, i:i + orientations] /= norms[:, :, np.newaxis]
 
     return descs
-
diff --git a/skimage/feature/tests/test_daisy.py b/skimage/feature/tests/test_daisy.py
index 9e4c663a..b0dd28de 100644
--- a/skimage/feature/tests/test_daisy.py
+++ b/skimage/feature/tests/test_daisy.py
@@ -17,14 +17,16 @@ def test_daisy_desc_dims():
     rings = 2
     histograms = 4
     orientations = 3
-    descs = daisy(img, rings=rings, histograms=histograms, orientations=orientations)
-    assert(descs.shape[2] == (rings*histograms + 1)*orientations)
+    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)
+    descs = daisy(img, rings=rings, histograms=histograms,
+                  orientations=orientations)
+    assert(descs.shape[2] == (rings * histograms + 1) * orientations)
 
 
 def test_descs_shape():
@@ -32,15 +34,15 @@ def test_descs_shape():
     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)))
+    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]
+    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)))
+    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():
@@ -56,14 +58,14 @@ def test_daisy_normalization():
     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)
+            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)
+            assert_almost_equal(sqrt(np.sum(descs[i, j, :] ** 2)), 1)
 
     orientations = 8
     descs = daisy(img, orientations=orientations, normalization='daisy')
@@ -72,13 +74,13 @@ def test_daisy_normalization():
         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)
+                    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_almost_equal(np.sum(descs[i, j, :]), 0)
 
     assert_raises(ValueError, daisy, img, normalization='does_not_exist')
 

From e5fb8a2a72023f7ea7bdae390e559a7bf6e8cb7c Mon Sep 17 00:00:00 2001
From: Anders Boesen Lindbo Larsen <abæll@imm.dtu.dk>
Date: Sun, 9 Dec 2012 12:53:52 +0100
Subject: [PATCH 03/10] Better Sphinx documentation.

---
 skimage/feature/_daisy.py | 69 +++++++++++++++++++++------------------
 1 file changed, 38 insertions(+), 31 deletions(-)

diff --git a/skimage/feature/_daisy.py b/skimage/feature/_daisy.py
index afe059be..b0cc6748 100644
--- a/skimage/feature/_daisy.py
+++ b/skimage/feature/_daisy.py
@@ -2,6 +2,7 @@ import numpy as np
 from scipy import sqrt, pi, arctan2, cos, sin, exp
 from scipy.ndimage import gaussian_filter
 from scipy.special import iv
+from skimage import img_as_float
 
 
 def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
@@ -12,15 +13,16 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
     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
+    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.
+
+      * 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
     ----------
@@ -36,52 +38,57 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
         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.
+    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
+        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
-            len(ring_radii) == len(sigmas)+1
-        since no radius is needed for the center histogram.
+        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``
+
 
     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
+          | ``P = ceil((M-radius*2)/step)``
+          | ``Q = ceil((N-radius*2)/step)``
+          | ``R = (rings*histograms + 1)*orientations``
 
     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
+    .. [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 == 'u':
-        img = img.astype(float)
-        img = img / 255.
+        img = img_as_float(img)
 
     # Validate parameters.
     if sigmas is not None and ring_radii is not None \

From 2849cbf8587de52d7d3090c3a8a96da2e51badb4 Mon Sep 17 00:00:00 2001
From: Anders Boesen Lindbo Larsen <abæll@imm.dtu.dk>
Date: Sat, 15 Dec 2012 18:43:33 +0100
Subject: [PATCH 04/10] Add circle perimeter to draw module.

---
 skimage/draw/__init__.py        |  2 +-
 skimage/draw/_draw.pyx          | 39 +++++++++++++++++++++++++++++++++
 skimage/draw/tests/test_draw.py | 33 +++++++++++++++++++++++++++-
 3 files changed, 72 insertions(+), 2 deletions(-)

diff --git a/skimage/draw/__init__.py b/skimage/draw/__init__.py
index 4eac9b63..b213f93f 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
 bresenham = line
diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx
index ca0c502a..b5b96a2b 100644
--- a/skimage/draw/_draw.pyx
+++ b/skimage/draw/_draw.pyx
@@ -187,3 +187,42 @@ 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
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')
 

From 6c19cf0693ca85eff0d43ff68693fcf3f2d83b48 Mon Sep 17 00:00:00 2001
From: Anders Boesen Lindbo Larsen <abæll@imm.dtu.dk>
Date: Sat, 15 Dec 2012 18:44:23 +0100
Subject: [PATCH 05/10] Add DAISY visualization + example plot.

---
 doc/examples/plot_daisy.py          | 28 ++++++++++++
 doc/examples/plot_shapes.py         |  8 +++-
 skimage/feature/_daisy.py           | 70 ++++++++++++++++++++++++++---
 skimage/feature/tests/test_daisy.py |  5 +++
 4 files changed, 103 insertions(+), 8 deletions(-)
 create mode 100644 doc/examples/plot_daisy.py

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/feature/_daisy.py b/skimage/feature/_daisy.py
index b0cc6748..85e259d0 100644
--- a/skimage/feature/_daisy.py
+++ b/skimage/feature/_daisy.py
@@ -6,7 +6,7 @@ from skimage import img_as_float
 
 
 def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
-          normalization='l1', sigmas=None, ring_radii=None):
+          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
@@ -65,6 +65,8 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
         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
@@ -75,10 +77,12 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
           | ``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-
+    .. [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
@@ -120,10 +124,11 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
     kappa = 1. / hist_sigma
     bessel = iv(0, kappa)
     hist = np.empty((orientations,) + img.shape, dtype=float)
-    for i in range(orientations):
-        mu = 2 * i * pi / orientations - pi
+    orientation_angles = [2 * o * pi / orientations - pi
+                          for o in range(orientations)]
+    for i, o in enumerate(orientation_angles):
         # Weigh bin contribution by the circular normal distribution
-        hist[i, :, :] = exp(kappa * cos(grad_ori - mu)) / (2 * pi * bessel)
+        hist[i, :, :] = exp(kappa * cos(grad_ori - o)) / (2 * pi * bessel)
         # Weigh bin contribution by the gradient magnitude
         hist[i, :, :] = np.multiply(hist[i, :, :], grad_mag)
 
@@ -169,4 +174,57 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
                                     axis=2))
                 descs[:, :, i:i + orientations] /= norms[:, :, np.newaxis]
 
-    return descs
+    if visualize:
+        from skimage import draw
+        from skimage import color
+        descs_img = 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])
+                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)
+                    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])
+                        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)
+                            set_color(descs_img, coords, color)
+        return descs, descs_img
+    else:
+        return descs
+
+
+def set_color(img, coords, color):
+    ''' Set pixel color at the given coordiantes in the image.'''
+    coords = filter(lambda (y, x): y >= 0 and y < img.shape[0] and x >= 0
+                    and x < img.shape[1], zip(*coords))
+    if coords:
+        rr, cc = zip(*coords)
+        img[rr, cc, :] = color
diff --git a/skimage/feature/tests/test_daisy.py b/skimage/feature/tests/test_daisy.py
index b0dd28de..40781a64 100644
--- a/skimage/feature/tests/test_daisy.py
+++ b/skimage/feature/tests/test_daisy.py
@@ -85,6 +85,11 @@ def test_daisy_normalization():
     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()

From 7aa943530144fb9478e0b09d57c94ea42efc5980 Mon Sep 17 00:00:00 2001
From: Anders Boesen Lindbo Larsen <abæll@imm.dtu.dk>
Date: Sun, 16 Dec 2012 13:47:35 +0100
Subject: [PATCH 06/10] Handle all non-float images.

---
 skimage/feature/_daisy.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/feature/_daisy.py b/skimage/feature/_daisy.py
index 85e259d0..d011d693 100644
--- a/skimage/feature/_daisy.py
+++ b/skimage/feature/_daisy.py
@@ -91,7 +91,7 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
     # Validate image format.
     if img.ndim > 2:
         raise ValueError('Only grey-level images are supported.')
-    if img.dtype.kind == 'u':
+    if img.dtype.kind != 'f':
         img = img_as_float(img)
 
     # Validate parameters.

From 097f25dbb5fe599491b067b5263f350990ada5d5 Mon Sep 17 00:00:00 2001
From: Anders Boesen Lindbo Larsen <abæll@imm.dtu.dk>
Date: Sun, 16 Dec 2012 16:44:56 +0100
Subject: [PATCH 07/10] Move set_color() to draw module.

---
 skimage/draw/__init__.py  |  2 +-
 skimage/draw/_draw.pyx    | 27 +++++++++++++++++++++++++++
 skimage/feature/_daisy.py | 24 +++++++-----------------
 3 files changed, 35 insertions(+), 18 deletions(-)

diff --git a/skimage/draw/__init__.py b/skimage/draw/__init__.py
index b213f93f..4b30ea18 100644
--- a/skimage/draw/__init__.py
+++ b/skimage/draw/__init__.py
@@ -1,2 +1,2 @@
-from ._draw import line, polygon, ellipse, circle, circle_perimeter
+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 b5b96a2b..509796bb 100644
--- a/skimage/draw/_draw.pyx
+++ b/skimage/draw/_draw.pyx
@@ -226,3 +226,30 @@ def circle_perimeter(int cy, int cx, int radius):
         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/feature/_daisy.py b/skimage/feature/_daisy.py
index d011d693..07caf22b 100644
--- a/skimage/feature/_daisy.py
+++ b/skimage/feature/_daisy.py
@@ -2,7 +2,8 @@ import numpy as np
 from scipy import sqrt, pi, arctan2, cos, sin, exp
 from scipy.ndimage import gaussian_filter
 from scipy.special import iv
-from skimage import img_as_float
+import skimage
+from skimage import img_as_float, draw
 
 
 def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
@@ -175,9 +176,7 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
                 descs[:, :, i:i + orientations] /= norms[:, :, np.newaxis]
 
     if visualize:
-        from skimage import draw
-        from skimage import color
-        descs_img = color.gray2rgb(img)
+        descs_img = skimage.color.gray2rgb(img)
         for i in range(descs.shape[0]):
             for j in range(descs.shape[1]):
                 # Draw center histogram sigma
@@ -185,7 +184,7 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
                 desc_y = i * step + radius
                 desc_x = j * step + radius
                 coords = draw.circle_perimeter(desc_y, desc_x, sigmas[0])
-                set_color(descs_img, coords, color)
+                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
@@ -194,7 +193,7 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
                     dx = sigmas[0] * bin_size * cos(o)
                     coords = draw.line(desc_y, desc_x, desc_y + dy,
                                        desc_x + dx)
-                    set_color(descs_img, coords, color)
+                    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)
@@ -204,7 +203,7 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
                         hist_x = desc_x + int(round(r * cos(t)))
                         coords = draw.circle_perimeter(hist_y, hist_x,
                                                        sigmas[r_num + 1])
-                        set_color(descs_img, coords, color)
+                        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 *
@@ -215,16 +214,7 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
                             dx = sigmas[r_num + 1] * bin_size * cos(o)
                             coords = draw.line(hist_y, hist_x, hist_y + dy,
                                                hist_x + dx)
-                            set_color(descs_img, coords, color)
+                            draw.set_color(descs_img, coords, color)
         return descs, descs_img
     else:
         return descs
-
-
-def set_color(img, coords, color):
-    ''' Set pixel color at the given coordiantes in the image.'''
-    coords = filter(lambda (y, x): y >= 0 and y < img.shape[0] and x >= 0
-                    and x < img.shape[1], zip(*coords))
-    if coords:
-        rr, cc = zip(*coords)
-        img[rr, cc, :] = color

From 3b43cf0dcfa96b453c7ab3f3b47d4ea9475ae567 Mon Sep 17 00:00:00 2001
From: Anders Boesen Lindbo Larsen <abæll@imm.dtu.dk>
Date: Sun, 16 Dec 2012 17:26:32 +0100
Subject: [PATCH 08/10] No need to normalize angular smoothing.

---
 skimage/feature/_daisy.py | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/skimage/feature/_daisy.py b/skimage/feature/_daisy.py
index 07caf22b..f5599dd5 100644
--- a/skimage/feature/_daisy.py
+++ b/skimage/feature/_daisy.py
@@ -1,7 +1,6 @@
 import numpy as np
 from scipy import sqrt, pi, arctan2, cos, sin, exp
 from scipy.ndimage import gaussian_filter
-from scipy.special import iv
 import skimage
 from skimage import img_as_float, draw
 
@@ -121,15 +120,13 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
     # to the histograms.
     grad_mag = sqrt(dx ** 2 + dy ** 2)
     grad_ori = arctan2(dy, dx)
-    hist_sigma = pi / orientations
-    kappa = 1. / hist_sigma
-    bessel = iv(0, kappa)
-    hist = np.empty((orientations,) + img.shape, dtype=float)
+    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(kappa * cos(grad_ori - o)) / (2 * pi * bessel)
+        hist[i, :, :] = exp(orientation_kappa * cos(grad_ori - o))
         # Weigh bin contribution by the gradient magnitude
         hist[i, :, :] = np.multiply(hist[i, :, :], grad_mag)
 

From c276ec06342b194fffb2cfd5efa4a87416bfcbaa Mon Sep 17 00:00:00 2001
From: Anders Boesen Lindbo Larsen <abæll@imm.dtu.dk>
Date: Mon, 17 Dec 2012 11:54:37 +0100
Subject: [PATCH 09/10] Fixed submodule import.

---
 skimage/feature/_daisy.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/feature/_daisy.py b/skimage/feature/_daisy.py
index f5599dd5..f83605b5 100644
--- a/skimage/feature/_daisy.py
+++ b/skimage/feature/_daisy.py
@@ -1,7 +1,7 @@
 import numpy as np
 from scipy import sqrt, pi, arctan2, cos, sin, exp
 from scipy.ndimage import gaussian_filter
-import skimage
+import skimage.color
 from skimage import img_as_float, draw
 
 

From 8b6b51ec04fda676de5fed7dd9b4c67809fe3aa2 Mon Sep 17 00:00:00 2001
From: Anders Boesen Lindbo Larsen <abæll@imm.dtu.dk>
Date: Thu, 27 Dec 2012 10:49:29 +0100
Subject: [PATCH 10/10] Add myself as contributor.

---
 CONTRIBUTORS.txt | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index 012fd53d..10f2a94a 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -123,3 +123,6 @@
 
 - Luis Pedro Coelho
   imread plugin
+
+- Anders Boesen Lindbo Larsen
+  Dense DAISY feature description, circle perimeter drawing.