Daisy features

menpo/external/skimage/_daisy.py

@@ -0,0 +1,171 @@
+import numpy as np
+from scipy import sqrt, pi, arctan2, cos, sin, exp
+from scipy.ndimage import gaussian_filter
+
+from menpo.feature import gradient
+
+
+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
+    '''
+    # Compute image derivatives.
+    # Get number of input image's channels
+    n_channels = img.shape[-1]
+
+    # Compute image gradient
+    grad = gradient(img)
+
+    # For each pixel, select gradient with highest magnitude
+    tmp_mag = np.zeros(img.shape)
+    tmp_ori = np.zeros(img.shape)
+    for c in range(n_channels):
+        tmp_mag[..., c] = sqrt(grad[..., 2*c] ** 2 + grad[..., 2*c+1] ** 2)
+        tmp_ori[..., c] = arctan2(grad[..., 2*c+1], grad[..., 2*c])
+    grad_mag_ind = np.argmax(tmp_mag, axis=2)
+
+    # Compute gradient orientation and magnitude and their contribution
+    # to the histograms.
+    b = np.concatenate([(grad_mag_ind == i)[..., None]
+                        for i in xrange(n_channels)], axis=-1)
+    grad_mag = tmp_mag[b].reshape(tmp_mag.shape[:2])
+    grad_ori = tmp_ori[b].reshape(tmp_ori.shape[:2])
+    orientation_kappa = orientations / pi
+    orientation_angles = [2 * o * pi / orientations - pi
+                          for o in range(orientations)]
+    hist = np.empty((orientations,) + img.shape[:2], 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]
+
+    # Change axes so that the channels go to the final axis
+    descs = np.ascontiguousarray(descs)
+
+    return descs

menpo/feature/__init__.py

@@ -1,3 +1,3 @@
-from .features import gradient, hog, lbp, es, igo, no_op, gaussian_filter
+from .features import gradient, hog, lbp, es, igo, no_op, gaussian_filter, daisy
 from .predefined import sparse_hog, double_igo
 from .base import ndfeature, imgfeature

menpo/feature/features.py

@@ -165,7 +165,6 @@ def hog(pixels, mode='dense', algorithm='dalaltriggs', num_bins=9,
         Vertical window step must be > 0
     ValueError
         Window step unit must be either pixels or cells
-
     """
     # Parse options
     if mode not in ['dense', 'sparse']:
@@ -356,6 +355,7 @@ def igo(pixels, double_angles=False, verbose=False):
     #                             'original_image_channels':
     #                                 self._image.pixels.shape[2]}
 
+
 @ndfeature
 def es(image_data, verbose=False):
     r"""
@@ -416,6 +416,124 @@ def es(image_data, verbose=False):
     #                               self._image.pixels.shape[2]}
 
 
+@ndfeature
+def daisy(pixels, step=1, radius=15, rings=2, histograms=2, orientations=8,
+          normalization='l1', sigmas=None, ring_radii=None, verbose=False):
+    r"""
+    Computes a 2-dimensional Daisy features image with N*C number of channels,
+    where N is the number of channels of the original image and C is the
+    feature channels determined by the input options. Specifically,
+    C = (rings * histograms + 1) * orientations.
+
+    Parameters
+    ----------
+    pixels :  ndarray
+        The pixel data for the image, where the last axis represents the
+        number of channels.
+
+    step : `int`, Optional
+        The sampling step that defines the density of the output image.
+
+    radius : `int`, Optional
+        The radius (in pixels) of the outermost ring.
+
+    rings : `int`, Optional
+        The number of rings to be used.
+
+    histograms : `int`, Optional
+        The number of histograms sampled per ring.
+
+    orientations : `int`, Optional
+        The number of orientations (bins) per histogram.
+
+    normalization : [ 'l1', 'l2', 'daisy', None ], Optional
+        It defines how to normalize the descriptors
+        If 'l1' then L1-normalization is applied at each descriptor.
+        If 'l2' then L2-normalization is applied at each descriptor.
+        If 'daisy' then L2-normalization is applied at individual histograms.
+        If None then no normalization is employed.
+
+    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``
+
+    verbose : `bool`
+        Flag to print Daisy related information.
+
+    Raises
+    -------
+    ValueError
+        `len(sigmas)-1 != len(ring_radii)`
+    ValueError
+        Invalid normalization method.
+    """
+    from menpo.external.skimage._daisy import _daisy
+
+    # Parse options
+    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 is None:
+        normalization = 'off'
+    if normalization not in ['l1', 'l2', 'daisy', 'off']:
+        raise ValueError('Invalid normalization method.')
+
+    # Compute daisy features
+    daisy_descriptor = _daisy(pixels, step=step, radius=radius, rings=rings,
+                              histograms=histograms, orientations=orientations,
+                              normalization=normalization, sigmas=sigmas,
+                              ring_radii=ring_radii)
+
+    # print information
+    if verbose:
+        info_str = "Daisy Features:\n"
+        info_str = "{}  - Input image is {}W x {}H with {} channels.\n".format(
+            info_str, pixels.shape[1], pixels.shape[0], pixels.shape[2])
+        info_str = "{}  - Sampling step is {}.\n".format(info_str, step)
+        info_str = "{}  - Radius of {} pixels, {} rings and {} histograms " \
+                   "with {} orientations.\n".format(
+                   info_str, radius, rings, histograms, orientations)
+        if not normalization == 'off':
+            info_str = "{}  - Using {} normalization.\n".format(info_str,
+                                                                normalization)
+        else:
+            info_str = "{}  - No normalization emplyed.\n".format(info_str)
+        info_str = "{}Output image size {}W x {}H x {}.".format(
+            info_str, daisy_descriptor.shape[1], daisy_descriptor.shape[0],
+            daisy_descriptor.shape[2])
+        print(info_str)
+
+    return daisy_descriptor
+
+
 @winitfeature
 def lbp(pixels, radius=None, samples=None, mapping_type='riu2',
         window_step_vertical=1, window_step_horizontal=1,
@@ -590,6 +708,7 @@ def lbp(pixels, radius=None, samples=None, mapping_type='riu2',
     #                             'original_image_channels':
     #                                 self._image.pixels.shape[2]}
 
+
 @ndfeature
 def no_op(image_data):
     r"""

menpo/feature/test_features.py

@@ -4,7 +4,7 @@
 import math
 
 from menpo.image import Image, MaskedImage
-from menpo.feature import hog, lbp, es, igo
+from menpo.feature import hog, lbp, es, igo, daisy
 import menpo.io as mio
 
 
@@ -189,6 +189,22 @@ def test_es_channels():
         assert_allclose(es_img.n_channels, 2 * channels[i, 0])
 
 
+def test_daisy_channels():
+    n_cases = 3
+    rings = np.random.randint(1, 3, [n_cases, 1])
+    orientations = np.random.randint(1, 7, [n_cases, 1])
+    histograms = np.random.randint(1, 6, [n_cases, 1])
+    channels = np.random.randint(1, 5, [n_cases, 1])
+    for i in range(n_cases):
+        image = Image(np.random.randn(40, 40, channels[i, 0]))
+        daisy_img = daisy(image, step=4, rings=rings[i, 0],
+                          orientations=orientations[i, 0],
+                          histograms=histograms[i, 0])
+        assert_allclose(daisy_img.shape, (3, 3))
+        assert_allclose(daisy_img.n_channels,
+                        ((rings[i, 0]*histograms[i, 0]+1)*orientations[i, 0]))
+
+
 def test_igo_values():
     image = Image([[1, 2], [2, 1]])
     igo_img = igo(image)
@@ -206,7 +222,6 @@ def test_igo_values():
 
 def test_es_values():
     image = Image([[1, 2], [2, 1]])
-    print image.shape
     es_img = es(image)
     k = 1 / (2 * (2**0.5))
     res = np.array([[[k, k], [-k, k]], [[k, -k], [-k, -k]]])
@@ -218,6 +233,16 @@ def test_es_values():
     assert_allclose(es_img.pixels, res)
 
 
+def test_daisy_values():
+    image = Image([[1, 2, 3, 4], [2, 1, 3, 4], [1, 2, 3, 4], [2, 1, 3, 4]])
+    daisy_img = daisy(image, step=1, rings=2, radius=1, orientations=8,
+                      histograms=8)
+    assert_allclose(np.around(daisy_img.pixels[0, 0, 10], 6), 0.000196)
+    assert_allclose(np.around(daisy_img.pixels[0, 1, 20], 6), 0.002842)
+    assert_allclose(np.around(daisy_img.pixels[1, 0, 30], 6), 0.006205)
+    assert_allclose(np.around(daisy_img.pixels[1, 1, 40], 6), 0.001946)
+
+
 def test_lbp_values():
     image = Image([[0., 6., 0.], [5., 18., 13.], [0., 20., 0.]])
     lbp_img = lbp(image, radius=1, samples=4, mapping_type='none',