Implement SIFT feature detector and descriptor

doc/examples/features_detection/plot_sift.py

@@ -0,0 +1,88 @@
+"""
+==============================================
+SIFT feature detector and descriptor extractor
+==============================================
+
+This example demonstrates the SIFT feature detection and its description
+algorithm.
+
+The scale-invariant feature transform (SIFT) [1]_ was published in 1999 and is
+still one of the most popular feature detectors available, as its promises to
+be "invariant to image scaling, translation, and rotation, and partially
+in-variant to illumination changes and affine or 3D projection" [2]_. Its
+biggest drawback is its runtime, that's said to be "at two orders of
+magnitude" [3]_ slower than ORB, which makes it unsuitable for real-time
+applications.
+
+References
+----------
+.. [1] https://en.wikipedia.org/wiki/Scale-invariant_feature_transform
+
+.. [2] D.G. Lowe. "Object recognition from local scale-invariant
+       features", Proceedings of the Seventh IEEE International
+       Conference on Computer Vision, 1999, vol.2, pp. 1150-1157.
+       :DOI:`10.1109/ICCV.1999.790410`
+
+.. [3] Ethan Rublee, Vincent Rabaud, Kurt Konolige and Gary Bradski
+       "ORB: An efficient alternative to SIFT and SURF"
+       http://www.gwylab.com/download/ORB_2012.pdf
+"""
+import matplotlib.pyplot as plt
+
+from skimage import data
+from skimage import transform
+from skimage.color import rgb2gray
+from skimage.feature import match_descriptors, plot_matches, SIFT
+
+img1 = rgb2gray(data.astronaut())
+img2 = transform.rotate(img1, 180)
+tform = transform.AffineTransform(scale=(1.3, 1.1), rotation=0.5,
+                                  translation=(0, -200))
+img3 = transform.warp(img1, tform)
+
+descriptor_extractor = SIFT()
+
+descriptor_extractor.detect_and_extract(img1)
+keypoints1 = descriptor_extractor.keypoints
+descriptors1 = descriptor_extractor.descriptors
+
+descriptor_extractor.detect_and_extract(img2)
+keypoints2 = descriptor_extractor.keypoints
+descriptors2 = descriptor_extractor.descriptors
+
+descriptor_extractor.detect_and_extract(img3)
+keypoints3 = descriptor_extractor.keypoints
+descriptors3 = descriptor_extractor.descriptors
+
+matches12 = match_descriptors(descriptors1, descriptors2, max_ratio=0.6,
+                              cross_check=True)
+matches13 = match_descriptors(descriptors1, descriptors3, max_ratio=0.6,
+                              cross_check=True)
+fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(11, 8))
+
+plt.gray()
+
+plot_matches(ax[0, 0], img1, img2, keypoints1, keypoints2, matches12)
+ax[0, 0].axis('off')
+ax[0, 0].set_title("Original Image vs. Flipped Image\n"
+                   "(all keypoints and matches)")
+
+plot_matches(ax[1, 0], img1, img3, keypoints1, keypoints3, matches13)
+ax[1, 0].axis('off')
+ax[1, 0].set_title("Original Image vs. Transformed Image\n"
+                   "(all keypoints and matches)")
+
+plot_matches(ax[0, 1], img1, img2, keypoints1, keypoints2, matches12[::15],
+             only_matches=True)
+ax[0, 1].axis('off')
+ax[0, 1].set_title("Original Image vs. Flipped Image\n"
+                   "(subset of matches for visibility)")
+
+plot_matches(ax[1, 1], img1, img3, keypoints1, keypoints3, matches13[::15],
+             only_matches=True)
+ax[1, 1].axis('off')
+ax[1, 1].set_title("Original Image vs. Transformed Image\n"
+                   "(subset of matches for visibility)")
+
+plt.tight_layout()
+plt.show()

skimage/feature/__init__.py

@@ -22,6 +22,7 @@
 from .brief import BRIEF
 from .censure import CENSURE
 from .orb import ORB
+from .sift import SIFT
 from .match import match_descriptors
 from .util import plot_matches
 from .blob import blob_dog, blob_log, blob_doh
@@ -75,6 +76,7 @@ def greycoprops(P, prop='contrast'):
            'BRIEF',
            'CENSURE',
            'ORB',
+           'SIFT',
            'match_descriptors',
            'plot_matches',
            'blob_dog',

skimage/feature/_sift.pyx

@@ -0,0 +1,163 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+
+import numpy as np
+cimport numpy as cnp
+from libc.math cimport M_PI, floor
+
+from .._shared.fused_numerics cimport np_floats
+
+
+cpdef _ori_distances(np_floats[::1] ori_bins,
+                     np_floats[::1] theta):
+    """Compute angular minima and their indices.
+
+    Parameters
+    ----------
+    ori_bins : ndarray
+        The orientation histogram bins.
+    theta : ndarray
+        The orientation computed for each pixel in the patch.
+
+    Returns
+    -------
+    near_t : ndarray
+        The histogram bin index corresponding to each pixel in the patch.
+    near_t_val : ndarray
+        The distance between the orientation of each pixel and the nearest bin.
+    """
+    cdef:
+        Py_ssize_t n_theta = theta.size
+        Py_ssize_t n_ori = ori_bins.size
+        Py_ssize_t i, j, idx_min
+        np_floats th, dist, dist_min
+        np_floats two_pi = 2 * M_PI
+        Py_ssize_t[::1] near_t = np.empty((n_theta, ), dtype=np.intp)
+
+    if np_floats == cnp.float32_t:
+        dtype = np.float32
+    else:
+        dtype = np.float64
+    cdef np_floats[::1] near_t_vals = np.empty((n_theta, ), dtype=dtype)
+
+    for i in range(n_theta):
+        dist_min = 100.
+        th = theta[i]
+        for j in range(n_ori):
+            dist = ori_bins[j] - th
+            if dist > two_pi:
+                dist -= two_pi
+            elif dist < 0:
+                dist += two_pi
+            if dist < dist_min:
+                idx_min = j
+                dist_min = dist
+        near_t[i] = idx_min
+        near_t_vals[i] = dist_min
+    return np.asarray(near_t), np.asarray(near_t_vals)
+
+
+cpdef _update_histogram(np_floats[:, :, ::1] histograms,
+                        const Py_ssize_t[::1] near_t,
+                        np_floats[::1] near_t_val,  # const
+                        np_floats[::1] magnitude,  # const
+                        np_floats[:, ::1] dist_r,  # const
+                        np_floats[:, ::1] dist_c,  # const
+                        np_floats rc_bin_spacing):
+    """Compute an array of orientation histograms (Eq. 28 of Otero et. al.)
+
+    Parameters
+    ----------
+    histograms : (n_hist, n_hist, n_ori) ndarray
+        An array of zeros that will contain the histogram output. `n_ori` is
+        the number of orientation bins and `n_hist` is the number of spatial
+        bins along each axis.
+    near_t : (n_patch,) ndarray
+        The orientation histogram bins obtained from `_ori_distances`.
+    near_t_val : (n_patch,) ndarray
+        The orientation histogram values obtained from `_ori_distances`.
+    magnitude : (n_patch,) ndarray
+        The magnitude weights based on spatial distance.
+    dist_r : (n_hist, n_patch) ndarray
+        Row distances between each point in the patch and the nearest row
+        histogram bin. Shape (n_hist, n_patch).
+    dist_c : (n_hist, n_patch) ndarray
+        Column distances between each point in the patch and the nearest column
+        histogram bin.  Shape (n_hist, n_patch).
+    rc_bin_spacing : float
+        The spacing between spatial histogram bins along a single axis.
+
+    """
+    cdef:
+        Py_ssize_t i, p, r, c, k_index, k_index2
+        Py_ssize_t n_patch = len(magnitude)
+        Py_ssize_t n_hist = histograms.shape[0]
+        Py_ssize_t n_ori = histograms.shape[2]
+        np_floats t_val, val_norm1, w0, w1, w2, inv_spacing
+    inv_spacing = 1 / rc_bin_spacing
+    val_norm1 = n_ori / (2 * M_PI)
+    for p in range(n_patch):
+        for r in range(n_hist):
+            if dist_r[r, p] > rc_bin_spacing:
+                continue
+            for c in range(n_hist):
+                if dist_c[c, p] > rc_bin_spacing:
+                    continue
+                w0 = ((1 - inv_spacing * dist_r[r, p])
+                      * (1 - inv_spacing * dist_c[c, p])
+                      * magnitude[p])
+                t_val = near_t_val[p]
+                w1 = w0 * val_norm1 * t_val
+                w2 = w0 - w1
+                k_index = near_t[p]
+                k_index2 = k_index + 1
+                if k_index2 == n_ori:
+                    k_index2 = 0
+                histograms[r, c, k_index] += w1
+                histograms[r, c, k_index2] += w2
+
+
+cpdef _local_max(np_floats[:, :, ::1] octave, double thresh):
+    cdef:
+        Py_ssize_t n_r = octave.shape[0]
+        Py_ssize_t n_c = octave.shape[1]
+        Py_ssize_t n_s = octave.shape[2]
+        Py_ssize_t r, c, s
+        np_floats* center_ptr
+        Py_ssize_t neighbor_offsets[26]
+        np_floats center_val, val
+        int n = 0
+
+    for r in range(-1, 2):
+        for c in range(-1, 2):
+            for s in range(-1, 2):
+                if (r != 0 or c != 0 or s != 0):
+                    neighbor_offsets[n] = (r * n_c + c) * n_s + s
+                    n += 1
+
+    maxima_coords = []
+    for r in range(1, n_r - 1):
+        for c in range(1, n_c - 1):
+            for s in range(1, n_s - 1):
+                center_ptr = &octave[r, c, s]
+                center_val = center_ptr[0]
+                if abs(center_val) < thresh:
+                    continue
+                is_local_min = True
+                is_local_max = False
+                for offset in neighbor_offsets:
+                    val = center_ptr[offset]
+                    if val <= center_val:
+                        is_local_max = True
+                        is_local_min = False
+                        for offset in neighbor_offsets:
+                            val = center_ptr[offset]
+                            if val >= center_val:
+                                is_local_max = False
+                                break
+                        break
+                if is_local_min or is_local_max:
+                    maxima_coords.append((r, c, s))
+    return np.asarray(maxima_coords, dtype=np.intp)

skimage/feature/setup.py

@@ -19,6 +19,7 @@ def configuration(parent_package='', top_path=None):
             '_texture.pyx',
             '_hessian_det_appx.pyx',
             '_hoghistogram.pyx',
+            '_sift.pyx',
             ], working_path=base_path)
     # _haar uses c++, so it must be cythonized separately
     cython(['_cascade.pyx',
@@ -44,6 +45,8 @@ def configuration(parent_package='', top_path=None):
     config.add_extension('_haar', sources=['_haar.cpp'],
                          include_dirs=[get_numpy_include_dirs(), '../_shared'],
                          language="c++")
+    config.add_extension('_sift', sources=['_sift.c'],
+                         include_dirs=[get_numpy_include_dirs(), '../_shared'])
 
     return config