diff --git a/bento.info b/bento.info index b1366173..9e3bd79b 100644 --- a/bento.info +++ b/bento.info @@ -90,6 +90,9 @@ Library: Extension: skimage.morphology._greyreconstruct Sources: skimage/morphology/_greyreconstruct.pyx + Extension: skimage.feature._brief_cy + Sources: + skimage/feature/_brief_cy.pyx Extension: skimage.feature.corner_cy Sources: skimage/feature/corner_cy.pyx diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 09ac4540..d0c51fb0 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -2,11 +2,13 @@ from ._daisy import daisy from ._hog import hog from .texture import greycomatrix, greycoprops, local_binary_pattern from .peak import peak_local_max -from .corner import (corner_kitchen_rosenfeld, corner_harris, corner_shi_tomasi, - corner_foerstner, corner_subpix, corner_peaks) +from .corner import (corner_kitchen_rosenfeld, corner_harris, + corner_shi_tomasi, corner_foerstner, corner_subpix, + corner_peaks) from .corner_cy import corner_moravec from .template import match_template - +from ._brief import brief, match_keypoints_brief +from .util import pairwise_hamming_distance __all__ = ['daisy', 'hog', @@ -21,4 +23,7 @@ __all__ = ['daisy', 'corner_subpix', 'corner_peaks', 'corner_moravec', - 'match_template'] + 'match_template', + 'brief', + 'pairwise_hamming_distance', + 'match_keypoints_brief'] diff --git a/skimage/feature/_brief.py b/skimage/feature/_brief.py new file mode 100644 index 00000000..22b8e75a --- /dev/null +++ b/skimage/feature/_brief.py @@ -0,0 +1,224 @@ +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 ._brief_cy import _brief_loop + + +def brief(image, keypoints, descriptor_size=256, mode='normal', patch_size=49, + sample_seed=1, variance=2): + """Extract BRIEF Descriptor about given keypoints for a given image. + + Parameters + ---------- + image : 2D ndarray + Input image. + keypoints : (P, 2) ndarray + Array of keypoint locations. + descriptor_size : int + Size of BRIEF descriptor about each keypoint. Sizes 128, 256 and 512 + preferred by the authors. Default is 256. + mode : string + Probability distribution for sampling location of decision pixel-pairs + around keypoints. Default is 'normal' otherwise uniform. + patch_size : int + Length of the two dimensional square patch sampling region around + the keypoints. Default is 49. + sample_seed : int + Seed for sampling the decision pixel-pairs. From a square window with + length patch_size, pixel pairs are sampled using the `mode` parameter + to build the descriptors using intensity comparison. The value of + `sample_seed` should be the same for the images to be matched while + building the descriptors. Default is 1. + variance : float + Variance of the Gaussian Low Pass filter applied on the image to + alleviate noise sensitivity. Default is 2. + + Returns + ------- + descriptors : (Q, `descriptor_size`) ndarray of dtype bool + 2D ndarray of binary descriptors of size `descriptor_size` about Q + keypoints after filtering out border keypoints with value at an index + (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. + + References + ---------- + .. [1] Michael Calonder, Vincent Lepetit, Christoph Strecha and Pascal Fua + "BRIEF : Binary robust independent elementary features", + http://cvlabwww.epfl.ch/~lepetit/papers/calonder_eccv10.pdf + + Examples + -------- + >>> import numpy as np + >>> from skimage.feature.corner import corner_peaks, corner_harris + >>> from skimage.feature import pairwise_hamming_distance, brief, match_keypoints_brief + >>> square1 = np.zeros([8, 8], dtype=np.int32) + >>> square1[2:6, 2:6] = 1 + >>> square1 + array([[0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 1, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0]], dtype=int32) + >>> keypoints1 = corner_peaks(corner_harris(square1), min_distance=1) + >>> keypoints1 + array([[2, 2], + [2, 5], + [5, 2], + [5, 5]]) + >>> descriptors1, keypoints1 = brief(square1, keypoints1, patch_size=5) + >>> keypoints1 + array([[2, 2], + [2, 5], + [5, 2], + [5, 5]]) + >>> square2 = np.zeros([9, 9], dtype=np.int32) + >>> square2[2:7, 2:7] = 1 + >>> square2 + array([[0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 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]], dtype=int32) + >>> keypoints2 = corner_peaks(corner_harris(square2), min_distance=1) + >>> keypoints2 + array([[2, 2], + [2, 6], + [6, 2], + [6, 6]]) + >>> descriptors2, keypoints2 = brief(square2, keypoints2, patch_size=5) + >>> keypoints2 + array([[2, 2], + [2, 6], + [6, 2], + [6, 6]]) + >>> pairwise_hamming_distance(descriptors1, descriptors2) + array([[ 0.03125 , 0.3203125, 0.3671875, 0.6171875], + [ 0.3203125, 0.03125 , 0.640625 , 0.375 ], + [ 0.375 , 0.6328125, 0.0390625, 0.328125 ], + [ 0.625 , 0.3671875, 0.34375 , 0.0234375]]) + >>> match_keypoints_brief(keypoints1, descriptors1, keypoints2, descriptors2) + array([[[ 2, 2], + [ 2, 2]], + + [[ 2, 5], + [ 2, 6]], + + [[ 5, 2], + [ 6, 2]], + + [[ 5, 5], + [ 6, 6]]]) + + """ + np.random.seed(sample_seed) + + image = np.squeeze(image) + if image.ndim != 2: + raise ValueError("Only 2-D gray-scale images supported.") + + image = img_as_float(image) + + # Gaussian Low pass filtering to alleviate noise + # sensitivity + image = gaussian_filter(image, variance) + + image = np.ascontiguousarray(image) + + keypoints = np.array(keypoints + 0.5, dtype=np.intp, order='C') + + # Removing keypoints that are within (patch_size / 2) distance from the + # image border + keypoints = _remove_border_keypoints(image, keypoints, patch_size // 2) + keypoints = np.ascontiguousarray(keypoints) + + descriptors = np.zeros((keypoints.shape[0], descriptor_size), dtype=bool, + order='C') + + # Sampling pairs of decision pixels in patch_size x patch_size window + if mode == 'normal': + + samples = (patch_size / 5.0) * np.random.randn(descriptor_size * 8) + samples = np.array(samples, dtype=np.int32) + samples = samples[(samples < (patch_size // 2)) + & (samples > - (patch_size - 2) // 2)] + + pos1 = samples[:descriptor_size * 2] + pos1 = pos1.reshape(descriptor_size, 2) + pos2 = samples[descriptor_size * 2:descriptor_size * 4] + pos2 = pos2.reshape(descriptor_size, 2) + + else: + + samples = np.random.randint(-(patch_size - 2) // 2, + (patch_size // 2) + 1, + (descriptor_size * 2, 2)) + pos1, pos2 = np.split(samples, 2) + + pos1 = np.ascontiguousarray(pos1) + pos2 = np.ascontiguousarray(pos2) + + _brief_loop(image, descriptors.view(np.uint8), keypoints, pos1, pos2) + + return descriptors, keypoints + + +def match_keypoints_brief(keypoints1, descriptors1, keypoints2, + descriptors2, threshold=0.15): + """Match keypoints described using BRIEF descriptors in one image to + those in second image. + + Parameters + ---------- + keypoints1 : (M, 2) ndarray + M Keypoints from the first image described using skimage.feature.brief + descriptors1 : (M, P) ndarray + BRIEF descriptors of size P about M keypoints in the first image. + keypoints2 : (N, 2) ndarray + N Keypoints from the second image described using skimage.feature.brief + descriptors2 : (N, P) ndarray + BRIEF descriptors of size P about N keypoints in the second image. + threshold : float in range [0, 1] + Maximum allowable hamming distance between descriptors of two keypoints + in separate images to be regarded as a match. Default is 0.15. + + Returns + ------- + match_keypoints_brief : (Q, 2, 2) ndarray + Location of Q matched keypoint pairs from two images. + + """ + if (keypoints1.shape[0] != descriptors1.shape[0] + or keypoints2.shape[0] != descriptors2.shape[0]): + raise ValueError("The number of keypoints and number of described " + "keypoints do not match. Make the optional parameter " + "return_keypoints True to get described keypoints.") + + if descriptors1.shape[1] != descriptors2.shape[1]: + raise ValueError("Descriptor sizes for matching keypoints in both " + "the images should be equal.") + + # Get hamming distances between keeypoints1 and keypoints2 + distance = pairwise_hamming_distance(descriptors1, descriptors2) + + temp = distance > threshold + row_check = np.any(~temp, axis=1) + matched_keypoints2 = keypoints2[np.argmin(distance, axis=1)] + matched_keypoint_pairs = np.zeros((np.sum(row_check), 2, 2), dtype=np.intp) + matched_keypoint_pairs[:, 0, :] = keypoints1[row_check] + matched_keypoint_pairs[:, 1, :] = matched_keypoints2[row_check] + + return matched_keypoint_pairs diff --git a/skimage/feature/_brief_cy.pyx b/skimage/feature/_brief_cy.pyx new file mode 100644 index 00000000..c53d85fc --- /dev/null +++ b/skimage/feature/_brief_cy.pyx @@ -0,0 +1,24 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +cimport numpy as cnp + + +def _brief_loop(double[:, ::1] image, char[:, ::1] descriptors, + Py_ssize_t[:, ::1] keypoints, + int[:, ::1] pos0, int[:, ::1] pos1): + + cdef Py_ssize_t k, d, kr, kc, pr0, pr1, pc0, pc1 + + for p in range(pos0.shape[0]): + pr0 = pos0[p, 0] + pc0 = pos0[p, 1] + pr1 = pos1[p, 0] + pc1 = pos1[p, 1] + for k in range(keypoints.shape[0]): + kr = keypoints[k, 0] + kc = keypoints[k, 1] + if image[kr + pr0, kc + pc0] < image[kr + pr1, kc + pc1]: + descriptors[k, p] = True diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index e769621d..f4c62957 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -13,11 +13,14 @@ def configuration(parent_package='', top_path=None): config.add_data_dir('tests') cython(['corner_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('_brief_cy', sources=['_brief_cy.c'], + include_dirs=[get_numpy_include_dirs()]) config.add_extension('_texture', sources=['_texture.c'], include_dirs=[get_numpy_include_dirs(), '../_shared']) config.add_extension('_template', sources=['_template.c'], diff --git a/skimage/feature/tests/test_brief.py b/skimage/feature/tests/test_brief.py new file mode 100644 index 00000000..b78270f7 --- /dev/null +++ b/skimage/feature/tests/test_brief.py @@ -0,0 +1,83 @@ +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 + + +def test_brief_color_image_unsupported_error(): + """Brief descriptors can be evaluated on gray-scale images only.""" + img = np.zeros((20, 20, 3)) + keypoints = [[7, 5], [11, 13]] + assert_raises(ValueError, brief, img, keypoints) + + +def test_match_keypoints_brief_lena_translation(): + """Test matched keypoints between lena image and its translated version.""" + img = data.lena() + img = rgb2gray(img) + img.shape + tform = tf.SimilarityTransform(scale=1, rotation=0, translation=(15, 20)) + translated_img = tf.warp(img, tform) + + keypoints1 = corner_peaks(corner_harris(img), min_distance=5) + descriptors1, keypoints1 = brief(img, keypoints1, descriptor_size=512) + + keypoints2 = corner_peaks(corner_harris(translated_img), min_distance=5) + descriptors2, keypoints2 = brief(translated_img, keypoints2, + descriptor_size=512) + + matched_keypoints = match_keypoints_brief(keypoints1, descriptors1, + keypoints2, descriptors2, + threshold=0.10) + + assert_array_equal(matched_keypoints[:, 0, :], matched_keypoints[:, 1, :] + + [20, 15]) + + +def test_match_keypoints_brief_lena_rotation(): + """Verify matched keypoints result between lena image and its rotated + version with the expected keypoint pairs.""" + img = data.lena() + img = rgb2gray(img) + img.shape + tform = tf.SimilarityTransform(scale=1, rotation=0.10, translation=(0, 0)) + rotated_img = tf.warp(img, tform) + + keypoints1 = corner_peaks(corner_harris(img), min_distance=5) + descriptors1, keypoints1 = brief(img, keypoints1, descriptor_size=512) + + keypoints2 = corner_peaks(corner_harris(rotated_img), min_distance=5) + descriptors2, keypoints2 = brief(rotated_img, keypoints2, + descriptor_size=512) + + matched_keypoints = match_keypoints_brief(keypoints1, descriptors1, + keypoints2, descriptors2, + threshold=0.07) + + expected = np.array([[[263, 272], + [234, 298]], + + [[271, 120], + [258, 146]], + + [[323, 164], + [305, 195]], + + [[414, 70], + [405, 111]], + + [[435, 181], + [415, 223]], + + [[454, 176], + [435, 221]]]) + + assert_array_equal(matched_keypoints, expected) + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite() diff --git a/skimage/feature/tests/test_util.py b/skimage/feature/tests/test_util.py new file mode 100644 index 00000000..6e25f51a --- /dev/null +++ b/skimage/feature/tests/test_util.py @@ -0,0 +1,32 @@ +import numpy as np +from numpy.testing import assert_array_equal +from skimage.feature.util import pairwise_hamming_distance + + +def test_pairwise_hamming_distance_range(): + """Values of all the pairwise hamming distances should be in the range + [0, 1].""" + a = np.random.random_sample((10, 50)) > 0.5 + b = np.random.random_sample((20, 50)) > 0.5 + dist = pairwise_hamming_distance(a, b) + assert np.all((0 <= dist) & (dist <= 1)) + + +def test_pairwise_hamming_distance_value(): + """The result of pairwise_hamming_distance of two fixed sets of boolean + vectors should be same as expected.""" + np.random.seed(10) + a = np.random.random_sample((4, 100)) > 0.5 + np.random.seed(20) + b = np.random.random_sample((3, 100)) > 0.5 + result = pairwise_hamming_distance(a, b) + expected = np.array([[0.5 , 0.49, 0.44], + [0.44, 0.53, 0.52], + [0.4 , 0.55, 0.5 ], + [0.47, 0.48, 0.57]]) + assert_array_equal(result, expected) + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite() diff --git a/skimage/feature/util.py b/skimage/feature/util.py new file mode 100644 index 00000000..aec4dfc8 --- /dev/null +++ b/skimage/feature/util.py @@ -0,0 +1,36 @@ + + +def _remove_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)] + + return keypoints + + +def pairwise_hamming_distance(array1, array2): + """Calculate hamming dissimilarity measure between two sets of + vectors. + + Parameters + ---------- + array1 : (P1, D) array + P1 vectors of size D. + array2 : (P2, D) array + P2 vectors of size D. + + Returns + ------- + distance : (P1, P2) array of dtype float + 2D ndarray with value at an index (i, j) representing the hamming + distance in the range [0, 1] between ith vector in array1 and jth + vector in array2. + + """ + distance = (array1[:, None] != array2[None]).mean(axis=2) + return distance