Expose Siftdescriptor

conda/meta.yaml

@@ -15,6 +15,7 @@ requirements:
     - python
     - vlfeat 0.9.*
     - numpy >=1.10
+    - scipy
 
 test:
   requires:

cyvlfeat/sift/cysift.pyx

@@ -142,6 +142,72 @@ cdef int korder(const void *a, const void *b) nogil:
     if x > 0: return +1
     return 0
 
+@cython.boundscheck(False)
+cpdef cy_siftdescriptor(np.ndarray[float, ndim=3, mode='c'] in_grad,
+                    np.ndarray[float, ndim=2, mode='c'] in_frames,
+                    int magnification, bint float_descriptors, float norm_threshold, bint verbose):
+
+
+    cdef:
+        int i = 0, j = 0
+        int m = in_grad.shape[0]
+        int n = in_grad.shape[1]
+        int o = in_grad.shape[2]
+        int num_frames = in_frames.shape[0]
+        double y, x, s, th
+
+        #create a filter to process the image
+        VlSiftFilt *filt = vl_sift_new(m, n, -1, -1, 0)
+
+
+        # Create empty 2D output array
+        np.ndarray[float, ndim=2, mode='c'] out_descriptors = np.empty(
+            (num_frames, 128), dtype=np.float32, order='C')
+        float *flat_descriptors = &out_descriptors[0, 0]
+
+        vl_sift_pix[128] single_descriptor_arr
+
+    out_descriptors = np.resize(out_descriptors,
+                                    (num_frames, 128))
+    flat_descriptors = &out_descriptors[0, 0]
+
+    if magnification >= 0:
+        vl_sift_set_magnif(filt, magnification)
+
+    if norm_threshold >= 0:
+        vl_sift_set_norm_thresh(filt, norm_threshold)
+
+    if verbose:
+        printf("vl_sift: filter settings:\n")
+        printf("vl_sift:   magnif                = %g\n", vl_sift_get_magnif(filt))
+        printf("vl_sift:   num of frames         = %d\n", num_frames)
+        printf("vl_sift:   float descriptor      = %d\n", float_descriptors)
+        printf("vl_sift:   norm thresh           = %g\n", vl_sift_get_norm_thresh(filt))
+
+       # MATLAB stores a two-dimensional matrix in memory as a one-
+       # dimensional array.  If the matrix is size MxN, then the
+       # first M elements of the one-dimensional array correspond to
+       # the first column of the matrix, and the next M elements
+       # correspond to the second column, etc.
+
+    for i in range(in_frames.shape[0]):
+        y = in_frames[i][0]-1
+        x = in_frames[i][1]-1
+        s = in_frames[i][2]-1
+        th = VL_PI/2 - in_frames [i][3]
+
+        vl_sift_calc_raw_descriptor(filt, &in_grad[0,0,0], single_descriptor_arr, m, n, x, y, s, th)
+
+        for j in range(128):
+            flat_descriptors[j] = min(512.0 * single_descriptor_arr[j], 255.0)
+
+    vl_sift_delete(filt)
+
+    if float_descriptors:
+        return out_descriptors
+    else:
+        return out_descriptors.astype(np.uint8)
+
 
 @cython.boundscheck(False)
 cpdef cy_sift(np.ndarray[float, ndim=2, mode='c'] data, int n_octaves,

cyvlfeat/sift/siftdescriptor.py

@@ -0,0 +1,104 @@
+import numpy as np
+from cyvlfeat.sift.cysift import cy_siftdescriptor
+
+
+def siftdescriptor(gradient_image, frames, magnification=3,
+                   float_descriptors=False, norm_thresh=None, verbose=False):
+    r"""
+    Calculates the SIFT descriptors of the keypoints ``frames`` on the
+     pre-processed image ``gradient_image``.
+
+    In order to match the standard SIFT descriptor, the gradient should be
+    calculated after mapping the image to the keypoint scale. This is obtained
+    by smoothing the image by a a Gaussian kernel of variance equal to the scale
+    of the keypoint. Additionally, SIFT assumes that the input image is
+    pre-smoothed at scale 0.5 (this roughly compensates for the effect of the
+    CCD integrators), so the amount of smoothing that needs to be applied is
+    slightly less.
+
+    Parameters
+    ----------
+    grad : [2, H, W] `float32` `ndarray`
+        gradient_image is a 2xMxN array. The first channel
+        ``gradient_image[0, :, :]`` contains the
+        modulus of gradient of the original image modulus. The second channel
+        `gradient_image[1, :, :]`` contains the gradient angle (measured in
+        radians, clockwise, starting from the X axis -- this assumes that the Y
+        axis points down).
+    frames : `[F, 4]` `float32` `ndarray`, optional
+        If specified, set the frames to use (bypass the detector). If frames are
+        not passed in order of increasing scale, they are re-orderded. A frame
+        is a vector of length 4 ``[Y, X, S, TH]``, representing a disk of center
+        frames[:2], scale frames[2] and orientation frames[3].
+    magnification : `int`, optional
+        Set the descriptor magnification factor. The scale of the keypoint is
+        multiplied by this factor to obtain the width (in pixels) of the spatial
+        bins. For instance, if there are there are 4 spatial bins along each
+        spatial direction, the ``side`` of the descriptor is approximately ``4 *
+        magnification``.
+    float_descriptors : `bool`, optional
+        If ``True``, the descriptor are returned in floating point rather than
+        integer format.
+    norm_thresh : `float`, optional
+        Set the minimum l2-norm of the descriptors before normalization.
+        Descriptors below the threshold are set to zero. If ``None``,
+        norm_thresh is ``-inf``.
+    verbose : `bool`, optional
+        If ``True``, be verbose.
+
+    Returns
+    -------
+    descriptors : `(F, 128)` `uint8` or `float32` `ndarray`, optional
+        ``F`` is the number of keypoints (frames) used. The 128 length vectors
+        per keypoint extracted. ``uint8`` by default.
+
+
+    Examples
+    --------
+    >>> import scipy.ndimage
+    >>> import numpy as np
+    >>> from cyvlfeat.sift import siftdescriptor
+    >>> from cyvlfeat.test_util import lena
+    >>> img = lena().astype(np.float32)
+    >>> # Create a single frame in the center of the image
+    >>> frames = np.array([[256, 256, 1.0, np.pi / 2]])
+    >>> sigma = np.sqrt(frames[0, 2] ** 2 - 0.25)  # 0.25 = 0.5^2
+    >>> img_smooth = scipy.ndimage.filters.gaussian_filter(img, sigma)
+    >>> y, x = np.gradient(img_smooth)
+    >>> mod = np.sqrt(x * x + y * y)
+    >>> ang = np.arctan2(y, x)
+    >>> gradient_image = np.stack((mod, ang), axis=0)
+    >>>
+    >>> d = siftdescriptor(gradient_image, frames, verbose=True)
+
+    Notes
+    -----
+    1. The above fragment generates results which are very close
+       but not identical to the output of ``sift`` as the latter
+       samples the scale space at finite steps.
+
+    2. For object categorization is sometimes useful to compute
+       SIFT descriptors without smoothing the image.
+    """
+
+    # Validate Gradient array size
+    if gradient_image.ndim != 3:
+        raise ValueError('Only 3D arrays are supported')
+
+    # Validate magnification
+    if magnification < 0:
+        raise ValueError('Magnification must be non-negative')
+
+    # Validate norm_thresh
+    if norm_thresh is not None and norm_thresh < 0:
+        raise ValueError('norm_thresh must be >= 0')
+    if norm_thresh is None:
+        norm_thresh = -1
+
+    # Ensure types are correct before passing to Cython
+    gradient_image = np.require(gradient_image, dtype=np.float32,
+                                requirements='C')
+    frames = np.require(frames, dtype=np.float32, requirements='C')
+
+    return cy_siftdescriptor(gradient_image, frames, magnification,
+                             float_descriptors, norm_thresh, verbose)

cyvlfeat/sift/tests/sift_test.py

@@ -1,6 +1,8 @@
 from __future__ import division
 from cyvlfeat.sift.dsift import dsift
 from cyvlfeat.sift.sift import sift
+from cyvlfeat.sift.siftdescriptor import siftdescriptor
+import scipy.ndimage as image
 import numpy as np
 from numpy.testing import assert_allclose
 from cyvlfeat.test_util import lena
@@ -9,6 +11,34 @@
 img = lena().astype(np.float32)
 half_img = img[:, :256]
 
+# Siftdescriptor test.
+result = sift(img, compute_descriptor=True)
+frame = result[0]
+sigma = np.sqrt(np.power(frame[3][2], 2) - np.power(0.5, 2))
+# smoothing
+img_smooth = image.filters.gaussian_filter(img, sigma)
+x, y = np.gradient(img_smooth)
+mod = np.sqrt(np.power(x, 2) + np.power(y, 2))
+ang = np.arctan2(y, x) * 180 / np.pi
+# Stack arrays in sequence depth wise (along third axis).
+grads = np.rollaxis(np.dstack((mod, ang)), 1)
+
+
+def test_siftdescriptor_non_float_descriptors():
+    descriptors = siftdescriptor(np.asfarray(grads, dtype='float'), frame, float_descriptors=False)
+    assert descriptors.dtype == np.uint8
+
+
+def test_siftdescriptor_float_descriptors():
+    descriptors = siftdescriptor(np.asfarray(grads, dtype='float'), frame, float_descriptors=True)
+    assert descriptors.dtype == np.float32
+
+
+def test_siftdescriptor_descriptors_shape():
+    descriptors = siftdescriptor(np.asfarray(grads, dtype='float'), frame)
+    assert descriptors.shape[0] == 730
+    assert descriptors.shape[1] == 128
+
 
 def test_dsift_non_float_descriptors():
     i = img.copy()