diff --git a/bento.info b/bento.info index 648445d6..acdb155d 100644 --- a/bento.info +++ b/bento.info @@ -90,6 +90,9 @@ Library: Extension: skimage.morphology._greyreconstruct Sources: skimage/morphology/_greyreconstruct.pyx + Extension: skimage.feature.censure_cy + Sources: + skimage/feature/censure_cy.pyx Extension: skimage.feature._brief_cy Sources: skimage/feature/_brief_cy.pyx diff --git a/doc/examples/plot_censure_keypoints.py b/doc/examples/plot_censure_keypoints.py new file mode 100644 index 00000000..fbba4e1b --- /dev/null +++ b/doc/examples/plot_censure_keypoints.py @@ -0,0 +1,54 @@ +""" +========================= +CenSurE Feature Detection +========================= + +In this example we detect and plot the CenSurE (Center Surround Extrema) +features at various scales using Difference of Boxes, Octagon and Star shaped +bi-level filters. + +""" + +from skimage.feature import keypoints_censure +from skimage.data import lena +from skimage.color import rgb2gray +import matplotlib.pyplot as plt + +# Initializing the parameters for Censure keypoints +img = lena() +gray_img = rgb2gray(img) +min_scale = 2 +max_scale = 6 +non_max_threshold = 0.15 +line_threshold = 10 + + +_, ax = plt.subplots(nrows=(max_scale - min_scale - 1), ncols=3, + figsize=(6, 6)) +plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.94, + bottom=0.02, left=0.06, right=0.98) + +# Detecting Censure keypoints for the following filters +for col, mode in enumerate(['dob', 'octagon', 'star']): + + ax[0, col].set_title(mode.upper(), fontsize=12) + + keypoints, scales = keypoints_censure(gray_img, min_scale, max_scale, + mode, non_max_threshold, + line_threshold) + + # Plotting Censure features at all the scales + for row, scale in enumerate(range(min_scale + 1, max_scale)): + mask = scales == scale + x = keypoints[mask, 1] + y = keypoints[mask, 0] + s = 0.5 * 2 ** (scale + min_scale + 1) + ax[row, col].imshow(img) + ax[row, col].scatter(x, y, s, facecolors='none', edgecolors='b') + ax[row, col].set_xticks([]) + ax[row, col].set_yticks([]) + ax[row, col].axis((0, img.shape[1], img.shape[0], 0)) + if col == 0: + ax[row, col].set_ylabel('Scale %d' % scale, fontsize=12) + +plt.show() diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index d0c51fb0..bde81e59 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -9,6 +9,7 @@ from .corner_cy import corner_moravec from .template import match_template from ._brief import brief, match_keypoints_brief from .util import pairwise_hamming_distance +from .censure import keypoints_censure __all__ = ['daisy', 'hog', @@ -26,4 +27,5 @@ __all__ = ['daisy', 'match_template', 'brief', 'pairwise_hamming_distance', - 'match_keypoints_brief'] + 'match_keypoints_brief', + 'keypoints_censure'] diff --git a/skimage/feature/_brief.py b/skimage/feature/_brief.py index 22b8e75a..27a8d085 100644 --- a/skimage/feature/_brief.py +++ b/skimage/feature/_brief.py @@ -2,7 +2,7 @@ import numpy as np from scipy.ndimage.filters import gaussian_filter from ..util import img_as_float -from .util import _remove_border_keypoints, pairwise_hamming_distance +from .util import _mask_border_keypoints, pairwise_hamming_distance from ._brief_cy import _brief_loop @@ -16,7 +16,7 @@ def brief(image, keypoints, descriptor_size=256, mode='normal', patch_size=49, image : 2D ndarray Input image. keypoints : (P, 2) ndarray - Array of keypoint locations. + Array of keypoint locations in the format (row, col). descriptor_size : int Size of BRIEF descriptor about each keypoint. Sizes 128, 256 and 512 preferred by the authors. Default is 256. @@ -44,8 +44,8 @@ def brief(image, keypoints, descriptor_size=256, mode='normal', patch_size=49, (i, j) either being True or False representing the outcome of Intensity comparison about ith keypoint on jth decision pixel-pair. keypoints : (Q, 2) ndarray - Keypoints after removing out those that are near border. - Returned only if return_keypoints is True. + Location i.e. (row, col) of keypoints after removing out those that + are near border. References ---------- @@ -142,7 +142,7 @@ def brief(image, keypoints, descriptor_size=256, mode='normal', patch_size=49, # Removing keypoints that are within (patch_size / 2) distance from the # image border - keypoints = _remove_border_keypoints(image, keypoints, patch_size // 2) + keypoints = keypoints[_mask_border_keypoints(image, keypoints, patch_size // 2)] keypoints = np.ascontiguousarray(keypoints) descriptors = np.zeros((keypoints.shape[0], descriptor_size), dtype=bool, diff --git a/skimage/feature/censure.py b/skimage/feature/censure.py new file mode 100644 index 00000000..00be16a8 --- /dev/null +++ b/skimage/feature/censure.py @@ -0,0 +1,233 @@ +import numpy as np +from scipy.ndimage.filters import maximum_filter, minimum_filter, convolve + +from skimage.transform import integral_image +from skimage.feature.corner import _compute_auto_correlation +from skimage.util import img_as_float +from skimage.morphology import octagon, star +from skimage.feature.util import _mask_border_keypoints + +from skimage.feature.censure_cy import _censure_dob_loop + + +# The paper(Reference [1]) mentions the sizes of the Octagon shaped filter +# kernel for the first seven scales only. The sizes of the later scales +# have been extrapolated based on the following statement in the paper. +# "These octagons scale linearly and were experimentally chosen to correspond +# to the seven DOBs described in the previous section." +OCTAGON_OUTER_SHAPE = [(5, 2), (5, 3), (7, 3), (9, 4), (9, 7), (13, 7), + (15, 10), (15, 11), (15, 12), (17, 13), (17, 14)] +OCTAGON_INNER_SHAPE = [(3, 0), (3, 1), (3, 2), (5, 2), (5, 3), (5, 4), (5, 5), + (7, 5), (7, 6), (9, 6), (9, 7)] + +# The sizes for the STAR shaped filter kernel for different scales have been +# taken from the OpenCV implementation. +STAR_SHAPE = [1, 2, 3, 4, 6, 8, 11, 12, 16, 22, 23, 32, 45, 46, 64, 90, 128] +STAR_FILTER_SHAPE = [(1, 0), (3, 1), (4, 2), (5, 3), (7, 4), (8, 5), + (9, 6), (11, 8), (13, 10), (14, 11), (15, 12), (16, 14)] + + +def _filter_image(image, min_scale, max_scale, mode): + + response = np.zeros((image.shape[0], image.shape[1], + max_scale - min_scale + 1), dtype=np.double) + + if mode == 'dob': + + # make response[:, :, i] contiguous memory block + item_size = response.itemsize + response.strides = (item_size * response.shape[0], item_size, + item_size * response.shape[0] * response.shape[1]) + + integral_img = integral_image(image) + + for i in range(max_scale - min_scale + 1): + n = min_scale + i + + # Constant multipliers for the outer region and the inner region + # of the bi-level filters with the constraint of keeping the + # DC bias 0. + inner_weight = (1.0 / (2 * n + 1)**2) + outer_weight = (1.0 / (12 * n**2 + 4 * n)) + + _censure_dob_loop(n, integral_img, response[:, :, i], + inner_weight, outer_weight) + + # NOTE : For the Octagon shaped filter, we implemented and evaluated the + # slanted integral image based image filtering but the performance was + # more or less equal to image filtering using + # scipy.ndimage.filters.convolve(). Hence we have decided to use the + # later for a much cleaner implementation. + elif mode == 'octagon': + # TODO : Decide the shapes of Octagon filters for scales > 7 + + for i in range(max_scale - min_scale + 1): + mo, no = OCTAGON_OUTER_SHAPE[min_scale + i - 1] + mi, ni = OCTAGON_INNER_SHAPE[min_scale + i - 1] + response[:, :, i] = convolve(image, + _octagon_filter_kernel(mo, no, mi, ni)) + + elif mode == 'star': + + for i in range(max_scale - min_scale + 1): + m = STAR_SHAPE[STAR_FILTER_SHAPE[min_scale + i - 1][0]] + n = STAR_SHAPE[STAR_FILTER_SHAPE[min_scale + i - 1][1]] + response[:, :, i] = convolve(image, _star_filter_kernel(m, n)) + + return response + + +def _octagon_filter_kernel(mo, no, mi, ni): + outer = (mo + 2 * no)**2 - 2 * no * (no + 1) + inner = (mi + 2 * ni)**2 - 2 * ni * (ni + 1) + outer_weight = 1.0 / (outer - inner) + inner_weight = 1.0 / inner + c = ((mo + 2 * no) - (mi + 2 * ni)) // 2 + outer_oct = octagon(mo, no) + inner_oct = np.zeros((mo + 2 * no, mo + 2 * no)) + inner_oct[c: -c, c: -c] = octagon(mi, ni) + bfilter = (outer_weight * outer_oct - + (outer_weight + inner_weight) * inner_oct) + return bfilter + + +def _star_filter_kernel(m, n): + c = m + m // 2 - n - n // 2 + outer_star = star(m) + inner_star = np.zeros_like(outer_star) + inner_star[c: -c, c: -c] = star(n) + outer_weight = 1.0 / (np.sum(outer_star - inner_star)) + inner_weight = 1.0 / np.sum(inner_star) + bfilter = (outer_weight * outer_star - + (outer_weight + inner_weight) * inner_star) + return bfilter + + +def _suppress_lines(feature_mask, image, sigma, line_threshold): + Axx, Axy, Ayy = _compute_auto_correlation(image, sigma) + feature_mask[(Axx + Ayy) * (Axx + Ayy) + > line_threshold * (Axx * Ayy - Axy * Axy)] = False + + +def keypoints_censure(image, min_scale=1, max_scale=7, mode='DoB', + non_max_threshold=0.15, line_threshold=10): + """ + Extracts CenSurE keypoints along with the corresponding scale using + either Difference of Boxes, Octagon or STAR bi-level filter. + + Parameters + ---------- + image : 2D ndarray + Input image. + min_scale : int + Minimum scale to extract keypoints from. + max_scale : int + Maximum scale to extract keypoints from. The keypoints will be + extracted from all the scales except the first and the last i.e. + from the scales in the range [min_scale + 1, max_scale - 1]. + mode : {'DoB', 'Octagon', 'STAR'} + Type of bi-level filter used to get the scales of the input image. + Possible values are 'DoB', 'Octagon' and 'STAR'. The three modes + represent the shape of the bi-level filters i.e. box(square), octagon + and star respectively. For instance, a bi-level octagon filter consists + of a smaller inner octagon and a larger outer octagon with the filter + weights being uniformly negative in both the inner octagon while + uniformly positive in the difference region. Use STAR and Octagon for + better features and DoB for better performance. + non_max_threshold : float + Threshold value used to suppress maximas and minimas with a weak + magnitude response obtained after Non-Maximal Suppression. + line_threshold : float + Threshold for rejecting interest points which have ratio of principal + curvatures greater than this value. + + Returns + ------- + keypoints : (N, 2) array + Location of the extracted keypoints in the ``(row, col)`` format. + scales : (N, 1) array + The corresponding scale of the N extracted keypoints. + + References + ---------- + .. [1] Motilal Agrawal, Kurt Konolige and Morten Rufus Blas + "CenSurE: Center Surround Extremas for Realtime Feature + Detection and Matching", + http://link.springer.com/content/pdf/10.1007%2F978-3-540-88693-8_8.pdf + + .. [2] Adam Schmidt, Marek Kraft, Michal Fularz and Zuzanna Domagala + "Comparative Assessment of Point Feature Detectors and + Descriptors in the Context of Robot Navigation" + http://www.jamris.org/01_2013/saveas.php?QUEST=JAMRIS_No01_2013_P_11-20.pdf + + """ + + # (1) First we generate the required scales on the input grayscale image + # using a bi-level filter and stack them up in `filter_response`. + # (2) We then perform Non-Maximal suppression in 3 x 3 x 3 window on the + # filter_response to suppress points that are neither minima or maxima in + # 3 x 3 x 3 neighbourhood. We obtain a boolean ndarray `feature_mask` + # containing all the minimas and maximas in `filter_response` as True. + # (3) Then we suppress all the points in the `feature_mask` for which the + # corresponding point in the image at a particular scale has the ratio of + # principal curvatures greater than `line_threshold`. + # (4) Finally, we remove the border keypoints and return the keypoints + # along with its corresponding scale. + + image = np.squeeze(image) + if image.ndim != 2: + raise ValueError("Only 2-D gray-scale images supported.") + + mode = mode.lower() + if mode not in ('dob', 'octagon', 'star'): + raise ValueError('Mode must be one of "DoB", "Octagon", "STAR".') + + if min_scale < 1 or max_scale < 1 or max_scale - min_scale < 2: + raise ValueError('The scales must be >= 1 and the number of scales ' + 'should be >= 3.') + + image = img_as_float(image) + image = np.ascontiguousarray(image) + + # Generating all the scales + filter_response = _filter_image(image, min_scale, max_scale, mode) + + # Suppressing points that are neither minima or maxima in their 3 x 3 x 3 + # neighbourhood to zero + minimas = minimum_filter(filter_response, (3, 3, 3)) == filter_response + maximas = maximum_filter(filter_response, (3, 3, 3)) == filter_response + + feature_mask = minimas | maximas + feature_mask[filter_response < non_max_threshold] = False + + for i in range(1, max_scale - min_scale): + # sigma = (window_size - 1) / 6.0, so the window covers > 99% of the + # kernel's distribution + # window_size = 7 + 2 * (min_scale - 1 + i) + # Hence sigma = 1 + (min_scale - 1 + i)/ 3.0 + _suppress_lines(feature_mask[:, :, i], image, + (1 + (min_scale + i - 1) / 3.0), line_threshold) + + rows, cols, scales = np.nonzero(feature_mask[..., 1:max_scale - min_scale]) + keypoints = np.column_stack([rows, cols]) + scales = scales + min_scale + 1 + + if mode == 'dob': + return keypoints, scales + + cumulative_mask = np.zeros(keypoints.shape[0], dtype=np.bool) + + if mode == 'octagon': + for i in range(min_scale + 1, max_scale): + c = (OCTAGON_OUTER_SHAPE[i - 1][0] - 1) // 2 \ + + OCTAGON_OUTER_SHAPE[i - 1][1] + cumulative_mask |= _mask_border_keypoints(image, keypoints, c) \ + & (scales == i) + elif mode == 'star': + for i in range(min_scale + 1, max_scale): + c = STAR_SHAPE[STAR_FILTER_SHAPE[i - 1][0]] \ + + STAR_SHAPE[STAR_FILTER_SHAPE[i - 1][0]] // 2 + cumulative_mask |= _mask_border_keypoints(image, keypoints, c) \ + & (scales == i) + + return keypoints[cumulative_mask], scales[cumulative_mask] diff --git a/skimage/feature/censure_cy.pyx b/skimage/feature/censure_cy.pyx new file mode 100644 index 00000000..1c352fde --- /dev/null +++ b/skimage/feature/censure_cy.pyx @@ -0,0 +1,72 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + + +def _censure_dob_loop(Py_ssize_t n, + double[:, ::1] integral_img, + double[:, ::1] filtered_image, + double inner_weight, double outer_weight): + # This function calculates the value in the DoB filtered image using + # integral images. If r = right. l = left, u = up, d = down, the sum of + # pixel values in the rectangle formed by (u, l), (u, r), (d, r), (d, l) + # is calculated as I(d, r) + I(u - 1, l - 1) - I(u - 1, r) - I(d, l - 1). + + cdef Py_ssize_t i, j + cdef double inner, outer + cdef Py_ssize_t n2 = 2 * n + cdef double total_weight = inner_weight + outer_weight + + # top-left pixel + inner = (integral_img[n2 + n, n2 + n] + + integral_img[n2 - n - 1, n2 - n - 1] + - integral_img[n2 + n, n2 - n - 1] + - integral_img[n2 - n - 1, n2 + n]) + + outer = integral_img[2 * n2, 2 * n2] + + filtered_image[n2, n2] = (outer_weight * outer + - total_weight * inner) + + # left column + for i in range(n2 + 1, integral_img.shape[0] - n2): + inner = (integral_img[i + n, n2 + n] + + integral_img[i - n - 1, n2 - n - 1] + - integral_img[i + n, n2 - n - 1] + - integral_img[i - n - 1, n2 + n]) + + outer = (integral_img[i + n2, 2 * n2] + - integral_img[i - n2 - 1, 2 * n2]) + + filtered_image[i, n2] = (outer_weight * outer + - total_weight * inner) + + # top row + for j in range(n2 + 1, integral_img.shape[1] - n2): + inner = (integral_img[n2 + n, j + n] + + integral_img[n2 - n - 1, j - n - 1] + - integral_img[n2 + n, j - n - 1] + - integral_img[n2 - n - 1, j + n]) + + outer = (integral_img[2 * n2, j + n2] + - integral_img[2 * n2, j - n2 - 1]) + + filtered_image[n2, j] = (outer_weight * outer + - total_weight * inner) + + # remaining block + for i in range(n2 + 1, integral_img.shape[0] - n2): + for j in range(n2 + 1, integral_img.shape[1] - n2): + inner = (integral_img[i + n, j + n] + + integral_img[i - n - 1, j - n - 1] + - integral_img[i + n, j - n - 1] + - integral_img[i - n - 1, j + n]) + + outer = (integral_img[i + n2, j + n2] + + integral_img[i - n2 - 1, j - n2 - 1] + - integral_img[i + n2, j - n2 - 1] + - integral_img[i - n2 - 1, j + n2]) + + filtered_image[i, j] = (outer_weight * outer + - total_weight * inner) diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index f4c62957..7df64c32 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -13,12 +13,15 @@ def configuration(parent_package='', top_path=None): config.add_data_dir('tests') cython(['corner_cy.pyx'], working_path=base_path) + cython(['censure_cy.pyx'], working_path=base_path) cython(['_brief_cy.pyx'], working_path=base_path) cython(['_texture.pyx'], working_path=base_path) cython(['_template.pyx'], working_path=base_path) config.add_extension('corner_cy', sources=['corner_cy.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('censure_cy', sources=['censure_cy.c'], + include_dirs=[get_numpy_include_dirs()]) config.add_extension('_brief_cy', sources=['_brief_cy.c'], include_dirs=[get_numpy_include_dirs()]) config.add_extension('_texture', sources=['_texture.c'], diff --git a/skimage/feature/tests/test_brief.py b/skimage/feature/tests/test_brief.py index b78270f7..1d26cbbd 100644 --- a/skimage/feature/tests/test_brief.py +++ b/skimage/feature/tests/test_brief.py @@ -2,9 +2,9 @@ import numpy as np from numpy.testing import assert_array_equal, assert_raises from skimage import data from skimage import transform as tf -from skimage.feature.corner import corner_peaks, corner_harris from skimage.color import rgb2gray -from skimage.feature import brief, match_keypoints_brief +from skimage.feature import (brief, match_keypoints_brief, corner_peaks, + corner_harris) def test_brief_color_image_unsupported_error(): diff --git a/skimage/feature/tests/test_censure.py b/skimage/feature/tests/test_censure.py new file mode 100644 index 00000000..4cd2ad68 --- /dev/null +++ b/skimage/feature/tests/test_censure.py @@ -0,0 +1,89 @@ +import numpy as np +from numpy.testing import assert_array_equal, assert_raises +from skimage.data import moon +from skimage.feature import keypoints_censure + + +def test_keypoints_censure_color_image_unsupported_error(): + """Censure keypoints can be extracted from gray-scale images only.""" + img = np.zeros((20, 20, 3)) + assert_raises(ValueError, keypoints_censure, img) + + +def test_keypoints_censure_mode_validity_error(): + """Mode argument in keypoints_censure can be either DoB, Octagon or + STAR.""" + img = np.zeros((20, 20)) + assert_raises(ValueError, keypoints_censure, img, mode='dummy') + + +def test_keypoints_censure_scale_range_error(): + """Difference between the the max_scale and min_scale parameters in + keypoints_censure should be greater than or equal to two.""" + img = np.zeros((20, 20)) + assert_raises(ValueError, keypoints_censure, img, min_scale=1, max_scale=2) + + +def test_keypoints_censure_moon_image_dob(): + """Verify the actual Censure keypoints and their corresponding scale with + the expected values for DoB filter.""" + img = moon() + actual_kp_dob, actual_scale = keypoints_censure(img, 1, 7, 'DoB', 0.15) + expected_kp_dob = np.array([[ 21, 497], + [ 36, 46], + [119, 350], + [185, 177], + [287, 250], + [357, 239], + [463, 116], + [464, 132], + [467, 260]]) + expected_scale = np.array([3, 4, 4, 2, 2, 3, 2, 2, 2]) + + assert_array_equal(expected_kp_dob, actual_kp_dob) + assert_array_equal(expected_scale, actual_scale) + + +def test_keypoints_censure_moon_image_octagon(): + """Verify the actual Censure keypoints and their corresponding scale with + the expected values for Octagon filter.""" + img = moon() + actual_kp_octagon, actual_scale = keypoints_censure(img, 1, 7, 'Octagon', + 0.15) + expected_kp_octagon = np.array([[ 21, 496], + [ 35, 46], + [287, 250], + [356, 239], + [463, 116]]) + + expected_scale = np.array([3, 4, 2, 2, 2]) + + assert_array_equal(expected_kp_octagon, actual_kp_octagon) + assert_array_equal(expected_scale, actual_scale) + + +def test_keypoints_censure_moon_image_star(): + """Verify the actual Censure keypoints and their corresponding scale with + the expected values for STAR filter.""" + img = moon() + actual_kp_star, actual_scale = keypoints_censure(img, 1, 7, 'STAR', 0.15) + expected_kp_star = np.array([[ 21, 497], + [ 36, 46], + [117, 356], + [185, 177], + [260, 227], + [287, 250], + [357, 239], + [451, 281], + [463, 116], + [467, 260]]) + + expected_scale = np.array([3, 3, 6, 2, 3, 2, 3, 5, 2, 2]) + + assert_array_equal(expected_kp_star, actual_kp_star) + assert_array_equal(expected_scale, actual_scale) + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite() diff --git a/skimage/feature/util.py b/skimage/feature/util.py index aec4dfc8..eb3817e8 100644 --- a/skimage/feature/util.py +++ b/skimage/feature/util.py @@ -1,16 +1,16 @@ -def _remove_border_keypoints(image, keypoints, dist): +def _mask_border_keypoints(image, keypoints, dist): """Removes keypoints that are within dist pixels from the image border.""" width = image.shape[0] height = image.shape[1] - keypoints = keypoints[(dist - 1 < keypoints[:, 0]) - & (keypoints[:, 0] < width - dist + 1) - & (dist - 1 < keypoints[:, 1]) - & (keypoints[:, 1] < height - dist + 1)] + keypoints_filtering_mask = ((dist - 1 < keypoints[:, 0]) & + (keypoints[:, 0] < width - dist + 1) & + (dist - 1 < keypoints[:, 1]) & + (keypoints[:, 1] < height - dist + 1)) - return keypoints + return keypoints_filtering_mask def pairwise_hamming_distance(array1, array2): diff --git a/skimage/morphology/selem.py b/skimage/morphology/selem.py index f85b3a95..4df52e94 100644 --- a/skimage/morphology/selem.py +++ b/skimage/morphology/selem.py @@ -280,10 +280,10 @@ def star(a, dtype=np.uint8): bfilter[:] = 1 return bfilter m = 2 * a + 1 - n = a / 2 + n = a // 2 selem_square = np.zeros((m + 2 * n, m + 2 * n)) selem_square[n: m + n, n: m + n] = 1 - c = (m + 2 * n - 1) / 2 + c = (m + 2 * n - 1) // 2 selem_rotated = np.zeros((m + 2 * n, m + 2 * n)) selem_rotated[0, c] = selem_rotated[-1, c] = selem_rotated[c, 0] = selem_rotated[c, -1] = 1 selem_rotated = convex_hull_image(selem_rotated).astype(int)